#!/usr/bin/perl
#
# Calculate the potential energy all over the plane in a 
#   rotating system.  Spit out PE with both (x, y) and (r, theta) coords.
#
# MWR 1/22/2004


$debug = 1;


# gaussian grav constant, for solar masses, AU and days
$k = 0.017202;
# grav constant, for mass in solar and dist in AU
$G = 8.95E8;
# softening radius for grav calcs, in units of AU
$soft = 1.0e-6;
#
$pi = 3.14159;
$DEGTORAD = ($pi/180.0);
# one AU (in meters)
$au = 1.49E11;
# one solar mass (in kg)
$solar_mass = 2.0E30;


# mass of big mass, in solar units
$m1 = 1.0;
# mass of smaller mass, in solar units
$m2 = 0.10;
# ratio of masses, or a really small number if m2 = 0
if ($m2 > 0) {
  $alpha = $m2/$m1;
}
else {
  $alpha = 1e-6;
}
# distance between big and small mass, in AU
$L = 5.0;
# period of planet's orbit, in years
$period = sqrt((4.0*$pi*$pi*$L*$L*$L)/(($m1 + $m2)*$k*$k));
$period = $period/365.25;
if ($debug > 1) {
  printf "period is %10.4e \n", $period;
}
# angular frequency, in radians per second
$circum = 2*$pi*$L*$au;
$time_in_sec = $period*86400*365.25;
$omega = (2.0*$pi)/($time_in_sec);
#$omega = 0.0;

# distance between each big mass and the center of mass
if ($m2 == 0) {
  $d1 = 0.0;
  $d2 = $L;
} else {
  $d1 = $L/(1.0 + (1.0/$alpha));
  $d2 = $L/(1.0 + $alpha);
}


# we walk in (x, y), calculating (r, theta) as we go

# we walk in x: positive x points from big mass to small mass
#$start_x = -2.0*$L;
#$end_x = 2.0*$L;
$start_x = 0.0;
$end_x = 5.0;
$d_x = ($end_x - $start_x)/200.0;
for ($pos_x = $start_x; $pos_x <= $end_x; $pos_x += $d_x) {

  # we walk in y
  #$start_y = -2.0*$L;
  #$end_y = 2.0*$L;
  $start_y = 1.0;
  $end_y = 6.0;
  $d_y = ($end_y - $start_y)/200.0;
  for ($pos_y = $start_y; $pos_y <= $end_y; $pos_y += $d_y) {

   # calculate (radius, theta) for this (x, y)
   $radius = sqrt($pos_x*$pos_x + $pos_y*$pos_y);
   $theta = atan2($pos_y, $pos_x);
   
     # GPE due to sun
     $dist = abs($radius) + $soft;
     $gpe_1 = -($G*$m1)/($dist);
   
     # GPE due to planet
     #   we do this in cartesian coords, planet is at (L, 0)
     $dx = $L - $pos_x;
     $dy = 0.0 - $pos_y;
     $dist = sqrt($dx*$dx + $dy*$dy + $soft*$soft);
     $gpe_2 = -($G*$m2)/($dist);
   
     # centrifugal potential energy
     #   we have to use the distance from the center of mass!
     $testmass_x = $pos_x;
     $testmass_y = $pos_y;
     $x_dist = $testmass_x - $d1;
     $y_dist = $testmass_y - 0.0;
     $cpe_dist = sqrt($x_dist*$x_dist + $y_dist*$y_dist);
     $cpe = -0.5*($cpe_dist*$au)*($cpe_dist*$au)*$omega*$omega;
   
     # total potential energy
     $tot = $gpe_1 + $gpe_2 + $cpe;
   
     if (1 == 1) {
       # print energies in linear units
       printf " %9.4e %9.4e %9.4e %9.4e  %9.4e  %9.4e  %9.4e  tot %15.10e \n",
             $pos_x, $pos_y, $radius, $theta, 
             $gpe_1, $gpe_2, $cpe, $tot;

     } else {

       # print energies in logarithmic units
       $log10 = log(10.0);
       printf " %9.4e %9.4e %9.4e %9.4e  %9.4e  %9.4e  %9.4e  tot %15.10e \n",
             $pos_x, $pos_y, $radius, $theta, 
             log(-$gpe_1)/$log10, log(-$gpe_2)/$log10, 
             log(-$cpe)/$log10, log(-$tot)/$log10;
     }
     
  }

  # print an empty line to make GNUplot happy
  printf "\n";
   
}

exit 0;

