#!/usr/bin/perl
#
# Create some data we can use to test the "collate" and "photom" programs.
# Then run those programs to make sure that they are working properly.
#
# Return with code 0 if all goes well,
#                  1 if an error(s) occurs -- and print message
#
# MWR 6/9/2001
use strict;

my($debug);
my($nstar);
my($nimages);
my($i_image);
my(@coeff_a_v, $coeff_b_v);
my(@coeff_a_i, $coeff_b_i);
my(@list_stem);
my($catalog_file);
my($retcode);

# if this is set to 1, then we print some diagnostic messages during
#   execution and don't delete temporary files
$debug = 0;


# number of stars
$nstar = 200;
# number of different pairs of images
$nimages = 5;
# zero-point offset coefficients
for ($i_image = 0; $i_image < $nimages; $i_image++) {
  $coeff_a_v[$i_image] = 1.0 + $i_image*0.1;
  $coeff_a_i[$i_image] = 1.1 + $i_image*0.1;
}
# color-term coefficients
$coeff_b_v = 0.05;
$coeff_b_i = 0.10;
# name of catalog file
$catalog_file = "catalog.dat";
# stem of file name for raw star lists
for ($i_image = 0; $i_image < $nimages; $i_image++) {
  $list_stem[$i_image] = "hxra202870$i_image.fits.ast";
}

$retcode = 0;

# generate a set of artificial stars, and place lists of their properties
#   into output file(s) which can be used as input to "collate"
if (generate_stars($nstar, $nimages, @coeff_a_v, @coeff_a_i, 
                       $coeff_b_v, $coeff_b_i, 
                       $catalog_file, @list_stem) != 0) {
  die("selftest.pl: failed to generate artificial star data");
}


# run the "collate" program to combine the two lists of raw photometry
#   into a single file which contains both V and I data
if (run_collate($nimages, @list_stem) != 0) {
  die("selftest.pl: collate failed to run properly");
}


# run the "photom" program to create a photometric solution for the 
#   night, using the collated output.  
if (run_photom($nimages, $catalog_file, @list_stem) != 0) {
  die("selftest.pl: run_photom failed to run properly");
}


# check the results of the "photom" program
if (check_results($nimages, @coeff_a_v, @coeff_a_i, $coeff_b_v, $coeff_b_i) != 0) {
  die("selftest.pl: check_results finds error in photometric solution");
}

# clean up
if (clean_up($nimages, $catalog_file, @list_stem) != 0) {
  die("selftest.pl: failed to clean up properly");
}

exit($retcode);



##############################################################################
# Generate data on artificial stars, and place it into raw star lists
#   which we can use as input to the "collate" program.
#
# We don't QUITE generate the magnitudes properly, so our color terms
#   are going to be a little bit different from the expected color terms --
#   but we should be in the ballpark.
#
# RETURNS:
#   0          if all goes well
#   1          if there's an error
#
sub generate_stars {
  my($nstar, $nimages);
  my(@coeff_a_v, @coeff_a_i, $coeff_b_v, $coeff_b_i, @list_stem);
  my($i_image);
  my($i, $id, $ra, $dec, $vmag, $imag, $color, $scatter);
  my($ra_min, $dec_min, $pos_spacing);
  my($mag_spacing, $vmag_sig, $imag_sig);
  my($output_file_v, $output_file_i);
  my($sky, $skysig, $quality);
  my($row, $col, $row_spacing);
  my($peak, $fwhm, $round, $sharp);
  my($catalog_v, $catalog_i, $catalog_file);
  my($this_v, $this_i);

  $nstar = $_[0];
  $nimages = $_[1];
  for ($i = 0; $i < $nimages; $i++) {
    $coeff_a_v[$i] = $_[2 + $i];
    $coeff_a_i[$i] = $_[2 + $nimages + $i];
    if ($debug > 0) {
      printf "generate_stars: coeff_a_v $i is $coeff_a_v[$i] \n";
      printf "generate_stars: coeff_a_i $i is $coeff_a_i[$i] \n";
    }
  }
  $coeff_b_v = $_[2 + 2*$nimages];
  $coeff_b_i = $_[3 + 2*$nimages];
  $catalog_file = $_[4 + 2*$nimages];
  for ($i = 0; $i < $nimages; $i++) {
    $list_stem[$i] = $_[5 + 2*$nimages + $i];
  }


  $ra_min = 180.0;
  $dec_min = 1.0;

 for ($i_image = 0; $i_image < $nimages; $i_image++) {

  # open the four output files: 2 for .coo, 2 for .ast
  # first, the .ast files
  $output_file_v = $list_stem[$i_image];
  $output_file_v =~ s/x/v/;
  $output_file_i = $list_stem[$i_image];
  $output_file_i =~ s/x/i/;
  open(OUTPUT_V_AST, ">$output_file_v") || 
           die("generate_stars: can't open file $output_file_v for writing");
  open(OUTPUT_I_AST, ">$output_file_i") || 
           die("generate_stars: can't open file $output_file_i for writing");
  # now, the .coo files
  $output_file_v =~ s/.ast/.coo/;
  $output_file_i =~ s/.ast/.coo/;
  open(OUTPUT_V_COO, ">$output_file_v") || 
           die("generate_stars: can't open file $output_file_v for writing");
  open(OUTPUT_I_COO, ">$output_file_i") || 
           die("generate_stars: can't open file $output_file_i for writing");

  # open the catalog file
  if ($i_image == 0) {
    open(CATALOG, ">$catalog_file") ||
           die("generate_stars: can't open file $catalog_file for writing");
  }



  # initialize properties of the first star -- we'll build from it
  $ra = $ra_min;
  $dec = $dec_min;
  $pos_spacing = 10.0/3600.0;
  $vmag = 10.0;
  $imag = 10.0;
  $mag_spacing = 0.01;
  $sky = 1000.0;
  $skysig = 1.0;
  $vmag_sig = 0.01;
  $imag_sig = 0.01;
  $quality = 0;
  $row = 10;
  $col = 10;
  $row_spacing = 10;
  $peak = 1000.0;
  $fwhm = 3.0;
  $round = 0.0;
  $sharp = 0.5;

  for ($i = 0; $i < $nstar; $i++) {
    $id = $i + 1;

    $color = ($vmag - $imag);
    $this_v = $vmag - $coeff_a_v[$i_image] - $coeff_b_v*($color);
    $this_i = $imag - $coeff_a_i[$i_image] - $coeff_b_i*($color);
       
    # print this star to the output files
    printf OUTPUT_V_AST "%5d %9.5f %9.5f  %4d %5.2f  %6.3f %6.3f  %2d\n",
            $id, $ra, $dec, $sky, $skysig, $this_v, $vmag_sig, $quality;
    printf OUTPUT_I_AST "%5d %9.5f %9.5f  %4d %5.2f  %6.3f %6.3f  %2d\n",
            $id, $ra, $dec, $sky, $skysig, $this_i, $imag_sig, $quality;

    printf OUTPUT_V_COO "%5d %7.2f %7.2f  %5d %6.3f  %7.3f %6.3f \n",
            $id, $row, $col, $peak, $fwhm, $round, $sharp;
    printf OUTPUT_I_COO "%5d %7.2f %7.2f  %5d %6.3f  %7.3f %6.3f \n",
            $id, $row, $col, $peak, $fwhm, $round, $sharp;

    # and, every now and then, put this star into the catalog file

    if (($i_image == 0) && ($i % 10 == 0)) {
      $catalog_v = $vmag;
      $catalog_i = $imag;
      printf CATALOG "%9.5f %9.5f  %6.3f  %6.3f \n", $ra, $dec, 
               $catalog_v, $catalog_i;
    } 

    # move onto the properties of the next star
    $ra += $pos_spacing;
    $dec += $pos_spacing;
    $vmag += $mag_spacing;
    $imag -= $mag_spacing;
    $row += $row_spacing;
    $col += $row_spacing;
  }

  if ($i_image == 0) {
    close(CATALOG);
  }

 }

  close(OUTPUT_V_AST);
  close(OUTPUT_I_AST);
  close(OUTPUT_V_COO);
  close(OUTPUT_I_COO);


  return(0);
}



#############################################################################
# Run the "collate" program on the test data files to produce a set
#   of merged data file, which may then be used as input for the
#   "photom" program.
#
# RETURNS
#    0          if all goes well
#    1          if there's an error
#
sub run_collate {
  my($nimages, $i_image);
  my(@list_stem);
  my($matchrad);
  my($jd);
  my($latitude);
  my($longitude);
  my($exptime);
  my($outfile);
  my($ast_file_v, $coo_file_v, $ast_file_i, $coo_file_i);
  my($pass_v, $pass_i);
  my($cmd_str);


  $nimages = $_[0];
  for ($i_image = 0; $i_image < $nimages; $i_image++) {
    $list_stem[$i_image] = $_[1 + $i_image];
  }

  # now, we generate one .clt file per pair of images
  for ($i_image = 0; $i_image < $nimages; $i_image++) {
  
    $outfile = get_collate_filename($list_stem[$i_image]);
    if ($debug > 0) {
      printf "run_collate: outfile will be ..$outfile..\n";
    }
    $pass_v = "V";
    $ast_file_v = $list_stem[$i_image];
    $ast_file_v =~ s/x/v/;
    $coo_file_v = get_coo_filename($ast_file_v);
    $pass_i = "I";
    $ast_file_i = $list_stem[$i_image];
    $ast_file_i =~ s/x/i/;
    $coo_file_i = get_coo_filename($ast_file_i);

    # match radius in arcsec
    $matchrad = 1.0;
    $jd = 2452029.700 + $i_image*0.01;
    $latitude = 45.0;
    $longitude = 90.0;
    $exptime = 1.0;

    # run the collate program
    $cmd_str = "./collate passband=$pass_v $coo_file_v $ast_file_v ";
    $cmd_str = $cmd_str . " passband=$pass_i $coo_file_i $ast_file_i ";
    $cmd_str = $cmd_str . " matchrad=$matchrad jd=$jd ";
    $cmd_str = $cmd_str . " lat=$latitude long=$longitude ";
    $cmd_str = $cmd_str . " exptime=$exptime outfile=$outfile ";
    if ($debug > 0) {
      printf "run_collate: cmd_str is ..$cmd_str..\n";
    }
  
    $retcode = system($cmd_str);
    if ($debug > 0) {
      printf "retcode is ..$retcode..\n";
    }
    if ($retcode != 0) {
      warn "run_collate: the collate program fails on test data";
      return(1);
    }

  }

  return(0);
}


##############################################################################
# Convert the filename stem to the appropriate .coo file name.
#
# RETURNS
#   the .coo output file name
#
sub get_coo_filename {
  my($list_stem, $outfile);

  $list_stem = $_[0];

  $outfile = $list_stem;
  $outfile =~ s/.fits.ast/.fits.coo/;

  return($outfile);
}




##############################################################################
# Convert the filename stem to the appropriate collate output file name.
#
# RETURNS
#   the collated output file name
#
sub get_collate_filename {
  my($list_stem, $outfile);

  $list_stem = $_[0];

  $outfile = $list_stem;
  $outfile =~ s/hx/Mh/;
  $outfile =~ s/.fits.ast/.clt/;

  return($outfile);
}



#############################################################################
# Delete temporary files.
#
# RETURNS
#    0          if all goes well
#    1          if there's an error
#
sub clean_up {
  my($nimages, $i_image);
  my($catalog_file);
  my(@list_stem);
  my($file);

  $nimages = $_[0];
  $catalog_file = $_[1];
  for ($i_image = 0; $i_image < $nimages; $i_image++) {
    $list_stem[$i_image] = $_[2 + $i_image];
  }

  if ($debug == 1) {
    printf "clean_up: doing nothing, since debug = 1\n";
    return(0);
  }

  # delete the catalog file
  unlink($catalog_file);

  # delete the temporary .ast and .coo files
  for ($i_image = 0; $i_image < $nimages; $i_image++) {
    $file = $list_stem[$i_image];
    $file =~ s/x/v/;
    unlink($file);
    $file = get_coo_filename($file);
    unlink($file);

    $file = $list_stem[$i_image];
    $file =~ s/x/i/;
    unlink($file);
    $file = get_coo_filename($file);
    unlink($file);

    # delete the collated output file
    $file = get_collate_filename($list_stem[$i_image]);
    unlink($file);
  }

  # delete all the files created by the photom program
  $file = "photom.coeff";
  unlink($file);
  $file = "photom_V.comp";
  unlink($file);
  $file = "photom_I.comp";
  unlink($file);
  $file = "photom.cal";
  unlink($file);

  return(0);
}


#############################################################################
# Run the "photom" program on collated output file(s);
#   if all goes well, this will produce a set of output files 
#   we'll check later
#
# RETURNS
#    0          if all goes well
#    1          if there's an error
#
sub run_photom {
  my($nimages, $i_image);
  my($catalog_file);
  my(@list_stem);
  my(@infile);
  my($cat_racol, $cat_deccol, $cat_nfilt, $cat_vcols, $cat_icols);
  my($det_racol, $det_deccol, $det_nfilt);
  my($det_jdcol, $det_aircol, $det_flagcol);
  my($det_vcols, $det_icols);
  my($matchrad);
  my($cmd_str);


  $nimages = $_[0];
  $catalog_file = $_[1];
  for ($i_image = 0; $i_image < $nimages; $i_image++) {
    $list_stem[$i_image] = $_[2 + $i_image];
  }
  

  for ($i_image = 0; $i_image < $nimages; $i_image++) {
    $infile[$i_image] = get_collate_filename($list_stem[$i_image]);
    if ($debug > 0) {
      printf "run_photom: collate input file is ..$infile[$i_image]..\n";
    }
  }

  $cat_racol = 0;
  $cat_deccol = 1;
  $cat_nfilt = 2;
  $cat_vcols = "V=2";
  $cat_icols = "I=3";

  $det_racol = 1;
  $det_deccol = 2;
  $det_jdcol = 3;
  $det_aircol = 4;
  $det_flagcol = 11;
  $det_nfilt = 2;
  $det_vcols = "5,6,7";
  $det_icols = "8,9,10";

  # match radius in arcsec
  $matchrad = 1.0;

  # run the photom program
  $cmd_str = "./photom $catalog_file $cat_racol $cat_deccol";
  $cmd_str = $cmd_str . " $cat_nfilt $cat_vcols $cat_icols";
  $cmd_str = $cmd_str . " $infile[0] ";
  $cmd_str = $cmd_str . " $det_racol $det_deccol $det_jdcol ";
  $cmd_str = $cmd_str . " $det_aircol $det_flagcol ";
  $cmd_str = $cmd_str . " $det_nfilt $det_vcols $det_icols";
  for ($i_image = 1; $i_image < $nimages; $i_image++) {
    $cmd_str = $cmd_str . " $infile[$i_image] ";
  }
  $cmd_str = $cmd_str . " radius=$matchrad > /dev/null";
  if ($debug > 0) {
    printf "run_photom: cmd_str is ..$cmd_str..\n";
  }

  $retcode = system($cmd_str);
  if ($debug > 0) {
    printf "retcode is ..$retcode..\n";
  }
  if ($retcode != 0) {
    warn "run_photom: the photom program fails on test data";
    return(1);
  }

  return(0);
}


##############################################################################
# Check the results of the photometric solution on the test data.
#   We look only at the values in the "photom.coeff" file, ignoring
#   all the other files we _could_ check.
#
# We pass the expected solution values to this routine:
#   arg 0       # of image pairs
#       1       zero-pt in V for first pair
#       2       zero-pt in V for second pair
#      ...
#       N       zero-pt in V for N'th pair
#      N+1      zero-pt in I for first pair
#      N+2      zero-pt in I for second pair
#      ...
#      2N       zero-pt in I for N'th pair
#     2N+1      color term for V (is same for all image pairs)
#     2N+2      color term for I (is same for all image pairs)
#
# RETURNS:
#    0            if all goes well
#    1            if there's an error
#
sub check_results {
  my($nimages, $i_image);
  my(@coeff_a_v, @coeff_a_i, $coeff_b_v, $coeff_b_i);
  my($coeff_file);
  my($nlines);
  my($line, @words);
  my($solution_a, $solution_b, $solution_rms);
  my($retcode);

  $nimages = $_[0];
  for ($i_image = 0; $i_image < $nimages; $i_image++) {
    $coeff_a_v[$i_image] = $_[1 + $i_image];
    $coeff_a_i[$i_image] = $_[1 + $nimages + $i_image];
  }
  $coeff_b_v = $_[2*$nimages + 1];
  $coeff_b_i = $_[2*$nimages + 2];
  if ($debug > 0) {
    printf "check_results: expect solutions with \n";
	 for ($i_image = 0; $i_image < $nimages; $i_image++) {
	   printf "V: pair %2d  a %6.3f  b %6.3f \n", 
                              $i_image, $coeff_a_v[$i_image], $coeff_b_v;
	   printf "I: pair %2d  a %6.3f  b %6.3f \n", 
                              $i_image, $coeff_a_i[$i_image], 0.0-$coeff_b_i;
    }
  }

  $retcode = 0;

  $coeff_file = "photom.coeff";

  # we read the results from the "photom.coeff" file
  if (open(COEFF, $coeff_file) == undef) {
    warn("check_results: can't open file $coeff_file with photometric results");
    return(1);
  }

  # there should be two lines in the COEFF file for each image:
  #   one for solution in V, the other for solution in I. 
  #   The lines are arranged like so:
  #        V sol for image pair 0
  #        V sol for image pair 1
  #           ....
  #        I sol for image pair 0
  #        I sol for image pair 1
  #           ....
  $nlines = 0;
  while (<COEFF>) {
    $line = $_;
    @words = split(/\s+/, $line);
    $solution_a = $words[5];
    $solution_b = $words[8];
    $solution_rms = $words[11];
  
    if ($debug > 0) {
      printf "soln %3d: a %6.3f b %6.3f  RMS %6.3f \n", 
                 $nlines, $solution_a, $solution_b, $solution_rms;
    }

    # Check results for V
    if ($nlines < $nimages) {
	   $i_image = $nlines;
      if (are_the_same($coeff_a_v[$i_image], $solution_a) == 0) {
        printf "check_results: pair %2d V-band zero point is %6.3f   should be %6.3f \n",
	                   $i_image, $solution_a, $coeff_a_v[$i_image];
        $retcode = 1;
      }
      if (are_the_same($coeff_b_v, $solution_b) == 0) {
        printf "check_results: pair %2d V-band color term is %6.3f   should be %6.3f \n",
	                   $i_image, $solution_b, $coeff_b_v;
        $retcode = 1;
      }
    }

    # Check results for I
    if ($nlines >= $nimages) {
	   $i_image = $nlines - $nimages;
      if (are_the_same($coeff_a_i[$i_image], $solution_a) == 0) {
        printf "check_results: pair %2d I-band zero point is %6.3f   should be %6.3f \n",
	                   $i_image, $solution_a, $coeff_a_i[$i_image];
        $retcode = 1;
      }
      # Note that due to the way we calculate colors, the I-band solution
      #   for color term should be equal to the NEGATIVE of the input
      #   coeff_b_i
      $solution_b = 0 - $solution_b;
      if (are_the_same($coeff_b_i, $solution_b) == 0) {
        printf "check_results: pair %2d I-band color term is %6.3f   should be %6.3f \n",
	                   $i_image, $solution_b, $coeff_b_i;
        $retcode = 1;
      }
    }

    # This should never happen
    if ($nlines >= 2*$nimages) {
      warn "check_results: too many lines in file $coeff_file?";
      $retcode = 1;
    }

    $nlines++;
  }

  return($retcode);
}


##############################################################################
# Compare two numbers, and decide if they are "close enough" 
#   to denote a successful execution.
#
# RETURNS
#     0                if the two are NOT the same
#     1                if the two are close enough
#
sub are_the_same {
  my($x, $y, $diff, $acceptable_diff);

  $x = $_[0];
  $y = $_[1];

  $acceptable_diff = 0.10;

  $diff = $x - $y;
  if (abs($diff) > $acceptable_diff) {
    return(0);
  }
  else {
    return(1);
  }

}