# MATLAB Simple Elliptical Orbits with J2 perturbation and drag

For some reason I need to build a Elliptical Orbits satellite dynamic model by myself, in here I am using MATALB.
The model its fine before (I think the J2 effect part works well) I add drag effect on it, the Altitude became ascent, it should descent by time right?
Here are some equation I am using:

And the simulator result:

Drag part code:

% Calculate the drag acceleration
v_rel_mag = norm(v_rel);
rho = atmospheric_density(norm(r_icrf) - RE); % call to your atmospheric density function
a_drag = -0.5 * Cd * A / m *  rho * v_rel_mag^2 * (v_rel / v_rel_mag)

% Total acceleration
a_total = a_gravity + a_J2 + a_drag;

% Output derivative of state vector
Ydot = [vvector; a_total];
end

function rho = atmospheric_density(altitude)
% Simple exponential atmosphere model (or use a more complex model)
rho0 = 3.614*10^-14; % kg/m^3 at sea level
H = 8500; % scale height in meters
rho = rho0 * exp(-((altitude-700) / 88.667));
end


[Full code][4] https://pastecode.io/s/pcw2mvjb
Can anyone help me? thank you!
• What is the integration scheme? Matlab's builtin ode45 or similar?
– AJN
Dec 15, 2023 at 11:12
• @AJN I am using ODE23s Dec 15, 2023 at 11:32
• I think there's a nearly-identical existing question on the site which makes it clearer that the problem is that the satellite is seeming to ascend because of atmospheric drag. I believe AJN's indication that integration precision issues can cause that was the resolution to that problem; I don't have time at the moment to figure out which question it was, but you'd be well-served looking for it on this site yourself. Dec 15, 2023 at 21:31

### When coding/scripting a new project, always independently test each of your functions. Don't just run the whole thing at once.

Also, the single scale height is useless at orbital altitudes. Above about 100 km the fall-off is MUCH slower.

"rho0 = 3.614*10^-14; % kg/m^3 at sea level" seems way off. It's about 1.3 mg/cm^3 and so about 1 .3 kg/m^3. And two lines below that you don't even use the scale height H (which looks good) but instead use 88.667 and have no idea what that number is.

function rho = atmospheric_density(altitude)
% Simple exponential atmosphere model (or use a more complex model)
rho0 = 3.614*10^-14; % kg/m^3 at sea level
H = 8500; % scale height in meters
rho = rho0 * exp(-((altitude-700) / 88.667));
end


My Python:

import numpy as np
import matplotlib.pyplot as plt

def atmospheric_density(altitude):
# Simple exponential atmosphere model (or use a more complex model)
rho0 = 3.614E-14 # kg/m^3 at sea level
H = 8500 # scale height in meters
rho = rho0 * np.exp(-((altitude - 700) / 88.667))
return rho

def atmospheric_density_fixed(altitude):
rho0 = 1.3 # kg/m^3 at sea level
H = 8500 # scale height in meters
rho = rho0 * np.exp(-((altitude - 700) / H))  # DIVIDE BY H
return rho

alti = np.linspace(0, 4E+05, 1000) # sea level to 400 km

density = atmospheric_density(alti)
density_fixed  = atmospheric_density_fixed(alti)

densities = density, density_fixed
labels = 'original', 'fixed'

fig, ax = plt.subplots(1, 1)

for den, label in zip(densities, labels):
ax.plot(alti/1000, den, label=label)
ax.set_yscale('log')
ax.set_xlabel('altitude (km)')
ax.set_ylabel('density (kg/m^3)')
ax.set_ylim(1E-30, 10)
ax.legend()
plt.show()


Result: