#!/usr/bin/perl
#
# Test of the ensemble photometry code.  This script should return 0
#   if all goes well, or 1 (plus print error messages) if not.
#
# Tests:
#     A:   ideal data, no variable stars: should find no differences
#     B:   ideal data, one (unknown) variable star: should find differences
#     C:   ideal data, one known variable star: should find no differences
#     D:   ideal data, one (unknown) bad image: should find differences
#     E:   ideal data, one known bad image: should find no differences
#
# MWR 9/12/2004

use POSIX;

# if set to 1, print out verbose messages as we go, plus leave all
#   intermediate files on disk for later inspection (rather than
#   deleting them all)
$debug = 0;


# by default, everything has passed so far ....
$retcode = 0;


# 
# Test A involves perfect data, no variable objects
#


#
# step 1A: create a set of artificial data
#

# first, we define the true properties for these stars
$mag_base = 10.0;
$mag_err = 0.05;
$max_x = 100.0;
$max_y = $max_x;
$nstar = 10;
$edge = $max_x/10.0;
for ($i = 0; $i < $nstar; $i++) {
  $input_x[$i] = $edge + ($max_x/(2.0*$nstar))*$i;
  $input_y[$i] = $edge + ((2.0*$max_y)/$nstar)*$i;
  if ($input_y[$i] > ($max_y - $edge)) {
		$input_y[$i] -= ($max_y - $edge);
  }
  $input_mag[$i] = $mag_base + $i*0.05;
  $input_err[$i] = $mag_err;
  if ($debug > 0) {
    printf "  input star %6d  %8.2f %8.2f  %7.4f %7.4f \n",
             $i, $input_x[$i], $input_y[$i], $input_mag[$i], $input_err[$i];
  }
}

# next, we define the number of images, and the date and shifts for each
$nimage = 10;
$max_shift = $edge;
$mag_offset = 0.10;
for ($i = 0; $i < $nimage; $i++) {
  $input_jd[$i] = 2450000.0 + 0.01*$i;
  $input_shift_x[$i] = rand($max_shift);
  $input_shift_y[$i] = rand($max_shift);
  $input_mag_offset[$i] = $i*$mag_offset;
  if ($debug > 0) {
    printf "  input image %6d  %12.5f %8.2f %8.2f  %7.4f \n",
             $i, $input_jd[$i],
             $input_shift_x[$i], $input_shift_y[$i], $input_mag_offset[$i];
  }
}

# now we can create files which have simulated observations:
#   one file per image, which contains the positions and mags of 
#   all the stars in that image.
$input_file_base = "./QQ_input";
for ($i_image = 0; $i_image < $nimage; $i_image++) {
  $filename = sprintf "%s.%03d\n", $input_file_base, $i_image;
  open(FILENAME, ">$filename") || die("can't open file $filename for writing");
  for ($i_star = 0; $i_star < $nstar; $i_star++) {
    $this_jd = $input_jd[$i_image];
    $this_x = $input_shift_x[$i_image] + $input_x[$i_star];
    $this_y = $input_shift_y[$i_image] + $input_y[$i_star];
    $this_mag = $input_mag_offset[$i_image] + $input_mag[$i_star];
    $this_err = $input_err[$i_star];
    printf FILENAME " %6d %12.5f %8.2f %8.2f %7.4f %7.4f \n", 
                    $i_star, $this_jd, $this_x, $this_y, $this_mag, $this_err;
  }
  close(FILENAME);
}
  

# 
# step 2A: run the "multipht" program to group data by star
#
$multiout = "./multipht.out";
$jdcol = 1;
$xcol = 2;
$ycol = 3;
$magcol = 4;
$errcol = 5;
# we need to run with auto flag 
@phtfiles = glob("${input_file_base}.*");
$cmd  = "./multipht jd=$jdcol x=$xcol y=$ycol mag=$magcol err=$errcol ";
$cmd .= " auto ";
$cmd .= " outfile=$multiout ";
foreach $phtfile (@phtfiles) {
  $cmd .= " $phtfile ";
}

$ret = run_cmd($cmd);


# 
# step 2A_check: check the multipht output file, to verify that it contains
#   the proper information.
#
if (multipht_is_okay($multiout, $nimage, $nstar) != 1) {
  printf STDOUT " multipht output file fails to pass test \n";
  $retcode++;
}
else {
  if ($debug > 0) {
    printf "multipht_is_okay passes all tests \n";
  }
}

#
# step 3A: run the "solvepht" program to find the ensemble solution
#
$solvepht_out = "./solvepht.out";
$solvepht_sig = "./solvepht.sig";
$solvepht_img = "./solvepht.img";
$cmd  = "./solvepht infile=$multiout  ";
$cmd .= " outfile=$solvepht_out ";
$cmd .= " sigfile=$solvepht_sig ";
$cmd .= " imfile=$solvepht_img ";

$ret = run_cmd($cmd);



#
# step 4A: compare ensemble solution output to input
#
# there are no variable stars in this dataset
$var_star_index = -1;
$var_image_index = -1;
$expect_to_fail = 0;
if (solvepht_is_okay($solvepht_out, $solvepht_img, $solvepht_sig, 
               $nimage, $nstar, $var_star_index, $var_image_index, 
               $expect_to_fail) != 1) {
  printf STDOUT " solvepht output file fails to pass Test A \n";
  $retcode++;
}
else {
  if ($debug > 0) {
    printf "solvepht_is_okay passes Test A \n";
  }
}


#
# Now, Tests B and C
#
#         We now tweak the magnitude values for one star in the set.
#         We will check to verify that we can detect the fact that
#         one star has varied relative to the rest.
#         Call this the dataset for "Tests B and C".
#
$var_star_index = 1;
for ($i_image = 0; $i_image < $nimage; $i_image++) {
  $var_mag_added[$i_image] = 0.20*sin($i_image/5.0);
}
# we now re-create files which have simulated observations:
#   one file per image, which contains the positions and mags of 
#   all the stars in that image.  But this time, we modify the
#   magnitude of the variable star slightly, so that it doesn't
#   match the expected sum of 
#       $input_mag_offset[$i_image] + $input_mag[$i_star];
# 
#          
$input_file_base = "./QQ_input";
for ($i_image = 0; $i_image < $nimage; $i_image++) {
  $filename = sprintf "%s.%03d\n", $input_file_base, $i_image;
  open(FILENAME, ">$filename") || die("can't open file $filename for writing");
  for ($i_star = 0; $i_star < $nstar; $i_star++) {
    $this_jd = $input_jd[$i_image];
    $this_x = $input_shift_x[$i_image] + $input_x[$i_star];
    $this_y = $input_shift_y[$i_image] + $input_y[$i_star];
    $this_mag = $input_mag_offset[$i_image] + $input_mag[$i_star];
    if ($i_star == $var_star_index) {
      $this_mag += $var_mag_added[$i_image];
    }
    $this_err = $input_err[$i_star];
    printf FILENAME " %6d %12.5f %8.2f %8.2f %7.4f %7.4f \n", 
                    $i_star, $this_jd, $this_x, $this_y, $this_mag, $this_err;
  }
  close(FILENAME);
}
  

# 
# step 2B: run the "multipht" program to group data by star
#
$multiout = "./multipht.out";
$jdcol = 1;
$xcol = 2;
$ycol = 3;
$magcol = 4;
$errcol = 5;
# we need to run with auto flag 
@phtfiles = glob("${input_file_base}.*");
$cmd  = "./multipht jd=$jdcol x=$xcol y=$ycol mag=$magcol err=$errcol ";
$cmd .= " auto ";
$cmd .= " outfile=$multiout ";
foreach $phtfile (@phtfiles) {
  $cmd .= " $phtfile ";
}

$ret = run_cmd($cmd);


# 
# step 2B_check: check the multipht output file, to verify that it contains
#   the proper information.
#
if (multipht_is_okay($multiout, $nimage, $nstar) != 1) {
  printf STDOUT " multipht output file fails to pass test \n";
  $retcode++;
}
else {
  if ($debug > 0) {
    printf "multipht_is_okay passes all tests \n";
  }
}

#
# step 3B: run the "solvepht" program to find the ensemble solution
#
$solvepht_out = "./solvepht.out";
$solvepht_sig = "./solvepht.sig";
$solvepht_img = "./solvepht.img";
$cmd  = "./solvepht infile=$multiout  ";
$cmd .= " outfile=$solvepht_out ";
$cmd .= " sigfile=$solvepht_sig ";
$cmd .= " imfile=$solvepht_img ";

$ret = run_cmd($cmd);



#
# step 4B: compare ensemble solution output to input
#
# there is one variable star this time (we set var_star_index above)
#   and we did NOT tell solvepht that it was a variable,
#   so it will throw off all the output magnitudes slightly.
#   This should cause our comparison of input to output to FAIL;
#   let us make sure it does FAIL
$expect_to_fail = 1;
if (solvepht_is_okay($solvepht_out, $solvepht_img, $solvepht_sig, 
                $nimage, $nstar, $var_star_index, $var_image_index,
                $expect_to_fail) == 1) {
  printf STDOUT " solvepht output file fails to pass Test B \n";
  $retcode++;
}
else {
  if ($debug > 0) {
    printf "solvepht_is_okay passes Test B \n";
  }
}



#
# step 3C: run the "solvepht" program to find the ensemble solution
#          but this time, we know there's a variable star in the
#          dataset.  So, we tell solvepht about it.
#
$solvepht_out = "./solvepht.out";
$solvepht_sig = "./solvepht.sig";
$solvepht_img = "./solvepht.img";
$cmd  = "./solvepht infile=$multiout  ";
$cmd .= " outfile=$solvepht_out ";
$cmd .= " sigfile=$solvepht_sig ";
$cmd .= " imfile=$solvepht_img ";
$cmd .= " varstar=$var_star_index ";

$ret = run_cmd($cmd);



#
# step 4C: compare ensemble solution output to input
#
# there is one variable star this time (we set var_star_index above)
#   and we DID tell solvepht that it was a variable,
#   so solvepht should not include it in the solutions.
#   This should cause our comparison of input to output to succeed,
#   as long as we don't compare the variable star.
$expect_to_fail = 0;
if (solvepht_is_okay($solvepht_out, $solvepht_img, $solvepht_sig, 
                $nimage, $nstar, $var_star_index, $var_image_index,
                $expect_to_fail) != 1) {
  printf STDOUT " solvepht output file fails to pass Test C \n";
  $retcode++;
}
else {
  if ($debug > 0) {
    printf "solvepht_is_okay passes Test C \n";
  }
}



#
# Now, Tests D and E
#
#         We now tweak the magnitude values for one image in the set.
#         We will check to verify that we can detect the fact that
#         one image has varied relative to the rest.
#         Call this the dataset for "Tests D and E".
#
$var_star_index = -1;
$var_image_index = $nimage-1;
for ($i_image = 0; $i_image < $nimage; $i_image++) {
  $var_mag_added[$i_image] = 0.10*sin($i_image/5.0);
}
# we now re-create files which have simulated observations:
#   one file per image, which contains the positions and mags of 
#   all the stars in that image.  But this time, we modify the
#   magnitudes of all stars in one image slightly, so that they doesn't
#   match the expected sum of 
#       $input_mag_offset[$i_image] + $input_mag[$i_star];
# 
#          
$input_file_base = "./QQ_input";
for ($i_image = 0; $i_image < $nimage; $i_image++) {
  $filename = sprintf "%s.%03d\n", $input_file_base, $i_image;
  open(FILENAME, ">$filename") || die("can't open file $filename for writing");
  for ($i_star = 0; $i_star < $nstar; $i_star++) {
    $this_jd = $input_jd[$i_image];
    $this_x = $input_shift_x[$i_image] + $input_x[$i_star];
    $this_y = $input_shift_y[$i_image] + $input_y[$i_star];
    $this_mag = $input_mag_offset[$i_image] + $input_mag[$i_star];
    if ($i_image == $var_image_index) {
      $this_mag += $var_mag_added[$i_image];
    }
    $this_err = $input_err[$i_star];
    printf FILENAME " %6d %12.5f %8.2f %8.2f %7.4f %7.4f \n", 
                    $i_star, $this_jd, $this_x, $this_y, $this_mag, $this_err;
  }
  close(FILENAME);
}
  

# 
# step 2D: run the "multipht" program to group data by star
#
$multiout = "./multipht.out";
$jdcol = 1;
$xcol = 2;
$ycol = 3;
$magcol = 4;
$errcol = 5;
# we need to run with auto flag 
@phtfiles = glob("${input_file_base}.*");
$cmd  = "./multipht jd=$jdcol x=$xcol y=$ycol mag=$magcol err=$errcol ";
$cmd .= " auto ";
$cmd .= " outfile=$multiout ";
foreach $phtfile (@phtfiles) {
  $cmd .= " $phtfile ";
}

$ret = run_cmd($cmd);


# 
# step 2D_check: check the multipht output file, to verify that it contains
#   the proper information.
#
if (multipht_is_okay($multiout, $nimage, $nstar) != 1) {
  printf STDOUT " multipht output file fails to pass test \n";
  $retcode++;
}
else {
  if ($debug > 0) {
    printf "multipht_is_okay passes all tests \n";
  }
}

#
# step 3D: run the "solvepht" program to find the ensemble solution
#
$solvepht_out = "./solvepht.out";
$solvepht_sig = "./solvepht.sig";
$solvepht_img = "./solvepht.img";
$cmd  = "./solvepht infile=$multiout  ";
$cmd .= " outfile=$solvepht_out ";
$cmd .= " sigfile=$solvepht_sig ";
$cmd .= " imfile=$solvepht_img ";

$ret = run_cmd($cmd);



#
# step 4D: compare ensemble solution output to input
#
# there is one bad image this time (we set var_image_index above)
#   and we did NOT tell solvepht that it was a bad image.
#   so it will throw off all the output magnitudes slightly.
#   This should cause our comparison of input to output to FAIL;
#   let us make sure it does FAIL
$expect_to_fail = 1;
if (solvepht_is_okay($solvepht_out, $solvepht_img, $solvepht_sig, 
                $nimage, $nstar, $var_star_index, $var_image_index, 
                $expect_to_fail) == 1) {
  printf STDOUT " solvepht output file fails to pass Test D \n";
  $retcode++;
}
else {
  if ($debug > 0) {
    printf "solvepht_is_okay passes Test D \n";
  }
}



#
# step 3E: run the "solvepht" program to find the ensemble solution
#          but this time, we know there's a bad image in the
#          dataset.  So, we tell solvepht about it.
#
$solvepht_out = "./solvepht.out";
$solvepht_sig = "./solvepht.sig";
$solvepht_img = "./solvepht.img";
$cmd  = "./solvepht infile=$multiout  ";
$cmd .= " outfile=$solvepht_out ";
$cmd .= " sigfile=$solvepht_sig ";
$cmd .= " imfile=$solvepht_img ";
$cmd .= " badim=$var_image_index ";

$ret = run_cmd($cmd);



#
# step 4E: compare ensemble solution output to input
#
# there is one bad image this time (we set var_image_index above)
#   and we DID tell solvepht that it was a bad image,
#   so solvepht should not include it in the solutions.
#   This should cause our comparison of input to output to succeed,
#   as long as we don't include the bad image.
$expect_to_fail = 0;
if (solvepht_is_okay($solvepht_out, $solvepht_img, $solvepht_sig, 
                $nimage, $nstar, $var_star_index, $var_image_index,
                $expect_to_fail) != 1) {
  printf STDOUT " solvepht output file fails to pass Test E \n";
  $retcode++;
}
else {
  if ($debug > 0) {
    printf "solvepht_is_okay passes Test E \n";
  }
}








#
# all done: clean up temporary files (if desired)
#
if ($debug == 0) {
  @files = glob("${input_file_base}*");
  foreach $file (@files) {
    unlink($file);
  }
  foreach $file ("multipht.out") {
    unlink($file);
  }
  foreach $file ("solvepht.out", "solvepht.img", "solvepht.sig") {
    unlink($file);
  }
}


if ($debug > 0) {
  printf "selftest.pl: final retcode is $retcode \n";
}

if ($retcode == 0) {
    printf "selftest.pl: passes all test okay \n";
}
else {
  printf "selftest.pl: fails some test(s) \n";
}


exit($retcode);



##########################################################################
# PROCEDURE: run_cmd
#
# DESCRIPTION: Run the given shell command and return its output.
#              If the global $debug variable is non-zero, then
#              print the command to stdout before running it,
#              and print the results of the command to stdout afterwards.
#
# RETURNS:
#      output of the command
#
sub run_cmd {

  $cmd = $_[0];

  if ($debug > 0) {
    printf " about to run ..$cmd.. \n";
  }
  $ret = `$cmd`;
  if ($debug > 0) {
    printf " result is ..$ret.. \n";
  }

  return($ret);
}


###########################################################################
# PROCEDURE: multipht_is_okay
#
# DESCRIPTION: Read the contents of the given multipht output file.
#              Verify that it contains information for all the
#              stars and all the images in the input set.
#
# RETURNS:
#     0        if the file is NOT okay
#     1        if it has the proper information and format
#

sub multipht_is_okay {

  my($multipht_file) = $_[0];
  my($nstar) = $_[1];
  my($nimage) = $_[2];
  
  if (open(MULTIPHT, $multipht_file) != 1) {
    printf STDOUT "multipht_is_okay: can't open file $multipht_file\n";
    return(0);
  }
  
  # the first line should look like this:
  #       N images     M stars
  # where N and M are the "nimage" and "nstar" values given as args
  $line = " " . <MULTIPHT>;
  @words = split(/\s+/, $line);
  if ($#words != 4) {
    printf STDOUT "multipht_is_okay: bad format initial line: \n";
    printf STDOUT "  ..$line.. \n";
    return(0);
  }
  if ($words[1] != $nimage) {
    printf STDOUT "multipht_is_okay: number of images should be $nimage, is $words[1] \n";
    printf STDOUT "  ..$line.. \n";
    return(0);
  }
  if ($words[3] != $nstar) {
    printf STDOUT "multipht_is_okay: number of stars should be $nstar, is $words[3] \n";
    printf STDOUT "  ..$line.. \n";
    return(0);
  }

  # now, we should next see "nimage" lines, one per image, which look 
  #    something like this each:
  #
  #        0          ./QQ_input.000    1.0      0.00000      9      9
  #
  for ($i = 0; $i < $nimage; $i++) {
    $line = " " . <MULTIPHT>;
    @words = split(/\s+/, $line);
    if ($#words != 6) {
      printf STDOUT "multipht_is_okay: bad image line: \n";
      printf STDOUT "  ..$line.. \n";
      return(0);
    }
    if ($words[1] != $i) {
      printf STDOUT "multipht_is_okay: image $i has ID $words[1] ?! \n";
      printf STDOUT "  ..$line.. \n";
      return(0);
    }
  }
   
  # after that comes information on individual stars.  We should see
  #   "nstar" stanzas.  Each stanza should have a header line,
  #   followed by lines for the magnitude in each image.
  #   The number of measurements is given by the final element of
  #   the header line (2 in example below).
  #   For example, 
  #
  #      6   29.86   54.07     2 
  #      6      0     29.86000     54.07000 0  10.1000   400.0
  #      6      1     29.86000     54.07000 0  10.2000   400.0
  #
  for ($i = 0; $i < $nstar; $i++) {
    $line = " " . <MULTIPHT>;
    @words = split(/\s+/, $line);
    if ($#words != 4) {
      printf STDOUT "multipht_is_okay: bad star header line: \n";
      printf STDOUT "  ..$line.. \n";
      return(0);
    }
    if ($words[1] != $i) {
      printf STDOUT "multipht_is_okay: star $i has ID $words[1] ?! \n";
      printf STDOUT "  ..$line.. \n";
      return(0);
    }
    # now we read in the measurements for this star
    $nmeasure = $words[4];
    for ($j = 0; $j < $nmeasure; $j++) {
      $line = " " . <MULTIPHT>;
      @words = split(/\s+/, $line);
      if ($words[1] != $i) {
        printf STDOUT "multipht_is_okay: bad measurement line for star $i ?!\n";
        printf STDOUT "  ..$line.. \n";
        return(0);
      }
    }
  }
   



  close(MULTIPHT);

  return(1);
}



###########################################################################
# PROCEDURE: solvepht_is_okay
#
# DESCRIPTION: Read the contents of the given solvepht output files.
#              Verify that they contain the proper information for all the
#              stars and all the images in the input set.
#
# RETURNS:
#     0        if the files are NOT okay
#     1        if they has the proper information and format
#

sub solvepht_is_okay {

  my($solvepht_out) = $_[0];
  my($solvepht_img) = $_[1];
  my($solvepht_sig) = $_[2];
  my($nstar) = $_[3];
  my($nimage) = $_[4];
  my($var_star_index) = $_[5];
  my($var_image_index) = $_[6];
  my($expect_to_fail) = $_[7];
  my($i);
  my($retcode) = 1;
  my(@words);
  my(@output_image_index);
  my(@output_image_name);
  my(@output_image_jd);
  my(@output_image_exptime);
  my(@output_image_zeropoint);
  my(@output_image_stdev);
  my(@output_image_stderrmean);
  my($input_delta, $output_delta);

  my(@output_mag_index);
  my(@output_mag_mean);
  my(@output_mag_stdev);
  my($num_output_mags);

  
  #
  # first, we check the "solvepht.img" output file.
  #
  if (open(SOLVEPHT_IMG, $solvepht_img) != 1) {
    printf STDOUT "solvepht_is_okay: can't open file $solvepht_img\n";
    return(0);
  }
  $num_output_images = 0;
  while (<SOLVEPHT_IMG>) {
    $line = " " . $_;
    @words = split(/\s+/, $line);
    # each line should have 7 words
    $num_words = 7;
    if ($#words != $num_words) {
      printf STDOUT "solvepht_is_okay: bad format in img file $solvepht_img\n";
      printf STDOUT " the following line has %d words, should have %d \n",
                       $#words, $num_words;
      return(0);
    }
      
    $output_image_index[$num_output_images] = $words[1];
    $output_name[$num_output_images] = $words[2];
    $output_image_exptime[$num_output_images] = $words[3];
    $output_image_jd[$num_output_images] = $words[4];
    $output_image_zeropoint[$num_output_images] = $words[5];
    $output_image_stdev[$num_output_images] = $words[6];
    $output_image_stderrmean[$num_output_images] = $words[7];

    $num_output_images++;
  }
  close(SOLVEPHT_IMG);

  # we should read one line for each (not-known-to-be-bad) image ...
  if (($var_image_index != -1) && ($expect_to_fail == 0)) {
    # in this case, we should have $nimage-1 good images in our set
    if ($num_output_images != ($nimage-1)) {
      printf STDOUT "solvepht_is_okay: num images %d should be %d \n",
                     $num_output_images, $nimage-1;
      return(0);
    }
  }
  else {
    # this is the ordinary case: expect $nimage images in the output
    if ($num_output_images != $nimage) {
      printf STDOUT "solvepht_is_okay: num images %d should be %d \n",
                   $num_output_images, $nimage;
      return(0);
    }
  }
  if ($debug > 0) {
    printf "solvepht_is_okay: read %d images in output \n", 
                   $num_output_images;
  }

  # and the properties of images on output should be the same as
  #   the properties between images on input
  for ($i = 0; $i < $num_output_images; $i++) {
    if ($output_image_jd[$i] != $input_jd[$i]) {
      printf STDOUT "solvepht_is_okay: img %d has jd %lf should be %lf \n",
                $i, $output_image_jd[$i], $input_jd[$i];
      $retcode = 0;
    }
    # need to look at the _differences_ between image zeropoints,
    #   as the absolute value may float
    $input_delta = $input_mag_offset[$i] - $input_mag_offset[0];
    $output_delta = $output_image_zeropoint[$i] - $output_image_zeropoint[0];
    if (is_equal($output_delta, $input_delta) != 1) {
      if ($expect_to_fail == 0) {
        printf STDOUT "solvepht_is_okay: img %d zeropoint is %lf should be %lf \n", 
                 $i, $output_delta, $input_delta;
      }
      else {
        if ($debug > 0) {
          printf STDOUT "solvepht_is_okay: img %d zeropoint is %lf not %lf okay \n",
                 $i, $output_delta, $input_delta;
        }
      }
      $retcode = 0;
    }
   
  }
  if ($debug > 0) {
    printf "solvepht_is_okay: finished with zeropoints of images \n";
  }




  #
  # next, we check the "solvepht.sig" output file.
  #
  if (open(SOLVEPHT_SIG, $solvepht_sig) != 1) {
    printf STDOUT "solvepht_is_okay: can't open file $solvepht_sig\n";
    return(0);
  }
  $num_output_mags = 0;
  while (<SOLVEPHT_SIG>) {
    $line = " " . $_;
    @words = split(/\s+/, $line);
    # each line should have 3 words
    $num_words = 3;
    if ($#words != $num_words) {
      printf STDOUT "solvepht_is_okay: bad format in sig file $solvepht_sig\n";
      printf STDOUT " the following line has %d words, should have %d \n",
                       $#words, $num_words;
      return(0);
    }
      
    $output_mag_index[$num_output_mags] = $words[1];
    $output_mag_mean[$num_output_mags] = $words[2];
    $output_mag_stdev[$num_output_mags] = $words[3];

    $num_output_mags++;
  }
  close(SOLVEPHT_SIG);
  # we should read one line for each star ...
  if ($num_output_mags != $nstar) {
    fprintf STDOUT "solvepht_is_okay: num stars %d should be %d \n",
                   $num_output_mags, $nstar;
    return(0);
  }
  if ($debug > 0) {
    printf "solvepht_is_okay: read %d output mags\n",
                   $num_output_mags;
  }
  # and the properties of stars on output should be the same as
  #   the properties of stars on input
  for ($i = 0; $i < $num_output_mags; $i++) {
    # need to look at the _differences_ between star mags,
    #   as the absolute value may float
    $input_delta = $input_mag[$i] - $input_mag[0];
    $output_delta = $output_mag_mean[$i] - $output_mag_mean[0];
    
    # now, if this star is the variable star, the values should not be equal
    #   The variable star (if it exists) has index value "$var_star_index".
    #   If $var_star_index == -1, there is no variable star.
    if ($i != $var_star_index) {
      # the magnitudes should match .... 
      if (is_equal($output_delta, $input_delta) != 1) {
        printf STDOUT "solvepht_is_okay: star %d mag is %lf should be %lf \n", 
                   $i, $output_delta, $input_delta;
        $retcode = 0;
      }
    }
    else {
      # this star is variable -- output should be different
      if (is_equal($output_delta, $input_delta) == 1) {
        printf STDOUT "solvepht_is_okay: star %d delta is %lf = %lf ?! \n", 
                   $i, $output_delta, $input_delta;
        $retcode = 0;
      }
    }

  }
  if ($debug > 0) {
    printf "solvepht_is_okay: finished with output mags \n";
  }



  return($retcode);
}



##########################################################################
# PROCEDURE: is_equal
#
# DESCRIPTION: Compare two given numbers, to see if they
#              are the same, or close enough.
#
# RETURNS:
#     0        if the values ARE equal (or close enough)
#     1        if the values are NOT equal
#

sub is_equal {

  my($val1) = $_[0];
  my($val2) = $_[1];
  my($tolerance) = 0.001;
  my($delta);

  $delta = fabs($val1 - $val2);
  if ($delta > $tolerance) {
    return(0);
  }
  else {
    return(1);
  }

}

