Suicide Burn Model

model SuicideBurn
  final constant Real pi=2*asin(1.0); // 3.14159265358979;
  final constant Real g0 = 9.81;

  type Mass = Real(unit="kg");
  type Thrust = Real(unit="N");
  type Length = Real(unit="m");
  type Velocity = Real(unit="m/s");

  parameter Real g = 1.62519 "Target body gravity";
  parameter Real isp = 240 "Rocket engine efficiency";
  parameter Mass dry_mass = 1.0e4 "Payload";
  parameter Thrust engine_thrust = 78.0e3 * 8 "Maximum thrust";

  /* Autopilot configuration */
  parameter Real burn_time(unit="s") = 120. "When to start the engine";
  parameter Velocity cutoff_x = 0.5 "Desired cutoff speed, x-axis";
  parameter Velocity cutoff_y = 3.0 "Desired cutoff speed, y-axis";
  parameter Real restart_factor = 1.5 "Restart when velocity > cutoff * restart_factor";

  /* Autopilot signals */
  Boolean stop = abs(vx) < cutoff_x and abs(vy) < cutoff_y;
  Boolean restart = ((abs(vx) > cutoff_x * restart_factor) or abs(vy) > cutoff_y * restart_factor)
                    and time > burn_time;
  Boolean start = time > burn_time;

  discrete Real throttle;
  discrete Real max_thrust(start=engine_thrust, fixed=true);

  /* Geometry */
  Length x,y;
  Real phi "Prograde direction";

  Velocity vx,vy;

  Mass m_p,m;

  /* Thrust */
  Thrust F,Fx,Fy;

  initial equation
  y = 15000;
  vx = 1000.0;
  m_p = 1400.0 * 8;
  vy = 0;

  equation

  /* Autopilot implementation */
  when start then
    throttle = 1.0;
  elsewhen stop then
    throttle = 0.0;
  elsewhen restart then
    throttle = 1.0;
  end when;

  /* Propellant limitations */
  when m_p < 0 then max_thrust = 0.0; end when;

  /* End-of-flight */
  when y < 0 and abs(vy) < (restart_factor * cutoff_y) and abs(vx) < (restart_factor * cutoff_x) then
    terminate("Landed");
  elsewhen y < 0 and abs(vy) >= (restart_factor * cutoff_y) then
    terminate("Crash: Vertical");
  elsewhen y < 0 and abs(vx) >= (restart_factor * cutoff_x) then
    terminate("Crash: Horizontal");
  end when;

  /* Flight physics */
  der(x) = vx;
  der(y) = vy;

  der(vx) = Fx / m;
  der(vy) = (Fy / m) - g;

  m = dry_mass + m_p;
  F = throttle * max_thrust;
  der(m_p) = -F / (g0 * isp);

  /* craft geometry, thrust points retrograde */
  phi = atan2(vx, vy);
  Fy = F*cos(pi + phi);
  Fx = F*sin(pi + phi);

  annotation(experiment(StartTime=0.0, StopTime=2000.0));

end SuicideBurn;