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