#!/usr/bin/perl
#
# Create a elliptical orbit, print out (x, y) positions along it
#    and (x, y) location of each focus.
#
# Then pick an axis which makes an angle "bigphi" relative to the
#    x-axis = semi-major axis of ellipse.
#
# Calculate positions (x', y') of ellipse in a coord system which is
#    rotated by "bigphi" relative to (x, y) coords.
#
# Then, project the (x', y') positions onto a plane which is tilted
#    by angle "phi" relative to the original plane.  Print out
#    the projected positions (x'', y'').
#
# MWR 1/15/2004
# 

use POSIX;

$debug = 1;

# eccentricity of original ellipse
$e = 0.8;
# semi-major axis of original ellipse 
$a = 10.0;
# number of pieces around ellipse for which we calc (x,y) pos
$npos = 100;
$pi = 3.14159265359;
$dtheta = (2.0*$pi)/$npos;
# angle of coord system relative to true semimajor axis (radians)
$bigphi = 80.0*($pi/180.0);
#
# You should set only one of the two tilts below to a non-zero value.
#
# angle by which orbital plane is tilted ALONG MAJOR AXIS relative to 
#    _normal_ to line of sight (i.e. 0 = face-on, pi/2 = edge-on)
$phi = 60*($pi/180.0);
# angle by which orbital plane is tilted ALONG MINOR AXIS relative to 
#    _normal_ to line of sight (i.e. 0 = face-on, pi/2 = edge-on)
$minor_phi = 0*($pi/180.0);


# Step 1: create (x, y) positions
#
# this is the distance between center and focus, along major axis
$foc_dist = $a*$e;
# as we calculate the (x, y) position, we add the foc_dist to 
#   each x-value in order to center the ellipse on the origin.
$theta = 0.0;
for ($i = 0; $i < $npos; $i++) {
  $r = ($a*(1.0 - $e*$e))/(1.0 + $e*cos($theta));
  $x[$i] = $r*cos($theta) + $foc_dist;
  $y[$i] = $r*sin($theta);
  $theta += $dtheta;
}
# and now add the position of the two foci
$x[$i++] = $foc_dist;
$y[$i] = 0.0;
$x[$i++] = 0.0 - $foc_dist;
$y[$i] = 0.0;


# Step 2: Calculate (x', y'), the same ellipse in a coord system
#   which is rotated by "bigphi" relative to the semimajor = x axis.
#
for ($i = 0; $i < $npos; $i++) {
  $r = sqrt($x[$i]*$x[$i] + $y[$i]*$y[$i]);
  $tt = atan2($y[$i], $x[$i]);
  $tt += $bigphi;
  $xp[$i] = $r*cos($tt);
  $yp[$i] = $r*sin($tt);
}
# and now rotate the position of the two foci
for ( ; $i < $npos + 2; $i++) {
  $r = sqrt($x[$i]*$x[$i] + $y[$i]*$y[$i]);
  $tt = atan2($y[$i], $x[$i]);
  $tt += $bigphi;
  $xp[$i] = $r*cos($tt);
  $yp[$i] = $r*sin($tt);
}


# Step 3: Calculate (x'', y''), the projected coords onto a plane
#   which is tilted around the x' axis by angle "phi"
#
for ($i = 0; $i < $npos; $i++) {
  $xpp[$i] = $xp[$i];
  $ypp[$i] = $yp[$i]*cos($phi);
}
# and now rotate the position of the two foci
for ( ; $i < $npos + 2; $i++) {
  $xpp[$i] = $xp[$i];
  $ypp[$i] = $yp[$i]*cos($phi);
}


# Step 4: Calculate (x''', y'''), the projected coords onto a plane
#   which is tilted around the y' axis by angle "minor_phi"
#   These start with (x', y'), not with (x'', y'').  So these
#   values include only a tilt around minor axis, not a tilt around major axis.
#
for ($i = 0; $i < $npos; $i++) {
  $xppp[$i] = $xp[$i]*cos($minor_phi);
  $yppp[$i] = $yp[$i];
}
# and now rotate the position of the two foci
for ( ; $i < $npos + 2; $i++) {
  $xppp[$i] = $xp[$i]*cos($minor_phi);
  $yppp[$i] = $yp[$i];
}



# Step 5: Calculate (new_x, new_y), in which we rotate the
#   rotated, tilted ellipse back into the original (x, y) coord system.
#   This un-does step 2, in essence.
#
for ($i = 0; $i < $npos; $i++) {
  $r = sqrt($xpp[$i]*$xpp[$i] + $ypp[$i]*$ypp[$i]);
  $tt = atan2($ypp[$i], $xpp[$i]);
  $tt -= $bigphi;
  $new_x[$i] = $r*cos($tt);
  $new_y[$i] = $r*sin($tt);
}
# and now rotate the position of the two foci
for ( ; $i < $npos + 2; $i++) {
  $r = sqrt($xpp[$i]*$xpp[$i] + $ypp[$i]*$ypp[$i]);
  $tt = atan2($ypp[$i], $xpp[$i]);
  $tt -= $bigphi;
  $new_x[$i] = $r*cos($tt);
  $new_y[$i] = $r*sin($tt);
}





# Final step: print out results, with (x, y), then (x', y'), then (x'', y''),
#   then (new_x, new_y),
#   all on the same line
#
for ($i = 0; $i < $npos + 2; $i++) {
  printf "%9.4e %9.4e   %9.4e %9.4e   %9.4e %9.4e  %9.4e %9.4e  %9.4e %9.4e\n", 
         $x[$i], $y[$i], $xp[$i], $yp[$i], $xpp[$i], $ypp[$i],
         $xppp[$i], $yppp[$i], $new_x[$i], $new_y[$i];
}


exit 0;
