import math
#Variable Declaration
u=3.986*(10**14) #Earth's Gravitational constant(m^3/sec^2)
#Calculation
n=(2*3.14)/(24*60*60) #Mean Motion(rad/sec)
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"
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
n0=(NN*2*3.142)/(24*60*60) #Mean Motion(rad/sec)
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"
import math
#Variable Declaration
R=6371 #Mean Earth's radius(km)
e=0.0011501 #Eccentricity
a=7192.3 #Semimajor axis(km)
#Calculation
ra=a*(1+e) #Radius Vector at apogee(km)
rp=a*(1-e) #Radius Vector at perigee(km)
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"
import math
#Variable Declaration
aE=6378.141 #Earth's equitorial radius(km)
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
n=(2*math.pi)/(12*60*60) #Mean Motion(rad/sec)
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"
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"
#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"
#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"
#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"
#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"
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"
#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"
#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"
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"
import math
#Variable Declaration
n=14.2171401*2*math.pi/86400 #Mean motion (rad/sec)
M=204.9779+0.001*180*5/math.pi #Mean anomaly(rad)
e=9.5981*10**-3 #Eccentricity
a=7194.9 #Semimajor axis(km)
#Calculation
v=(M*math.pi/180)+2*e*math.sin(M*math.pi/180)+(5*e**2*math.sin(2*M*math.pi)/(4*180)) #True Anomaly (radians)
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"
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"
import math
#Variable Declaration
T=1.009638376 #Time in Julian centuries from Example 2.11
UT=13 #Universal time(hours)
#Calculation
GST=(99.6910+36000.7689*T+0.004*T**2)*3.142/180 #GST(radians)
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"
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=(GST+EL)*math.pi/180 #LST(radians)
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"
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"
#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"
#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)
a=math.pi #Guess value for LST(radians)
b=math.atan(rK/rI) #Guess Value for latitude(radians)
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"