#!/usr/bin/perl
#
# Given the output file of the 'nbody' program, we read through
#   it to analyze the results.  The file has a format like this:
#
#       timeN    body0 mass  x y z   vx vy vz
#       timeN    body1 mass  x y z   vx vy vz
#       timeN    body2 mass  x y z   vx vy vz
#         ...
#       timeN+1  body0 mass  x y z   vx vy vz
#       timeN+1  body1 mass  x y z   vx vy vz
#       timeN+1  body2 mass  x y z   vx vy vz
#
# Our job is to 
#          read in all the data for a given time
#          compute KE, GPE, totE, angmom for that time
#          print out one line per time with the information
#
# We send output to stdout
#
# Usage:  analyze.pl   nbody.out   [every=N]
#
# MWR 12/24/2008
#
# Added "every=N" argument, so that we only print output
#    every N lines, if desired
# Also, now print out one line at end with change in total energy 
#    from start to end.
# 
# MWR 1/1/2009

$debug = 0;

$G = 6.67E-11;

# print output once for every "every'th" line of input
$every = 1;

# we will save initial and final total energy, and report on the
#   change at the end of the analysis
$initial_tote = 0.0;
$final_tote = 0.0;


if ($#ARGV < 0) {
  printf STDERR "usage: analyze.pl nbody.out [every=] \n";
  exit(1);
}

$infile = $ARGV[0];
if ($debug > 0) {
  printf " reading input from file ..%s.. \n", $infile;
}
if ($#ARGV > 0) {
  if ($ARGV[1] =~ /every=/) {
    $every = $ARGV[1];
    $every =~ s/every=//;
    if ($debug > 0) {
      printf "  every is ..%d.. \n", $every;
    }
  }
} 

$cur_time = -1.0;
$nbody = 0;
$num_lines = 0;
open(INFILE, "$infile") || die("can't open file $infile");
while (<INFILE>) {
  $line = " " . $_;
  chomp($line);
  @words = split(/\s+/, $line);
  if ($#words != 9) {
    printf STDERR " bad line doesn't have 9 words ..%s.. \n", $line;
  }
  $this_time = $words[1];
  $this_index = $words[2];
  $this_mass = $words[3];
  $this_px = $words[4];
  $this_py = $words[5];
  $this_pz = $words[6];
  $this_vx = $words[7];
  $this_vy = $words[8];
  $this_vz = $words[9];

  if ($this_time == $cur_time) {
    # just read information on another body 
    $mass_array[$nbody] = $this_mass;
    $px_array[$nbody] = $this_px;
    $py_array[$nbody] = $this_py;
    $pz_array[$nbody] = $this_pz;
    $vx_array[$nbody] = $this_vx;
    $vy_array[$nbody] = $this_vy;
    $vz_array[$nbody] = $this_vz;
    $nbody++;
  }
  else {
	  # ah, this represents a new timestep.  We should go back
     #   at this point and compute quantities based on
     #   the information for the previous timestep
     if ($cur_time >= 0) {

        if ($debug > 0) {
          printf " total of %5d bodies for time %12.5e \n",
              $nbody, $cur_time;
        }

        # compute stuff here
        #  
        # total kinetic energy
        $ke = 0.0;
        for ($i = 0; $i < $nbody; $i++) {
          $v_squared = ($vx_array[$i]*$vx_array[$i] +
                        $vy_array[$i]*$vy_array[$i] +
                        $vz_array[$i]*$vz_array[$i]);
          $ke += 0.5*$mass_array[$i]*$v_squared;
        }

        # total potential energy
        $gpe = 0.0;
        for ($i = 0; $i < $nbody - 1; $i++) {
          for ($j = $i + 1; $j < $nbody; $j++) {
            $dx = $px_array[$i] - $px_array[$j];
            $dy = $py_array[$i] - $py_array[$j];
            $dz = $pz_array[$i] - $pz_array[$j];
            $dist = sqrt($dx*$dx + $dy*$dy + $dz*$dz);
            if ($dist <= 0) {
              printf STDERR "distance <= 0 ?! time %12.5e bodies %d %d \n",
                      $cur_time, $i, $j;
              exit(1);
            }
            $this_gpe = 0.0 - $G*$mass_array[$i]*$mass_array[$j];
            $this_gpe /= $dist;
          
            $gpe += $this_gpe;
          }
        }
            
        # calc total energy
        $tot_e = $ke + $gpe;

        # total angular momentum around center of mass
        #   first, we need center of mass location
        $com_x = 0.0;
        $com_y = 0.0;
        $com_z = 0.0;
        $total_mass = 0.0;
        for ($i = 0; $i < $nbody; $i++) {
          $com_x += $mass_array[$i]*$px_array[$i];
          $com_y += $mass_array[$i]*$py_array[$i];
          $com_z += $mass_array[$i]*$pz_array[$i];
          $total_mass += $mass_array[$i];
        }
        $com_x /= $total_mass;
        $com_y /= $total_mass;
        $com_z /= $total_mass;
        #
        # now we can compute vector angular momentum of 
        #   each body around center of mass
        $ang_mom_x = 0.0;
        $ang_mom_y = 0.0;
        $ang_mom_z = 0.0;
        for ($i = 0; $i < $nbody; $i++) {
          $this_x = $py_array[$i]*$vz_array[$i] - $pz_array[$i]*$vy_array[$i];
          $this_y = $pz_array[$i]*$vx_array[$i] - $px_array[$i]*$vz_array[$i];
          $this_z = $px_array[$i]*$vy_array[$i] - $py_array[$i]*$vx_array[$i];
          $ang_mom_x += $this_x * $mass_array[$i];
          $ang_mom_y += $this_y * $mass_array[$i];
          $ang_mom_z += $this_z * $mass_array[$i];
        }
        # find magnitude of total angular momentum
        $mag_total_ang_mom = sqrt($ang_mom_x*$ang_mom_x +
                                  $ang_mom_y*$ang_mom_y + 
                                  $ang_mom_z*$ang_mom_z);
            
        # calc total energy
        $tot_e = $ke + $gpe;

        # save initial total energy, if this is the first line
        if ($num_lines == 0) {
          $initial_tote = $tot_e;
        }

        # print summary line to stdout (if this is one of
        #    the every'th lines)
        if ($num_lines % $every == 0) {
          printf "%15.9e  %15.9e %15.9e  %15.9e  %15.9e \n", 
              $cur_time, $ke, $gpe, $tot_e, $mag_total_ang_mom;
        }
        $num_lines++;

        # reset counter of number of bodies
        $nbody = 0;

     }

     # and read in the data for the first body in new timestep
     $cur_time = $this_time;
     $mass_array[$nbody] = $this_mass;
     $px_array[$nbody] = $this_px;
     $py_array[$nbody] = $this_py;
     $pz_array[$nbody] = $this_pz;
     $vx_array[$nbody] = $this_vx;
     $vy_array[$nbody] = $this_vy;
     $vz_array[$nbody] = $this_vz;

     $nbody++;
  } 


}
close(INFILE);


# Now, when we get here, we won't have computed and printed out
#   quantities for the final timestep.  So, we'll do that now.
#   We assume that the positions and velocities of all bodies
#   for the final time are sitting in the arrays

        if ($debug > 0) {
          printf " total of %5d bodies for time %12.5e \n",
              $nbody, $cur_time;
        }

        # compute stuff here
        #  
        # total kinetic energy
        $ke = 0.0;
        for ($i = 0; $i < $nbody; $i++) {
          $v_squared = ($vx_array[$i]*$vx_array[$i] +
                        $vy_array[$i]*$vy_array[$i] +
                        $vz_array[$i]*$vz_array[$i]);
          $ke += 0.5*$mass_array[$i]*$v_squared;
        }

        # total potential energy
        $gpe = 0.0;
        for ($i = 0; $i < $nbody - 1; $i++) {
          for ($j = $i + 1; $j < $nbody; $j++) {
            $dx = $px_array[$i] - $px_array[$j];
            $dy = $py_array[$i] - $py_array[$j];
            $dz = $pz_array[$i] - $pz_array[$j];
            $dist = sqrt($dx*$dx + $dy*$dy + $dz*$dz);
            if ($dist <= 0) {
              printf STDERR "distance <= 0 ?! time %12.5e bodies %d %d \n",
                      $cur_time, $i, $j;
              exit(1);
            }
            $this_gpe = 0.0 - $G*$mass_array[$i]*$mass_array[$j];
            $this_gpe /= $dist;
          
            $gpe += $this_gpe;
          }
        }
            
        # calc total energy
        $tot_e = $ke + $gpe;

        # total angular momentum around center of mass
        #   first, we need center of mass location
        $com_x = 0.0;
        $com_y = 0.0;
        $com_z = 0.0;
        $total_mass = 0.0;
        for ($i = 0; $i < $nbody; $i++) {
          $com_x += $mass_array[$i]*$px_array[$i];
          $com_y += $mass_array[$i]*$py_array[$i];
          $com_z += $mass_array[$i]*$pz_array[$i];
          $total_mass += $mass_array[$i];
        }
        $com_x /= $total_mass;
        $com_y /= $total_mass;
        $com_z /= $total_mass;
        #
        # now we can compute vector angular momentum of 
        #   each body around center of mass
        $ang_mom_x = 0.0;
        $ang_mom_y = 0.0;
        $ang_mom_z = 0.0;
        for ($i = 0; $i < $nbody; $i++) {
          $this_x = $py_array[$i]*$vz_array[$i] - $pz_array[$i]*$vy_array[$i];
          $this_y = $pz_array[$i]*$vx_array[$i] - $px_array[$i]*$vz_array[$i];
          $this_z = $px_array[$i]*$vy_array[$i] - $py_array[$i]*$vx_array[$i];
          $ang_mom_x += $this_x * $mass_array[$i];
          $ang_mom_y += $this_y * $mass_array[$i];
          $ang_mom_z += $this_z * $mass_array[$i];
        }
        # find magnitude of total angular momentum
        $mag_total_ang_mom = sqrt($ang_mom_x*$ang_mom_x +
                                  $ang_mom_y*$ang_mom_y + 
                                  $ang_mom_z*$ang_mom_z);
            
        # calc total energy
        $tot_e = $ke + $gpe;

        # save final total energy, since this is the end
        $final_tote = $tot_e;

        # print summary line to stdout, even if this is not
        #    one of the every'th lines
        printf "%15.9e  %15.9e %15.9e  %15.9e  %15.9e \n", 
            $cur_time, $ke, $gpe, $tot_e, $mag_total_ang_mom;

 
        # compute change in total energy, and print it out
        #   in a line which starts with comment char
        if ($initial_tote == 0.0) {
          $absolute_change = $final_tote - $initial_tote;
          printf "# total_energy init %9.4e  final %9.4e  abs  change %9.4e \n",
                     $initial_tote, $final_tote, $absolute_change;
        }
        else {
          $fractional_change = ($final_tote - $initial_tote) / ($initial_tote);
          printf "# total_energy init %9.4e  final %9.4e  frac change %9.4e \n",
                     $initial_tote, $final_tote, $fractional_change;
        }





exit 0;


