"""
Sasha Bakker, HW#3, 11/23/2018, Computational Physics, Shubha Tewari
WITH DRAG:
this program plots the orbit and altitude of a satellite with drag
"""
# Import packages
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
# Initialize constants
G = 6.674*10**(-11) # Graviational constant
M = 5.9722*10**(24) # Earth mass
m = 10000 # Satellite mass
R = 6.3781*10**(6) # Earth radius
h = 500*1000 # Satellite altitude
r = R + h # Distance btwn Earth & Satellite
# Calculate the orbital speed
v = np.sqrt((G*M)/r) # Derived from a = v**2/r
# Initialize variables
ri = r # initial radius
vi = v # initial velocity
bi = (2*np.pi*ri)/vi # time for 1 loop
tDay = 24*60**2 # seconds in day
nLoops = tDay/bi # loops for 1 day
print("Input %.5f to simulate altitude lost over 1 day"%(nLoops))
o = float(input("Number of orbits is"))
N = 10000 # number of bins
a = 0 # initial time
b = (2*o*np.pi*ri)/vi # final time
h = (b-a)/N # time change per bin
# air resistance values
cv = 2
p = 1.7*10**(-11)
A = 50
# Create arrays
rVals = [[0 for x in range(N+1)] for y in range(2)]
vVals = [[0 for x in range(N+1)] for y in range(2)]
rVals[0][0] = ri
vVals[1][0] = vi
# thrust calculation
def thrust(v):
return 0.5*cv*p*A*v**2
# work from a to b
def work(force,xi,yi,xf,yf):
return force*(np.sqrt((xf-xi)**2+(yf-yi)**2))
# linear fit
def f(x,p0,p1):
return p0+p1*x
# time for an amplitude loss
def lossTime(km,p0,p1):
return (500*1000-km*1000-p0)/p1
# radius midpoint
def rMid(r,v):
return r+(h/2)*v
# velocity midpoint
def vMid(r,v,x,y,vX,vY):
value = v-(h/2)*G*M*r*(1/np.sqrt(x**2+y**2)**3)
value -= (h/2)*(1/(2*m))*cv*p*A*v*(np.sqrt(vX**2+vY**2))
return value
# new velocity
def vNew(v,rMid,x,y,vMid,vX,vY):
value = v - h*G*M*rMid*(1/(np.sqrt(x**2+y**2))**3)
value -= h*(1/(2*m))*cv*p*A*vMid*(np.sqrt(vX**2+vY**2))
return value
# new radius
def rNew(r,vMid):
return r + h*vMid
for j in range(N):
x = rVals[0][j]
y = rVals[1][j]
vX = vVals[0][j]
vY = vVals[1][j]
rMidx = rMid(x,vX)
rMidy = rMid(y,vY)
vMidx = vMid(x,vVals[0][j],x,y,vX,vY)
vMidy = vMid(y,vVals[1][j],x,y,vX,vY)
vX = vNew(vVals[0][j],rMidx,rMidx,rMidy,vMidx,vMidx,vMidy)
vY = vNew(vVals[1][j],rMidy,rMidx,rMidy,vMidy,vMidx,vMidy)
rX = rNew(x,vMidx)
rY = rNew(y,vMidy)
rVals[0][j+1] = rX
rVals[1][j+1] = rY
vVals[0][j+1] = vX
vVals[1][j+1] = vY
# Calculate thrust force
useV = np.sqrt(vVals[0][0]**2+vVals[1][0]**2)
thrusti = thrust(useV)
# print thrust force
print("\n\n*The thrust force to cancel the air drag is %.5f N"%(thrusti))
# Calculate work done
w = work(thrusti,rVals[0][0],rVals[1][0],rVals[0][-1],rVals[1][-1])
print("*Energy the thrust requires over the orbits is %.5f J"%(w))
# Calclate the angles for the polar plot
theta = np.arctan2(rVals[1][:],rVals[0][:])
# Calculate the radii for the polar plot
m = list()
for i in range(N+1):
m.append(np.sqrt(rVals[0][i]**2 + rVals[1][i]**2))
m = np.array(m) - R
# state altitude loss
minVal = min(m)
maxVal = max(m)
loss = m[0]-m[-1]
# plot polar
print("*The altitude loss for the orbits is %.5f meters\n"%(loss))
plt.polar(theta,m,'.')
plt.ylim(minVal,maxVal)
plt.title("Polar Satellite Trajectory")
plt.show()
# Make the altitudes plot
altitudes = m
times = np.linspace(a,b,N)
plt.figure()
plt.plot(times,altitudes[:N])
plt.title("Altitude (m) vs. Time (s)")
# Best fit the altitude plot
guesses = (500200.0,-0.01)
(p0,p1),cc = curve_fit(f,times,altitudes[:N],p0=guesses)
yfit = f(times,p0,p1)
plt.plot(times,yfit)
plt.show()
# Request input
k = input("By how many kilometers should the altitude be decreased (from 500km)?")
km = float(k)
# Calculate time (s) to lose altitude
loss = (lossTime(km,p0,p1)/60**2)/24
print("It takes %.1f days to lose %.1f km altitude."%(loss,km))
# Message
print("\nThe only thing I fixed is the altitude loss part.")
print("The thrust force should be equal in magnitude to the drag force.")