RocketLab-Electron-AscentProfile

#!/usr/bin/env python3

"""
* Planned mission: burn time not exact, but good start
* http://www.spacelaunchreport.com/electron.html: rocket parameters (wikiepdia is outdated/wrong)
* Live webcast: audio info + estimated info from the video
* Lots of fidling with Cd and Pitch over maneuver
"""

import numpy as np
from matplotlib import pyplot as plt

#
# Simple math helper
#

class Vector(object):
    def __init__(self, x, y):
        self.x = x
        self.y = y

    @staticmethod
    def init_from_magangle(magnitude, angle):  # angle: deg
        angle *= np.pi/180
        return Vector(magnitude*np.cos(angle), magnitude*np.sin(angle))

    def magnitude(self):
        return np.sqrt(self.x**2 + self.y**2)

    def __add__(self, other):
        return Vector(self.x+other.x, self.y+other.y)

    def __sub__(self, other):
        return Vector(self.x-other.x, self.y-other.y)

    def __repr__(self):
        return 'Vector(%.2f %.2f)' % (self.x, self.y)

#
# Rocket related objects
#

m_payload = 160
m_s2_dry_mass = 250
m_s2_fuel = 2500
m_s1_dry_mass = 950
m_s1_fuel = 9250

print('Rocker mass: %.3f t' % (sum((
    m_payload,
    m_s2_dry_mass,
    m_s2_fuel,
    m_s1_dry_mass,
    m_s1_fuel,
))/1000))

class StageS1(object):
    """
    Variable thrust motor (due to varying atmospheric density)
    """
    def __init__(self):
        self.fuel = m_s1_fuel  # Kg
        self.dry_mass = m_s1_dry_mass + m_s2_fuel + m_s2_dry_mass + m_payload

    def get_thrust(self, h):
        sl_thrust = 15.65*9.81*1000  # N
        vac_thrust = 18.82*9.81*1000  # N

        rho = 1.3 * np.e ** ( -h / 7000 )
        d_thrust = (vac_thrust-sl_thrust)/1.3  # kN
        thrust = vac_thrust - d_thrust*rho
        return thrust  # N

    def get_fc(self, h):
        return m_s1_fuel/153

        #rho = 1.3 * np.e ** ( -h / 7000 )
        #d_fc = (64.59-54.50)/1.3
        #fc = 64.59 - d_fc*rho
        #return fc  # Kg/s

    def jettison(self):
        return self.fuel < 0


class StageSeparationS12(object):
    """
    0 N thrust motor to simulate a 6 sec ballistic flight between MECO and 2nd stage ignition
    """
    def __init__(self):
        self.fuel = 6/10**6
        self.dry_mass = m_s2_fuel + m_s2_dry_mass + m_payload

    def get_thrust(self, h):
        return 0

    def get_fc(self, h):
        return 1/10**6

    def jettison(self):
        return self.fuel < 0


class StageS2(object):
    def __init__(self):
        self.fuel = m_s2_fuel  # Kg
        self.dry_mass = m_s2_dry_mass + m_payload  # Kg

    def get_thrust(self, h):
        return 2.27*9.81*1000  # 22 269 N

    def get_fc(self, h):
        return m_s2_fuel/380  # Kg/s

    def jettison(self):
        return self.fuel < 0


#
# Rocket definition
# TODO: refactor this
#

# None = follow ideal angle
ANGLES = (
    [90] * 10
    + [65] * (160-10)
    # stage sep
    + [0] * 140
    # 200 km
    + [None] * 240
)
max_angle = 20

cd = 0.3

#
# Main functions
#

def fly_rocket():
    s1 = StageS1()
    s12 = StageSeparationS12()
    s2 = StageS2()

    stages = [s1, s12, s2]

    s1_dv = 9.81*303*np.log((s1.dry_mass+s1.fuel)/s1.dry_mass)
    s2_dv = 9.81*333*np.log((s2.dry_mass+s2.fuel)/s2.dry_mass)
    print('Rocket calculated dv: %d + %d = %d m/s' % (s1_dv, s2_dv, s1_dv+s2_dv))

    angle = 90

    ts = [0] # time
    accs = [Vector(0, 0)]
    vels = [Vector(0, 0)]
    poss = [Vector(0, 0)]
    dps_p = [0]  # dynamic pressure - pascal
    dps_n = [0]  # dyn pre - newton
    drag_losses = [0]
    g_losses = [0]

    _thrust = 0
    _cent = 0

    #
    # Simulate flight
    #

    for i, stage in enumerate(stages):
        print('=== Stage %d (t=%d) ===' % (i+1, ts[-1]))

        while not stage.jettison():
            thrust = stage.get_thrust(poss[-1].magnitude())  # N = Kg*m/s**2
            mass = stage.fuel + stage.dry_mass  # Kg

            t = ts[-1] + 1

            # compute acceleration

            am_thrust = thrust/mass  # m/s**2
            am_drag = -dps_n[-1]/mass  # m/s**2
            am_g = -9.81  # m/s**2
            am_c = vels[-1].x**2 / (6371000 + poss[-1].y)

            # Pitch over


            angle = ANGLES[ts[-1]]
            if angle is None:  # ideal angle
                _a = 180/np.pi*np.arcsin(- (am_g + am_c) / am_thrust)
                angle = min(_a, max_angle)

            av_thrust = Vector.init_from_magangle(am_thrust, angle)
            av_drag = Vector.init_from_magangle(am_drag, angle)
            av_g = Vector.init_from_magangle(am_g, 90)
            av_centrifugal = Vector.init_from_magangle(am_c, 90)

            _thrust += am_thrust
            _cent += am_c

            a = av_thrust + av_drag + av_g + av_centrifugal

            # compute other params

            v = vels[-1] + a  # m/s
            p = poss[-1] + v  # m

            rho = 1.3 * np.e ** ( -p.magnitude() / 7000 )  # kg/m**3
            dp_p = rho * v.magnitude()**2 / 2  # pascal = kg/m**3 * (m/s)**2
            dp_n = dp_p * np.pi*0.6**2 * cd

            t_spent = av_thrust.magnitude()
            t_actual = (av_thrust + av_g).magnitude()
            g_loss = t_actual - t_spent

            # save for next iteration and for plotting

            ts.append(t)
            accs.append(a)
            vels.append(v)
            poss.append(p)
            dps_p.append(dp_p)
            dps_n.append(dp_n)
            drag_losses.append(am_drag)
            g_losses.append(g_loss)

            # dont forget to actually burn fuel

            stage.fuel -= stage.get_fc(poss[-1].magnitude())

            #print('%3d %5.1f %6.1f %8.1f' % (t, a.magnitude(), v.magnitude(), p.magnitude()))

    print('=== Burn out (t=%d)' % ts[-1])

    assert ts[-1] == len(ANGLES), '%d == %d is False' % (ts[-1], len(ANGLES))

    # Print info

    print('Simulation results:')

    print('\tAltitude at t=100 sec:      %.3f km' % (poss[100].y/1000))
    print('\tTraveled downrange          %6d km' % (poss[-1].x/1000))
    print('\tFinal altitude              %6d km' % (poss[-1].y/1000))
    print('\tFinal speed                 %6d m/s' % (vels[-1].magnitude()))

    print()

    print('\tDrag losses                 %6d m/s' % sum(drag_losses))
    print('\tGravity losses              %6d m/s' % sum(g_losses))

    print()

    print('Rocket dv: %d m/s' % _thrust)
    print('Rocket simulated dv: %d - %d - %d - %d = %d m/s' % (
        vels[-1].magnitude(), sum(drag_losses), sum(g_losses), _cent,
        (vels[-1].magnitude() - sum(drag_losses) - sum(g_losses) - _cent),
    ))

    Er = 6371000
    mu = 3.986004418*10**14
    e = v.magnitude()**2/2 - mu/(Er+p.y)
    print('Energy: %10.3f*10**6 (final e = [-29.55 ; -29.66]*10**6)' % (e/10**6))

    a = 1/( 2/(Er+poss[-1].y) - vels[-1].magnitude()**2/mu )

    #print(a-Er)

    T = 2*np.pi*np.sqrt(a**3/mu)
    n = 2*np.pi / T  # mean motion

    v = vels[-1].magnitude()
    dr = vels[-1].y
    r = Er + poss[-1].y

    for true_anomaly in np.linspace(0, np.pi, T/2):  # we are rising (dr > 0) <=> after periapsis and before apoapsis <=> v in [0 ; pi]
        # Solve for e
        # e**2 + r*cos(true_anomaly)/a * e + r/a - 1 = 0
        _b = r/a*np.cos(true_anomaly)
        _c = r/a - 1

        e_squared = (_b/2)*(_b/2 - 1) - _c
        if e_squared < 0:
            continue

        e = np.sqrt(e_squared)

        r0 = a * (1-e**2) / (1*e*np.cos(true_anomaly))
        r1 = a * (1-e**2) / (1*e*np.cos(true_anomaly + 2*np.pi/T))

        if abs(dr - (r1-r0)) > 50:
            continue

        rp = a * (1-e**2) / (1+e*np.cos(0))
        ra = a * (1-e**2) / (1+e*np.cos(np.pi))

        #if rp < Er+100*1000 or ra < Er+100*1000:
        #    continue

        #print('%.2f %.2f' % (true_anomaly, e))
        print('%.2f %.2f -> %.2f %.2f' % (true_anomaly, e, dr, r1-r0))
        print('\t\t\t%d %d' % ((rp-Er)/1000, (ra-Er)/1000))




    """
Ariane A-44L:       Gravity Loss: 1576 m/s Drag Loss: 135 m/s
Atlas I:            Gravity Loss: 1395 m/s Drag Loss: 110 m/s
Delta 7925:         Gravity Loss: 1150 m/s Drag Loss: 136 m/s
Shuttle:            Gravity Loss: 1222 m/s Drag Loss: 107 m/s
Saturn V:           Gravity Loss: 1534 m/s Drag Loss:  40 m/s (!!)
Titan IV/Centaur:   Gravity Loss: 1442 m/s Drag Loss: 156 m/s

    """


    #
    # Plot data
    #


    def iterx(array):
        return [a.x for a in array]
    def itery(array):
        return [a.y for a in array]
    def itermag(array):
        return [a.magnitude() for a in array]

    plt.subplot(2, 2, 1)
    plt.plot(ts, iterx(accs))
    plt.plot(ts, itery(accs))
    plt.plot(ts, itermag(accs))
    plt.hlines([9.81*i for i in range(6)], ts[0], ts[-1])
    plt.legend(['Horizontal', 'Vertical', 'Magnitude'])
    plt.title('Acceleration (m/s**2)')

    plt.subplot(2, 2, 2)
    plt.plot(ts, iterx(vels))
    plt.plot(ts, itery(vels))
    plt.plot(ts, itermag(vels))
    plt.legend(['Horizontal', 'Vertical', 'Magnitude'])
    plt.title('Velocity (m/s)')

    plt.subplot(2, 2, 3)
    plt.plot(ts, iterx(poss))
    plt.plot(ts, itery(poss))
    plt.legend(['Downrange (horiz)', 'Altitude (vert)'])
    plt.title('Position (m)')

    plt.subplot(2, 2, 4)
    plt.plot(ts, dps_p)
    plt.plot(ts, dps_n)
    plt.plot(ts, drag_losses)
    plt.legend(['Pressure (Pascal = kg/(m*s**2))', 'Force (Newton = Kg/(m*s**2))', 'Loss (m/s**2)'])
    plt.title('Dynamic Pressure (Aerodynamic Drag)')
    #plt.xlim((45, 110))
    #plt.ylim((10000, 25000))

    plt.suptitle("Rocket Lab's Electron launch 2018-11-11")

    # optimal pitch over angle
    #plt.show()
    #plt.subplot(2, 2, 1)
    #plt.plot(ts, itermag(accs))
    #plt.subplot(2, 2, 2)
    #plt.plot(ts, 180/np.pi*np.arcsin(1.2 / (1+np.array(itermag(accs)))))


def print_info():
    def calc_speed(h):  # altitude (km)
        Er = 6371000
        mu = 3.986004418 * 10**14
        r = Er + h*1000
        a = (Er+h*1000 + Er+500*1000)/2
        return np.sqrt(mu*(2/r - 1/a))

    print('Known or estimated info on the launch:')  # from press kit / mission or live info / webcast

    # based on Rocket Lab webcast and Gmaps, estimated distance, adding ~2 mn at
    # end of burn traveling at orbital speed
    # 9 cm on map / 200 km scale * 1.5 cm scale + 120 sec*8km/s
    downrange = 9*200/1.5 + 120*8

    print('\tAltitude at t=100 sec:      20.000 km')

    """
         50 sec  dyn pres = 50 % * max
1mn20 =  80 sec  max dyn press +/- 3 sec
1mn35 =  95 sec  altitude 20 km
1mn45 = 105 sec  dy press = 50 %
3mn30 = 210 sec  speed 2.5 km/s
4mn50 = 290 sec  speed 3 km/sec
5mn   = 300 sec  altitude 200 km


    """


    print('\tTraveled downrange:           %d km' % downrange)

    print()

    print('\tFinal periapsis:       [200 ; 250] km')
    print('\tSpeed at periapsis:  [%d ; %d] m/s' % (calc_speed(200), calc_speed(250)))
    print('\tFinal apoapsis:                500 km')


def main():
    print_info()
    fly_rocket()
    plt.show()


main()