Calculations for FESS 1 Problem Set #2, Problem #1

pset2-1.jl

#!/usr/bin/env julia
################################################################################
# pset2-1.jl 
#
# Josh Kaplan
# _jk@jhu.edu
#
# Computations for Problem Set #2, Problem #1.
################################################################################

#----------( Constants )----------#

R_E = 6378.1        # Radius of Earth   [km]
µ_E = 3.986e5       # µ of Earth        [km^3/s^2] 


#----------( Functions )----------#

"""
Computes the period in seconds given the semi-major axis. It assumes the 
semi-major axis is given in km.

Source: Space Mission Engineering. Appendix C, Eq. C-19 (p. 964)
"""
function period(a)
    return 2*π * sqrt(a^3 / µ_E)
end


"""
Computes the velocity for a given semi-major axis and radius. It assumes the
semi-major axis and radius are given km. The returned velocity is in km/s.
 
Source: Space Mission Engineering. Appendix C, Eq. C-64 (p. 967)
"""
function velocity(a, r)
    return sqrt(2*µ_E/r - µ_E/a)
end


"""
Computes the maximum angular rate in deg/s as seen from the ground, given the 
semi-major axis (in km).

Source: Space Mission Engineering. Appendix C, Eq. C-32 (p.964)
"""
function max_angular_rate(a)
    return a * 360 / (period(a) * (a - R_E))
end


"""
Computes the maximum angular rate as seen from the ground for an elliptical 
orbit, given the the velocity at apogee (km/s, radius at apogee (km), and a 
radius. Alternatively, velocity at perigee and radius at perigee can be used 
instead of those at apogee.

Source: Space Mission Engineering. Appendix C, Eq. C-82 (p.968)
"""
function max_angular_rate_ellip(va, ra, r)
    nu_dot = (360 / (2*π)) * true_anomaly_rate(va, ra, r)
    return nu_dot * r / (r - R_E)
end


"""
Computes the rate of change of true anomaly in radians per second. Assumes the 
inputs are in compatible units.

Source: Space Mission Engineering. Appendix C, Eq. C-62 (p.966)
"""
function true_anomaly_rate(va, ra, r)
    nu_dot = va*ra / r^2
    return nu_dot
end


"""
Computes the maximum eclipse time (in seconds) for a given semi-major 
axis (in km).

Source: Space Mission Engineering. Appendix C, Eq. C-37 (p.965)
"""
function max_eclipse(a)
    P = period(a)
    return P * asind(R_E / a) / 180
end


"""
Computes the maximum eclipse time (in seconds) for a given velocity at 
apogee (km/s), radius at apogee (km), and a radius (km).

Source: Space Mission Engineering. Appendix C, Eq. C-88 (p.968)
"""
function max_eclipse_ellip(va, ra, r)
    nu_dot = true_anomaly_rate(va, ra, r)
    return (2/nu_dot) * asin(R_E / r)
end


################################################################################
##                               Problem 1                                    ##
################################################################################

th = "Orbit    Period   Max Eclipse   Velocity   Max Angular Rate\n"
un = "         [min]       [min]       [km/s]        [deg/s]     \n"
sp = "-----   -------   -----------   --------   ----------------\n"
write(STDOUT, "\n")
write(STDOUT, th) # Print table headers
write(STDOUT, un) # Print units
write(STDOUT, sp) # Print separator


##############################
##  Orbit #1 - Planet Labs  ##
##############################

# Given:
a = R_E + 410           # Semi-major axis (circular)    [km]
i = 51.7                # Inclination                   [deg]
P = 92.77*60            # Period                        [s]
v = 7.663               # Velocity                      [km/s]

# Compute remaining table items
eclipse = max_eclipse(a)
angular = max_angular_rate(a)

# Print results
@printf("#1")
@printf("%12.2f", P/60)
@printf("%12.2f", eclipse / 60)
@printf("%14.3f ", v)
@printf("%14.3f\n", angular)


#############################
##  Orbit #2 - Skybox      ##
#############################
a = R_E + 600           # Semi-major axis (circular)    [km]
i = 97.6                # Inclination                   [deg]
P = 96.69*60            # Period                        [s]       
v = 7.558               # Velocity                      [km/s]

# Compute remaining table items
eclipse = max_eclipse(a)
angular = max_angular_rate(a)

# Print Results
@printf("#2")
@printf("%12.2f", P/60)
@printf("%12.2f", eclipse / 60)
@printf("%14.3f ", v)
@printf("%14.3f\n", angular)


###################################
##  Orbit #3 - Van Allen Probes  ##
###################################
rp = R_E + 600          # Radius at apogee          [km]
ra = R_E + 30000        # Radius at perigee         [km]
a = (rp + ra)/2         # Semi-major axis           [km]
i = 10                  # Inclination               [deg]
P = period(a)           # Period                    [s]       
vp = velocity(a, rp)    # Velocity at perigee       [km/s]
va = velocity(a, ra)    # Velocity at apogee        [km/s]

# Compute eclipse times
eclipse_per = max_eclipse_ellip(va, ra, rp)
eclipse_apo  = max_eclipse_ellip(va, ra, ra)

# Compute angular rates
angular_per = max_angular_rate_ellip(va, ra, rp)
angular_apo  = max_angular_rate_ellip(va, ra, ra)

# Print Results
@printf("#3")
@printf("%12.2f\n", P/60)
@printf(" per.")
@printf("  %19.2f", eclipse_per / 60)
@printf("  %12.3f", vp)
@printf("  %13.3f\n", angular_per)
@printf(" apo.")
@printf("  %19.2f", eclipse_per / 60)  
@printf("  %12.3f", va)
@printf("  %13.3f\n", angular_per)


##############################
##  Orbit #4 - GEO Telecom  ##
##############################
a = R_E + 35786         # Semi-major axis (circular)    [km]
i = 0                   # Inclination                   [deg]
P = period(a)           # Period                        [s]       
v = velocity(a, a)      # Velocity                      [km/s]

# Compute remaining table items
eclipse = max_eclipse(a)
angular = max_angular_rate(a)

@printf("#4")
@printf("%12.2f", period(a)/60)
@printf("%12.2f", eclipse / 60)
@printf("%14.3f ", v)
@printf("%14.3f\n", angular)
println("")