MATLAB Simple Elliptical Orbits with J2 perturbation and drag

Question

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.61410^-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!

Answer 1

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.

• see this answer to What orbit would a space station need to stay in orbit for N years?
• Why does Earth's atmospheric density have a big "knee" around 100 km? Is there a good analytical approximation?
• Would the >100 km "knee" in Earth's atmosphere still be where MFP exceeded scale height if it were pure Ne or Ar?
• How many solar system bodies have "knees" in their atmospheres?

"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.

Your Matlab:

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: