#!/usr/bin/perl
#
# Compute the (x,y) coordinates of both a circle
#   and an ellipse.
#
# MWR 9/9/2018
use POSIX;

$debug = 1;

# we will compute the coordinates of a circle centered at the
#   origin, with radius 1, 
#   and an ellipse centered on the origin, with radius 1 but
#   eccentricity "e"
#
# This eccentricty corresponds to that of Saturn's orbit
$e = 0.057;
$e = 0.50;
# Radius of circle
$radius = 1.0;
# Parameters of the ellipse
$a = $radius;
$b = $a*sqrt(1.0 - $e*$e);
# tolerance for our check on the shape of the ellipse
$tolerance = 1.0E-6;

$xmin = -1.0*$radius;
$xmax = 1.0*$radius;
$dx = 0.001*$radius;

for ($x = $xmin; $x <= $xmax; $x += $dx) {

  # first, the circle.  We compute both pos and neg y values
  $circ_pos_y = $radius*sqrt(1.0 - $x*$x);
  $circ_neg_y = 0.0 - $circ_pos_y;

  # now, the ellipse.  
  $ell_pos_y = $b*sqrt(1.0 - (($x*$x)/($a*$a)));
  $ell_neg_y = 0.0 - $ell_pos_y;

  # sanity check
  $val = ($x*$x)/($a*$a) + ($ell_pos_y*$ell_pos_y)/($b*$b);
  $diff = $val - 1.0;
  if (abs($diff) > $tolerance) {
    printf "error: x %10.5f y %10.5f  val %10.4e \n", $x, $y, $val;
    exit(1);
  }

  printf " %10.5f %10.5f %10.5f   %10.5f %10.5f \n", 
      $x, $circ_pos_y, $circ_neg_y, $ell_pos_y, $ell_neg_y;


}

# now do the special case of x = $xmax, since we might not 
#   have reached $xmax exactly in the loop
$x = $xmax;
  # first, the circle.  We compute both pos and neg y values
  $circ_pos_y = $radius*sqrt(1.0 - $x*$x);
  $circ_neg_y = 0.0 - $circ_pos_y;

  # now, the ellipse.  
  $ell_pos_y = $b*sqrt(1.0 - (($x*$x)/($a*$a)));
  $ell_neg_y = 0.0 - $ell_pos_y;

  # sanity check
  $val = ($x*$x)/($a*$a) + ($ell_pos_y*$ell_pos_y)/($b*$b);
  $diff = $val - 1.0;
  if (abs($diff) > $tolerance) {
    printf "error: x %10.5f y %10.5f  val %10.4e \n", $x, $y, $val;
    exit(1);
  }

  printf " %10.5f %10.5f %10.5f   %10.5f %10.5f \n", 
      $x, $circ_pos_y, $circ_neg_y, $ell_pos_y, $ell_neg_y;







exit 0;

