Orbital propagation: atmospheric drag

Question

I'm developing the numerical orbital propagator. I've already implemented the Newtonian gravity (Earth, Sun, Moon), Earth harmonics, SRP and relativistic effect.

The software is compared with GMAT with the same configuration: for 70 days the difference now is just 30m!

Now, I'm trying to implement the atmospheric drag. I've found the equation from Vallado: $a_d=-\frac{1}{2}\rho \frac{C_dA}{m_s}|v_{rel}|v_{rel}$

I'm thinking in this way: transform the frame from ECI to ECEF. The velocity in ECEF will be $v_{rel}$. Get the density using any model (Jaccia-Roberts, etc) and calculate.

Is it correct and are there any pitfalls in the implementation? The Mathlab realization is a little different.

UPDATE

I've implemented this, something is wrong: the difference with GMAT became 100km.

The code in Julia (I tested the ECI_ECEF, ECEF_ECI and LLA functions, they are correct):

function drag(dy,y,date)

    y_ecef=ECI_ECEF([y;[0,0,0]],date)[1:6]
    lla=LLA(ECEF(y_ecef[1:3]),wgs84)

    density=nrlmsise00(date,                   # Julian Day
                       lla.alt,                # Altitude [m]
                       lla.lat*pi/180,         # Latitude [rad]
                       lla.lon*pi/180,         # Longitude [rad]
                       flux_avg,               # 81 day average F10.7
                       flux_daily,             # Daily F10.7
                       kp_index).den_Total     # KP index

    acc=-0.5*C_d*(A_d/sat_mass)*density*norm(y_ecef[4:6])*y_ecef[4:6]

    dy= ECEF_ECI([y_ecef;acc],date)[7:9]
end

Answer 1

Yes, this equation is correct, with two possible minor additions.

1) Depending on the altitude of perigee and the accuracy required, you might need to consider atmospheric winds. At altitudes typical of many LEO objects like ISS, or higher, this won't matter much. But winds farther down are typically ~100 m/s and can be as much as 800 m/s. This can induce forces parallel to and perpendicular to the velocity vector.

2) You might need to consider aerodynamic lift. For most spacecraft shapes this won't matter much, but for some shapes it can generate a non-negligible force perpendicular to the velocity vector. The magnitude of the force is given by the same equation you use for drag (with the velocity magnitude substituted for the veolcity vector), multiplied by the lift-to-drag ratio L/D. However, it is really difficult to know the hypersonic L/D accurately.

I hope this helps!

Answer 2

The problem is more simple than I expected. The equation and its implementation are correct.

The problem is, the nrlmsise00 algorithm doesn't use KP index, it uses AP index. However, in GMAT the density model configuration requires KP index.