"""
Sasha Bakker, HW#3, 11/9/2018, Computational Physics, Shubha Tewari
WITHOUT DRAG:
this program plots the orbit and altitude of a satellite without drag
"""

#Import packages
import numpy as np
import matplotlib.pyplot as plt

# 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
sA = 500*1000 # Satellite altitude
r = R + sA # Distance btwn Earth & Satellite
v = np.sqrt((G*M)/r) # orbital speed
ri = r # initial radius
vi = v # initial velocity

# Initialize variables
o = float(input("Number of orbits is"))
print("Number of orbits is",o)
N = 100000 # number of bins
a = 0 # initial time
b = (2*o*np.pi*ri)/vi # final time
h = (b-a)/N # time change per bin
m = list() # R's of the polar plot
alt = list() # altitudes

# Create arrays
rVals = [[0 for x in range(N+1)] for y in range(2)] #radii
vVals = [[0 for x in range(N+1)] for y in range(2)] #velocities
rVals[0][0] = ri
vVals[1][0] = vi

# define a function to calculate the radius midpoint
# this value inputs r = r of n, v = v of n
def rMid(r,v):
    return r+(h/2)*v

# define a function to calculate the velocity midpoint
# this value inputs r = r of n, v = v of n
def vMid(r,v,x,y):
    return v-(h/2)*G*M*r*(1/np.sqrt(x**2+y**2)**3)

# calculates new velocity: v of (n+1)
def vNew(v,rMid,x,y):
    return v - h*G*M*rMid*(1/(np.sqrt(x**2+y**2))**3)

# calculates new radius: r of (n+1)
def rNew(r,vMid):
    return r + h*vMid


# Fill the rVals and vVals arrays
for j in range(N):
    
    # Define the current r of n values for x,y
    x = rVals[0][j]
    y = rVals[1][j]
    
    # Calculate the radius midpoints for x,y
    rMidx = rMid(x,vVals[0][j])
    rMidy = rMid(y,vVals[1][j])
    
    # Calculate the velocity midpoints for x,y
    vMidx = vMid(x,vVals[0][j],x,y)
    vMidy = vMid(y,vVals[1][j],x,y)
    
    # Calculate the velocity, v of n+1 for x,y
    vX = vNew(vVals[0][j],rMidx,rMidx,rMidy)
    vY = vNew(vVals[1][j],rMidy,rMidx,rMidy) 
    
    # Calculate the radius, r of n+1 for x,y
    rX = rNew(x,vMidx)
    rY = rNew(y,vMidy)
    
    # Save the new radius and velocity values for x,y
    rVals[0][j+1] = rX
    rVals[1][j+1] = rY
    vVals[0][j+1] = vX
    vVals[1][j+1] = vY
    
# Calclate the angles for the polar plot
theta = np.arctan2(rVals[1][:],rVals[0][:])

# Calculate the radii for the polar plot
for i in range(N+1):
    val = np.sqrt(rVals[0][i]**2 + rVals[1][i]**2)
    m.append(val) # polar plot radius
    alt.append(val-R) # altitudes

# Plot the polar
plt.polar(theta,m,'.')
plt.title("Polar Satellite Trajectory")
plt.show()

    
# Make the altitude plot
altitudes = np.array(alt)
times = np.linspace(a,b,N) # create range of times 
plt.figure()
plt.plot(times,altitudes[0:N])
plt.title("Altitude (m) vs. Time (s)")
plt.show()

Number of orbits is 10.0

"""
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.")

Input 15.21919 to simulate altitude lost over 1 day

*The thrust force to cancel the air drag is 0.04926 N
*Energy the thrust requires over the orbits is 431153.25600 J
*The altitude loss for the orbits is 497.24198 meters

It takes 14.4 days to lose 10.0 km altitude.

The only thing I fixed is the altitude loss part.
The thrust force should be equal in magnitude to the drag force.