#!/usr/bin/perl
#
# Run a query on SIMBAD to lookup the coordinates of an object.
#
# MWR 6/22/2004

use strict;
use POSIX;
use HTTP::Request::Common qw(POST);
use LWP::UserAgent;


# global variables
my($debug);
my($name);
my($retval, $ra, $dec, $epoch_jd, $period);
my($num_markiv_lines, @markiv_lines);


# should we be verbose?
$debug = 0;

# flush output immediately
$| = 1;



if (1 == 1) {
  ($retval, $name, $ra, $dec, $epoch_jd, $period) = get_variable_info("CQ Cep");
  if ($retval != 0) {
    printf STDERR "an error occurred ... \n";
    exit(1);
  }
  else {
    printf "$name: ra $ra  dec $dec  epoch_jd $epoch_jd  period $period \n";
  }
}
else {
  $name = "RR Lyr";
  $ra = 291.36630375;
  $dec = 42.78436;
  $epoch_jd = 2442923.419300;
  $period = 0.5668677600;
}


# now grab data from the Mark IV database (if it exists)
($num_markiv_lines, @markiv_lines) = get_markiv_data($ra, $dec);
if ($num_markiv_lines <= 0) {
  printf STDERR "get_markiv_data returns no data \n";
  exit(1);
}

# make a nice graph showing the measurements phased with the known period
if (make_phased_graph($name, $ra, $dec, $epoch_jd, $period,
                         $num_markiv_lines, @markiv_lines) != 0) {
  printf STDERR "make_phased_graph returns with error \n";
  exit(1);
}

exit(0);


###########################################################################
# PROCEDURE: get_markiv_data
# 
# DESCRIPTION: go to the Mark IV database and request measurements 
#              of the star at a given (RA, Dec).  Place the
#              resulting measurements into an array and return them.
#
# RETURNS:
#              (num_lines, @lines)
#                    num_lines         number of lines of data
#                    lines             array of data, one V,I pair per line
#
sub get_markiv_data {

  my($i);
  my($ra, $dec);
  my($ua, $url, $req, $reply);
  my($match_radius, $match_scale, $radec_string);
  my($retval, $rah, $ram, $ras, $decsign, $decd, $decm, $decs);
  my($num_reply_lines, @reply_lines);
  my($num_data_lines, @data_lines);

  $ra = $_[0];
  $dec = $_[1];

  $ua = LWP::UserAgent->new;

  $num_data_lines = 0;
  @data_lines = "";

  $match_radius = 5;
  $match_scale = "Seconds";
  # we must convert RA and Dec to sexigesimal notation in order to
  #   make use of full accuracy of positions (sigh)
  ($rah, $ram, $ras, $decsign, $decd, $decm, $decs) = 
                convert_deg_to_baby($ra, $dec);
  if ($debug > 0) {
    printf "get_markiv_data: rah $rah ram $ram ras $ras \n";
    printf "get_markiv_data: decsign $decsign decd $decd decm $decm decs $decs\n";
  }
  # now make a nice string that the Mark IV database engine will accept
  $radec_string = sprintf "%02d%02d%04.1f%s%02d%02d%04.1f",
                $rah, $ram, $ras, $decsign, $decd, $decm, $decs;
  if ($debug > 0) {  
    printf "get_markiv_data: radec_string: ..%s..\n", $radec_string;
  }

  # make a query to the Mark IV database
  $url = "http://sallman.tass-survey.org/servlet/markiv/action/DataDownload";
  $req = POST $url, [ "header" => 'false',
                      "compress" => 'false',
                      "position" => $radec_string,
                      "radius" => $match_radius,
                      "scale" => $match_scale
                    ];
  $reply = $ua->request($req)->as_string;
  if ($debug > 0) {
    printf "get_markiv_data: reply from Mark IV database follows .. \n";
    printf "..%s..", $reply;
  }

  # pick out just the lines with measurements 
  # first, skip lines until we reach a "<pre>"
  ($num_reply_lines, @reply_lines) = split_reply($reply);
  for ($i = 0; $i < $num_reply_lines; $i++) {
    if ($reply_lines[$i] =~ /<pre>/) {
      last;
    }
  }
  if ($i == $num_reply_lines) {
    # there wasn't any real data 
    if ($debug > 0) {
      printf "get_markiv_data: no real data ... \n";
    }
  }
  else {
    # now, we copy the real data into the "@data_lines" array
    for ($i++; $i < $num_reply_lines; $i++) {
      if ($reply_lines[$i] =~ /Sorry/) {
        last;
      }
      if ($reply_lines[$i] =~ /<\/pre>/) {
        last;
      }
      $data_lines[$num_data_lines] = " " . $reply_lines[$i];
      chomp($data_lines[$num_data_lines]);
      if ($debug > 0) {
        printf "get_markiv_data: data line %5d is ..%s..\n",
                     $num_data_lines, $data_lines[$num_data_lines];
      }
      $num_data_lines++;
    }
  }
  # there is a single blank line just before the </pre>, so get rid of it
  if ($num_data_lines > 0) {
    $num_data_lines--;
  }

  return($num_data_lines, @data_lines);
}


###########################################################################
# PROCEDURE: convert_deg_to_baby
#
# DESCRIPTION: convert (RA, Dec) from decimal degrees to babylonian
#              notation: HH MM SS.sss and "sign" DD MM SS.ss.
#
#              Note that we return values which are ALWAYS positive
#              for decd, decm, decs; only the separate "sign" value
#              indicates the sign of the Declination.  The "sign"
#              will be either "+" or "-".
#
# RETURNS:
#             (rah, ram, ras, decsign, decd, decm, decs)
#
sub convert_deg_to_baby {

  my($ra, $dec);
  my($rah, $ram, $ras, $decsign, $decd, $decm, $decs);

  $ra = $_[0];
  $dec = $_[1];

  if ($debug > 0) {
    printf "convert_deg_to_baby: ra ..%s..  dec ..%s.. \n";
  }

  # pick apart the RA first
  $rah = int($ra/15.0);
  $ram = int(($ra - 15*$rah)*4.0);
  $ras = $ra - ($rah*15) - ($ram/4.0);
  $ras *= 240.0;
  
  # now figure out the Dec value
  if ($dec < 0) {
    $decsign = "-";
    $dec = 0.0 - $dec;
  } else {
    $decsign = "+";
  }
  $decd = int($dec);
  $decm = int(($dec - $decd)*60.0);
  $decs = $dec - ($decd) - ($decm/60);
  $decs *= 3600.0;

  if ($debug > 0) {
    printf "convert_deg_to_baby: RA  ..%s.. ..%s.. ..%s.. \n", $rah, $ram, $ras;
    printf "convert_deg_to_baby: Dec ..%s.. ..%s.. ..%s.. ..%s.. \n", 
                     $decsign, $decd, $decm, $decs;
  }

  return($rah, $ram, $ras, $decsign, $decd, $decm, $decs);

}


###########################################################################
# PROCEDURE: get_variable_info 
#
# DESCRIPTION: Given a string which consists of either 
#                     a name:         "TT Boo"
#              or
#                     a position:     "12.2343 +23.8873"
#                       (J2000)       "12 33 43.7 -02 12 58"
#
#              (which we assume is the position of a variable star)
#              use SIMBAD to look up information on the star.
#              We get 5 quantities:
#
#                       name              official name of variable
#                       ra                Right Ascension (J2000)
#                       dec               Declination (J2000)
#                       epoch             time of max or min (JD)
#                       period            period of variability (days)
#
# RETURNS:
#              a list of 6 quantities:
#                       proc_retval       0 if all okay, 1 if error occurs
#                       name
#                       ra
#                       dec
#                       epoch
#                       period
#
#                  
sub get_variable_info {

  my($input_name_pos);
  my($name);
  my($proc_retval, $ra, $dec, $epoch_jd, $period);
  my($ua);
  my($line);
  my($num_reply_lines, @reply_lines);
  my($search_string);
  
  $input_name_pos = $_[0];

  # by default, we fail
  $proc_retval = 1;   

  # initialize some return values to nonsensical values 
  $ra = -99.0;
  $dec = -99.0;
  $epoch_jd = -99.0;
  $period = -99.0;

  $ua = LWP::UserAgent->new;

  # these are the parameters we'll need for the lookup request to SIMBAD
  my $Protocol = "html";
  my $Ident = $input_name_pos;
  my $NbIdent = 1;
  my $Radius = "5";
  my $Radius_unit = "arc sec";
  my $CooFrame = "FK5";
  my $CooEpoch = 2000;
  my $CooEquinox = 2000;
  my $Output_max = "all";
  my $Output_mesdisp = "N";
  my $Bibyear1 = 1983;
  my $Bibyear2 = 2004;
  my $Frame1 = "FK5";
  my $Equi1 = "2000.0";
  my $Epoch1 = "2000.0";
  my $Frame2 = "none";
  my $Equi2 = "2000.0";
  my $Epoch2 = "2000.0";
  my $Frame3 = "none";
  my $Equi3 = "2000.0";
  my $Epoch3 = "2000.0";
  
  my $url = "http://simbad.u-strasbg.fr/sim-id.pl";
  my $req = POST $url, [ 
                         # protocol => "ascii",
                         "Ident" => $Ident,
                         "NbIdent" => $NbIdent,
                         "Radius" => $Radius,
                         "Radius.unit" => $Radius_unit,
                         "CooFrame"  => $CooFrame,
                         "CooEpoch"  => $CooEpoch,
                         "CooEquinox"  => $CooEquinox,
                         "output.max" => $Output_max,
                         "output.mesdisp" => $Output_mesdisp,
                         "Bibyear1" => $Bibyear1,        
                         "Bibyear2" => $Bibyear2,        
                         "Frame1" => $Frame1,
                         "Equi1" => $Equi1,
                         "Frame2" => $Frame2,
                         "Equi2" => $Equi2,
                         "Frame3" => $Frame3,
                         "Equi3" => $Equi3
                       ];
  my $reply_stuff = $ua->request($req)->as_string;
  
  if ($debug > 0) {
    printf "reply_stuff is ..$reply_stuff.. \n";
  }
  
  #   split the reply into individual lines
  ($num_reply_lines, @reply_lines) = split_reply($reply_stuff);
  if ($num_reply_lines == 0) {
      printf STDERR "SIMBAD ID query returns empty?! \n";
      return($proc_retval, $name, $ra, $dec, $epoch_jd, $period);
  }
  # make the special search through the data for the line
  #   containing the RA and Dec; if found, we'll get 
  #   back RA and Dec in decimal degrees (J2000)
  ($retval, $ra, $dec) = get_radec($num_reply_lines, @reply_lines);
  if ($retval != 0) {
    printf STDERR "get_radec fails ?! \n";
    return($proc_retval, $name, $ra, $dec, $epoch_jd, $period);
  }
 
  
  #        look for the entry for this object in the GCVS, 
  #        to get the name, period and time of max/min
  ($retval, $name, $epoch_jd, $period) = get_epoch_period($num_reply_lines,
                                     @reply_lines);
  if ($retval != 0) {
    printf STDERR "get_variable_info: get_epoch_period fails !? \n";
    return($proc_retval, $name, $ra, $dec, $epoch_jd, $period);
  }


  return($retval, $name, $ra, $dec, $epoch_jd, $period);
}



############################################################################
# PROCEDURE: get_epoch_period
# 
# DESCRIPTION: Given an array of lines containing the reply to our
#              request to SIMBAD, search through it for the 
#              line which contains a URL to _another_ SIMBAD 
#              bit of info: the entry for this star in the GCVS.
#              Post that URL and look in _its_ reply for the
#              official GCVS name, period and epoch of maximum/minimum light.
#
# RETURNS:
#              (proc_retval, name, epoch_jd, period)
#                where     proc_retval    = 0 if all OK, = 1 if error occurs
#                          name           official GCVS name of variable
#                          epoch_jd       JD of max/min light
#                          period         period of variable (days)
#              
#              we use values of -99 for epoch_jd and period by default
#                 to indicate no real data available
#              
#     
sub get_epoch_period {
  
  my($proc_retval);
  my($i);
  my($num_reply_lines);
  my($name);
  my($gcvs_num_lines, @gcvs_lines);
  my($gcvs_url, $gcvs_reply);
  my($retval, $search_string, $line);
  my($epoch_jd, $period);
  

  $num_reply_lines = $_[0];

  $proc_retval = 1;
  $name = "no_name";
  $epoch_jd = -99.0;
  $period = -99.0;
 
  for ($i = 0;  $i < $num_reply_lines; $i++) {
    # skip ahead to the "Catalogue information" line
    if ($_[1 + $i] =~ /^\<LI\>Catalogue information/) {
      if ($debug > 0) {
        printf "get_epoch_period: found Catalogue info in line ..%s..", 
                       $_[1 + $i];
      }
      last;
    }
  }
  if ($i == $num_reply_lines) {
    # there wasn't any "Catalogue information" line
    if ($debug > 0) {
      printf "get_epoch_period: no Catalogue info ... \n";
    }
    return($proc_retval, $name, $epoch_jd, $period);
  }

  # get the RA and Dec, which occur a few lines later ...
  for ( ;  $i < $num_reply_lines; $i++) {
    if ($_[1 + $i] =~ /\V\*/) {
      if ($debug > 0) {
        printf "get_epoch_period: found V* in line ..%s..", $_[1 + $i];
      }
      last;
    }
    else {
      if ($debug > 0) {
        printf "get_epoch_period: skipping line $i is ..%s..\n", $_[1 + $i];
      }
    }
  }
  if ($i == $num_reply_lines) {
    # there wasn't any link to GCVS entry 
    if ($debug > 0) {
      printf "get_epoch_period: no GCVS entry ... \n";
    }
    return($proc_retval, $name, $epoch_jd, $period);
  }

    $line = $_[1 + $i];
    chomp($line);
    if ($debug > 0) {
      printf "get_epoch_perod: look for GCVS entry in ..%s.. \n", $line;
    }
    if ($line !~ /.*HREF="(.*)"/) {
      # rats!  No link to the GCVS query?!
      if ($debug > 0) {
        printf "get_epoch_period: no HREF inside line ..%s.. ?! \n", $line;
        return($proc_retval, $name, $epoch_jd, $period);
      }
    }
    # the line did have the right format -- here's the URL for GCVS query
    $gcvs_url = $1;
    if ($debug > 0) {
      printf "get_epoch_period: gcvs_url is ..%s..\n", $gcvs_url;
    }

    # and now we run the GCVS query to get information on the variable star
    my $gcvs_ua = LWP::UserAgent->new;
    my $gcvs_req = new HTTP::Request ('GET', => $gcvs_url);
    $gcvs_reply = $gcvs_ua->request($gcvs_req)->as_string;
    if ($debug > 0) {
      printf "get_epoch_period: gcvs_reply is ..$gcvs_reply.. \n";
    }

    # now we have to scan through the output of the gcvs reply to find
    #   the values for "Period" and "Epoch"
    ($gcvs_num_lines, @gcvs_lines) = split_reply($gcvs_reply);
    if ($gcvs_num_lines == 0) {
      printf STDERR "get_epoch_period: gcvs query returns empty?! \n";
      return($proc_retval, $name, $epoch_jd, $period);
    }

    # look for the one line which contains the GCVS Name of the variable
    $search_string = "Variable star designation";
    ($retval, $line) = find_reply_line($search_string,
                           $gcvs_num_lines, @gcvs_lines);
    if ($retval != 0) {
      printf STDERR "get_epoch_period: can't find string $search_string\n";
      return($proc_retval, $name, $epoch_jd, $period);
    }
    # strip out the name from this line
    if ($line !~ /.*Vname=.*?">(.*?)<\/A/) {
      printf STDERR "get_epoch_period: can't match pattern for name in line ..%s.. ?! \n",
                            $line;
      return($proc_retval, $name, $epoch_jd, $period);
    }
    else {
      $name = $1;
      # and get rid of extra white space
      $name =~ s/\s\s+/ /g;
      if ($debug > 0) {
        printf "get_epoch_period: name is ..%s.. \n", $name;
      }  
    }
    

    # look for the one line which contains the Epoch of max/min light
    $search_string = "Epoch for maximum light";
    ($retval, $line) = find_reply_line($search_string,
                           $gcvs_num_lines, @gcvs_lines);
    if ($retval != 0) {
      printf STDERR "get_epoch_period: can't find string $search_string\n";
      return($proc_retval, $name, $epoch_jd, $period);
    }
    # strip out the Julian Date of max/min light from this line
    if ($line !~ /.*<B>\s*([0-9].*)<\/B>/) {
      printf STDERR "get_epoch_period: can't match pattern for JD in Epoch line ..%s.. ?! \n",
                            $line;
      return($proc_retval, $name, $epoch_jd, $period);
    }
    else {
      $epoch_jd = $1;
      if ($debug > 0) {
        printf "get_epoch_period: epoch_jd is ..%s.. \n", $epoch_jd;
      }  
    }
    
    # look for the one line which contains the Period (in days)
    $search_string = "Period of the variable star";
    ($retval, $line) = find_reply_line($search_string,
                           $gcvs_num_lines, @gcvs_lines);
    if ($retval != 0) {
      printf STDERR "get_epoch_period: can't find string $search_string\n";
      return($proc_retval, $name, $epoch_jd, $period);
    }
    # strip out the Period (we assume in days) from this line
    if ($line !~ /.*<B>\s*([0-9].*)<\/B>/) {
      printf STDERR "get_epoch_period: can't match pattern for period in Period line ..%s.. ?! \n",
                            $line;
      return($proc_retval, $name, $epoch_jd, $period);
      return(1);
    }
    else {
      $period = $1;
      if ($debug > 0) {
        printf "get_epoch_period: period is ..%s.. \n", $period;
      }  
    }
    
  return($retval, $name, $epoch_jd, $period);

}


###########################################################################
# PROCEDURE: get_radec
#
# DESCRIPTION: We are given an array with lines from the reply to 
#              our HTTP request to SIMBAD.  We scan through them
#              for the line with the RA and Dec.  After we find it,
#              we read the RA and Dec in sexigesimal form and convert
#              to decimal degrees.
#
# RETURNS:
#              (retval, ra, dec)
#                where     retval  =       0 if all OK, 1 if error
#                          ra              RA (decimal degrees J2000)
#                          dec             Dec (decimal degrees J2000)
#
sub get_radec {
  my($num_lines, $line);
  my($i);
  my($search_string);
  my($retval, $ra, $dec);
  my(@words);
  my($rah, $ram, $ras, $decd, $decm, $decs, $sign);

  $num_lines = $_[0];
 
  $retval = 1;
  $ra = -99.0;
  $dec = -99.0;
  $search_string = "ICRS 2000.0 coordinates";
  
  # first, skip lines until we reach a line which starts
  #    "ICRS 2000.0 coordinates"
  for ($i = 0; $i < $num_lines; $i++) {
    if ($_[1 + $i] =~ /^$search_string/) {
      last;
    }
  }
  if ($i == $num_lines) {
    # there wasn't any real data 
    if ($debug > 0) {
      printf "get_radec: couldn't find string $search_string \n";
    }
    return($retval, $ra, $dec);
  }

  # get the RA and Dec, which occur a few lines later ...
  $i++;
  # skip </TD> line 
  $i++;
  # skip <TD COLSPAN=2> line
  $i++;
  # now we should see a line that starts like this
  # <b>06 45 08.9173 -16 42 58.017
  # strip out the <b>, then pick out RA and Dec
  $line = $_[1 + $i];
  $line =~ s/\<b\>//;
  @words = split(/\s+/, $line);
  $rah = $words[0];
  $ram = $words[1];
  $ras = $words[2];
  $decd = $words[3];
  $decm = $words[4];
  $decs = $words[5];
  if ($debug > 0) {
    printf "words are rah ..%s.. ram ..%s.. ras ..%s.. \n", $rah, $ram, $ras;
    printf "words are decd ..%s.. decm ..%s.. decs ..%s.. \n", 
                                      $decd, $decm, $decs;
  }
  if (substr($decd, 0, 1) eq "-") {
    if ($debug > 0) {
      printf "sign is negative \n";
    }
    $sign = -1;
    # make the decd value positive, so we can convert to degrees properly
    if ($decd < 0) {
      $decd = 0.0 - $decd;
    }
  }
  else {
    if ($debug > 0) {
      printf "sign is positive \n";
    }
    $sign = 1;
  }

  if ($debug > 0) {
    printf "get_radec: RA is %02d %02d %04.1f  %+02d %02d %04.1f \n",
         $rah, $ram, $ras, $decd, $decm, $decs;
  }

  # now convert RA and Dec to decimal degrees
  $ra = ($rah*15.0) + ($ram/4.0) + ($ras/240.0);
  $dec = abs($decd) + abs($decm/60.0) + abs($decs/3600.0);
  $dec *= $sign;
  if ($debug > 0) {
    printf "get_radec: RA is %12.5f  Dec is %12.5f \n", $ra, $dec;
  }


  # and we're done
  $retval = 0;
  return($retval, $ra, $dec);

}



#############################################################################
# PROCEDURE: split_reply
#
# DESCRIPTION: Given the entire reply from an HTTP request as a single
#              string, break it up into individual lines.  Count the lines,
#              and create an array with one line per element.
#
# RETURNS:
#              (num_of_lines, array_of_lines)
#
sub split_reply {
  my($num_reply_lines);
  my(@reply_lines);
  my($big_reply_string);
  my($ll);
  
  $big_reply_string = $_[0];

  @reply_lines = split(/\n/, $big_reply_string);
  $num_reply_lines = 0;
  foreach $ll (@reply_lines) {
    $num_reply_lines++;
    if ($debug > 0) {
      printf "next ll is ..$ll.. \n";
    }
  }

  return($num_reply_lines, @reply_lines);
}


#########################################################################
# PROCEDURE: find_reply_line
#
# DESCRIPTION: Given an array of many lines containing the reply
#              from an HTTP request, search them for the (first)
#              line containing the given string.  
#
# RETURNS:
#              (retval, line)
#                    where     retval  =  0       if we find it
#                                         1       if we don't find it
#
#                              line    = the entire line containing the string
#
sub find_reply_line {

  my($search_string);
  my($num_lines);
  my($i);
  my($retval, $line);

  $search_string = $_[0];
  $num_lines = $_[1];
  # args 2 - N are the lines we'll search

  # initialize to "unsuccessful search"
  $retval = 1;

  for ($i = 0; $i < $num_lines; $i++) {
    if ($_[2 + $i] =~ /$search_string/) {
      $line = $_[2 + $i];
      $retval = 0;
      if ($debug > 0) {
        printf "find_reply_line: found ..$search_string.. in line \n";
        printf "    ..%s.. \n", $line;
      }
      return($retval, $line);
    }
  } 

  # if we get here, there was no match
  return($retval, "");
}


###########################################################################
# PROCEDURE: make_phased_graph
# 
# DESCRIPTION: Given information about a variable star, including its
#              period, and given a bunch of lines with Mark IV measurements,
#              create a nice plot of the phased light curve.
#
# RETURNS:
#              (proc_retval, plotfile)
#                    proc_retval       = 0 if all goes well, 1 if error
#                    plotfile          name of a file containing the plot
#
sub make_phased_graph {

  my($i);
  my($proc_retval, $plotfile);
  my($name, $ra, $dec, $epoch_jd, $period);
  my($num_data_lines, @data_lines);
  my($min_v_mag, $max_v_mag);
  my($min_i_mag, $max_i_mag);
  my($minx, $maxx, $miny, $maxy);
  my($plot_cmd_file, $plot_data_file);
  my($retval, $cmd);

  $proc_retval = 1;
  $plotfile = "./make_phased.png";
  $plot_cmd_file = "./make_phased.gnu";
  $plot_data_file = "./make_phased.dat";

  $name = $_[0];
  $ra = $_[1];
  $dec = $_[2];
  $epoch_jd = $_[3];
  $period = $_[4];
  $num_data_lines = $_[5];
  for ($i = 0; $i < $num_data_lines; $i++) {
    $data_lines[$i] = $_[6 + $i];
  }

  # if there are no data lines, we return with no graph
  if ($num_data_lines <= 0) {
    printf STDERR "make_phased_graph: no data to plot\n";
    return($proc_retval, $plotfile);
  }

  # if the epoch_jd value is invalid, just pick an arbitrary date for the
  #    phase = 0.  This can happen in the GCVS for some irregular variables
  if ($epoch_jd < 0) {
    $epoch_jd = 2452000.0;
  }
  # if the PERIOD is less than zero, then we can't make a phased graph;
  #    so return with error
  if ($period < 0) {
    return($proc_retval, $plotfile);
  }

  # calculate the phase for each observed magnitude, and put data into
  #   a file with lines that look like this:
  #
  #          JD   phase1 phase2  vmag1  vmag2  imag1 imag2
  #
  if (create_phased_data($plot_data_file, $epoch_jd, $period,
                                $num_data_lines, @data_lines) != 0) {
    printf STDERR "make_phased_graph: create_phased_data returns with error";
    return($proc_retval, $plotfile);
  }


  # figure out the min, max values of V-band magnitude
  ($min_v_mag, $max_v_mag) = find_extreme_mags("V", 
                              $num_data_lines, @data_lines);
  ($min_i_mag, $max_i_mag) = find_extreme_mags("I", 
                              $num_data_lines, @data_lines);
  if ($debug > 0) {
    printf "min_v_mag is %6.2f max_v_mag is %6.2f \n", $min_v_mag, $max_v_mag;
    printf "min_i_mag is %6.2f max_i_mag is %6.2f \n", $min_i_mag, $max_i_mag;
  }

  # figure out the min, max magnitude values on the graph
  ($minx, $maxx, $miny, $maxy) = find_graph_limits($min_v_mag,
                                  $max_v_mag, $min_i_mag, $max_i_mag);
  if ($debug > 0) {
    printf "graph limits: %7.2f %7.2f    %7.2f %7.2f  \n",
             $minx, $maxx, $miny, $maxy;
  }

  # open a file to hold GNUPLOT commands
  #   and put into it all the command necessary to make a nice plot
  open(CMD_FILE, ">$plot_cmd_file") || 
           die("make_phased_graph: can't open file $plot_cmd_file");
  printf CMD_FILE "set output '%s' \n", $plotfile;
  printf CMD_FILE "set term png color \n";
  printf CMD_FILE "set grid \n";
  printf CMD_FILE "set xlabel 'Phase (two periods shown for clarity)' \n";
  printf CMD_FILE "set ylabel 'Magnitude' \n";
  printf CMD_FILE "set title 'TASS data on %s using known period %.5f \n",
                          $name, $period;
  # here comes the big command, in which we plot all the data in both bands
  #   we build this very long line one bit at a time
  printf CMD_FILE "plot [%f:%f][%f:%f] '%s' using 2:4 t 'V good' lt 1 pt 1 ",
                         $minx, $maxx, $miny, $maxy, $plot_data_file;
  printf CMD_FILE " ,  '%s' using 3:4 t '' lt 1 pt 1 ",
                            $plot_data_file;
  printf CMD_FILE " , '%s' using 2:5 t 'V flag' lt 1 pt 4 ps 0.8 " ;
                            $plot_data_file;
  printf CMD_FILE " , '%s' using 3:5 t '' lt 1 pt 4 ps 0.8 " ;
                            $plot_data_file;
  printf CMD_FILE " ,  '%s' using 2:6 t 'I good' lt 2 pt 1 ",
                            $plot_data_file;
  printf CMD_FILE " ,  '%s' using 3:6 t '' lt 2 pt 1 ",
                            $plot_data_file;
  printf CMD_FILE " ,  '%s' using 2:7 t 'I flag' lt 2 pt 4 ps 0.8 ",
                            $plot_data_file;
  printf CMD_FILE " ,  '%s' using 3:7 t '' lt 2 pt 4 ps 0.8 ",
                            $plot_data_file;
  printf CMD_FILE "\n";
  close(CMD_FILE);
                               
  # execute the "gnuplot" program on the commands to create
  #    a file holding the graph
  $cmd = "gnuplot $plot_cmd_file";
  $retval = `$cmd`;
  if ($? != 0) {
    printf STDERR "make_phased_graph: gnuplot returns with error \n";
    return($proc_retval, $plotfile);
  }

  # all is done
  $proc_retval = 0;

  return($proc_retval, $plotfile);
}


#############################################################################
# PROCEDURE: find_extreme_mags
# 
# DESCRIPTION: Given a passband and an array of data lines from 
#              Mark IV database, find the min and max values of the
#              magnitude in the given band.  We have some default
#              values which are "sensible" in case there are 
#              no valid values.
#
# RETURNS:
#              (min_mag, max_mag)
#
sub find_extreme_mags {

  my($i);
  my($default_min_mag, $default_max_mag);
  my($bad_min_mag, $bad_max_mag);
  my($min_mag, $max_mag); 
  my($passband);
  my($mag, $mag_column);
  my($line, @words);
  my($num_lines);

  $passband = $_[0];
  $num_lines = $_[1];

  if ($debug > 0) {
    printf "find_extreme_mags: passband ..%s.. num_lines ..%s..  \n",
                 $passband, $num_lines;
  }

  if ($passband eq "V") {
    $mag_column = 7;
  } 
  elsif ($passband eq "I") {
    $mag_column = 11;
  }
  else {
    printf STDERR "find_extreme_mags: given bad passband $passband -- abort\n";
    exit(1);
  }

  # these will signal no real data if they persist
  $bad_min_mag = 100;
  $bad_max_mag = -100;

  # these are sensible defaults
  $default_min_mag = 6.0;
  $default_max_mag = 15.0;

  # find the min, max values in the data
  $min_mag = $bad_min_mag;
  $max_mag = $bad_max_mag;
  for ($i = 0; $i < $num_lines; $i++) {
    $line = " " . $_[2 + $i];
    @words = split(/\s+/, $line);
    $mag = $words[$mag_column];
    if ($debug > 0) {
      printf "find_extreme_mags: passband $passband line %3d mag %7.3f \n",
                       $i, $mag;
    }
    # we don't include values of 99 
    if ($mag > 90) {
      next;
    }
    if ($mag < $min_mag) {
      $min_mag = $mag;
    }
    if ($mag > $max_mag) {  
      $max_mag = $mag;
    }
  }
  if ($min_mag == $bad_min_mag) {
    $min_mag = $default_min_mag;
  }
  if ($max_mag == $bad_max_mag) {
    $max_mag = $default_max_mag;
  }

  return($min_mag, $max_mag);
}


#############################################################################
# PROCEDURE: find_graph_limits
# 
# DESCRIPTION: Given the min,max V and min,max I magnitudes of star,
#              figure out good limits for a phased light curve.
#              The "x" values are phase, and "y" values are magnitudes
#              (with min at top and max at bottom)
#
# RETURNS:
#              (minx, maxx, miny, maxy)
#
sub find_graph_limits {

  my($extra);
  my($min_v, $max_v, $min_i, $max_i);
  my($min_mag, $max_mag, $delta_mag);
  my($minx, $maxx, $miny, $maxy);

  $min_v = $_[0];
  $max_v = $_[1];
  $min_i = $_[2];
  $max_i = $_[3];
  if ($min_v < $min_i) {
    $min_mag = $min_v;
  } 
  else {
    $min_mag = $min_i;
  }
  if ($max_v > $max_i) {
    $max_mag = $max_v;
  } 
  else {
    $max_mag = $max_i;
  }

  # the limits in phase are simple
  $minx = -0.1;
  $maxx = 2.1;
  
  # the limits in magnitude are harder -- we allow a little extra
  #   room at top and bottom
  $delta_mag = $max_mag - $min_mag;
  $extra = $delta_mag*0.10;
  $miny = $max_mag + $extra;
  $maxy = $min_mag - 2.0*$extra;
 
  if ($debug > 0) {
    printf "find_graph_limits: X %6.2f %6.2f   Y %6.2f %6.2f \n",
                 $minx, $maxx, $miny, $maxy;
  }

  return($minx, $maxx, $miny, $maxy);
}



#############################################################################
# PROCEDURE: create_phased_data
# 
# DESCRIPTION: Given a filename, epoch, period, an array of data on 
#              a particular star, calculate the phase for each
#              time of measurement.  Create a datafile with format
#
#                    JD   phase1 phase+1  vmag1  vmag2  imag1 imag2
#
#              where      vmag1      is the V-band mag, if no warning flags
#                         vmag2      is the V-band mag, if warning flag set
#                         imag1      is the I-band mag, if no warning flags
#                         imag2      is the I-band mag, if warning flag set
#
#              We use a value of 99.0 for any mag which we're not going
#              to plot.
#  
#              Write the datafile to disk.
#
# RETURNS:
#              0        if all goes well
#              1        if an error occurs
#
sub create_phased_data {

  my($i);
  my($retval);
  my($datafile_name);
  my($epoch_jd, $period);
  my($num_lines);
  my($line, @words);
  my($jd, $vmag, $vsig, $vflag, $imag, $isig, $iflag);
  my($phase, $x);
  my($vmag1, $vmag2, $imag1, $imag2);

  # by default, we haven't succeeded
  $retval = 1;

  $datafile_name = $_[0];
  $epoch_jd = $_[1];
  $period = $_[2];
  $num_lines = $_[3];

  if ($debug > 0) {
    printf "create_phased_data: file ..%s.. epoch ..%12.5f.. period %12.5f \n",
                 $datafile_name, $epoch_jd, $period;
  }
  # sanity checks
  if ($epoch_jd < 0) {
    printf STDERR "create_phased_data: bad epoch %12.5f \n", $epoch_jd;
    return($retval);
  }
  if ($period <= 0) {
    printf STDERR "create_phased_data: bad period 12.5f \n", $period;
    return($retval);
  }
  
  open(PHASED_DATA, ">$datafile_name") || 
              die("create_phased_data: can't open $datafile_name");

  for ($i = 0; $i < $num_lines; $i++) {
    $line = " " . $_[4 + $i];
    @words = split(/\s+/, $line);
    $jd = $words[5];
    $vmag = $words[7];
    $vsig = $words[8];
    $vflag = $words[9];
    $imag = $words[11];  
    $isig = $words[12];
    $iflag = $words[13];

    if ($debug > 0) { 
       printf "  jd ..%s..  vmag ..%s..  imag ..%s.. \n", $jd, $vmag, $imag;
    }

    $x = ($jd - $epoch_jd)/$period;
    $phase = $x - floor($x);
    if ($debug > 0) {
       printf " phase %8.3f \n", $phase;
    }

    # look at flags, set the quantities we'll print
    if ($vflag == 0) {
      $vmag1 = $vmag;
      $vmag2 = 99.0;
    }
    else {
      $vmag1 = 99.0;
      $vmag2 = $vmag;
    }
    if ($iflag == 0) {
      $imag1 = $imag;
      $imag2 = 99.0;
    }
    else {
      $imag1 = 99.0;
      $imag2 = $imag;
    }
    
    # now we print a line into the output datafile
    printf PHASED_DATA " %12.5f %7.3f %7.3f  %7.3f %7.3f %7.3f %7.3f \n",
               $jd, $phase, $phase+1.0, $vmag1, $vmag2, $imag1, $imag2;

  }
  close(PHASED_DATA);

  # if we get here, all is well!
  $retval = 0;
  return($retval);
}