triple y_axis = triple( 0, 1, 0 );

real g_acc = -9.81 m/s2 * 1/6; /* acceleration due to gravity */
gravity = g_acc * booster.mass;

component booster {
  shape body_shape = shape "booster_shape.txt" * ( 4.8, 1, 4.8 );
  
  triple position = triple( 0, 12500, 0 ) m;
  triple velocity = -triple( 0, 50, 0 ) m/s;
  mass = 180 kg + fuel;

  altitude = dot( y_axis, booster.position ) - 1 m, clock coriolis;
  fuel = fuel - ( fuel_rate( old total_thrust ) * delta_time ) : 145 kg, 
	clock coriolis;
}

/* target velocity for booster */
real vtarget = if( booster.altitude > 2 kilometer, -150 m/s, 
	           booster.altitude > 1 kilometer, 
		   -( ( booster.altitude - 1 kilometer ) * 0.1 / 1 s + 50 m/s ),
		   booster.altitude > 5 m,
		  -( booster.altitude - 5 m ) / 20 s - .25 m/s,
		  -0.25 m/s ), clock coriolis;

/* if dropping too quickly then delta_v negative */
real delta_v = dot( booster.velocity, y_axis ) - vtarget, clock coriolis;

/* dv1 >= 0 */
dv1 = if( delta_v >= 0 m/s, 0 m/s, -delta_v * 1.2 ), clock coriolis;
/* error < 1 */
error = if (delta_v >= 0 m/s, 0 m/s, error * .7 + dv1 * .1) : 0 m/s, clock coriolis;

/* desired_thrust >= 0? */
thrust = ( dv1 * 1.2 s/m + error * 1.4 s/m ) * -gravity, clock coriolis;

isp = 220 sec;
fuel_rate( thrust ) = thrust / isp / 9.81 m/s2;

/* real thrust values in newtons */
w_thrust = max( min( max_thrust, thrust / 6 ), 0 newton ), clock coriolis;
e_thrust = max( min( max_thrust, thrust / 6 ), 0 newton ), clock coriolis;
nw_thrust = max( min( max_thrust, thrust / 6 ), 0 newton ), clock coriolis;
ne_thrust = max( min( max_thrust, thrust / 6 ), 0 newton ), clock coriolis;
sw_thrust = max( min( max_thrust, thrust / 6 ), 0 newton ), clock coriolis;
se_thrust = max( min( max_thrust, thrust / 6 ), 0 newton ), clock coriolis;
total_thrust = sum( key "*_thrust" of . ) : 0 newton, clock coriolis;

/* gravity */
link booster : force gravity * y_axis, triple ( 0, 0, 0 );

/* constraints */
constraint out_of_fuel = booster.fuel >= 0 kg, clock coriolis;

constraint check_landed = booster.altitude > .75 meter, clock coriolis;

constraint no_crash = check_landed || mag( booster.velocity ) < 5 m/s, 
  clock coriolis;

constraint thrusting = thrust == 0 newton, clock coriolis;

/* THRUSTER STUFF */
transparent max_thrust = 3300 newton;

/* thrusters -- two sets of three.  center thruster lies in the xy plane, at a 
   10 degree angle from the negative y axis.
   the side two are 10 degrees from the negative y axis in the xy plane and 15 
   degrees from the negative y axis in the yz plane */
x_thruster_angle = 10 degree;
z_thruster_angle = 15 degree;

x_tan = tan( x_thruster_angle );
z_tan = tan( z_thruster_angle );

y_side = 1 / sqrt( x_tan * x_tan + 1 + z_tan * z_tan );
x_side = x_tan * y_side;
z_side = z_tan * y_side;

y_main = 1 / sqrt( x_tan * x_tan + 1 );
x_main = x_tan * y_main;

/* thrusters are labeled east_wast if in the pos/neg x direction, or 
   north_south for the pos/neg z direction */
triple thrust_west = gvec( x_main, y_main, 0 );
triple thrust_east = gvec(-x_main, y_main, 0 );

triple thrust_northwest = gvec( x_side, y_side, -z_side );
triple thrust_northeast = gvec(-x_side, y_side, -z_side );
triple thrust_southwest = gvec( x_side, y_side, z_side );
triple thrust_southeast = gvec(-x_side, y_side, z_side );

triple wthrust_bp = triple( -2.4, -.30, 0 ); /* -.50, -.30 scaled */
triple ethrust_bp = triple( 2.4, -.30, 0 );  /* .50, -.30 scaled */

/* forces due to thrusters */
link booster : force w_thrust * thrust_west, wthrust_bp;
link booster : force e_thrust * thrust_east, ethrust_bp;
link booster : force nw_thrust * thrust_northwest, wthrust_bp;
link booster : force ne_thrust * thrust_northeast, ethrust_bp;
link booster : force sw_thrust * thrust_southwest, wthrust_bp;
link booster : force se_thrust * thrust_southeast, ethrust_bp;

default zzz = 10 kg;
default zyz = 10 sec/m;

component a {
  zzz = 1;
  bool zyz = false;
  x = x + 1 : 0;
}
