1

I have to calculate 2 body gravity system (Earth & Moon) without vpython.

Purpose of this code is practicing numerical calculation and there is my code.

import numpy as np
from math import *

from astropy.constants import *

import matplotlib.pyplot as plt

import time

start_time = time.time()

"""
G = Gravitational constant
g0 = Standard acceleration of gravity ( 9.8 m/s2)
M_sun = Solar mass
M_earth = Earth mass
R_sun = Solar darius
R_earth = Earth equatorial radius
au = Astronomical unit


"""

M_moon = 7.342E22
R_moon = 1.737E6
# Mean radius of moon.


M_earth = M_earth.value
R_earth = R_earth.value
G = G.value

perigee, apogee = 3.626E8, 4.054E8

position_E = np.array([0,0])
position_M = np.array([(perigee+apogee)/2.,0])
position_com = (M_earth*position_E+M_moon*position_M)/(M_earth+M_moon)

rel_pE = position_E - position_com
rel_pM = position_M - position_com


p_E = {"x":rel_pE[0], "y":rel_pE[1],"v_x":0, "v_y":10}
p_M = {"x":rel_pM[0], "y":rel_pM[1],"v_x":0, "v_y":-100}

t = range(0,365)

data_E , data_M = [0]*len(t), [0]*len(t)

def s(initial_velocity, acceleration, time):
    result = initial_velocity*time + 0.5*acceleration*time**2
    return result

def v(initial_velocity, acceleration, time):
    result = initial_velocity + acceleration*time
    return result

dist = float(sqrt((p_E["x"]-p_M['x'])**2 + (p_E["y"]-p_M["y"])**2))

# position data of Earth and Moon. make new list to make easy to draw plot

xE=[]
yE=[]
xM=[]
yM=[]

for i in t:
    dist = float(sqrt((p_E["x"]-p_M["x"])**2 + (p_E["y"]-p_M["y"])**2))


    a_Ex = -G*M_moon*p_E["x"]/(dist**2)
    a_Ey = -G*M_moon*p_E["y"]/(dist**2)

    data_E[i] = p_E

    p_E["x"] += s(p_E['v_x'], a_Ex, 24*3600)
    p_E["v_x"] += v(p_E['v_x'], a_Ex, 24*3600)
    p_E["y"] += s(p_E['v_y'], a_Ey, 24*3600)
    p_E["v_y"] += v(p_E['v_y'], a_Ey, 24*3600)

    xE += [p_E["x"]]
    yE += [p_E["y"]]

    a_Mx = -G*M_earth*p_M["x"]/(dist**2)
    a_My = -G*M_earth*p_M["y"]/(dist**2)

    data_M[i] = p_M

    p_M["x"] += s(p_M['v_x'], a_Mx, 24*3600)
    p_M["v_x"] += v(p_M['v_x'], a_Mx, 24*3600)
    p_M["y"] += s(p_M['v_y'], a_My, 24*3600)
    p_M["v_y"] += v(p_M['v_y'], a_My, 24*3600)

    xM += [p_M["x"]]
    yM += [p_M["y"]]

print("\n Run time \n --- %d seconds ---" %(time.time()-start_time))

But this code makes both of x&y increase.

It doesn't show a elliptical orbit. How can I fix my code or write a new code describe perfect elliptical orbit about Earth & Moon gravity system.

Thanks for your help!

Gangil Seo
  • 195
  • 1
  • 1
  • 13
  • I want to know how can i calculate the orbit by the gravitational force. So ellipse function is not what I want to know. – Gangil Seo Jun 03 '17 at 08:36
  • 1
    Just a thought but have you tried reformulating your x, y positions in polar coordinates? If you iterate theta from 0 to 2 pi, the magnitude of r should change as a function of every x, y pair; then you could do a polar plot. Even if this doesn't solve your problem, you can still use print statements to see where the orbit isn't but should be elliptical. Will look more into it when by a computer (on mobile currently). –  Jun 03 '17 at 09:10
  • Why are you doing `+= v(` ? Your velocity calculation already includes the initial velocity. However, using those simple formulas for velocity and distance is not a good idea. You need a better integration technique. For orbit calculations I recommend [Synchronised Leapfrog](https://en.wikipedia.org/wiki/Leapfrog_integration) because it conserves energy. BTW, it's not a good idea to use `time` as a variable name when you're already using it to refer to the `time` module. – PM 2Ring Jun 03 '17 at 09:19
  • @GangilSeo using Newton D'Lambert and want to change circular orbit to elliptical? simple just change the initial speed leaving position as is the gravity will do the rest. However to obtain realistic data this is far from enough as Earth-Moon-Sun is very specific system biased not only by gravity but by liquid nature of are planet too ... Also you might want to look at [Is it possible to make realistic n-body solar system simulation in matter of size and mass?](https://stackoverflow.com/a/28020934/2521214) – Spektre Jun 07 '17 at 08:16
  • I don't have any idea to solve this problem. My code given increase value not circular value. How can I do this – Gangil Seo Jun 10 '17 at 17:01
  • @GangilSeo use Newton D'Lambert physics... so in each `dt` time step compute the gravitational forces for each body `F`, transform it to actual acceleration vector `a=F/m` update actual speed `v += a*dt` update actual position `p += v*dt` and render ... What is the problem ? – Spektre Jun 12 '17 at 07:02

0 Answers0