GTO deficit delta-v calculator

gto.py

#!/usr/bin/env python 
import sys
import os
import math

RE = 6371                   # radius earth
MU = 3.986005 * 10**14      # G*M for earth
ALT_GEO = 35786             # GEO altitude
GEO_V = 3075                # speed at GEO


#
# Based on LouScheffer's algorithm @
# https://forum.nasaspaceflight.com/index.php?topic=36954.0
#

def semi_major_axis(perigee, apogee):
    """
    return semi-major axis in meters from a given perigee and apogee
    """
    return ((perigee + apogee)/2.0 + (RE * 1000))

def orbital_speed(distance, sma):
    """
    returns the orbital speed in m/s of an object at distance meters
    in an orbit with semi-major axis sma 
    """
    distance = distance + (RE * 1000)
    return math.sqrt(MU * ((2.0/distance) - (1.0/sma)))

def orbital_period(sma):
    """
    return the orbital period from the given semi-major axis.
    the elements returned are (revolutions_per_day, revolution_hours and
    revolution_minutes)
    """
    rpd = 8681663.653/((sma/1000) ** (3.0/2.0))    # revolutions per day
    rev_hours = (24 * (1/rpd))
    rev_min = round(60 * (rev_hours - math.floor(rev_hours)), 0)
    rev_hours = math.floor(rev_hours)
    return (rpd, rev_hours, rev_min)

def gto_circularization(perigee, apogee, inclination):
    sub_sync_cost = 0
    # calculate initial GTO insertion orbit
    sma_initial = semi_major_axis(perigee, apogee)
    speed_perigee_initial = orbital_speed(perigee, sma_initial)
    speed_apogee_initial = orbital_speed(apogee, sma_initial)
    rev_per_day, rev_hours, rev_mins = orbital_period(sma_initial) 
    
    print "Initial GTO orbit: %.2f x %.2f x %.2f" % (perigee/1000.0, apogee/1000.0, inclination)
    print "  speed at perigee: %.3f m/s" % (speed_perigee_initial) 
    print "  speed at apogee: %.3f m/s" % (speed_apogee_initial) 
    print "  revs/day: %.3f" % (rev_per_day)
    print "  days/rev: %.3f" % (1/rev_per_day)
    print "  rev time: %02uh%02um" % (rev_hours, rev_mins)

    if apogee < ALT_GEO * 1000:
        # if sub-sync, raise apogee at first perigee
        print "sub-sync injection, raising apogee to GEO at first perigee"
        apogee_raised = ALT_GEO * 1000
        sma_raised = semi_major_axis(perigee, apogee_raised)
        speed_from_raised_apogee = orbital_speed(perigee, sma_raised)
        print "  speed at perigee after raising apogee: %.3f m/s" % (speed_from_raised_apogee)
        sub_sync_cost = speed_from_raised_apogee - speed_perigee_initial
        print "  sub-sync apogee raise costs %.3f m/s" % (sub_sync_cost)
        # set new apogee to GEO altitude
        apogee = apogee_raised 
        # set new initial speeds for calculations below
        speed_perigee_initial = speed_from_raised_apogee
        speed_apogee_initial = orbital_speed(apogee, sma_raised)
        print "  orbit now %.2f x %.2f x %.2f" % (perigee/1000.0, apogee/1000.0, inclination)

    # now raise the perigee to GEO altitude
    perigee_raised = ALT_GEO * 1000
    sma_raised = semi_major_axis(perigee_raised, apogee)
    speed_apogee_raised = orbital_speed(apogee, sma_raised)
    speed_perigee_raised = orbital_speed(perigee_raised, sma_raised)
    print "after perigee raise:"
    print "  speed at perigee: %.3f m/s" % (speed_perigee_raised)
    print "  speed at apogee: %.3f m/s" % (speed_apogee_raised)


    # velocity as right angles to equator
    speed_cross = speed_apogee_initial * math.sin(inclination/180 * math.pi)
    # velocity along the equator
    print "delta-v requirements:"
    speed_along = speed_apogee_initial * math.cos(inclination/180 * math.pi)
    print "  cross speed at apogee: %.3f m/s" % (speed_cross)
    print "  along track: %.3f m/s" % (speed_along)
    # how much extra we need to raise perigee
    need_along = speed_apogee_raised - speed_along
    # removing cross component and raising perigee
    dv_top = math.sqrt(speed_cross**2 + need_along**2)
    dv_bot = speed_perigee_raised - GEO_V
    # use absolute value here in order to ensure that
    # sub synchronous injections are properly handled
    dv_total = dv_top + math.fabs(dv_bot) + sub_sync_cost
    print "  dv at top: %.3f" % (dv_top)
    print "  dv at bottom: %.3f" % (dv_bot)
    print "  dv total: %.3f" % (dv_total)

     
def usage(p):
    print "usage: %s <perigee km> <apogee km> <inclination degrees>" % p
    sys.exit(1)

if __name__ == "__main__":
    if len(sys.argv) != 4:
        usage(sys.argv[0])
    perigee = float(sys.argv[1]) * 1000 # convert to meters
    apogee = float(sys.argv[2]) * 1000 # convert to meters
    inclination = float(sys.argv[3])
    gto_circularization(perigee, apogee, inclination)