#!/usr/bin/perl
#
# Graphs to illustrate astrometry
#

$cmdfile = "./makeplot.gnu";



#####################################################################
# Create a graph showing the "supjup" (Super Jupiter) system
#    face-on: both objects, in initial positions only
#

  $infile = "nbody.supjup.out";
  $outfile = "supjup_a.ps";

  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set nokey \n";
  printf CMDFILE "set xlabel 'X position (m)' \n";
  printf CMDFILE "set ylabel 'Y position (m)' \n";
  
  $cmd = "plot [-1.5e11:1.5e11][-1.5e11:1.5e11] '$infile' ";
  $cmd .= "      index 0 using 4:5   ";
  $cmd .= "      with points pt 6 ps 2 lc 1 lw 5 t '' ";
  $cmd .= " , ";
  $cmd .= "                                     '$infile' ";
  $cmd .= "      index 1 using 4:5   ";
  $cmd .= "      with points pt 6 ps 2 lc 3 lw 5 t '' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;



#####################################################################
# Create a graph showing the "supjup" (Super Jupiter) system
#    face-on: both objects, initial positions and motion over
#    a small fraction of one revolution
#

  $infile = "nbody.supjup.out";
  $outfile = "supjup_b.ps";

  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set nokey \n";
  printf CMDFILE "set xlabel 'X position (m)' \n";
  printf CMDFILE "set ylabel 'Y position (m)' \n";
  
  $cmd = "plot [-1.5e11:1.5e11][-1.5e11:1.5e11] '$infile' ";
  $cmd .= "      index 0 using 4:5   ";
  $cmd .= "      with points pt 6 ps 2 lc 1 lw 5 t '' ";
  $cmd .= " , ";
  $cmd .= "                                     '$infile' ";
  $cmd .= "      index 1 using 4:5   ";
  $cmd .= "      with points pt 6 ps 2 lc 3 lw 5 t '' ";
  $cmd .= " , ";
  $cmd .= "                                     '$infile' ";
  $cmd .= "      index 2 every 4::0::70 using 4:5   ";
  $cmd .= "      with points pt 6 ps 0.1 lc 1 lw 2 t '' ";
  $cmd .= " , ";
  $cmd .= "                                     '$infile' ";
  $cmd .= "      index 2 every 4::1::70 using 4:5   ";
  $cmd .= "      with points pt 6 ps 0.1 lc 3 lw 2 t '' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;



#####################################################################
# Create a graph showing the "supjup" (Super Jupiter) system
#    face-on: both objects, initial positions and motion over
#    the entirety of one revolution
#

  $infile = "nbody.supjup.out";
  $outfile = "supjup_c.ps";

  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set nokey \n";
  printf CMDFILE "set xlabel 'X position (m)' \n";
  printf CMDFILE "set ylabel 'Y position (m)' \n";
  
  $cmd = "plot [-1.5e11:1.5e11][-1.5e11:1.5e11] '$infile' ";
  $cmd .= "      index 0 using 4:5   ";
  $cmd .= "      with points pt 6 ps 2 lc 1 lw 5 t '' ";
  $cmd .= " , ";
  $cmd .= "                                     '$infile' ";
  $cmd .= "      index 1 using 4:5   ";
  $cmd .= "      with points pt 6 ps 2 lc 3 lw 5 t '' ";
  $cmd .= " , ";
  $cmd .= "                                     '$infile' ";
  $cmd .= "      index 2 every 4::0 using 4:5   ";
  $cmd .= "      with points pt 6 ps 0.1 lc 1 lw 2 t '' ";
  $cmd .= " , ";
  $cmd .= "                                     '$infile' ";
  $cmd .= "      index 2 every 4::1 using 4:5   ";
  $cmd .= "      with points pt 6 ps 0.1 lc 3 lw 2 t '' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;




#####################################################################
# Create a graph showing the "supjup" (Super Jupiter) system
#    face-on: both objects, initial positions and motion over
#    the entirety of one revolution
# 
# This version makes a big graph, which we need for zooming in
#    to show center of mass
#

  $infile = "nbody.supjup.out";
  $outfile = "supjup_c1.ps";

  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' size 14,10 \n";
  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set nokey \n";
  printf CMDFILE "set xlabel 'X position (m)' \n";
  printf CMDFILE "set ylabel 'Y position (m)' \n";
  
  $cmd = "plot [-1.5e11:1.5e11][-1.5e11:1.5e11] '$infile' ";
  $cmd .= "      index 0 using 4:5   ";
  $cmd .= "      with points pt 6 ps 2 lc 1 lw 5 t '' ";
  $cmd .= " , ";
  $cmd .= "                                     '$infile' ";
  $cmd .= "      index 1 using 4:5   ";
  $cmd .= "      with points pt 6 ps 2 lc 3 lw 5 t '' ";
  $cmd .= " , ";
  $cmd .= "                                     '$infile' ";
  $cmd .= "      index 2 every 4::0 using 4:5   ";
  $cmd .= "      with points pt 6 ps 0.1 lc 1 lw 2 t '' ";
  $cmd .= " , ";
  $cmd .= "                                     '$infile' ";
  $cmd .= "      index 2 every 4::1 using 4:5   ";
  $cmd .= "      with points pt 6 ps 0.1 lc 3 lw 2 t '' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;



#####################################################################
# Create a graph showing the "supjup" (Super Jupiter) system
#    face-on: only the star, 
#    for the entirety of one revolution
#

  $infile = "nbody.supjup.out";
  $outfile = "supjup_d.ps";

  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set nokey \n";
  printf CMDFILE "set xlabel 'X position (m)' \n";
  printf CMDFILE "set ylabel 'Y position (m)' \n";
  
  $cmd = "plot [-1.5e11:1.5e11][-1.5e11:1.5e11] '$infile' ";
  $cmd .= "      index 0 using 4:5   ";
  $cmd .= "      with points pt 6 ps 2 lc 1 lw 5 t '' ";
  $cmd .= " , ";
#  $cmd .= "                                     '$infile' ";
#  $cmd .= "      index 1 using 4:5   ";
#  $cmd .= "      with points pt 6 ps 2 lc 3 lw 5 t '' ";
#  $cmd .= " , ";
  $cmd .= "                                     '$infile' ";
  $cmd .= "      index 2 every 4::0 using 4:5   ";
  $cmd .= "      with points pt 6 ps 0.1 lc 1 lw 2 t '' ";
#  $cmd .= " , ";
#  $cmd .= "                                     '$infile' ";
#  $cmd .= "      index 2 every 4::1 using 4:5   ";
#  $cmd .= "      with points pt 6 ps 0.1 lc 3 lw 2 t '' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;





#####################################################################
# Create a graph showing the "supjup" (Super Jupiter) system
#    face-on: only the star, 
#    for the entirety of one revolution
# 
#    This version makes a 'big' version
#

  $infile = "nbody.supjup.out";
  $outfile = "supjup_d1.ps";

  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' size 14,10  \n";
  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set nokey \n";
  printf CMDFILE "set xlabel 'X position (m)' \n";
  printf CMDFILE "set ylabel 'Y position (m)' \n";
  
  $cmd = "plot [-1.5e11:1.5e11][-1.5e11:1.5e11] '$infile' ";
  $cmd .= "      index 0 using 4:5   ";
  $cmd .= "      with points pt 6 ps 2 lc 1 lw 5 t '' ";
  $cmd .= " , ";
#  $cmd .= "                                     '$infile' ";
#  $cmd .= "      index 1 using 4:5   ";
#  $cmd .= "      with points pt 6 ps 2 lc 3 lw 5 t '' ";
#  $cmd .= " , ";
  $cmd .= "                                     '$infile' ";
  $cmd .= "      index 2 every 4::0 using 4:5   ";
  $cmd .= "      with points pt 6 ps 0.1 lc 1 lw 2 t '' ";
#  $cmd .= " , ";
#  $cmd .= "                                     '$infile' ";
#  $cmd .= "      index 2 every 4::1 using 4:5   ";
#  $cmd .= "      with points pt 6 ps 0.1 lc 3 lw 2 t '' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;



#####################################################################
# Create a graph illustrating the nature of the "one-body
#    Kepler problem."  
#

  $outfile = "gamma_func_a.ps";

  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
  # printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set nokey \n";
  printf CMDFILE "set samples 1000 \n";
  printf CMDFILE "set xlabel 'mass ratio {/Symbol g} = (m/M) ' \n";
  printf CMDFILE "set ylabel 'function (1 + 1/{/Symbol g})^3 / (1 + {/Symbol g}) ' \n";
  
  $cmd = "plot [0:1][] ((1 + 1/x)**3)/(1+x) ";
  $cmd .= "      with lines  lc 1 lw 3 t '' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;



#####################################################################
# Create a graph illustrating the nature of the "one-body
#    Kepler problem."  
#    This is log-log version
#

  $outfile = "gamma_func_b.ps";

  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
  # printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set nokey \n";
  printf CMDFILE "set logscale x \n";
  printf CMDFILE "set logscale y \n";
  printf CMDFILE "set samples 1000 \n";
  printf CMDFILE "set xlabel 'mass ratio {/Symbol g} = (m/M) ' \n";
  printf CMDFILE "set ylabel 'function (1 + 1/{/Symbol g})^3 / (1 + {/Symbol g}) ' \n";
  
  $cmd = "plot [1e-5:1][] ((1 + 1/x)**3)/(1+x) ";
  $cmd .= "      with lines  lc 1 lw 3 t '' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;




#####################################################################
# Create a graph illustrating the radial velocities in
#    two systems, if BOTH WERE EDGE-ON.
#

  $outfile = "rv_compare_a.ps";

  $TWOPI = 6.28319;
  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
  # printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set key bot right \n";
#  printf CMDFILE "set logscale x \n";
#  printf CMDFILE "set logscale y \n";
  printf CMDFILE "set samples 1000 \n";
  printf CMDFILE "set title 'Two systems, both edge-on ' \n";
  printf CMDFILE "set xlabel 'Time (years) ' \n";
  printf CMDFILE "set ylabel 'Stellar radial velocity (m/s) ' \n";
  
  $cmd = "plot [0:2.5][-75:75] 29.76*sin($TWOPI*x/1.0)  ";
  $cmd .= "      with lines  lc 3 lw 3 t '1 M_J planet' ";
  $cmd .= " , ";
#  $cmd .= "                     59.45*sin($TWOPI*x/1.0)  ";
#  $cmd .= "      with lines  lc 1 lw 3 t '2 M_J planet' ";
  $cmd .= "     'make_rv.out'  using 1:(\$2*2.0)         ";
  $cmd .= "      with points ps 1.5 pt 6  lc 1 lw 3 t '2 M_J planet' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;



#####################################################################
# Create a graph illustrating the radial velocities in
#    two systems, if THE 2-JUP-MASS SYSTEM IS TILTED.
#

  $outfile = "rv_compare_b.ps";

  $TWOPI = 6.28319;
  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
  # printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set key bot right \n";
#  printf CMDFILE "set logscale x \n";
#  printf CMDFILE "set logscale y \n";
  printf CMDFILE "set samples 1000 \n";
  printf CMDFILE "set title 'Two systems, one tilted by 60 degrees ' \n";
  printf CMDFILE "set xlabel 'Time (years) ' \n";
  printf CMDFILE "set ylabel 'Stellar radial velocity (m/s) ' \n";
  
  $cmd = "plot [0:2.5][-75:75] 29.76*sin($TWOPI*x/1.0)  ";
  $cmd .= "      with lines lc 3 lw 3 t '1 M_J planet' ";
  $cmd .= " , ";
  $cmd .= "      'make_rv.out' using 1:2   ";
  $cmd .= "      with points ps 1.5 pt 6 lc 1 lw 2 t '2 M_J planet' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;




#####################################################################
# Create a graph illustrating a tilted circular orbit and a
#    face-on elliptical orbit.
#

  $outfile = "time_compare_a.ps";

  $TWOPI = 6.28319;
  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set key bot right \n";
#  printf CMDFILE "set logscale x \n";
#  printf CMDFILE "set logscale y \n";
#  printf CMDFILE "set samples 1000 \n";
  printf CMDFILE "set title 'Two systems of unknown orientation ' \n";
  printf CMDFILE "set xlabel 'Observed X-position (m) ' \n";
  printf CMDFILE "set ylabel 'Observed Y-position (m) ' \n";
  
  $cmd = "plot [-1e8:1e8][-1e8:1e8] 'nbody.circ.center.out' ";
  $cmd .= "      every 97 using (\$3*0.5405):(\$4*0.2644) ";
  $cmd .= "      with points pt 6 ps 0.6 lc 1 lw 5 t 'A' ";
  $cmd .= " , ";
  $cmd .= "                   'nbody.big_eccen.center.out' ";
  $cmd .= "      every 47 using (\$3*0.8500):(\$4*0.8500) ";
  $cmd .= "      with points pt 5 ps 0.8 lc 3 lw 5 t 'B' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;



#####################################################################
# Create a graph illustrating a face-on circular orbit and a
#    face-on elliptical orbit -- same orbits are in previous graph,
#    but this time seen face-on.
#

  $outfile = "time_compare_b.ps";

  $TWOPI = 6.28319;
  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set key bot right \n";
#  printf CMDFILE "set logscale x \n";
#  printf CMDFILE "set logscale y \n";
#  printf CMDFILE "set samples 1000 \n";
  printf CMDFILE "set title 'Two systems, both face-on' \n";
  printf CMDFILE "set xlabel 'Observed X-position (m) ' \n";
  printf CMDFILE "set ylabel 'Observed Y-position (m) ' \n";
  
  $cmd = "plot [-1e8:1e8][-1e8:1e8] 'nbody.circ.center.out' ";
  $cmd .= "      every 97 using (\$3*0.5405):(\$4*0.5405) ";
  $cmd .= "      with points pt 6 ps 0.6 lc 1 lw 5 t 'A' ";
  $cmd .= " , ";
  $cmd .= "                   'nbody.big_eccen.center.out' ";
  $cmd .= "      every 67 using (\$3*0.8500):(\$4*0.8500) ";
  $cmd .= "      with points pt 5 ps 0.6 lc 3 lw 5 t 'B' ";
  $cmd .= " , ";
  $cmd .= "      'time_compare.dat' index 0 using 1:2  ";
  $cmd .= "      with points pt 6 ps 1.0 lc 1 lt 1 lw 5 t '' ";
  $cmd .= " , ";
  $cmd .= "      'time_compare.dat' index 1 using (\$1*0.85):(\$2*0.85)  ";
  $cmd .= "      with points pt 6 ps 1.0 lc 3 lt 1 lw 5 t '' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;



#####################################################################
# Create a graph on which we can draw Keplerian triangles
#    to illustrate Kepler's Second Law -- face-on, circular orbit
#

  $outfile = "kepler_triangles_a.ps";

  $TWOPI = 6.28319;
  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set key bot right \n";
#  printf CMDFILE "set logscale x \n";
#  printf CMDFILE "set logscale y \n";
#  printf CMDFILE "set samples 1000 \n";
#  printf CMDFILE "set title 'Two systems, both face-on' \n";
  printf CMDFILE "set xlabel 'Observed X-position (m) ' \n";
  printf CMDFILE "set ylabel 'Observed Y-position (m) ' \n";
  
  $cmd = "plot [-1e8:1e8][-1e8:1e8] 'nbody.circ.center.out' ";
  $cmd .= "      every 970 using (\$3*0.5405):(\$4*0.5405) ";
  $cmd .= "      with points pt 6 ps 0.6 lc 1 lw 5 t '' ";
  $cmd .= " , ";
  $cmd .= "                          'nbody.circ.center.out' ";
  $cmd .= "      using (\$3*0.5405):(\$4*0.5405) ";
  $cmd .= "      with lines lc 1 lt 1 lw 2 t '' ";
#  $cmd .= " , ";
#  $cmd .= "                   'nbody.big_eccen.center.out' ";
#  $cmd .= "      every 67 using (\$3*0.8500):(\$4*0.8500) ";
#  $cmd .= "      with points pt 5 ps 0.6 lc 3 lw 5 t 'B' ";
  $cmd .= " , ";
  $cmd .= "      'time_compare.dat' index 0 using 1:2  ";
  $cmd .= "      with points pt 6 ps 1.0 lc 1 lt 1 lw 5 t '' ";
#  $cmd .= " , ";
#  $cmd .= "      'time_compare.dat' index 1 using (\$1*0.85):(\$2*0.85)  ";
#  $cmd .= "      with points pt 6 ps 1.0 lc 3 lt 1 lw 5 t '' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;




#####################################################################
# Create a graph on which we can draw Keplerian triangles
#    to illustrate Kepler's Second Law -- tilted, circular orbit
#    (so in this case, the areas will not be equal)
#

  $outfile = "kepler_triangles_b.ps";

  $TWOPI = 6.28319;
  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set key bot right \n";
#  printf CMDFILE "set logscale x \n";
#  printf CMDFILE "set logscale y \n";
#  printf CMDFILE "set samples 1000 \n";
#  printf CMDFILE "set title 'Two systems, both face-on' \n";
  printf CMDFILE "set xlabel 'Observed X-position (m) ' \n";
  printf CMDFILE "set ylabel 'Observed Y-position (m) ' \n";
  
  $cmd = "plot [-1e8:1e8][-1e8:1e8] 'nbody.circ.center.out' ";
  $cmd .= "      every 437 using (\$3*0.5405):(\$4*0.2644) ";
  $cmd .= "      with points pt 6 ps 0.6 lc 1 lw 5 t '' ";
  $cmd .= " , ";
  $cmd .= "                          'nbody.circ.center.out' ";
  $cmd .= "      using (\$3*0.5405):(\$4*0.2644) ";
  $cmd .= "      with lines lt 1 lc 1 lw 2 t '' ";
#  $cmd .= " , ";
#  $cmd .= "                   'nbody.big_eccen.center.out' ";
#  $cmd .= "      every 67 using (\$3*0.8500):(\$4*0.8500) ";
#  $cmd .= "      with points pt 5 ps 0.6 lc 3 lw 5 t 'B' ";
  $cmd .= " , ";
  $cmd .= "      'time_compare.dat' index 0 using 1:2  ";
  $cmd .= "      with points pt 6 ps 1.0 lc 1 lt 1 lw 5 t '' ";
#  $cmd .= " , ";
#  $cmd .= "      'time_compare.dat' index 1 using (\$1*0.85):(\$2*0.85)  ";
#  $cmd .= "      with points pt 6 ps 1.0 lc 3 lt 1 lw 5 t '' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;



#####################################################################
# Create a graph on which we can draw Keplerian triangles
#    to illustrate Kepler's Second Law -- face-on, elliptical orbit
#

  $outfile = "kepler_triangles_c.ps";

  $TWOPI = 6.28319;
  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set key bot right \n";
#  printf CMDFILE "set logscale x \n";
#  printf CMDFILE "set logscale y \n";
#  printf CMDFILE "set samples 1000 \n";
#  printf CMDFILE "set title 'Two systems, both face-on' \n";
  printf CMDFILE "set xlabel 'Observed X-position (m) ' \n";
  printf CMDFILE "set ylabel 'Observed Y-position (m) ' \n";
  
  $cmd = "plot [-1e8:1e8][-1e8:1e8]  ";
#  $cmd .= "                        'nbody.circ.center.out' ";
#  $cmd .= "      every 437 using (\$3*0.5405):(\$4*0.2644) ";
#  $cmd .= "      with points pt 6 ps 0.6 lc 1 lw 5 t '' ";
#  $cmd .= " , ";
#  $cmd .= "                          'nbody.circ.center.out' ";
#  $cmd .= "      using (\$3*0.5405):(\$4*0.2644) ";
#  $cmd .= "      with lines lc 1 lw 2 t '' ";
#  $cmd .= " , ";
  $cmd .= "                   'nbody.big_eccen.center.out' ";
  $cmd .= "      every 87 using (\$3*0.8500):(\$4*0.8500) ";
  $cmd .= "      with points pt 5 ps 0.6 lc 3 lw 5 t 'B' ";
  $cmd .= " , ";
  $cmd .= "                   'nbody.big_eccen.center.out' ";
  $cmd .= "      every 1 using (\$3*0.8500):(\$4*0.8500) ";
  $cmd .= "      with lines lt 1 lc 3 lw 2 t '' ";
#  $cmd .= " , ";
#  $cmd .= "      'time_compare.dat' index 0 using 1:2  ";
#  $cmd .= "      with points pt 6 ps 1.0 lc 1 lt 1 lw 5 t '' ";
  $cmd .= " , ";
  $cmd .= "      'time_compare.dat' index 1 using (\$1*0.85):(\$2*0.85)  ";
  $cmd .= "      with points pt 6 ps 1.0 lc 3 lt 1 lw 5 t '' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;


#####################################################################
# Create a graph on which we'll show to find the focus
#

  $outfile = "ellipse_geometry_a.ps";

  $TWOPI = 6.28319;
  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' size 9,9 \n";
  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set key bot right \n";
#  printf CMDFILE "set logscale x \n";
#  printf CMDFILE "set logscale y \n";
#  printf CMDFILE "set samples 1000 \n";
#  printf CMDFILE "set title 'Two systems, both face-on' \n";
  printf CMDFILE "set xlabel 'Observed X-position (m) ' \n";
  printf CMDFILE "set ylabel 'Observed Y-position (m) ' \n";
  
  $cmd = "plot [-1.5e8:1.5e8][-1.5e8:1.5e8]  ";
  $cmd .= "                        'nbody.med_eccen.center.out' ";
  $cmd .= "      every 119 using (\$4*1.0000):(\$5*1.0000) ";
  $cmd .= "      with points pt 6 ps 0.6 lc 1 lw 5 t '' ";
#  $cmd .= " , ";
#  $cmd .= "                          'nbody.circ.center.out' ";
#  $cmd .= "      using (\$3*0.5405):(\$4*0.2644) ";
#  $cmd .= "      with lines lc 1 lw 2 t '' ";
#  $cmd .= " , ";
#  $cmd .= "                   'nbody.big_eccen.center.out' ";
#  $cmd .= "      every 87 using (\$3*0.8500):(\$4*0.8500) ";
#  $cmd .= "      with points pt 5 ps 0.6 lc 3 lw 5 t 'B' ";
#  $cmd .= " , ";
#  $cmd .= "                   'nbody.big_eccen.center.out' ";
#  $cmd .= "      every 1 using (\$3*0.8500):(\$4*0.8500) ";
#  $cmd .= "      with lines lt 1 lc 3 lw 2 t '' ";
#  $cmd .= " , ";
#  $cmd .= "      'time_compare.dat' index 0 using 1:2  ";
#  $cmd .= "      with points pt 6 ps 1.0 lc 1 lt 1 lw 5 t '' ";
#  $cmd .= " , ";
#  $cmd .= "      'time_compare.dat' index 1 using (\$1*0.85):(\$2*0.85)  ";
#  $cmd .= "      with points pt 6 ps 1.0 lc 3 lt 1 lw 5 t '' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/\.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;


#####################################################################
# Create a graph with the foci marked, for convenience
#

  $outfile = "ellipse_geometry_b.ps";

  $TWOPI = 6.28319;
  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' size 9,9 \n";
  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set key bot right \n";
#  printf CMDFILE "set logscale x \n";
#  printf CMDFILE "set logscale y \n";
#  printf CMDFILE "set samples 1000 \n";
#  printf CMDFILE "set title 'Two systems, both face-on' \n";
  printf CMDFILE "set xlabel 'Observed X-position (m) ' \n";
  printf CMDFILE "set ylabel 'Observed Y-position (m) ' \n";
  
  $cmd = "plot [-1.5e8:1.5e8][-1.5e8:1.5e8]  ";
  $cmd .= "                        'nbody.med_eccen.center.out' ";
  $cmd .= "      every 119 using (\$4*1.0000):(\$5*1.0000) ";
  $cmd .= "      with points pt 6 ps 0.6 lc 1 lw 5 t '' ";
#  $cmd .= " , ";
#  $cmd .= "                          'nbody.circ.center.out' ";
#  $cmd .= "      using (\$3*0.5405):(\$4*0.2644) ";
#  $cmd .= "      with lines lc 1 lw 2 t '' ";
#  $cmd .= " , ";
#  $cmd .= "                   'nbody.big_eccen.center.out' ";
#  $cmd .= "      every 87 using (\$3*0.8500):(\$4*0.8500) ";
#  $cmd .= "      with points pt 5 ps 0.6 lc 3 lw 5 t 'B' ";
#  $cmd .= " , ";
#  $cmd .= "                   'nbody.big_eccen.center.out' ";
#  $cmd .= "      every 1 using (\$3*0.8500):(\$4*0.8500) ";
#  $cmd .= "      with lines lt 1 lc 3 lw 2 t '' ";
  $cmd .= " , ";
  $cmd .= "      'ellipse_geometry.dat' using 1:2  ";
  $cmd .= "      with points pt 3 ps 0.5 lc -1 lt 1 lw 3 t '' ";
#  $cmd .= " , ";
#  $cmd .= "      'time_compare.dat' index 1 using (\$1*0.85):(\$2*0.85)  ";
#  $cmd .= "      with points pt 6 ps 1.0 lc 3 lt 1 lw 5 t '' ";

  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/\.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;


#####################################################################
# Create a graph showing area swept out by radius vector
#     for two foci
#

  $outfile = "area_swept_a.ps";

  $TWOPI = 6.28319;
  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
#  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set key top right \n";
#  printf CMDFILE "set logscale x \n";
#  printf CMDFILE "set logscale y \n";
#  printf CMDFILE "set samples 1000 \n";
  printf CMDFILE "set title 'Area swept per hour, measured from two foci' \n";
  printf CMDFILE "set xlabel 'Starting hour in orbit' \n";
  printf CMDFILE "set ylabel 'Area of triangle (sq. meters)' \n";
  
  $cmd = "plot [][1.6e12:1.9e13]  ";
  $cmd .= "                        'sweep_left.out' ";
  $cmd .= "      using 2:4 ";
  $cmd .= "      with points pt 6 ps 0.5 lc 1 lw 2 t 'left-hand focus' ";
  $cmd .= " , ";
  $cmd .= "                        'sweep_right.out' ";
  $cmd .= "      using 2:4 ";
  $cmd .= "      with points pt 6 ps 0.5 lc 3 lw 2 t 'right-hand focus' ";


  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/\.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;



#####################################################################
# Create a graph showing area swept out by radius vector
#     for two foci -- for an ellipse which has been projected,
#     so that neither provides a constant areal velocity
#

  $outfile = "area_swept_b.ps";

  $TWOPI = 6.28319;
  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
#  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set key top right \n";
#  printf CMDFILE "set logscale x \n";
#  printf CMDFILE "set logscale y \n";
#  printf CMDFILE "set samples 1000 \n";
  printf CMDFILE "set title 'Area swept per hour, measured from two foci' \n";
  printf CMDFILE "set xlabel 'Starting hour in orbit' \n";
  printf CMDFILE "set ylabel 'Area of triangle (sq. meters)' \n";
  
  #$cmd = "plot [][1.6e12:1.9e13]  ";
  $cmd = "plot [][]  ";
  $cmd .= "                        'sweep_proj_left.out' ";
  $cmd .= "      using 2:4 ";
  $cmd .= "      with points pt 6 ps 0.5 lc 1 lw 2 t 'left-hand focus' ";
  $cmd .= " , ";
  $cmd .= "                        'sweep_proj_right.out' ";
  $cmd .= "      using 2:4 ";
  $cmd .= "      with points pt 6 ps 0.5 lc 3 lw 2 t 'right-hand focus' ";


  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/\.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;



#####################################################################
# Create a graph showing area swept out by radius vector
#     for a set of ellipses which have been de-projected by
#     various angles, so that we can see the change with angle
#

  $outfile = "best_deproject_a.ps";

  $TWOPI = 6.28319;
  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
#  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set key top right \n";
#  printf CMDFILE "set logscale x \n";
#  printf CMDFILE "set logscale y \n";
#  printf CMDFILE "set samples 1000 \n";
  printf CMDFILE "set title 'Area swept per hour, for choices of de-projection ' \n";
  printf CMDFILE "set xlabel 'Starting hour in orbit' \n";
  printf CMDFILE "set ylabel 'Area of triangle (sq. meters)' \n";
  
  $cmd = "plot [][3.9e12:1.05e13]  ";
  #$cmd = "plot [][]  ";
  $cmd .= "                        'nbody.med_eccen.proj.fit' ";
  $cmd .= "      index 1 using 2:4 ";
  $cmd .= "      with points pt 6 ps 0.5 lc 1 lw 2 t 'observed orbit' ";
  $cmd .= " , ";
  $cmd .= "                        'nbody.med_eccen.deprojA.fit' ";
  $cmd .= "      index 1 using 2:4 ";
  $cmd .= "      with points pt 6 ps 0.5 lc 2 lw 2 t 'de-projected 10 deg' ";
  $cmd .= " , ";
  $cmd .= "                        'nbody.med_eccen.deprojB.fit' ";
  $cmd .= "      index 1 using 2:4 ";
  $cmd .= "      with points pt 6 ps 0.5 lc 3 lw 2 t 'de-projected 20 deg' ";
  $cmd .= " , ";
  $cmd .= "                        'nbody.med_eccen.deprojC.fit' ";
  $cmd .= "      index 1 using 2:4 ";
  $cmd .= "      with points pt 6 ps 0.5 lc 4 lw 2 t 'de-projected 30 deg' ";
  $cmd .= " , ";
  $cmd .= "                        'nbody.med_eccen.deprojE.fit' ";
  $cmd .= "      index 1 using 2:4 ";
  $cmd .= "      with points pt 6 ps 0.5 lc -1 lw 2 t 'de-projected 35 deg' ";




  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/\.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;


#####################################################################
# Create a graph showing the ellipses created as we de-project
#     an observed orbit by various angles
#

  $outfile = "best_deproject_b.ps";

  $TWOPI = 6.28319;
  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' size 11,11 \n";
  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set key at 0,-2e7\n";
#  printf CMDFILE "set logscale x \n";
#  printf CMDFILE "set logscale y \n";
#  printf CMDFILE "set samples 1000 \n";
  printf CMDFILE "set xtics 5e7 \n";
  printf CMDFILE "set title '' \n";
  printf CMDFILE "set xlabel 'X position (meters) ' \n";
  printf CMDFILE "set ylabel 'Y position (meters) ' \n";
  
  $cmd = "plot [-1.02e8:1.02e8][-1.02e8:1.02e8]  ";
  #$cmd = "plot [][]  ";
  $cmd .= "                        'nbody.med_eccen.proj.out' ";
  $cmd .= "      using 3:4 ";
  $cmd .= "      with points pt 6 ps 0.4 lc 1 lw 2 t 'observed' ";
  $cmd .= " , ";
  $cmd .= "                        'nbody.med_eccen.deprojA.out' ";
  $cmd .= "      using 4:5 ";
  $cmd .= "      with points pt 6 ps 0.2 lc -1 lw 1 t 'de-proj 10 deg' ";
  $cmd .= " , ";
  $cmd .= "                        'nbody.med_eccen.deprojB.out' ";
  $cmd .= "      using 4:5 ";
  $cmd .= "      with points pt 6 ps 0.5 lc 3 lw 2 t 'de-proj 20 deg' ";
  $cmd .= " , ";
  $cmd .= "                        'nbody.med_eccen.deprojC.out' ";
  $cmd .= "      using 4:5 ";
  $cmd .= "      with points pt 6 ps 0.5 lc 4 lw 2 t 'de-proj 30 deg' ";
  $cmd .= " , ";
  $cmd .= "                        'nbody.med_eccen.deprojE.out' ";
  $cmd .= "      using 4:5 ";
  $cmd .= "      with points pt 6 ps 0.5 lc 2 lw 2 t 'de-proj 35 deg' ";
  $cmd .= " , ";
  $cmd .= "                        'deproj_foci.dat' ";
  $cmd .= "      index 0 using 1:2 ";
  $cmd .= "      with points pt 6 ps 0.5 lc 1 lw 2 t '' ";
  $cmd .= " , ";
  $cmd .= "                        'deproj_foci.dat' ";
  $cmd .= "      index 1 using 1:2 ";
  $cmd .= "      with points pt 6 ps 0.5 lc -1 lw 2 t '' ";
  $cmd .= " , ";
  $cmd .= "                        'deproj_foci.dat' ";
  $cmd .= "      index 2 using 1:2 ";
  $cmd .= "      with points pt 6 ps 0.8 lc 3 lw 4 t '' ";
  $cmd .= " , ";
  $cmd .= "                        'deproj_foci.dat' ";
  $cmd .= "      index 3 using 1:2 ";
  $cmd .= "      with points pt 6 ps 0.8 lc 4 lw 3 t '' ";
  $cmd .= " , ";
  $cmd .= "                        'deproj_foci.dat' ";
  $cmd .= "      index 4 using 1:2 ";
  $cmd .= "      with points pt 6 ps 0.8 lc 2 lw 4 t '' ";




  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/\.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;


#####################################################################
# Create a graph showing an observed orbit of host star 
#     around center of mass -- serves as basis for the 
#     de-projection discussion
#

  $outfile = "best_deproject_c.ps";

  $TWOPI = 6.28319;
  
  open (CMDFILE, ">$cmdfile") || die("can't open $cmdfile for writing");
  printf CMDFILE "set output '$outfile' \n";
  printf CMDFILE "set term postscript color enhanced 'Helvetica,18' \n";
  printf CMDFILE "set size square \n";
  printf CMDFILE "set grid \n";
  printf CMDFILE "set key at 0,-2e7\n";
#  printf CMDFILE "set logscale x \n";
#  printf CMDFILE "set logscale y \n";
#  printf CMDFILE "set samples 1000 \n";
  printf CMDFILE "set xtics 5e7 \n";
  printf CMDFILE "set title 'An observed orbit and its foci' \n";
  printf CMDFILE "set xlabel 'X position (meters) ' \n";
  printf CMDFILE "set ylabel 'Y position (meters) ' \n";
  
  $cmd = "plot [-1.1e8:1.1e8][-1.1e8:1.1e8]  ";
  #$cmd = "plot [][]  ";
  $cmd .= "                        'nbody.med_eccen.proj.out' ";
  $cmd .= "      every 59 using 3:4 ";
  $cmd .= "      with points pt 6 ps 0.4 lc 1 lw 2 t '' ";
#  $cmd .= " , ";
#  $cmd .= "                        'nbody.med_eccen.deprojA.out' ";
#  $cmd .= "      using 4:5 ";
#  $cmd .= "      with points pt 6 ps 0.2 lc -1 lw 1 t 'de-proj 10 deg' ";
#  $cmd .= " , ";
#  $cmd .= "                        'nbody.med_eccen.deprojB.out' ";
#  $cmd .= "      using 4:5 ";
#  $cmd .= "      with points pt 6 ps 0.5 lc 3 lw 2 t 'de-proj 20 deg' ";
#  $cmd .= " , ";
#  $cmd .= "                        'nbody.med_eccen.deprojC.out' ";
#  $cmd .= "      using 4:5 ";
#  $cmd .= "      with points pt 6 ps 0.5 lc 4 lw 2 t 'de-proj 30 deg' ";
#  $cmd .= " , ";
#  $cmd .= "                        'nbody.med_eccen.deprojE.out' ";
#  $cmd .= "      using 4:5 ";
#  $cmd .= "      with points pt 6 ps 0.5 lc 2 lw 2 t 'de-proj 35 deg' ";
  $cmd .= " , ";
  $cmd .= "                        'deproj_foci.dat' ";
  $cmd .= "      index 0 using 1:2 ";
  $cmd .= "      with points pt 6 ps 1.5 lc 1 lw 3 t '' ";
#  $cmd .= " , ";
#  $cmd .= "                        'deproj_foci.dat' ";
#  $cmd .= "      index 1 using 1:2 ";
#  $cmd .= "      with points pt 6 ps 0.5 lc -1 lw 2 t '' ";
#  $cmd .= " , ";
#  $cmd .= "                        'deproj_foci.dat' ";
#  $cmd .= "      index 2 using 1:2 ";
#  $cmd .= "      with points pt 6 ps 0.8 lc 3 lw 4 t '' ";
#  $cmd .= " , ";
#  $cmd .= "                        'deproj_foci.dat' ";
#  $cmd .= "      index 3 using 1:2 ";
#  $cmd .= "      with points pt 6 ps 0.8 lc 4 lw 3 t '' ";
#  $cmd .= " , ";
#  $cmd .= "                        'deproj_foci.dat' ";
#  $cmd .= "      index 4 using 1:2 ";
#  $cmd .= "      with points pt 6 ps 0.8 lc 2 lw 4 t '' ";




  printf CMDFILE "$cmd \n";
  printf CMDFILE "set output \n";
  printf CMDFILE "quit \n";
  close(CMDFILE) ;
  
  $retval = `gnuplot < $cmdfile`;
  printf "retval is ..%s..\n", $retval;

  # create a PNG copy of any Postscript file
    $png_file = $outfile;
    $png_file =~ s/\.ps/.png/;
    $cmd = "convert -flatten -background white -rotate 90 $outfile $png_file ";
    printf "cmd is ..$cmd.. \n";
    $retval = `$cmd`;












exit 0;



