;Edward Burin des Roziers
;ASEN5050

;function SAT.PRO calculates the satellite's position.

function sat, t

;Earth Constants
mu_E = 3.986d5                                ;km^3/sec^2
Re = 6378.14d                                ;km

;Sat Constants
rp_Sat = 30.0d * Re                        ;km
ra_Sat = 2.56d * Re                        ;km
a_Sat = (rp_Sat + ra_Sat)/2                ;km
inc_Sat = 14.1d / !radeg                ;rad
w_Sat = 0.0                                ;rad
Omega_Sat = 0.0                                ;rad
e_Sat = 1 - rp_Sat/a_Sat
p_Sat = a_Sat*(1-e_Sat^2)                ;km

tp_Sat = 0.0                                ;sec
n_Sat = sqrt(mu_E/a_Sat^3)

;Function PMA2 used to iterate to find eccentric anomaly
ma_Sat = (tp_Sat+t)*n_Sat + 1d-10
maae_Sat = [ma_Sat,a_Sat,e_Sat]
eata_Sat = pma2(maae_Sat)
ta_Sat = eata_Sat[1]

;Calculate Satellite position and velocity
peiOmegawnu_Sat = [p_Sat,e_Sat,inc_Sat,Omega_sat,w_Sat,ta_Sat]
posvel_Sat = k2c(peiOmegawnu_Sat)
posijk_Sat = posvel_Sat[0:2]
velijk_Sat = posvel_Sat[3:5]

return, posvel_Sat

end