While propagating the satellite motion faced the strange effect. Without the atmosphere results (x,y,z coordinates in meters and vx,vy,vz velocities in meters per second) are [3.30971e5, -6.55755e6, 2.57717e6, -1261.46, -2762.24, -6871.88] but including the atmosphere drag effect i obtain [-2.90003e5, -7.04404e6, -8.97772e5, -1268.31, 996.889, -7305.6]. Notice that difference between satelite positions is 300+ km. Can the atmospheric drag effect be so big or it is the implementation and modelling issues? The formula for atmospheric drag was taken from GMAT documentation:$$a=\frac{1}{2}\rho v^2_{rel}\frac{C_dh}{m_s}\hat{v_{rel}}$$ and the barometric formula(from Wikipedia) to calculate the density: $$\rho=\rho_d\exp\left[ \frac{-g_0M(h-h_0)}{ RT_b}\right ].$$ The $h_0,T_b$ and $\rho_d$ constants were taken from http://www.braeunig.us/space/atmos.htm. Here is the code:

using DifferentialEquations
jd  = Dates.datetime2julian(DateTime(2100,01,01,0,0,0)) * 86400
jd2 = Dates.datetime2julian(DateTime(2100,01,11,0,0,0)) * 86400
y = [-976.3107644649057E+03, -4835.627052558522E+03, -5031.728586125443E+03, -0.7944031487816871E+03, 5.474532271429767E+03, -5.094496750907486E+03]
GMe = BigFloat(398600.4415E+9)
req = BigFloat(6378136.3) #m
mass = 850
A = 15.0 #m^2
rho_b  = 6.65E-14 #kg/m^3,
g   = 9.80665 #m/s^2,
M = 0.0289644 #kg/mol,
h_0 = 650000 #m,
R = 8.3144598 #N·m/(mol·K),
T_b = 1011.5365 #K
omega = [0.0, 0.0, 7.292115078468551e-5]
v_wind = [0.0,0.0,0.0]
C_d = 0.47#from GMAT

#× vector product
function atmospheric_drag(y)
    v = [y[4],y[5],y[6]]
    r = [y[1],y[2],y[3]]
    h = sqrt(y[1]^2+y[2]^2+y[3]^2) - req
    rho = rho_b*exp((-g*M*(h-h_0))/(R*T_b))
    #println("rho = $rho")
    v_rel = v - omega×r + v_wind
    v_rel2 = v_rel[1]^2+v_rel[2]^2+v_rel[3]^2
    a = -0.5*rho*(C_d*A/mass)*v_rel2*v_rel


function f2(dy,y,p,t)

    atm_dr = atmospheric_drag(y)
    dy[1] = y[4]
    dy[2] = y[5]
    dy[3] = y[6]
    dy[4] = - GMe*y[1]/re3 - atm_dr[1]
    dy[5] = - GMe*y[2]/re3 - atm_dr[2]
    dy[6] = - GMe*y[3]/re3 - atm_dr[3]

prob = ODEProblem(f2,y,(jd,jd2))
solution = solve(prob,Vern9(),abstol=1e-13,reltol=1e-13)
sol = solution[end]
#order 0 degree 0 no Sun, no Moon
GMATresult = [330.97011851316,-6557.5472439174,2577.1624246640,-1.2614569607713,-2.7622388770690,-6.8718847034529]

a = 0
for i = (1:6)
  b = sol[i] - GMATresult[i]*1000
  println("$b m")
	if i < 4
	 a = a + b^2
println("Distance error = $(sqrt(a)) m")

I've noticed that even little values like 1e-6 being added to f2 affect significantly (70+ km) on the result. The values of the acceleration provided by atmospheric_drag function are of 1e-5 order. Maybe that is the reason?

  • $\begingroup$ What is the time interval (or number of orbits) used when propagating the satellite's motion? Does the difference between the calculation with and without atmospheric effects increase over time? At what height do you start? $\endgroup$
    – Uwe
    Aug 2, 2018 at 8:58
  • $\begingroup$ @Uwe, the time interval is 10 days.The difference acts strange: it gave me 3.563345952606408e6 meters on the interval of 10 days, on 20 - 3.1575382220368516e6, 30 - 6.027211848129026e6, 60 days - 350452.0664774459, 90 - 6.883028246992346e6. I can not see any sense or law in it. The height i start from is 668477.373579517 meters( in program as value of h variable). $\endgroup$ Aug 2, 2018 at 10:17
  • 1
    $\begingroup$ Just write 668 km instead of 668477.373579517 meters. We do not need nanometers here for the height of a satellite. $\endgroup$
    – Uwe
    Aug 2, 2018 at 10:23

1 Answer 1


Back of the spherical cow's envelope:

It's moving at $v$ = 7500 m/s. After ten days it's off by 3E+05 meters.

$$x = \frac{1}{2} a t^2$$

gives an acceleration of 8E-08 m/s^2. With a mass $m$ = 850 kg that's a force of 7E-04 Newtons.

Using a simple drag equation for $A$ = 10 m^2 and $C_D$ = 1, plugging that force into

$$F_D = \frac{1}{2} \rho v^2 C_D A$$

and solving for $\rho$ gives 2.5E-12 kg/m^3.

Looking that up in this answer or the reference there gives an altitude of only 450 km, which means it's ballparkish reasonable.

However, I don't like your simple scale height-ish atmosphere expression for density versus height. The exponential factor $ R T / (g M) $ is about 30 km with a temperature of 1000 Kelvin and I'm not sure why that works even as well as it does because you are using an $h_0$ of sea level. With such a low slope, I think you are meant to use a $\rho_d$ and $h_0$ that's already very high in the atmosphere.

Oh! If you are just using that from Wikipedia, then it's just plain wrong. You'll need to find a model for drag that's more appropriate for space!

So instead you could of course just interpolate density from your own link http://www.braeunig.us/space/atmos.htm or even just use their scale height values at say 600 km:

altitude:      600       km
scale height:  74.8      km
rho_0          9.89E-14  kg/m^2

and use $$ \rho = 9.89 \times 10^{-14} \exp\left(\frac{r - (6378.137 + 600) km}{74.8 km}\right). $$

That should give you a roughly factor of 10 smaller difference in drag vs no-drag, maybe a 10 to 50 km instead of 300 km.

In a 0th order, spherical cowish kind of way, you could even use an interpolator on the data I've shown there plus the drag equation plus $a = F/m$ to get a better acceleration than you're using now at least.

Good luck!

enter image description here

below: ugly sketch from here to illustrate that you can't use a single, simple scale height formula from sea level to LEO.

enter image description here

  • $\begingroup$ but in the first formula $x=\frac{1}{2}at^2$ you don't take the velocity into consideration. Is your formula correct? I tried both of your advices, but none of them gave the better answer. I think, the problem is that the acceleration is too high. $\endgroup$ Aug 3, 2018 at 7:06
  • $\begingroup$ @AlexJohnson It's a back of the envelope approximation. Acceleration due to a very small force causes a lag between the two positions. Of course they are both moving with almost the exact same speed, but since they start at the same spot, only the acceleration causes the difference in path length. However, I don't know what you think "better answer" is. It's a hard problem, you should compare to a known-to-be-correct propagator, not just judging what the "right answer" is by eye. $\endgroup$
    – uhoh
    Aug 3, 2018 at 7:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.