#!/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.04;
# 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 > 0) {
  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 radius outwards from the big mass.
$start_rad = 0.1;
$end_rad = 2.0*$L;
$d_rad = ($end_rad - $start_rad)/20.0;
for ($radius = $start_rad; $radius <= $end_rad; $radius += $d_rad) {

   
   # we walk in angle around the Sun.  Zero is direction of planet.
   #
   $start_theta = 0.0;
   $end_theta = 360.0;
   $dtheta = 2;
   for ($theta = $start_theta; $theta <= $end_theta; $theta += $dtheta) {
   
     # GPE due to sun
     $dist = abs($radius) + $soft;
     $gpe_1 = -($G*$m1)/($dist);
   
     # GPE due to planet
     $alpha = 0.5*(180.0 - $theta);
     $alpha *= $DEGTORAD;
     $dist = abs(2.0*$radius*cos($alpha)) + $soft;
     $gpe_2 = -($G*$m2)/($dist);
   
     # centrifugal potential energy
     #   we have to use the distance from the center of mass!
     $testmass_x = $radius*cos($theta*$DEGTORAD);
     $testmass_y = $radius*sin($theta*$DEGTORAD);
     $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 ($debug > 0) {
       printf " %9.4e %9.4e %9.4e %9.4e  %9.4e  %9.4e  %9.4e  tot %15.10e \n",
             $testmass_x, $testmass_y, $radius, $theta, 
             $gpe_1, $gpe_2, $cpe, $tot;
     }
     
  }
   
}

exit 0;

