... and here is a translation of the VBA code to the Prime
... TODO: add a graphical user interface
I think, sometimes it is better to start over from scratch instead of trying to convert "old spagheti" to "cuisine nouvelle".
Code:
//
// Simulation of a vertical descent to a planet
//
// Based on analytical integration of equations
// of motion for variable mass
//
// Martin Hepperle
// 2015
// forward declarations
ShowState (t,h,mf,v,Thrust);
EXPORT Lander()
BEGIN
// metric units
// no resistance of air
// no lateral motion
// no graphics
// current time and time step
LOCAL t, dt;
// current altitude
LOCAL h, dh;
// current speed
LOCAL v, dv;
// current total mass
LOCAL m, dm;
// fuel mass
LOCAL mf;
// specific fuel consumption of engine
LOCAL SFC;
// gravity acceleration
LOCAL g;
// engine thrust (commanded by pilots input)
LOCAL Thrust;
// count
LOCAL i, xx;
// ===== set initial conditions
// time
t := 0;
// time step [s]
dt := 0.25;
// thrust (commanded by pilot)
Thrust := 0;
// initial altitude (+=up)
h := 100;
// initial velocity of vehicle [kg] (+=up)
v := 0;
// initial mass of vehicle [kg]
m := 1000;
// fuel mass (part of m)
mf := 46;
// specific fuel consumption of propulsion unit [kg/N/s]
// SFC is the fuel mass per second burnt for one Newton of thrust
SFC := 10 / 3600;
// gravity acceleration [m/s^2]
g := 1.225;
RECT_P (0,0,400,200,0,0);
// output initial state
TEXTOUT_P( "t [s]" + " " + "H [m]" + " " + "mf [kg]" + " " + "v [m/s]"+ " " + "T [N]",
10,40,4,#FFFFFF,400,#000000);
ShowState (t,h,mf,v,Thrust);
FOR i FROM 1 TO 30000 DO
// handle one time step with constant thrust
WAIT(0.1);
IF mf > 0.0 THEN
IF ISKEYDOWN(2) THEN
Thrust:=Thrust+250;
ELSE
IF ISKEYDOWN(12) THEN
Thrust:=Thrust-250;
END;
END;
ELSE
Thrust:=0;
END;
// here we could implement a control law (autopilot)
IF h < 0.5 THEN
// feet sensors (0.5 meters below the lander)
// are touching the ground
Thrust := 0
END;
//
// what happens during this time step dt?
//
// change of mass (negative: burning fuel)
dm := -SFC * Thrust * dt;
// change of velocity (momentum equation for variable mass
// analytically integrated over time)
dv := -(g * dt) - 1 / SFC * LN(1 - Thrust * SFC * dt / m);
// change of height (momentum equation for variable mass
// analytically integrated over time)
IF Thrust < 0.000000000001 THEN
// limit case: Thrust=0
dh := (v - g * dt / 2) * dt;
ELSE
// general case: Thrust>0
xx := LN(1 - SFC * Thrust * dt / m);
dh := ((1 - xx) / SFC + v - g / 2 * dt) * dt + xx * m /
(Thrust * SFC * SFC);
END;
// check height and fuel
IF h < -dh THEN
dh := -h;
END;
IF mf < -dm THEN
dm := -mf;
END;
// update state
m := m + dm;
mf := mf + dm;
v := v + dv;
h := h + dh;
t := t + dt;
// output current state
ShowState (t,h,mf,v,Thrust);
IF h <= 0.01 THEN
IF v > -5 THEN
TEXTOUT_P("soft landing",
10,125,4,#DDFFDD,400,#000000);
ELSE
TEXTOUT_P("you broke: landing gear, your spine",
10,125,4,#FFFFFF,400,#000000);
END;
TEXTOUT_P("fuel remaining: " + STRING(mf,2,1) + " kg",
10,90,4,#FFFFFF,400,#000000);
BREAK; // FOR
END;
END; // FOR
// output final state
ShowState (t,h,mf,v,Thrust);
FREEZE;
return 0;
END;
// show the current state vector of the vehicle motion
ShowState (t,h,mf,v,Thrust)
BEGIN
RECT_P (10,70,300,90,0,0);
TEXTOUT_P(STRING(t,2,2), 10,70,4,#FFFFFF,500,#000000);
TEXTOUT_P(STRING(h,2,1), 70,70,4,#FFFFFF,500,#000000);
TEXTOUT_P(STRING(mf,2,1), 130,70,4,#FFFFFF,500,#000000);
TEXTOUT_P(STRING(v,2,1), 190,70,4,#FFFFFF,500,#000000);
TEXTOUT_P(STRING(Thrust,2,0),250,70,4,#FFFFFF,500,#000000);
// have a look at the fuel gauge
IF mf < 1.0 THEN
TEXTOUT_P(" FUEL! ", 10,150,5,#FFFF00,100,#FF0000);
END;
END;