orbitx/physics/ode_solver.py
def simulation_differential_function(
    t: float, y_1d: np.ndarray,
    pass_through_state: PhysicalState,
    masses: np.ndarray, artificials: np.ndarray
) -> np.ndarray:
    """
    Given the state of a simulation, calculate its instantaneous rate-of-change.

    This function is probably one of the most important ones in OrbitX :)
    Everything from "an object with a velocity will change position" to
    "the habitat reactor heats up when it's being used" is encoded here.

    This function is used by scipy.solve_ivp, which unfortunately requires
    that the y_1d input of this function is a single 1-dimensional array.
    We also have the same requirement of this function's output. This means
    that we have to encode almost the entire state of the OrbitX simulation
    in the y_1d input array, and every single value in that array has to have
    a corresponding output value in the output array. This is a big array!

    This is roughly the ordering of fields in the input array:

    y_1d =
     [X, Y, VX, VY, Heading, Spin, Fuel, Throttle, LandedOn, Broken] +
     SRB_time_left + time_acc (these are both single values) +
     [ComponentData, CoolantLoopData, RadiatorData]
    returns the derivative of y_1d, i.e.
    [VX, VY, AX, AY, Spin, 0, Fuel consumption, 0, 0, 0] + -constant + 0
    (zeroed-out fields are changed elsewhere)

    The ordering of fields is determined by the ordering in orbitx.proto.
    If the structure of the y_1d array changes, make sure to update this
    function and anywhere annotated with #Y_VECTOR_CHANGESITE!

    !!!!!!!!!!! IMPORTANT !!!!!!!!!!!
    This function should return a DERIVATIVE. The numpy.solve_ivp function
    will do the rest of the work of the simulation, this function just
    describes how things _move_.
    At its most basic level, this function takes in the _position_ of
    everything (plus some random stuff), and returns the _velocity_ of
    everything (plus some random stuff).
    Essentially, numpy.solve_ivp does this calculation:
    new_positions_of_system = t_delta * _derive(
                                            current_t_of_system,
                                            current_y_of_system)

    Any arguments after `t` and `y_1d` are just extra constants and
    pass-through state, which should remain constant during normal simulation.
    They are passed in to speed up computation, since this function is the most
    performance-sensitive part of the orbitx codebase(!!!)
    """

    # Note: we create this y as a PhysicsState for convenience, but if you
    # set any values of y, the changes will be discarded! The only way they
    # will be propagated out of this function is by numpy using the return
    # value of this function as a derivative, as explained above.
    # If you want to set values in y, look at _reconcile_entity_dynamics.
    y = PhysicsState(y_1d, pass_through_state)
    acc_matrix = calc.grav_acc(y.X, y.Y, masses, y.Fuel)
    zeros = np.zeros(y._n)
    fuel_cons = np.zeros(y._n)

    # Engine thrust and fuel consumption
    for artif_index in artificials:
        if y[artif_index].fuel > 0 and y[artif_index].throttle > 0:
            # We have fuel remaining, calculate thrust
            entity = y[artif_index]
            capability = common.craft_capabilities[entity.name]

            fuel_cons[artif_index] = \
                -abs(capability.fuel_cons * entity.throttle)
            eng_thrust = (capability.thrust * entity.throttle
                          * calc.heading_vector(entity.heading))
            mass = entity.mass + entity.fuel

            if entity.name == AYSE and \
                    y[HABITAT].landed_on == AYSE:
                # It's bad that this is hardcoded, but it's also the only
                # place that this comes up so IMO it's not too bad.
                hab = y[HABITAT]
                mass += hab.mass + hab.fuel

            eng_acc = eng_thrust / mass
            acc_matrix[artif_index] += eng_acc

    # And SRB thrust
    srb_usage = 0
    try:
        if y.srb_time >= 0:
            hab_index = y._name_to_index(HABITAT)
            hab = y[hab_index]
            srb_acc = common.SRB_THRUST / (hab.mass + hab.fuel)
            srb_acc_vector = srb_acc * calc.heading_vector(hab.heading)
            acc_matrix[hab_index] += srb_acc_vector
            srb_usage = -1
    except PhysicsState.NoEntityError:
        # The Habitat doesn't exist.
        pass

    # Drag effects
    craft = y.craft
    if craft is not None:
        craft_index = y._name_to_index(y.craft)
        drag_acc = calc.drag(y)
        acc_matrix[craft_index] -= drag_acc

    # Centripetal acceleration to keep landed entities glued to each other.
    landed_on = y.LandedOn
    for landed_i in landed_on:
        lander = y[landed_i]
        ground = y[landed_on[landed_i]]

        centripetal_acc = (lander.pos - ground.pos) * ground.spin ** 2
        acc_matrix[landed_i] = \
            acc_matrix[landed_on[landed_i]] - centripetal_acc

    # Time for the Engineering section.
    component_resistances = electrofunctions.component_resistances(y.engineering)
    electrical_bus = electrofunctions.bus_electricity(component_resistances)
    bus_power = electrical_bus.voltage * electrical_bus.current

    # Just lifted this from DM's enghabw.bas code. Probably has deeper significance.
    fuel_usage_rate = .05 * bus_power / 50_000_000
    # TODO: Replace the existing fuel consumption calculations with this calculation.

    component_heating_rate_array = electrofunctions.component_heating_rate(
        y.engineering,
        component_resistances,
        electrical_bus.voltage
    )

    # Sets velocity and spin of a couple more entities.
    # If you want to set the acceleration of an entity, do it above and
    # keep that logic in _derive. If you want to set the velocity and spin
    # or any other fields that an Entity has, you should put that logic in
    # this _reconcile_entity_dynamics helper.
    y = helpers._reconcile_entity_dynamics(y)

    # If the ordering of fields in orbitx.proto changes, change here too. #Y_VECTOR_CHANGESITE
    y_differential = np.concatenate((
        y.VX, y.VY, np.hsplit(acc_matrix, 2), y.Spin,
        zeros, fuel_cons, zeros, zeros, zeros, np.array([srb_usage, 0]),
        np.zeros(N_COMPONENTS), np.zeros(N_COMPONENTS),  # connected and capacity fields
        component_heating_rate_array,
        np.zeros(N_COMPONENTS), np.zeros(N_COMPONENTS), np.zeros(N_COMPONENTS),  # coolant connection states
        np.zeros(N_COOLANT_LOOPS * _N_COOLANT_FIELDS),
        np.zeros(N_RADIATORS * _N_RADIATOR_FIELDS)
    ), axis=None)

    return y_differential

    Given the state of a simulation, calculate its instantaneous rate-of-change.

    This function is probably one of the most important ones in OrbitX :)
    Everything from "an object with a velocity will change position" to
    "the habitat reactor heats up when it's being used" is encoded here.

    This function is used by scipy.solve_ivp, which unfortunately requires
    that the y_1d input of this function is a single 1-dimensional array.
    We also have the same requirement of this function's output. This means
    that we have to encode almost the entire state of the OrbitX simulation
    in the y_1d input array, and every single value in that array has to have
    a corresponding output value in the output array. This is a big array!

    This is roughly the ordering of fields in the input array:

    y_1d =
     [X, Y, VX, VY, Heading, Spin, Fuel, Throttle, LandedOn, Broken] +
     SRB_time_left + time_acc (these are both single values) +
     [ComponentData, CoolantLoopData, RadiatorData]
    returns the derivative of y_1d, i.e.
    [VX, VY, AX, AY, Spin, 0, Fuel consumption, 0, 0, 0] + -constant + 0
    (zeroed-out fields are changed elsewhere)

    The ordering of fields is determined by the ordering in orbitx.proto.
    If the structure of the y_1d array changes, make sure to update this
    function and anywhere annotated with #Y_VECTOR_CHANGESITE!

    !!!!!!!!!!! IMPORTANT !!!!!!!!!!!
    This function should return a DERIVATIVE. The numpy.solve_ivp function
    will do the rest of the work of the simulation, this function just
    describes how things _move_.
    At its most basic level, this function takes in the _position_ of
    everything (plus some random stuff), and returns the _velocity_ of
    everything (plus some random stuff).
    Essentially, numpy.solve_ivp does this calculation:
    new_positions_of_system = t_delta * _derive(
                                            current_t_of_system,
                                            current_y_of_system)

    Any arguments after `t` and `y_1d` are just extra constants and
    pass-through state, which should remain constant during normal simulation.
    They are passed in to speed up computation, since this function is the most
    performance-sensitive part of the orbitx codebase(!!!)
    
Given the state of a simulation, calculate its instantaneous rate-of-change.
This function is probably one of the most important ones in OrbitX :)
Everything from "an object with a velocity will change position" to
"the habitat reactor heats up when it's being used" is encoded here.

This function is used by scipy.solve_ivp, which unfortunately requires
that the y_1d input of this function is a single 1-dimensional array.
We also have the same requirement of this function's output. This means
that we have to encode almost the entire state of the OrbitX simulation
in the y_1d input array, and every single value in that array has to have
a corresponding output value in the output array. This is a big array!

This is roughly the ordering of fields in the input array.

y_1d
[X, Y, VX, VY, Heading, Spin, Fuel, Throttle, LandedOn, Broken] +
SRB_time_left + time_acc (these are both single values) +
[ComponentData, CoolantLoopData, RadiatorData]
returns the derivative of y_1d, i.e.
[VX, VY, AX, AY, Spin, 0, Fuel consumption, 0, 0, 0] + -constant + 0
(zeroed-out fields are changed elsewhere)

The ordering of fields is determined by the ordering in orbitx.proto.
If the structure of the y_1d array changes, make sure to update this
function and anywhere annotated with #Y_VECTOR_CHANGESITE!

IMPORTANT
This function should return a DERIVATIVE. The numpy.solve_ivp function
will do the rest of the work of the simulation, this function just
describes how things _move_.
At its most basic level, this function takes in the _position_ of
everything (plus some random stuff), and returns the _velocity_ of
everything (plus some random stuff).

Any arguments after `t` and `y_1d` are just extra constants and
pass-through state, which should remain constant during normal simulation.
They are passed in to speed up computation, since this function is the most
performance-sensitive part of the orbitx codebase(!!!)
Given the state of a simulation calculate its instantaneous rate - of - change . This function is probably one of the most important ones in OrbitX : ) Everything from " an object with a velocity will change position " to " the habitat reactor heats up when it ' s being used " is encoded here . This function is used by scipy . solve_ivp which unfortunately requires that the y_1d input of this function is a single 1 - dimensional array . We also have the same requirement of this function ' s output . This means that we have to encode almost the entire state of the OrbitX simulation in the y_1d input array and every single value in that array has to have a corresponding output value in the output array . This is a big array! This is roughly the ordering of fields in the input array . y_1d [ X Y VX VY Heading Spin Fuel Throttle LandedOn Broken ] + SRB_time_left + time_acc ( these are both single values ) + [ ComponentData CoolantLoopData RadiatorData ] returns the derivative of y_1d i . e . [ VX VY AX AY Spin 0 Fuel consumption 0 0 0 ] + - constant + 0 ( zeroed - out fields are changed elsewhere ) The ordering of fields is determined by the ordering in orbitx . proto . If the structure of the y_1d array changes make sure to update this function and anywhere annotated with #Y_VECTOR_CHANGESITE! IMPORTANT This function should return a DERIVATIVE . The numpy . solve_ivp function will do the rest of the work of the simulation this function just describes how things _move_ . At its most basic level this function takes in the _position_ of everything ( plus some random stuff ) and returns the _velocity_ of everything ( plus some random stuff ) . Any arguments after ` t ` and ` y_1d ` are just extra constants and pass - through state which should remain constant during normal simulation . They are passed in to speed up computation since this function is the most performance - sensitive part of the orbitx codebase ( !!! )
def simulation_differential_function(
    t: float, y_1d: np.ndarray,
    pass_through_state: PhysicalState,
    masses: np.ndarray, artificials: np.ndarray
) -> np.ndarray:
    y = PhysicsState(y_1d, pass_through_state)
    acc_matrix = calc.grav_acc(y.X, y.Y, masses, y.Fuel)
    zeros = np.zeros(y._n)
    fuel_cons = np.zeros(y._n)
    for artif_index in artificials:
        if y[artif_index].fuel > 0 and y[artif_index].throttle > 0:
            entity = y[artif_index]
            capability = common.craft_capabilities[entity.name]
            fuel_cons[artif_index] = \
                -abs(capability.fuel_cons * entity.throttle)
            eng_thrust = (capability.thrust * entity.throttle
                          * calc.heading_vector(entity.heading))
            mass = entity.mass + entity.fuel
            if entity.name == AYSE and \
                    y[HABITAT].landed_on == AYSE:
                hab = y[HABITAT]
                mass += hab.mass + hab.fuel
            eng_acc = eng_thrust / mass
            acc_matrix[artif_index] += eng_acc
    srb_usage = 0
    try:
        if y.srb_time >= 0:
            hab_index = y._name_to_index(HABITAT)
            hab = y[hab_index]
            srb_acc = common.SRB_THRUST / (hab.mass + hab.fuel)
            srb_acc_vector = srb_acc * calc.heading_vector(hab.heading)
            acc_matrix[hab_index] += srb_acc_vector
            srb_usage = -1
    except PhysicsState.NoEntityError:
        pass
    craft = y.craft
    if craft is not None:
        craft_index = y._name_to_index(y.craft)
        drag_acc = calc.drag(y)
        acc_matrix[craft_index] -= drag_acc
    landed_on = y.LandedOn
    for landed_i in landed_on:
        lander = y[landed_i]
        ground = y[landed_on[landed_i]]
        centripetal_acc = (lander.pos - ground.pos) * ground.spin ** 2
        acc_matrix[landed_i] = \
            acc_matrix[landed_on[landed_i]] - centripetal_acc
    component_resistances = electrofunctions.component_resistances(y.engineering)
    electrical_bus = electrofunctions.bus_electricity(component_resistances)
    bus_power = electrical_bus.voltage * electrical_bus.current
    fuel_usage_rate = .05 * bus_power / 50_000_000
    component_heating_rate_array = electrofunctions.component_heating_rate(
        y.engineering,
        component_resistances,
        electrical_bus.voltage
    )
    y = helpers._reconcile_entity_dynamics(y)
    y_differential = np.concatenate((
        y.VX, y.VY, np.hsplit(acc_matrix, 2), y.Spin,
        zeros, fuel_cons, zeros, zeros, zeros, np.array([srb_usage, 0]),
        np.zeros(N_COMPONENTS), np.zeros(N_COMPONENTS),  
        component_heating_rate_array,
        np.zeros(N_COMPONENTS), np.zeros(N_COMPONENTS), np.zeros(N_COMPONENTS),  
        np.zeros(N_COOLANT_LOOPS * _N_COOLANT_FIELDS),
        np.zeros(N_RADIATORS * _N_RADIATOR_FIELDS)
    ), axis=None)
    return y_differential
def simulation_differential_function ( t : float , y_1d : np . ndarray , pass_through_state : PhysicalState , masses : np . ndarray , artificials : np . ndarray ) -> np . ndarray : y = PhysicsState ( y_1d , pass_through_state ) acc_matrix = calc . grav_acc ( y . X , y . Y , masses , y . Fuel ) zeros = np . zeros ( y . _n ) fuel_cons = np . zeros ( y . _n ) for artif_index in artificials : if y [ artif_index ] . fuel > 0 and y [ artif_index ] . throttle > 0 : entity = y [ artif_index ] capability = common . craft_capabilities [ entity . name ] fuel_cons [ artif_index ] = - abs ( capability . fuel_cons * entity . throttle ) eng_thrust = ( capability . thrust * entity . throttle * calc . heading_vector ( entity . heading ) ) mass = entity . mass + entity . fuel if entity . name == AYSE and y [ HABITAT ] . landed_on == AYSE : hab = y [ HABITAT ] mass += hab . mass + hab . fuel eng_acc = eng_thrust / mass acc_matrix [ artif_index ] += eng_acc srb_usage = 0 try : if y . srb_time >= 0 : hab_index = y . _name_to_index ( HABITAT ) hab = y [ hab_index ] srb_acc = common . SRB_THRUST / ( hab . mass + hab . fuel ) srb_acc_vector = srb_acc * calc . heading_vector ( hab . heading ) acc_matrix [ hab_index ] += srb_acc_vector srb_usage = - 1 except PhysicsState . NoEntityError : pass craft = y . craft if craft is not None : craft_index = y . _name_to_index ( y . craft ) drag_acc = calc . drag ( y ) acc_matrix [ craft_index ] -= drag_acc landed_on = y . LandedOn for landed_i in landed_on : lander = y [ landed_i ] ground = y [ landed_on [ landed_i ] ] centripetal_acc = ( lander . pos - ground . pos ) * ground . spin ** 2 acc_matrix [ landed_i ] = acc_matrix [ landed_on [ landed_i ] ] - centripetal_acc component_resistances = electrofunctions . component_resistances ( y . engineering ) electrical_bus = electrofunctions . bus_electricity ( component_resistances ) bus_power = electrical_bus . voltage * electrical_bus . current fuel_usage_rate = .05 * bus_power / 50_000_000 component_heating_rate_array = electrofunctions . component_heating_rate ( y . engineering , component_resistances , electrical_bus . voltage ) y = helpers . _reconcile_entity_dynamics ( y ) y_differential = np . concatenate ( ( y . VX , y . VY , np . hsplit ( acc_matrix , 2 ) , y . Spin , zeros , fuel_cons , zeros , zeros , zeros , np . array ( [ srb_usage , 0 ] ) , np . zeros ( N_COMPONENTS ) , np . zeros ( N_COMPONENTS ) , component_heating_rate_array , np . zeros ( N_COMPONENTS ) , np . zeros ( N_COMPONENTS ) , np . zeros ( N_COMPONENTS ) , np . zeros ( N_COOLANT_LOOPS * _N_COOLANT_FIELDS ) , np . zeros ( N_RADIATORS * _N_RADIATOR_FIELDS ) ) , axis = None ) return y_differential
Given the state of a simulation, calculate its instantaneous rate-of-change.
Given the state of a simulation calculate its instantaneous rate - of - change .
"""
    Given the state of a simulation, calculate its instantaneous rate-of-change.

    This function is probably one of the most important ones in OrbitX :)
    Everything from "an object with a velocity will change position" to
    "the habitat reactor heats up when it's being used" is encoded here.

    This function is used by scipy.solve_ivp, which unfortunately requires
    that the y_1d input of this function is a single 1-dimensional array.
    We also have the same requirement of this function's output. This means
    that we have to encode almost the entire state of the OrbitX simulation
    in the y_1d input array, and every single value in that array has to have
    a corresponding output value in the output array. This is a big array!

    This is roughly the ordering of fields in the input array:

    y_1d =
     [X, Y, VX, VY, Heading, Spin, Fuel, Throttle, LandedOn, Broken] +
     SRB_time_left + time_acc (these are both single values) +
     [ComponentData, CoolantLoopData, RadiatorData]
    returns the derivative of y_1d, i.e.
    [VX, VY, AX, AY, Spin, 0, Fuel consumption, 0, 0, 0] + -constant + 0
    (zeroed-out fields are changed elsewhere)

    The ordering of fields is determined by the ordering in orbitx.proto.
    If the structure of the y_1d array changes, make sure to update this
    function and anywhere annotated with #Y_VECTOR_CHANGESITE!

    !!!!!!!!!!! IMPORTANT !!!!!!!!!!!
    This function should return a DERIVATIVE. The numpy.solve_ivp function
    will do the rest of the work of the simulation, this function just
    describes how things _move_.
    At its most basic level, this function takes in the _position_ of
    everything (plus some random stuff), and returns the _velocity_ of
    everything (plus some random stuff).
    Essentially, numpy.solve_ivp does this calculation:
    new_positions_of_system = t_delta * _derive(
                                            current_t_of_system,
                                            current_y_of_system)

    Any arguments after `t` and `y_1d` are just extra constants and
    pass-through state, which should remain constant during normal simulation.
    They are passed in to speed up computation, since this function is the most
    performance-sensitive part of the orbitx codebase(!!!)
    """ # Note: we create this y as a PhysicsState for convenience, but if you # set any values of y, the changes will be discarded! The only way they # will be propagated out of this function is by numpy using the return # value of this function as a derivative, as explained above. # If you want to set values in y, look at _reconcile_entity_dynamics. # Engine thrust and fuel consumption # We have fuel remaining, calculate thrust # It's bad that this is hardcoded, but it's also the only # place that this comes up so IMO it's not too bad. # And SRB thrust # The Habitat doesn't exist. # Drag effects # Centripetal acceleration to keep landed entities glued to each other. # Time for the Engineering section. # Just lifted this from DM's enghabw.bas code. Probably has deeper significance. # TODO: Replace the existing fuel consumption calculations with this calculation. # Sets velocity and spin of a couple more entities. # If you want to set the acceleration of an entity, do it above and # keep that logic in _derive. If you want to set the velocity and spin # or any other fields that an Entity has, you should put that logic in # this _reconcile_entity_dynamics helper. # If the ordering of fields in orbitx.proto changes, change here too. #Y_VECTOR_CHANGESITE # connected and capacity fields # coolant connection states
t
float
y_1d
np.ndarray
pass_through_state
PhysicalState
masses
np.ndarray
artificials
np.ndarray