function cannonball_b(start_velocity, timestep, method, output_file)
//
// Follow the motion of a projectile launched from a big cannon
//   on the surface of the Earth.  Allow for the force of
//   air resistance on the projectile, in a simple way,
//   but do NOT include any changing gravitational force
//
// The "method" argument may be either 
//
//        "euler"             use Euler's method
//  or
//        "heun"              use Heun's method
//
// We compute KE, GPE and total E at each timestep and
//    print values only to screen for debugging 
//
// MWR 4/29/2007

  //
  verbose = 1;

  euler_method = 1;
  heun_method = 2;

  // mass of projectile (kg)
  mass = 10000;
  // radius of projectile (m)
  radius = 2.0;


  // sanity checks
  if (start_velocity <= 0)
    mprintf("cannonball_b: start_velocity %f must be positive \n", start_velocity);
    error("cannonball_b: quitting now \n");
  end
  if (timestep <= 0)
    mprintf("cannonball_b: timestep %f must be positive \n", timestep);
    error("cannonball_b: quitting now \n");
  end

  // parse the method -- allow for initial capital
  if ((method == "Euler") | (method == "euler"))
    use_method = euler_method;
  elseif ((method == "Heun") | (method == "heun")) 
    use_method = heun_method;
  else 
    mprintf("cannonball_b: invalid method %s  \n", method);
    error("cannonball_b: quitting now \n");
  end
  

  // we'll write output to this file
  fid = mopen(output_file, "w");
  if (fid < 0)
    mprintf("cannonball_b: could not open file %s \n", output_file);
    error("cannonball_b: quitting now \n");
  end

  // some physical constants
  dt = timestep;

  // initial properties
  time = 0;
  cur_v = start_velocity;
  cur_height = 0;
  iter = 0;
  

  // main loop, in which we update position and velocity
  //   of the cannonball
  while ((iter == 0) | (cur_height > 0))

    // compute energy of projectile
    [ke, gpe, tot_e] = compute_energy(mass, cur_height, cur_v);

    if (verbose > 0)
      mprintf(" iter %6d  time %9.2f  height %9.2f  vel %9.2f  %9.4e %9.4e %9.4e \n", ...
               iter, time, cur_height, cur_v, ke, gpe, tot_e);
    end
    mfprintf(fid, " %9.2f  %9.2f  %9.2f  %9.4e %9.4e %9.4e \n", ...
             time, cur_height, cur_v, ke, gpe, tot_e);

    // compute current acceleration
    accel = acceleration(mass, radius, cur_height, cur_v)

    if (use_method == euler_method)
      // use Euler's method to compute future vel, height

      // compute the new velocity
      new_v = cur_v + accel*dt;

      // compute the new height 
      new_height = cur_height + cur_v*dt;

    else
      // we are using Heun's method, then

      // first, we predict (poorly) the future
      //   velocity, height, and acceleration
      predict_v = cur_v + accel*dt;
      predict_height = cur_height + predict_v*dt;
      predict_accel = acceleration(mass, radius, predict_height, predict_v);

      // use current and future accel to compute a decent
      //    average acceleration over this timestep
      avg_accel = 0.5*(accel + predict_accel);

      // use average accel to calculate good new velocity
      new_v = cur_v + avg_accel*dt;

      // now compute an average velocity over this timestep,
      //   and use it to get the new height
      avg_v = 0.5*(cur_v + new_v);
      new_height = cur_height + avg_v*dt;

    end


    // prepare for next time through loop
    iter = iter + 1;
    time = time + dt;
    cur_height = new_height;
    cur_v = new_v;

  end

  // print final values, just after the ball strikes
  //    the surface
  [ke, gpe, tot_e] = compute_energy(mass, cur_height, cur_v);

  if (verbose > 0)
    mprintf(" iter %6d  time %9.2f  height %9.2f  vel %9.2f  %9.4e %9.4e %9.4e \n", ...
             iter, time, cur_height, cur_v, ke, gpe, tot_e);
  end
  mfprintf(fid, " %9.2f  %9.2f  %9.2f  %9.4e %9.4e %9.4e \n", ...
             time, cur_height, cur_v, ke, gpe, tot_e);


  
  // close the file
  mclose(fid);



endfunction



/////////////////////////////////////////////////////////////
function accel = acceleration(mass, radius, height, velocity)
//
// Given the mass of a projectile, its height above the
//   surface of the Earth, and its velocity, calculate
//   the acceleration on it.
//
// This could be easy, if assume no air resistance,
//   or complicated, if there's air resistance and
//   perhaps multiple bodies.
//

  // physical constants
  G = 6.67E-11;
  mass_earth = 5.98E24;
  radius_earth = 6.37E6;
  pi = 3.14159265358979312;


  // air resistance will be computed like so:
  //
  //    Force  =  C_air * density * area * velocity^2
  //
  C_air = 0.02;
  area = pi*radius*radius;
  density = air_density(height);
  force_air = C_air * area * density * (velocity*velocity);
  // make sure that the force of air resistance
  //   points opposite to velocity
  if (velocity > 0)
    force_air = 0.0 - force_air;
  end
  
  // this version has a constant gravitational force
  g = 9.8;
  force_grav = 0.0 - (mass * g);

  // acceleration is total force divided by mass
  tot_force = force_air + force_grav;
  accel = tot_force / mass;

endfunction

  

  
/////////////////////////////////////////////////////////////
function [ke, gpe, tot_e] = compute_energy(mass, height, velocity)
//
// Given the mass of a projectile, its height above the
//   surface of the Earth, and its velocity, calculate
//   its current KE, GPE and total energy.

  // physical constants
  G = 6.67E-11;
  mass_earth = 5.98E24;
  radius_earth = 6.37E6;

  // this version has no air resistance, just uses
  //   the Law of Universal Gravitation with Earth
  //   and projectile alone
  tot_rad = radius_earth + height;

  ke = 0.5*mass*(velocity*velocity);
  gpe = 0.0 - ((G*mass_earth*mass)/(tot_rad));
  tot_e = ke + gpe;

endfunction

  

/////////////////////////////////////////////////////////////
function density = air_density(height)
//
// Given the height above the surface of the Earth,
//   compute the current air density (kg/m^2).
//
// If given height is less than zero, use value for height zero.

  // air density at sea level (kg/m^3)
  rho_0 = 1.21;
  // scale height of atmosphere (m)
  scale_height = 8000.0;

  if (height <= 0)
    density = rho_0;
  else
    density = rho_0 * exp(0.0 - (height/scale_height));
  end

endfunction