function cannonball_a(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
//   gravity on the projectile to change with altitude.
//   Do NOT include any air resistance.
//
// 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;


  // sanity checks
  if (start_velocity <= 0)
    mprintf("cannonball_a: start_velocity %f must be positive \n", start_velocity);
    error("cannonball_a: quitting now \n");
  end
  if (timestep <= 0)
    mprintf("cannonball_a: timestep %f must be positive \n", timestep);
    error("cannonball_a: 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_a: invalid method %s  \n", method);
    error("cannonball_a: quitting now \n");
  end
  

  // we'll write output to this file
  fid = mopen(output_file, "w");
  if (fid < 0)
    mprintf("cannonball_a: could not open file %s \n", output_file);
    error("cannonball_a: 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, 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, 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, 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;

  // this version has no air resistance, just uses
  //   the Law of Universal Gravitation with Earth
  //   and projectile alone
  tot_rad = radius_earth + height;
  accel = -(G*mass_earth)/(tot_rad*tot_rad);

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