#!/usr/bin/perl
#
# Calculate the potential energy along the line joining big and small
#   mass in a rotating system.
#
# 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;
# 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.00;
# 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 along a line between sun and planet.  Zero is at Sun,
#    negative is away from planet, positive is towards planet.
$start_x = -20.0;
$end_x = 20.0;
$dx = 0.01*$L;
for ($x = $start_x; $x <= $end_x; $x += $dx) {

  # GPE due to sun
  $dist = abs($x) + $soft;
  $gpe_1 = -($G*$m1)/($dist);

  # GPE due to planet
  $dist = abs(($x - $L)) + $soft;
  $gpe_2 = -($G*$m2)/($dist);

  # centrifugal potential energy
  #   we have to use the distance from the center of mass!
  $cpe_dist = ($x - $d1);
  $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  tot %9.4e \n",
          $x, $gpe_1, $gpe_2, $cpe, $tot;
  }
  



}

exit 0;

