N

#### Nzjeux

##### Guest

My main problem at the moment is that the orbiting satellite is oscaliting around the target planet when it should be in a stable 320km circular orbit.

I have made four different functions for four different integrations. Euler, Leapfrog, Verlet and RK4.I'd expect Euler & leapfrog to have some inaccuracy but not RK4 or Verlet. My calculus knowledge is limited so i'd like some extra eyes to check.

They are as follows:

Code:

```
def leapfrog_integration(satellite, planet, dt): #most accurate under 5 orbits
# Update velocity by half-step
satellite.velocity += 0.5 * satellite.acceleration * dt
# Update position
satellite.position += (satellite.velocity / 1000)
# Calculate new acceleration
satellite.acceleration = satellite.acceleration_due_to_gravity(planet)
# Update velocity by another half-step
satellite.velocity += 0.5 * satellite.acceleration * dt
def euler_integration(satellite, planet, dt):
# satellite.accleration = satellite.calculate_gravity(planet)
satellite.acceleration = satellite.acceleration_due_to_gravity(planet)
satellite.velocity += satellite.acceleration * dt
satellite.position += (satellite.velocity / 1000) #convert into kilometers
def verlet_integration(satellite, planet, dt):
acc_c = (satellite.acceleration_due_to_gravity(planet) / 1000)#convert to km/s
satellite.velocity = (satellite.position - satellite.previous_position)
new_pos = 2 * satellite.position - satellite.previous_position + (acc_c * dt)
satellite.previous_position = satellite.position #km
satellite.position = new_pos #km
satellite.velocity = (satellite.position - satellite.previous_position)
def rk4_integration(satellite, planet, dt):# need to resolve the conversion to km for position. If i remove the DT from the kx_r then it's the excat same as Verlet and Euler
def get_acceleration(position, velocity):
temp_mass = PointMass(position,satellite.mass,satellite.radius,satellite.colour,(velocity), np.zeros_like(satellite.acceleration))
return temp_mass.acceleration_due_to_gravity(planet)
k1_v = dt * get_acceleration(satellite.position, (satellite.velocity))
k1_r = (satellite.velocity / 1000)
k2_v = dt * get_acceleration(satellite.position + 0.5 * k1_r, (satellite.velocity) + 0.5 * k1_v)
k2_r = (satellite.velocity + 0.5 * k1_v) / 1000
k3_v = dt * get_acceleration(satellite.position + 0.5 * k2_r, (satellite.velocity) + 0.5 * k2_v)
k3_r = (satellite.velocity + 0.5 * k2_v) / 1000
k4_v = dt * get_acceleration(satellite.position + 0.5 * k3_r, (satellite.velocity) + 0.5 * k3_v)
k4_r = (satellite.velocity + 0.5 * k3_v) / 1000
satellite.position +=(k1_r + 2*k2_r + 2*k3_r + k4_r) / 6
satellite.velocity +=(k1_v + 2*k2_v + 2*k3_v + k4_v) / 6
```

To give you the class for the pointmasses:

Code:

```
class PointMass:
def __init__(self, position, mass, radius, colour, velocity, acceleration):
self.position = np.array(position) #in KM
self.mass = mass #in Kilograms
self.radius = radius #in meters
self.colour = colour
self.velocity = np.array(velocity) #in m/s
self.acceleration = np.array(acceleration) #in m/s per second
self.gForce = None #This is in Newtons
self.previous_position = self.position - (self.velocity / 1000) # Initialize previous position for Verlet integration
def acceleration_due_to_gravity(self,other):
dVector = self.position - other.position # distance vector from self to the other point mass in pixels(km)
distance_km = np.linalg.norm(dVector) #Compute the Euclidean distance in km
distance_m = distance_km * 1000 + other.radius #the distance including the radius to the centre in meters
unit_vector = (dVector / distance_km) #the unit vector for the direction of the force
acceleration_magnitude = -Constants.mu_earth / distance_m**2
return acceleration_magnitude * (unit_vector * 1000) #Return the acceleration vector by multiplying the magnitude with the unit vector(converted to meters)
#the returned acceleration vector is in m/s
```

With the intial conditions being:

Code:

```
planet = PointMass(
position=[600.0,400.0,0.0],
mass=5.9722e24,
radius=6.371e6,
velocity=[0.0,0.0,0.0], #in m/s
acceleration=[0.0,0.0,0.0], #in m/s^2
colour=[125,100,100]
)
satellite = PointMass(
position=[280.0,400.0,0.0],
mass=100,
radius=1,
velocity=[0.0,7718.0,0.0], #need an intial velocity or else it'll jsut fall to the central mass
acceleration=[1.0,0.0,0.0],#added an non-zero acceleration jsut to make sure there's no issues with the integrations.
colour=[255,255,255]
)
```

The code is a bit messy since I'm constantly trying new things and haven't really formatted it.

Now I'm not 100% the problem lays with the integrations but I'd thought i would start there and go form that. It could easily be my gravity force function, the delta time being applied wrong, unit conversion causing problems or how it's being placed on screen that is the problem.

Github Link

Above is the repo link just to keep things in here tidy, since it could be anything. Screen shots of tests etc.

Now it could be just expect behaviour from python with loosing accuracy with floating point numbers but I'd hit X on that as it's more likely my poor calculus knowledge.

After extensive testing, I got a DT of 0.021 as the best result for delta time. Leapfrog seems the most accurate for either low orbit counts or high orbit counts.

Euler and Verlet seem about even under 100 orbits where Euler goes a little buts and RK4 seems unstable in the fact it slowly slows down and thus the orbits get larger and larger.

I have the conversion from Meters to KM in different places for each integration since they wouldn't work right depending on where i put it.

I had to remove the DT from some of the integration parts since if i left them in the object would just spiral slightly and fall to the central mass.