#!/usr/bin/perl
#
# Run the "photom" program to find the photometric solution
#    for a dataset and calibrate the stars in it.
#
# MWR 10/31/2009
use POSIX;


$debug  = 1;
$DEGTORAD = (3.14159/180.0);
$match_radius = 3.0;
$temp = "./run_photom.tmp";


# this is for the Henden sequence files I grabbed long ago
if (0 == 1) {
  $catalog = "./henden_bv_cut.catalog";
  $cat_ra_col = 1;
  $cat_dec_col = 2;
  $cat_nfilt = 2;
  $cat_filt1 = "V=3";
  $cat_filt2 = "I=4";
}
# this is for the AAVSO sequences
if (1 == 1) {
  # create a catalog by concatenating several reference files
  @refs = ("tass_gidel.cat", "tass_mrdel.cat", "tass_hodel.cat", "aavso_hodel.cat");
  $catalog = "./input.cat";
  $cmd = "/bin/rm -f $catalog";
  $ret = exec_cmd($cmd);
  foreach $ref (@refs) {
    $cmd = "cat $ref >> $catalog ";
    $ret = exec_cmd($cmd);
  }
  $cat_ra_col = 1;
  $cat_dec_col = 2;
  $cat_nfilt = 2;
  $cat_filt1 = "V=5";   
  $cat_filt2 = "I=6";  

  # remove stars with V > 13 from reference set
  $cmd = "awk '{ if (\$6 < 12.5) { print \$0 }}' $catalog > $temp ";
  $ret = exec_cmd($cmd);
  $cmd = "mv $temp $catalog ";
  $ret = exec_cmd($cmd);
}



$detected_file = "./collate_VI/collate_VI.out";
$det_ra_col = 1;
$det_dec_col = 2;
$det_jd_col = 3;
$det_air_col = 4;
$det_flag_col = 11;
$det_nfilt = 2;
$det_filt1 = "5,6,7";
$det_filt2 = "8,9,10";
$mode = "extinct";
$radius = 3.0;
$outfile = "photom_vi_";

$cmd  = "~/old/tass/photom/photom ";
$cmd .= " $catalog $cat_ra_col $cat_dec_col ";
$cmd .= " $cat_nfilt $cat_filt1 $cat_filt2 ";
$cmd .= " $detected_file $det_ra_col $det_dec_col ";
$cmd .= " $det_jd_col $det_air_col $det_flag_col ";
$cmd .= " $det_nfilt $det_filt1 $det_filt2 ";
$cmd .= " mode=$mode ";
$cmd .= " radius=$radius ";
$cmd .= " outfile=$outfile ";

$ret = exec_cmd($cmd);


# Now pick out the measurements of MR Del field, comp star "A"
$mag_col = 6;
$num = 0;
$sum = 0;
$sumsq = 0;
$want_ra = 307.8507;
$want_dec = 5.21389;
$calfile = $outfile . ".cal";
open(CALFILE, "$calfile") || die("can't open file $calfile");
while (<CALFILE>) {
  $line = $_;
  if ($line =~ /^#/) {
    next;
  }
  $line = " " . $line;
  chomp($line);
  if ($debug > 1) {
    printf " next output line is ..%s.. \n", $line;
  }
  @words = split(/\s+/, $line);
  $ra = $words[2];
  $dec = $words[3];
  
  $cosdec = cos($dec*$DEGTORAD);
  $dra = ($ra - $want_ra)*$cosdec;
  $ddec = ($dec - $want_dec);
  $dist = sqrt($dra*$dra + $ddec*$ddec);
  $dist_arcsec = $dist*3600.0;
  if ($dist_arcsec < $match_radius) {
    $num++;
    $sum += $words[$mag_col];
    $sumsq += $words[$mag_col]*$words[$mag_col];
  }
}
close(CALFILE);
if ($num < 1) {
  printf " not found ?! \n";
}
else {
  $mean = $sum / $num;
  if ($num == 1) {
    $stdev = 0.0;
  }
  else {
    $stdev = sqrt(($sumsq - $num*$mean*$mean) / ($num - 1.0));
  }
  printf " A N %5d  mean %6.3f  stdev %6.3f \n", $num, $mean, $stdev;
}


# Now pick out the measurements of MR Del field, comp star "B"
$num = 0;
$sum = 0;
$sumsq = 0;
$want_ra = 307.8481;
$want_dec = 5.30086;
$calfile = $outfile . ".cal";
open(CALFILE, "$calfile") || die("can't open file $calfile");
while (<CALFILE>) {
  $line = $_;
  if ($line =~ /^#/) {
    next;
  }
  $line = " " . $line;
  chomp($line);
  if ($debug > 1) {
    printf " next output line is ..%s.. \n", $line;
  }
  @words = split(/\s+/, $line);
  $ra = $words[2];
  $dec = $words[3];
  
  $cosdec = cos($dec*$DEGTORAD);
  $dra = ($ra - $want_ra)*$cosdec;
  $ddec = ($dec - $want_dec);
  $dist = sqrt($dra*$dra + $ddec*$ddec);
  $dist_arcsec = $dist*3600.0;
  if ($dist_arcsec < $match_radius) {
    $num++;
    $sum += $words[$mag_col];
    $sumsq += $words[$mag_col]*$words[$mag_col];
  }
}
close(CALFILE);
if ($num < 1) {
  printf " not found ?! \n";
}
else {
  $mean = $sum / $num;
  if ($num == 1) {
    $stdev = 0.0;
  }
  else {
    $stdev = sqrt(($sumsq - $num*$mean*$mean) / ($num - 1.0));
  }
  printf " B N %5d  mean %6.3f  stdev %6.3f \n", $num, $mean, $stdev;
}




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);
}