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

/* booster cylinder 2.3 d, .5 h
   on tripod 1m high, top dist from center ?*/

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

y_vel = dot( y_axis, booster.velocity );

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

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

altitude = dot( y_axis, booster.position ) - 3 m;

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


max_thrust = 3300 newton;

/* desired_thrust >= 0? */
triple desired_thrust = ( dv1 * 1.2 s/m + error * 1.4 s/m /*+.99*/ ) * -gravity;

/* thrust => 0 */
triple thrust = y_axis * max( min( max_thrust, dot( y_axis, desired_thrust ) ), 0 newton );
	            
constraint out_of_fuel = fuel >= 0 kg;
		     
real g_acc = -9.81 m/s2 * 1/6; /* acceleration due to gravity */
triple gravity = y_axis * g_acc * booster.mass;

link booster : force thrust, ( 0, -1, 0 ), ( 0, 0, 0 );
link booster : force gravity, ( 0, 0, 0 ), ( 0, 0, 0 ); // gravity

constraint check_high = altitude > 10 kilometer;
constraint check_low = altitude > 2.5 kilometer;
constraint check_landed = altitude > .75 meter;

constraint no_crash = check_landed || mag( velocity ) < 5 m/s;

triple act_vel = ( velocity - old velocity ) * ( 1 / delta_time );

isp = 220 sec;
fuel_rate( triple thrust ) = mag( thrust ) / isp / 9.81 m/s2;
fuel = fuel - ( fuel_rate( old thrust ) * delta_time ) : 145 kg;

/* 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 );
