"""
Satellite life expectancy simulation
@jap 2008 having fun with Quark!
"""
from math import *
from earth import *
from msis import msise90 as rho
m = 13000. # starting mass of satellite
d = 3.0 # diameter of satellite (m)
S = pi*(d/2)**2 # cross-section of satellite (m^2)
Ca = 0.8 # aerodinamic coeficient for airdrag formula (dimensionless)
alpha = m/(Ca*S) # ballistic coefficient of satellite
sday = 24*60*60 #seconds per day
syear = sday*365 #second per year
f = open('fall.dat','w')
def T(h):
"Period of earth satellite at height h"
r = RT + h
return 2*pi*sqrt(r**3/(G*MT))
def v(h):
"Velocity of satellite at height h"
r = RT+h
return 2*pi*r/T(h)
def Ec(h):
"Kinetic energy of satellite at height h"
return 0.5*m*v(h)**2
def Ep(h):
"Potential energy of satellite at height h"
r = RT + h
return -G*MT*m/r
def Em(h):
"Mechanical energy of satellite at height h"
return Ec(h) + Ep(h)
def height(E):
return -G*MT*m/(2.*E)-RT
if __name__ == "__main__":
h = input("Initial height of satellite (m): ")
t = 0
E = Em(h)
beta = pi*G*MT*m/alpha
n = 0 # number of revolutions
while h > 0:
print >>f,t/syear,h/1E3,v(h)*3.6,Em(h) # t in years, h in km, v in km/h
t += T(h)
E -= beta*rho(h)
h = height(E)
n += 1
days, years = t/sday, t/syear
print "Satellite fate:"
print "Expected number of revolutions before falling on earth", n
print "Lifetime: %.2f days or %.5f years" % (days,years)
f.close()