approximate flight parameters for an iron man-like jet suit

jet.py

from bokeh.plotting import figure, output_file, show
from numpy import sqrt, sin, arctan2

import scipy.integrate as integrate
import numpy as np

import pint

"""
This script approximates the flight range and speed of a jet-powered Iron Man-
style suit according to various parameters.
"""

ureg = pint.UnitRegistry()

# using value for kerosene, which should be close enough
FUEL_DENSITY = 0.819 * ureg.g / ureg.mL
# this value seemed reasonable given starting speed. could be tuned.
MAX_FUEL_VOLUME = 30 * ureg.L

# these parameters found here: http://www.pbsvb.com/customer-industries/aerospace/aircraft-engines/tj-20-turbojet-engine
# http://www.pbsvb.com/getattachment/Zakaznicka-odvetvi/Letectvi/Aircraft-UAV-engines/proudovy-motor-tj-20/Turbojet-engine_TJ20U.pdf.aspx
SPECIFIC_FUEL_CONSUMPTION = 0.165 * ureg.kg / (ureg.N * ureg.hour)
# note: this figure is for takeoff ONLY, so we should expect worse
THRUST_MAX = 210 * ureg.N
ENGINE_MASS = 2.1 * ureg.kg
ENGINE_COUNT = 6

HUMAN_WIDTH = 2 * ureg.foot
HUMAN_HEIGHT = 6 * ureg.foot
HUMAN_MASS = 175 * ureg.lb
HUMAN_DRAG_COEFFICIENT = 1

THRUST_TOTAL = THRUST_MAX * ENGINE_COUNT
ENGINE_MASS_TOTAL = ENGINE_MASS * ENGINE_COUNT

FUEL_MASS_FLOW_RATE = ENGINE_COUNT * THRUST_MAX * SPECIFIC_FUEL_CONSUMPTION
FUEL_VOLUME_FLOW_RATE = FUEL_MASS_FLOW_RATE / FUEL_DENSITY
FLIGHT_TIME = (MAX_FUEL_VOLUME / FUEL_VOLUME_FLOW_RATE).to(ureg.second)

AIR_DENSITY = 0.001275 * ureg.g / ureg.mL
G = 9.8 * ureg.meter / (ureg.second ** 2)
MPH = ureg.mile / ureg.hour

def fuel_volume(t):
    """Remaining fuel volume."""
    return MAX_FUEL_VOLUME - (t * FUEL_VOLUME_FLOW_RATE)

def total_mass(t):
    return HUMAN_MASS + fuel_volume(t) * FUEL_DENSITY + ENGINE_MASS_TOTAL

def lift_force(t):
    """The required lift force to keep our ramshackle aircraft aloft."""
    return total_mass(t) * G

def free_thrust(t):
    return THRUST_TOTAL - lift_force(t)

"""
Assume the pilot forms a rigid, flat body with thrusters oriented axially.
I.e. thrust is being provided in one precise direction in order to maintain
altitude, while using all remaining thrust to contribute to forward motion.
"""
def angle_of_attack(t):
    return arctan2(lift_force(t), free_thrust(t))

def human_area(t):
    return HUMAN_WIDTH * HUMAN_HEIGHT * sin(angle_of_attack(t))

def terminal_velocity(t):
    """
    This expression is derived from the terminal velocity equation for a
    falling body, found here: https://www.wikiwand.com/en/Terminal_velocity,
    where the 'free thrust' term is substituted for the force of gravity on the
    object.
    """
    return sqrt(2 * free_thrust(t) / (AIR_DENSITY * human_area(t) * HUMAN_DRAG_COEFFICIENT))

def terminal_dimensionless(t):
    """
    Compute the terminal velocity as above, but without requiring Pint
    dimensions.

    Primarily for use with scipy and numpy functions that don't play nicely.
    """
    if not hasattr(t, 'magnitude'):
        t = t * ureg.second

    return terminal_velocity(t).to(ureg.meter / ureg.second).magnitude

range_ = integrate.quad(terminal_dimensionless, 0., FLIGHT_TIME.magnitude)[0] * ureg.meter
print("range is {:.2f}, min speed {:.1f}, max speed {:.1f}".format(range_.to(ureg.mile), terminal_velocity(0 * ureg.second).to(MPH), terminal_velocity(FLIGHT_TIME).to(MPH)))

# Produce a graph of the maximum speed of the aircraft over time.
output_file("jet.html")

N = 250
x = np.linspace(0, FLIGHT_TIME.magnitude, N)
v = np.vectorize(terminal_dimensionless)(x) * ureg.meter / ureg.second

fig1 = figure(title='Iron Man Maximum Speed by Flight Time', x_axis_label='Time (minutes)', y_axis_label='Speed (mph)')
fig1.circle(x/60, v.to(MPH).magnitude)

show(fig1)