# Orbital propagation: atmospheric drag

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]
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

• You have an extra v term at the end of your equation. – djr Oct 16 '18 at 18:48
• @djr Yes, I was thinking about it just now. Now it should be ok – Leeloo Oct 16 '18 at 18:50
• I'm not sure I understand what the different elements of the vector are, but in norm(y[4:6])*y[4:6] should they be y_ECEF? – djr Oct 16 '18 at 19:35
• @djr You're absolutely right. I forgot to modify it while simplifying the code to copy here. In my code I make these transformations outside the function. – Leeloo Oct 16 '18 at 19:37
• The (updated) implementation looks correct. You might have a difficult time getting your results to match GMAT exactly though - you'll need to be certain you're using the exact same source for your space weather data. – od-guy Oct 16 '18 at 20:30

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!

• Thank you! However, I'm implementing this for attitude > 600m. And use a simple SC shape, comparing with GMAT. Yes, the equation is correct, but the implementation is not – Leeloo Oct 16 '18 at 19:34
• Please, look at the update – Leeloo Oct 16 '18 at 19:48

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.