#!/usr/bin/perl
#
# Prepare a subset of the data for a run through the "photom"
#    photometric calibration program.  We need to pick two filters
#    and match up stars detected in those two filters from pairs
#    of images taken at nearly the same time.
#
# MWR 10/31/2009
use POSIX;

$debug = 1;
$DEGTORAD = (180.0/3.14159);

$filt1 = "V";
$filt2 = "I";
# matching radius, in arcseconds
$match_radius_arcsec = 3.0;
# log file with a list of all images and exposure times
$logfile = "../addheader.log.all";


# columns holding data from the .radec files.  All are one-indexed
$index_col = 1;
$ra_col = 2;
$dec_col = 3;
$jd_col = 13;
$airmass_col = 14;
# which magnitude do we use?  choose zero-indexed column index here
$mag_col = 8;
$magerr_col = $mag_col + 1;

$filt1_lc = lc($filt1);
$filt2_lc = lc($filt2);

# create a subdirectory for the collated data
$subdir = sprintf "collate_%s%s", $filt1, $filt2;
if (!-e $subdir) {
  mkdir($subdir);
}

# create an output file in that subdirectory
$output_file = sprintf "%s/collate_%s%s.out", $subdir, $filt1, $filt2;
open(OUTPUT_FILE, ">$output_file") || die("can't open file $output_file");

$total_matches = 0;

$globstr = sprintf "*[0-9][0-9][0-9]%s.pht.radec", $filt1;
@filt1_list = glob("$globstr");
foreach $filt1_file (@filt1_list) {
  if ($debug > 0) {
    printf " file for filt1 %s is %s \n", $filt1, $filt1_file;
  }

  # Step 1: look for a matching .radec file in filter 2
  $filt2_file = $filt1_file;
  for ($i = 1; $i < length($filt2_file); $i++) {
    if (substr($filt2_file, $i - 1, 1) eq "_") {
      if (substr($filt2_file, $i, 1) eq "$filt1_lc") {
        if (substr($filt2_file, $i + 1, 1) eq "-") {
          substr($filt2_file, $i, 1) = "$filt2_lc";
        }
      }
      else {
        if ($debug > 1) {
          printf " failed, next is %s not %s\n", 
            substr($filt2_file, $i, 1), $filt1_lc;
        }
      }
    }
  }
  if ($debug > 1) {
    printf " filt1_file %s filt2_file %s \n", $filt1_file, $filt2_file;
  }
  $filt2_file =~ s/${filt1}\.pht/${filt2}\.pht/;
  if ($debug > 1) {
    printf " filt1_file %s filt2_file %s \n", $filt1_file, $filt2_file;
  }

  if (!-e $filt2_file) {
    # failed to find a matching .radec file in the other filter
    if ($debug > 0) {
      printf " failed to find matching file for $filt1_file \n";
      next;
    }
  }
  else {
    # okay, we can move forward 
    if ($debug > 0) { 
      printf "   yes, matching file: %s %s \n", $filt1_file, $filt2_file;
    }
  }

  # Step 1.5: read in all information from each file into memory
  $num_filt1 = 0;
  if ($debug > 0) {
    printf "   reading data from file $filt1_file \n";
  }
  open(INFILE, "$filt1_file") || die("can't open file $filt1_file");
  while (<INFILE>) {
    $line = $_;
    if ($line =~ /^#/) {
      next;
    }
    chomp($line);
    $line = " " . $line;
    @words = split(/\s+/, $line);
    $filt1_index[$num_filt1] = $words[$index_col];
    $filt1_ra[$num_filt1] = $words[$ra_col];
    $filt1_dec[$num_filt1] = $words[$dec_col];
    $filt1_mag[$num_filt1] = $words[$mag_col];
    $filt1_magerr[$num_filt1] = $words[$magerr_col];
    $filt1_jd[$num_filt1] = $words[$jd_col];
    $filt1_airmass[$num_filt1] = $words[$airmass_col];
    if ($debug > 1) {
     printf " %5d %12.5f %12.5f  %6.3f %6.3f %12.5f %6.2f \n",
          $filt1_index[$num_filt1], 
          $filt1_ra[$num_filt1], 
          $filt1_dec[$num_filt1], 
          $filt1_mag[$num_filt1], 
          $filt1_magerr[$num_filt1], 
          $filt1_jd[$num_filt1], 
          $filt1_airmass[$num_filt1];
    }

    $num_filt1++;
  }
  close(INFILE);
  if ($debug > 0) {
    printf "     read %5d stars from file %s \n", $num_filt1, $filt1_file;
  }
  # figure out the exposure time for this image
  $filt1_exptime = -1;
  $base = $filt1_file;
  $base =~ s/.pht//;
  $base =~ s/.radec//;
  $cmd = "grep $base $logfile ";
$old_debug = $debug; $debug = 0;
  $ret = exec_cmd($cmd);
$debug = $old_debug;
  @words = split(/\s+/, $ret);
  $filt1_exptime = $words[4];
  if ($filt1_exptime <= 0) {
    printf "error: can't get exptime for filt1 file $filt1_file \n";
    exit(1);
  }
  if ($debug > 0) {
    printf "    $filt1_file has exptime %6.1f \n", $filt1_exptime;
  }
  

  $num_filt2 = 0;
  if ($debug > 0) {
    printf "   reading data from file $filt2_file \n";
  }
  open(INFILE, "$filt2_file") || die("can't open file $filt2_file");
  while (<INFILE>) {
    $line = $_;
    if ($line =~ /^#/) {
      next;
    }
    chomp($line);
    $line = " " . $line;
    @words = split(/\s+/, $line);
    $filt2_index[$num_filt2] = $words[$index_col];
    $filt2_ra[$num_filt2] = $words[$ra_col];
    $filt2_dec[$num_filt2] = $words[$dec_col];
    $filt2_mag[$num_filt2] = $words[$mag_col];
    $filt2_magerr[$num_filt2] = $words[$magerr_col];
    $filt2_jd[$num_filt2] = $words[$jd_col];
    $filt2_airmass[$num_filt2] = $words[$airmass_col];
    if ($debug > 1) {
     printf " %5d %12.5f %12.5f  %6.3f %6.3f %12.5f %6.2f \n",
          $filt2_index[$num_filt2], 
          $filt2_ra[$num_filt2], 
          $filt2_dec[$num_filt2], 
          $filt2_mag[$num_filt2], 
          $filt2_magerr[$num_filt2], 
          $filt2_jd[$num_filt2], 
          $filt2_airmass[$num_filt2];
    }

    $num_filt2++;
  }
  close(INFILE);
  if ($debug > 0) {
    printf "    read %5d stars from file %s \n", $num_filt2, $filt2_file;
  }
  # figure out the exposure time for this image
  $filt2_exptime = -1;
  $base = $filt2_file;
  $base =~ s/.pht//;
  $base =~ s/.radec//;
  $cmd = "grep $base $logfile ";
$old_debug = $debug; $debug = 0;
  $ret = exec_cmd($cmd);
$debug = $old_debug;
  @words = split(/\s+/, $ret);
  $filt2_exptime = $words[4];
  if ($filt2_exptime <= 0) {
    printf "error: can't get exptime for filt2 file $filt2_file \n";
    exit(1);
  }
  if ($debug > 0) {
    printf "    $filt2_file has exptime %6.1f \n", $filt2_exptime;
  }
  


    

  # Step 2: walk through the filt1 radec file.  For each star,
  #    we look for a matching star in the filt2 file.
  #    Use simple linear search, and if we find a match, stop, don't bother
  #        to look for closest match
  $num_match_in_file = 0;
  for ($i1 = 0; $i1 < $num_filt1; $i1++) {
    if ($debug > 1) {
      printf " looking for match for star1 %12.5f %12.5f \n", 
          $filt1_ra[$i1], $filt1_dec[$i1];
    }
    
    $found_match = 0;
    for ($i2 = 0; $i2 < $num_filt2; $i2++) {
      $dra = $filt1_ra[$i1] - $filt2_ra[$i2];
      $ddec = $filt1_dec[$i1] - $filt2_dec[$i2];
      $cosdec = cos(($filt1_dec[$i1]/$DEGTORAD));
      $dra *= $cosdec;
      $dist = sqrt($dra*$dra + $ddec*$ddec);
      $dist_arcsec = $dist*3600.0;
      if ($dist_arcsec <= $match_radius_arcsec) {
        if ($debug > 1) {
          printf " found match star2 %12.5f %12.5f  dist %6.2f \n", 
              $filt2_ra[$i2], $filt2_dec[$i2], $dist_arcsec;
        }
        $found_match = 1;
        last;
      }
   
    }

    if ($found_match == 0) {
      next;
    }

    # if we get here, we found a match.  Now we can print a single line
    #    to the output file with information on the star in the 
    #    two filters.
    
    $avg_ra = ($filt1_ra[$i1] + $filt2_ra[$i2]) / 2.0;
    $avg_dec = ($filt1_dec[$i1] + $filt2_dec[$i2]) / 2.0;
    $avg_jd = ($filt1_jd[$i1] + $filt2_jd[$i2]) / 2.0;
    $avg_air = ($filt1_airmass[$i1] + $filt2_airmass[$i2]) / 2.0;

    # scale the magnitudes in each filter to correspond to an exposure
    #    time of 1 second
    $mag1_adj = 2.5*log10($filt1_exptime);
    $mag1 = $filt1_mag[$i1] + $mag1_adj;
    if ($debug > 1) {
      printf " filt1 mag = %6.3f + %6.3f = %6.3f \n",
             $filt1_mag[$i1], $mag1_adj, $mag1;
    }
    $mag2_adj = 2.5*log10($filt2_exptime);
    $mag2 = $filt2_mag[$i2] + $mag2_adj;
    if ($debug > 1) {
      printf " filt2 mag = %6.3f + %6.3f = %6.3f \n",
             $filt2_mag[$i2], $mag2_adj, $mag2;
    }


    $quality_flag = 0;
    printf OUTPUT_FILE "%6d %12.5f %12.5f %12.5f %7.3f %s %6.3f %5.3f %s %6.3f %5.3f %2d \n",
        $total_matches,
        $avg_ra,
        $avg_dec,
        $avg_jd,
        $avg_air,
        $filt1,
        $mag1, 
        $filt1_magerr[$i1],
        $filt2,
        $mag2, 
        $filt2_magerr[$i2],
        $quality_flag;
        
    $num_match_in_file++;
    $total_matches++;
  }
  if ($debug > 0) {
    printf "    found %5d matches out of %5d %5d possible \n",
          $num_match_in_file, $num_filt1, $num_filt2;
  }





}



close(OUTPUT_FILE);





exit 0;





##############################################################################
# PROCEDURE: exec_cmd
#
# DESCRIPTION: Execute the given shell command line.  If the $debug flag
#              is set, we print to stdout the command line, and
#              also print out the result string.
#
# RETURNS:
#              the result of the command
#
#
sub exec_cmd {

  my($cmd, $ret);

  $cmd = $_[0];

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

  return($ret);
}