Launch to Orbit

To encourage fellow coders to try this task, below is a python-based implementation that allows you to execute a rocket flight. You still need to write and tune your own flight control program - the example used is the one provided in the text. Graphical output is provided via pyplot. This might make it easier to experiment with different approaches at home.

from math import pi, sin, cos, atan
import matplotlib.pyplot as plt

# constants
R_EARTH = 6_371_000  # radius earth (m)
G = -9.81  # gravity of earth (m/s^2)

M0 = 272_000  # default mass of rocket (kg)
FTH = 3_000_000  # default thrust force (N = kg*m/s^2)
VM = -800  # default fuel burning rate (kg/s)
DT = 1 / 4  # deafult time step (s)

H_TOLERANCE = 20  # height tolerenace (km)
V_TOLERANCE = 200  # velocity tolerance (m/s)
A_TOLERANCE = 1.0  # angle tolerance (degree)


class Rocket:
    ON_GROUND, IN_FLIGHT, CRASH, SHUTDOWN = "on_ground", "in_flight", "crash", "shutdown"

    def __init__(self, mass_init, force_thrust, fuel_burning_rate, time_step, h_target, v_target):
        self.h_target = h_target  # target height (m)
        self.v_target = v_target  # target velocity (m/s)
        self.m0, self.m = mass_init, mass_init  # initial and current mass (kg)
        self.fth = force_thrust  # thrust force (N = kg*m/s^2)
        self.vm = fuel_burning_rate  # fuel burning rate (kg/s)
        self.x, self.y = [0.0], [R_EARTH]  # coordinates (m)
        self.vx, self.vy = 0.0, 0.0  # velocity (m/s)
        self.att = 0.0  # attitude / angle (radiant)
        self.va = 0.0  # attitude velocity
        self.fa = 0.0  # attitude turning force
        self.t, self.dt = 0.0, time_step  # time (s)
        self.stage = 0
        self.off = 0
        self.flight_status = Rocket.ON_GROUND
        self.show_output = False

    def height(self):
        return (self.x[-1] ** 2 + self.y[-1] ** 2) ** 0.5 - R_EARTH

    def velocity(self):
        return (self.vx ** 2 + self.vy ** 2) ** 0.5

    def height_error(self):
        return self.height() - self.h_target

    def velocity_error(self):
        return self.velocity() - self.v_target

    def angle_error(self):
        at1 = atan(self.x[-1] / self.y[-1])
        at2 = atan(self.vx / self.vy) if self.vx < self.vy else atan(self.vy / - self.vx) + pi / 2
        return (at1 + pi / 2 - at2) * 180 / pi  # corrected

    def set_show_output(self, yn):
        self.show_output = yn

    def set_flight_status(self, status):
        self.flight_status = status
        if self.show_output:
            print(f"t={self.t}: {self.flight_status}")

    def fly(self, program=None):
        self.set_flight_status(Rocket.IN_FLIGHT)
        while True:
            self.exec_move()
            if self.flight_status in [Rocket.CRASH, Rocket.SHUTDOWN]:
                break
            self.pilot(program)
            if self.stage:
                self.exec_stage()
            if self.off:
                self.set_flight_status(Rocket.SHUTDOWN)
                break

    def exec_move(self):
        r = (self.x[-1] ** 2 + self.y[-1] ** 2) ** 0.5
        ag = G * (R_EARTH / r) ** 2
        ath = self.fth / self.m
        ax = ag * self.x[-1] / r + ath * sin(self.att)
        ay = ag * self.y[-1] / r + ath * cos(self.att)
        da = 0.001 * (self.m0 / self.m) * self.fa
        self.x += [self.x[-1] + self.vx * self.dt]
        self.y += [self.y[-1] + self.vy * self.dt]
        self.vx += ax * self.dt
        self.vy += ay * self.dt
        self.att += self.va * self.dt
        self.va += da * self.dt
        self.m += self.vm * self.dt
        self.t += self.dt
        if self.show_output and self.t == int(self.t):
            print(f"t={self.t:,.0f}  "
                  f"x={self.x[-1]:,.0f}  y={self.y[-1]:,.0f}  h={self.height():,.0f}  "
                  f"vx={self.vx:,.0f}  vy={self.vy:,.0f}  v={self.velocity():,.0f}  "
                  f"att={self.att:,.2f}  fa={self.fa:,.4f}  m={self.m:,.0f}  "
                  f"h_error={self.height_error():,.0f}  "
                  f"v_error={self.velocity_error():,.0f}  "
                  f"a_error={self.angle_error():,.1f}  ")
        if self.att > pi or self.att < - pi / 4 or self.height() < 0:
            self.set_flight_status(Rocket.CRASH)
        elif self.m <= 10_000:
            self.set_flight_status(Rocket.SHUTDOWN)

    def exec_stage(self):
        if self.vm <= -800:
            self.vm /= 4.0
            self.fth /= 4.0
        if self.show_output:
            print(f"t={self.t}: stage")

    def pilot(self, program=None):
        if program is not None:
            self.stage, self.off = 0, 0
            for i in range(0, len(program), 2):
                tc, cmd = program[i], program[i + 1]
                if tc == self.t:
                    if cmd == "stage":
                        self.stage = 1
                    elif cmd == "off":
                        self.off = 1
                    else:
                        self.fa = min(1.0, max(-1.0, cmd))
                    return

Continued:

    def display_flight(self):
        plt.plot([R_EARTH * cos(2 * pi / 1_000 * i) for i in range(1_000)],
                 [R_EARTH * sin(2 * pi / 1_000 * i) for i in range(1_000)],
                 "g-", label="earth")
        plt.plot([(R_EARTH + self.h_target) * cos(2 * pi / 1_000 * i) for i in range(1_000)],
                 [(R_EARTH + self.h_target) * sin(2 * pi / 1_000 * i) for i in range(1_000)],
                 "b--", label="orbit")
        plt.plot(self.x, self.y, "r-", label="rocket")
        plt.xlim([-500_000, 2_500_000])
        plt.ylim([R_EARTH - 300_000, R_EARTH + 1_000_000])
        plt.title("rocket flight")
        plt.legend()
        plt.show()


# inputs
h_target, v_target = 530, 8200

# ***CHANGE HERE*** manual flight program
program = [10, 1, 30, -0.98, 50, 0, 200, "stage", 210, 0.005, 700, "off"]

# simulate rocket path
rocket = Rocket(M0, FTH, VM, DT, h_target * 1_000, v_target)
rocket.set_show_output(True)
rocket.fly(program)

# results
print()
print(f"flight status: {rocket.flight_status}")
if rocket.flight_status == Rocket.SHUTDOWN:
    h_error = rocket.height_error() / 1_000
    v_error = rocket.velocity_error()
    a_error = rocket.angle_error()
    print(f"h_error: {h_error:.1f} km   {'(ok)' if abs(h_error) <= H_TOLERANCE else '(too large)'}")
    print(f"v_error: {v_error:.1f} m/s   {'(ok)' if abs(v_error) <= V_TOLERANCE else '(too large)'}")
    print(f"a_error: {a_error:.1f} degree   {'(ok)' if abs(a_error) <= A_TOLERANCE else '(too large)'}")
    print("flight program:", *program)

rocket.display_flight()