# Chapter 2: Orbits and Launching Methods¶

## Example 2.1, Page 23¶

In [4]:
import math

#Variable Declaration

u=3.986*(10**14)         #Earth's Gravitational constant(m^3/sec^2)

#Calculation

a=((u/n**2)**(0.33333))/1000      #Radius of the orbit by kepler's 3rd law(km)
a=round(a,2)

#Result

print "The Radius of the circular orbit with 1 day period is:", a,"km"

The Radius of the circular orbit with 1 day period is: 42247.94 km


## Example 2.2, Page 29¶

In [5]:
import math

#Variable Declaration
NN=14.22296917        #Mean Motion (1/day)
u=3.986005*(10**14)   #Earth's Gravitational COnstant(m^3/sec^2)

#Calculation

a=((u/n0**2)**(0.33333))/1000   #Radius of the orbit by kepler's 3rd law(km)
a=round(a,2)
#Result

print "The Semimajor axis for given satellite parameters is:", a,"km"

The Semimajor axis for given satellite parameters is: 7193.97 km


## Example 2.3, Page 29¶

In [6]:
import math

#Variable Declaration

e=0.0011501  #Eccentricity
a=7192.3     #Semimajor axis(km)

#Calculation

ha=round(ra-R,2)      #Apogee height(km)
hp=round(rp-R,2)      #Perigee height(km)

#Result

print "The Apogee height for given orbital parameters is:", ha,"km"
print "The Apogee height for given orbital parameters is:", hp,"km"

The Apogee height for given orbital parameters is: 829.57 km
The Apogee height for given orbital parameters is: 813.03 km


## Example 2.4, Page 31¶

In [7]:
import math

#Variable Declaration

e=0.002  #Eccentricity
p=12     #period from perigee to perigee (hours)
K1=66063.1704  #Constant (km^2)
u=3.986005*(10**14)         #Earth's Gravitational constant(m^3/sec^2)

#Calculation

anp=((u/n**2)**(0.3333))/1000      #Radius of the orbit by kepler's 3rd law(km)
k2=(1-e**2)**1.5
# Solving for perturbed value of semimajor axis
import scipy
import scipy.optimize
def f(a):
y=(n-((u/a**3)**0.5)*(1+K1/a**2*k2))
return y
a=scipy.optimize.fsolve(f,2)
a=a/1000  #Converting a into km

#Result

print "The nonperturbed value of semimajor axis is", anp,"km"
print "The perturbed value of semimajor axis is", a,"km"

The nonperturbed value of semimajor axis is 26564.7679852 km
The perturbed value of semimajor axis is [ 26610.22410209] km


## Example 2.5, Page 34¶

In [8]:
import math

#Variable Declaration

i=98.6328      #Angle(degrees)
e=0.0011501    #eccentricity
n=14.23304826  #Mean Motion(1/day)
a=7192.3       #Semimajor axis(km)
K1=66063.1704  #Known constant(km^2)

#Calculation

n0=(2*180*n)  #Mean Motion (deg/sec)
K=(n0*K1)/((a**2)*((1-e**2)**2))  #Constant (deg/day)
w=-K*math.cos(i*3.142/180)    #Rate of regression of nodes(deg/day)
W=K*(2-2.5*(math.sin(i*3.142/180))**2)   #Rate of rotation of line of apsides(deg/day)
w=round(w,3)
W=round(W,3)
#Results

print "The rate of regression of nodes is:", w,"deg/day"
print "The rate of rotation of line of apsides is:", W,"deg/day"

The rate of regression of nodes is: 0.984 deg/day
The rate of rotation of line of apsides is: -2.902 deg/day


## Example 2.6, Page 34¶

In [9]:
#Variable Declaration

w=0.982  #rate of regression of nodes from Example 2.5(deg/day)
W=-2.903 #rate of rotation of line of apsides from Example 2.5)deg/day)
n=14.23304826  #Mean Motion(1/day)
W0=113.5534   #Argument of perigee(deg)
w0=251.5324   #Right ascension of the ascending node(deg)

#Calculation

PA=1/n   #Period
w=round(w0+w*PA,2)  #New value of rate of regression of nodes(deg)
W=round(W0+W*PA,2)  #New Value of rate of rotation of line of apsides(deg)

#Result

print "New value of rate of regression of nodes is:", w,"deg"
print "New value of rate of rotation of line of apsides is:", W,"deg"

New value of rate of regression of nodes is: 251.6 deg
New value of rate of rotation of line of apsides is: 113.35 deg


## Example 2.7, Page 38¶

In [10]:
#Calculation

ndays=400*365    #Nominal number of days in 400years
nleapyrs=400/4   #Nominal number of leap years
gregoriandays=ndays+nleapyrs-3  #number of days in 400 years of Gregorian calendar
gregavg=round(gregoriandays/float(400),2)    #number of days in 400 years of Gregorian calendar

#Result
print gregoriandays
print "The average length of the civil year in gregorian calender is:", gregavg,"days"

146097
The average length of the civil year in gregorian calender is: 365.24 days


## Example 2.8, Page 38¶

In [11]:
#Calculation and Results

if 1987%4==0:
print "1987 is a leap year"
else:
print "1987 is not a leap year"

if 1988%4==0:
print "1988 is a leap year"
else:
print "1988 is not a leap year"

if 2000%400==0:
print "2000 is a leap year"
else:
print "2000 is not a leap year"

if 2100%400==0:
print "2100 is a leap year"
else:
print "2100 is not a leap year"

1987 is not a leap year
1988 is a leap year
2000 is a leap year
2100 is not a leap year


## Example 2.9, Page 38¶

In [12]:
#Calculation
days=324   #Number of days
hours=math.floor(24*0.95616765)  # Number of hours
decimalfraction1=24*0.95616765-hours
minutes=math.floor(60*decimalfraction1)  # Number of minutes
decimalfraction2=60*decimalfraction1-minutes
seconds=round(60*decimalfraction2,2)  # Number of seconds

#Result

print decimalfraction1,decimalfraction2
print "An Epoch day has", days,"days",hours,"hours",minutes,"minutes",seconds,"seconds"

0.9480236 0.881416
An Epoch day has 324 days 22.0 hours 56.0 minutes 52.88 seconds


## Example 2.10, Page 39¶

In [13]:
import math

#Variable Declaration

y=2000      #year
mon=12      #month
dy=18       #day
hours=13    #hours of the day
minutes=0   #Minutes of the day
seconds=0   #Seconds of the day

#Calculation
d=dy+(hours/24)+(minutes/(24*60))+seconds  #Days in December
if mon<=2:
y=y-1
mon=mon+12
else:
y=y
mon=mon

A=math.floor(y/100)      #Converting years to days
B=2-A+math.floor(A/4)    #Converting years to days
C=math.floor(365.25*y)     #rounding the days
D=math.floor(30.6001*(mon+1)) #Converting months to days
JD=B+C+D+d+1720994.5     #Adding reeference to number of days

#Result

print "The Julian day of given day is:", JD,"Days"

The Julian day of given day is: 2451896.5 Days


## Example 2.11, Page 41¶

In [14]:
#Variable Declaration

JDref=2415020   #Reference Julian days
JC=36525
JD=2451897.0417 #Julian days with reference from Example 2.10

#Calculation

T=round((JD-JDref)/JC,4)  #Time in julian Centuries

#Result

print "The time for given date is:", T,"Julian Centuries"

The time for given date is: 1.0096 Julian Centuries


## Example 2.12, Page 43¶

In [15]:
#Variable Declaration

n=14.23304826  #Mean Motion (rev/day)
M0=246.6853    #Mean Anomaly (degrees)
t0=223.79688452 #Time of anomaly

#Calculation

T=round(t0-(M0/(n*360)),3)  #Time of perigee passage

#Result

print "The time of perigee passage for NASA elements is:", T,"days"

The time of perigee passage for NASA elements is: 223.749 days


## Example 2.13, Page 44¶

In [16]:
import math
#Variable Declaration

M=205  #Mean anomaly(degrees)
e=0.0025  #Eccentricity
E=math.pi #Initial guess value for eccentric anomaly

#Calculation

import scipy
import scipy.optimize
def f(E):
y=M-E+e*math.sin(E)
return y
E=scipy.optimize.fsolve(f,3.142)
E=round(E,3)

print "The Eccentric anomaly is:",E,"degrees"

The Eccentric anomaly is: 204.998 degrees


## Example 2.14, Page 45¶

In [17]:
import math
#Variable Declaration

e=9.5981*10**-3   #Eccentricity
a=7194.9   #Semimajor axis(km)

#Calculation

v=v*180/math.pi  #True anomaly(degrees)
r=a*(1-e**2)/(1+e*math.cos(v))  #Magnitude of radius vector after 5s(km)
v=round(v,2)
r=round(r,2)

#Results

print "The true anomaly is:",v,"degrees"
print "The magnitude of radius vector 5s after epoch is:",r,"km"

The true anomaly is: 204.79 degrees
The magnitude of radius vector 5s after epoch is: 7252.02 km


## Example 2.15, Page 46¶

In [18]:
import math
#Variable Declaration

v=204.81   #True anomaly(degrees) from Example 2.14
r=7257     #Magnitude of radius vector(km) from Example 2.14

#Calculation

rP=r*math.cos(v*math.pi/180)  #P coordinate of radius vector(km)
rQ=r*math.sin(v*math.pi/180)  #Q coordinate of radius vector(km)
rP=round(rP,2)
rQ=round(rQ,2)
#Result

print "r in the perifocal coordinate system is", rP,"Pkm", rQ,"Qkm"

r in the perifocal coordinate system is -6587.21 Pkm -3045.11 Qkm


## Exanple 2.17, Page 50¶

In [19]:
import math
#Variable Declaration

T=1.009638376  #Time in Julian centuries from Example 2.11
UT=13  #Universal time(hours)

#Calculation

UT=2*math.pi*UT/24  #Universal time converted to fraction of earth rotation (radians)

GST=GST+UT

GST=(math.fmod(GST,2*math.pi))*180/math.pi#using fmod multiplr revolutions are removed (degrees)
GST=round(GST,2)

#Result

print "The GST for given date and time is ", GST,"degrees"

The GST for given date and time is  287.18 degrees


## Example 2.18, Page 51¶

In [20]:
import math
#Variable Declararion

WL=-89.26  #Expressing the longitude in degrees west
GST=282.449  #GST from Example 2.17 (degrees)

#Calculation

EL=2*math.pi+WL  #Longitude in degrees East
LST=(math.fmod(LST,2*math.pi))*180/math.pi   #fmod removes multiple revolutions(Degrees)
LST=round(LST,2)

#Results

print "LST for Thunder Bay on given day is:" ,LST,"Degrees"

LST for Thunder Bay on given day is: 199.47 Degrees


## Example 2.19, Page 52¶

In [21]:
import math
#Variable Declaration

LST=167.475 #LST(degrees)
LE=48.42    #Latitude at thunder bay(degrees)
H=200       #Height above sea level(metres)
aE=6378.1414 #Semimajor axis(km)
eE=0.08182   #Eccentricity

#Calculation

l=(aE/math.sqrt(1-eE**2*math.sin(LE*3.142/180)**2)+H/1000)*math.cos(LE*3.142/180)

z=((aE*(1-eE**2))/math.sqrt(1-eE**2*math.sin(LE*3.142/180)**2)+H/1000)*math.sin(LE*3.142/180)

RI=round(l*math.cos(LST*3.142/180),2) #I component of radius vector at thunder bay(km)

RJ=round(l*math.sin(LST*3.142/180),2) #J component of radius vector at thunder bay(km)

RK=round(z,2)      #Z component of radius vector at thunder bay(km)

R=math.sqrt(RI**2+RJ**2+RK**2)
R=round(R,2)

#Results

print "The Radius vector components are" ,RI,"ikm+",RJ,"jkm+",RK,"kkm"
print "The Magnitude of radius component is" ,R,"km"

The Radius vector components are -4139.81 ikm+ 918.02 jkm+ 4748.46 kkm
The Magnitude of radius component is 6366.21 km


## Example 2.20, Page 55¶

In [22]:
#Variable Declaration

PI=-1280  #I component of range vector for a satellite(km)
PJ=-1278  #J component of range vector for a satellite(km)
PK=66     #K component of range vector for a satellite(km)
GST=240   #GST(degrees)
LE=48.42  #Latitude(Degrees)
PE=-89.26 #Longitude(Degrees)
H=200     #Height above mean sea level(metres)
aE=6378.1414 #Semimajor axis(km)
eE=0.08182   #Eccentricity

import math
#Calculation

l=(aE/math.sqrt(1-eE**2*math.sin(LE*3.142/180)**2)+H/1000)*math.cos(LE*3.142/180)

z=((aE*(1-eE**2))/math.sqrt(1-eE**2*math.sin(LE*3.142/180)**2)+H/1000)*math.sin(LE*3.142/180)

SE=(math.atan(z/l))*180/3.142 #Geocentric latitude angle (degrees)
LST=240+PE

from numpy import matrix
from numpy import linalg
a=math.sin(SE*3.142/180)*math.cos(LST*3.142/180)
b=math.sin(SE*3.142/180)*math.sin(LST*3.142/180)
c=-math.cos(SE*3.142/180)
d=-math.sin(LST*3.142/180)
e=math.cos(LST*3.142/180)
f=0
g=math.cos(SE*3.142/180)*math.cos(LST*3.142/180)
h=math.cos(SE*3.142/180)*math.sin(LST*3.142/180)
i=math.sin(SE*3.142/180)

D = matrix( [[a,b,c],[d,e,f],[g,h,i]])

P=matrix( [[PI],[PJ],[PK]] )

R=D*P  #Components of range of earth station
Ro=math.sqrt(R[0,0]**2+R[1,0]**2+R[2,0]**2) #Magnitude of range of earth station(km)
Ro=round(Ro,2)
El=math.asin(R[2,0]/Ro) #Antenna elevation angle for the earth station(radians)
El=round(El*180/3.142,2) #Converting El to degrees
alpha=(math.atan(R[1,0]/R[2,0]))*180/3.142
if (R[0,0]<0 and R[1,0]>0):
Aza=alpha
else:
Aza=0
if (R[0,0]>0 and R[1,0]>0):
Azb=180-alpha
else:
Azb=0
if (R[0,0]>0 and R[1,0]<0):
Azc=180+alpha
else:
Azc=0
if (R[0,0]<0 and R[1,0]<0):
Azd=360-alpha
else:
Azd=0
Az=round(Aza+Azb+Azc+Azd,2)  #Azimuth angle (degrees)

print "The magnitude of range of earth station is" ,Ro,"km"
print "The antenna elevation angle for the earth station are",El,"degrees"
print "The Azimuth angle for the earth station is",Az,"degrees"

The magnitude of range of earth station is 1809.98 km
The antenna elevation angle for the earth station are 12.03 degrees
The Azimuth angle for the earth station is 102.24 degrees


## Example 2.21, Page 58¶

In [23]:
#Variable Declaration

rI=-4685.3  #I component of radius vector from Example 2.16(km)
rJ=5047.7   #J component of radius vector from Example 2.16(km)
rK=-3289.1  #K component of radius vector from Example 2.16(km)
aE=6378.1414 #Semimajor axis (km)
eE=0.08182   #Eccentricity

#Calculation

r=math.sqrt(rI**2+rJ**2+rK**2)
c=r-aE   #Guess value for height(km)

from scipy.optimize import fsolve

def equations(p):
L,h,LST = p
return (rI-((aE/math.sqrt(1-eE**2*math.sin(L)**2))+h)*math.cos(L)*math.cos(LST),rJ-((aE/math.sqrt(1-eE**2*math.sin(L)**2))+h)*math.cos(L)*math.sin(LST),rK-((aE*(1-eE**2)/math.sqrt(1-eE**2*math.sin(L)**2))+h)*math.sin(L))

L,h,LST =  fsolve(equations, (b,c,a))
L=round(L*180/3.142,2)  #Converting L into degrees
h=round(h)
LST=round(LST*180/3.142,1) #Converting LST into degrees

print "The latitude of subsatellite is",L,"degrees"
print "The height of subsatellite is",h,"km"
print "The LST of subsatellite is",LST,"degrees"

The latitude of subsatellite is -25.65 degrees
The height of subsatellite is 1258.0 km
The LST of subsatellite is 132.9 degrees