I am trying to build a function that will calculate Newton's law of gravitation (G*M/r^2) of an object for a given timestep while being able to change the position and velocity of that point.
As an example I will try to show what I have done to set the initial velocity and starting position of a point (I will however at some point need to be able to display the point in real-time and move it around while the loops are running).
I have written two examples that try to demonstrate the trouble I am having.
If you look closely at the plot you might notice that those two lines don't overlap (even though all starting params were the same).
This plot shows when all starting params = 0
printed output:
Concept 1 (calculated intervals = 1 second):
v:9.82 pos:9.82
v:19.64 pos:29.46
v:29.46 pos:58.93
v:39.28 pos:98.22
Concept 2 (calculated intervals = 1 second):
v:0.0 pos:0.0
v:9.82 pos:4.91
v:19.64 pos:19.64
v:29.46 pos:44.20
According to this online calculator it indicates Concept1 is incorrect, and Concept2 is very close to actual.
The code below plots both examples. And I have added two variables that allow for setting the initial_velocity and initial_position (please see the additional two plot images below for details).
import numpy as np
import matplotlib.pyplot as plt
G = 6.67408 * 10 ** -11 # Newton meters^2/kg^2
Me = 5.9736 * 10 ** 24 # mass of earth in kg
Re = 6.371 * 10 ** 6 # radius of earth in meters
def get_g(mass, distance, seconds):
_g = G * mass / distance ** 2
_d = (_g * (seconds ** 2)) / 2
_v = _g * seconds
return _v, _d, _g
def Concept1(seconds=1, timestep=0.1):
plt_x = np.linspace(0, seconds, (seconds/timestep))
plt_y = np.zeros_like(plt_x)
"""
Simple concept demonstrating values that are very close to actual.
Time = 0.1 vel = 0.9822280588290989 pos = 0.049111402941454954 _g = 9.822280588290988
Time = 0.2 vel = 1.9644561176581978 pos = 0.19644561176581982 _g = 9.822280588290988
Time = 0.3 vel = 2.946684176487297 pos = 0.44200262647309463 _g = 9.822280588290988
Time = 0.4 vel = 3.9289122353163957 pos = 0.7857824470632793 _g = 9.822280588290988
"""
for _i in range(len(plt_y)):
t = (_i+1)
_v, _d, _g = get_g(Me, Re, t*timestep)
plt_y[_i] = _d
#print("Time = ", t*timestep, " vel = ", _v, " pos = ", _d, " _g = ", _g)
return plt_x,plt_y
def Concept2(seconds=1, timestep=0.1, initial_velocity=0, initial_position=0):
plt_x = np.linspace(0, seconds, seconds/timestep)
plt_y = np.zeros_like(plt_x, dtype=np.float32)
"""
Dynamic concept showing variable velocity and starting position (output is not accurate)
Time = 0.1 vel = 0.9822280588290989 pos = 0.0982228058829099 _g = 9.822280588290988
Time = 0.1 vel = 1.9644561176581978 pos = 0.2946684176487297 _g = 9.822280588290988
Time = 0.1 vel = 2.9466841764872966 pos = 0.5893368352974594 _g = 9.822280588290988
Time = 0.1 vel = 3.9289122353163957 pos = 0.982228058829099 _g = 9.822280588290988
"""
vel = initial_velocity * timestep
pos = initial_position
for _i in range(len(plt_y)):
t = timestep
_v, _d, _g = get_g(Me, Re, t)
vel += _g * timestep
pos += vel * timestep
plt_y[_i] = pos
#print("Time = ", t, " vel = ", vel, " pos = ", pos, " _g = ", _g)
return plt_x, plt_y
def Concept3(seconds=1, timestep=0.1, initial_velocity=0, initial_position=0):
plt_x = np.linspace(0, seconds, seconds/timestep)
plt_y = np.zeros_like(plt_x, dtype=np.float32)
"""
Dynamic concept showing variable velocity and starting position (output is accurate)
Time = 0.1 vel = 0.982228058829099 pos = 0.049111402941454954 _g = 9.822280588290988
Time = 0.1 vel = 1.964456117658198 pos = 0.19644561176581982 _g = 9.822280588290988
Time = 0.1 vel = 2.9466841764872975 pos = 0.44200262647309463 _g = 9.822280588290988
Time = 0.1 vel = 3.928912235316396 pos = 0.7857824470632793 _g = 9.822280588290988
"""
vel = initial_velocity * timestep
pos = initial_position
gforce = 0
for _i in range(len(plt_y)):
t = timestep
_v, _d, _g = get_g(Me, Re, t)
vel += (gforce + _g)*timestep**2/2
pos += vel
gforce = _g
plt_y[_i] = pos
#print("Time = ", t, " vel = ", (vel + _d)/timestep, " pos = ", pos, " _g = ", gforce)
return plt_x, plt_y
seconds=1
timestep=0.1
initial_velocity = 0
initial_position = 0
x1, y1 = Concept1(seconds=seconds, timestep=timestep)
x2, y2 = Concept2(seconds=seconds, timestep=timestep, initial_velocity=initial_velocity, initial_position=initial_position)
x3, y3 = Concept3(seconds=seconds, timestep=timestep, initial_velocity=initial_velocity, initial_position=initial_position)
plt.xlabel("Time (Seconds)")
plt.ylabel("Distance")
plt.plot(x1, y1,'g', label='Concept 1 (gauge)')
plt.plot(x2, y2,'b', label='Concept 2')
plt.plot(x3, y3,'r', label='Concept 3')
plt.legend(loc='upper left')
plt.show()
I am relatively new to physics (and this is a learning experience). I have been wondering if possibly my approach is incorrect. And if so I really would like to learn more about how to accomplish this. All suggestions will be gladly noted. And thank you!
Plot showing initial_velocity = 2000
Plot showing initial_position = 100000
[UPDATE]
Added a Concept3 to the code.
I "think" it achieves what I needed by storing the current gforce value in a variable for each loop iteration.
I then add it to the current gforce value and divide that by 2. After which I update the current velocity.
Here's the output after making the alterations (seconds=1 timestep=0.1):
