#!/usr/bin/perl

# function orbit(timestep, period, output_file)
#
# Calculate the motion of the Earth in its orbit
# around the Sun.  Use Euler's method with the given
# timestep.  Place output lines of the form
#
#        time   xpos  ypos    xvel  yvel
#
# for the Earth into the given 'output_file'.  The units
# are days, meters, and meter/second.
#
# Arguments:
#
#        timestep    (input)     timestep of integration (days)
#
#        period      (input)     follow motion for this long
#                                        (days)
#
#        output_file (input)     name of file into which we 
#                                     write the output data
#
# MWR 4/29/2002
#
# Modified so that we use the NEW velocity to calculate new position
#   in the main loop.  It is _much_ more accurate than using the
#   old velocity.  Thanks to Melanie Day for pointing this out.
#
# MWR 5/4/2003
#    
#
#

use POSIX;

# print values every this number of days
$printstep = 20;

$timestep = 1;
$period = 1000;
$output_file = "orbit.out";


# set to 1 to see diagnostic messsages
$debug = 1;

# sanity check of input values
if ($timestep <= 0) {
   die("timestep must be greater than zero");
}
if ($period <= 0) {
   die("period must be greater than zero");
}


# convert timestep from days to seconds for internal use
$timestep_sec = $timestep*86400;
# number of bodies to follow
$nbody = 1;

# properties of each body, including initial position and 
# velocity.  Units are kg, m, m/s.
#
# Earth
#$mass[1] = 5.98e24;
#$xpos[1] = 1.50e11;
#$e = 0.0167;
#
# Mars
$mass[1] = 0.642e24;
$xpos[1] = 1.50e11*1.5237;
$e = 0.0934;
$xpos[1] = $xpos[1]*(1 + $e);
#xpos(1) = 0.75e11;
$ypos[1] = 0.0;
$xvel[1] = 0.0;
$yvel[1] = 29747;
$factor = sqrt( 1/(1+$e) );
$factor2 = sqrt( 2 * ( ((1/(1-$e)) - (1/(1+$e))) / ( (((1+$e)*(1+$e))/((1-$e)*(1-$e))) - 1 ) ) );
# Earth
#$yvel[1] = 29747*$factor2;
# Mars
$yvel[1] = 29747*(1.0/sqrt(1.5237))*$factor2;
if ($debug > 1) {
  printf " factor is %9.4f  factor2 is %9.4f  yvel is %9.1f \n", 
     $factor, $factor2, $yvel[1];
}

# open file for output
open(OUTFILE, ">$output_file") || die("can't open output file $output_file");


# calculate initial total mechanical energy
$initial_energy = tot_energy($mass[1], $xpos[1], $ypos[1], $xvel[1], $yvel[1]);

# prepare to enter main loop
$time = 0;

# main loop
while ($time < $period) {
   
   for ($i = 1; $i < 2; $i++) {
	   # print out current values
      if ($debug > 0) {
         printf "%12.4f %12.7e %12.7e   %12.7e %12.7e \n",
               $time, $xpos[$i], $ypos[$i], $xvel[$i], $yvel[$i];
      }
      if ($time % $printstep == 0) {
        printf OUTFILE "%12.4f %12.7e %12.7e   %12.7e %12.7e \n", 
               $time, $xpos[$i], $ypos[$i], $xvel[$i], $yvel[$i];
      }
      
      #% calculate new values, using simple Euler's method
      #% consider only forces between Sun and body 
      ($xacc, $yacc) = calc_accel($mass[$i], $xpos[$i], $ypos[$i]);   
      $new_xvel[$i] = $xvel[$i] + $timestep_sec*$xacc;
      $new_yvel[$i] = $yvel[$i] + $timestep_sec*$yacc;
      # it works MUCH better to use new velocity here ...
      $new_xpos[$i] = $xpos[$i] + $timestep_sec*$new_xvel[$i];
      $new_ypos[$i] = $ypos[$i] + $timestep_sec*$new_yvel[$i];
      
      #% check total mechanical energy
      $final_energy = tot_energy($mass[1], $xpos[1], $ypos[1], $xvel[1], $yvel[1]);
      
   }
   
   #% update quantities for next time through loop
   for ($i = 1; $i < 2; $i++) {
      $xpos[$i] = $new_xpos[$i];
      $ypos[$i] = $new_ypos[$i];
      $xvel[$i] = $new_xvel[$i];
      $yvel[$i] = $new_yvel[$i];
   }
   $time = $time + $timestep;
   
}
# end of main loop

# check total mechanical energy
$final_energy = tot_energy($mass[1], $xpos[1], $ypos[1], $xvel[1], $yvel[1]);
$frac_change = ($final_energy - $initial_energy)/abs($initial_energy);
printf "energy: initial %12.5e  final %12.5e  frac_change %12.5e \n",
       $initial_energy, $final_energy, $frac_change;


# close up files and clean up
close(OUTFILE);

# end of main program

exit(0);


#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# calculate the acceleration of a body
# simple idea: include only gravitational force due to Sun
#function [xacc, yacc] = calc_accel(mass, xpos, ypos)
#
# Arguments:
#      mass      (input)        mass of body (kg)
#      
#      xpos      (input)        x-position of body (m)
#
#      ypos      (input)        y-position of body (m)
#
#      xacc     (output)        x-acceleration of body (m/s^2)
#
#      yacc     (output)        y-acceleration of body (m/s^2)
#

sub calc_accel {

  my($mass, $xpos, $ypos);

  $mass = $_[0];
  $xpos = $_[1];
  $ypos = $_[2];

# some physical constants
$G = 6.67e-11;
# the central body (the Sun) is fixed in place
# units are kg, m
$solar_mass = 1.99e30;
$solar_xpos = 0.0;
$solar_ypos = 0.0;

# distance between body and Sun
$r = sqrt($xpos*$xpos + $ypos*$ypos);
$accel = ($G*$solar_mass)/($r*$r);
$xacc = $accel*(-$xpos/$r);
$yacc = $accel*(-$ypos/$r);

return($xacc, $yacc);

}
# end of function 'calc_accel'


#%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
# calculate the total mechanical energy of the given body
# function energy = tot_energy(mass, xpos, ypos, xvel, yvel)
#
# Arguments:
#      mass      (input)        mass of body (kg)
#      
#      xpos      (input)        x-position of body (m)
#
#      ypos      (input)        y-position of body (m)
#
#      xvel      (input)        x-velocity of body (m/s)
#
#      yvel      (input)        y-velocity of body (m/s)
#
#      energy   (output)        total mechanical (KE+GPE) of body (J)

sub tot_energy {

  my($mass, $xpos, $ypos, $xvel, $yvel);

  $mass = $_[0];
  $xpos = $_[1];
  $ypos = $_[2];
  $xvel = $_[3];
  $yvel = $_[4];

# some physical constants
$G = 6.67e-11;
# the central body (the Sun) is fixed in place
# units are kg, m
$solar_mass = 1.99e30;
$solar_xpos = 0.0;
$solar_ypos = 0.0;

# distance between body and Sun
$r = sqrt($xpos*$xpos + $ypos*$ypos);

# gravitational potential energy
$gpe = -($G*$solar_mass*$mass)/$r;

# kinetic energy
$ke = 0.5*$mass*($xvel*$xvel + $yvel*$yvel);

$energy = $gpe + $ke;

  if ($debug > 0) {
    printf "    gpe %12.6e  ke %12.6e  tote %12.6e \n",
           $gpe, $ke, $energy;
  }

return($energy);

}
# end of function 'tot_energy'

