#!/usr/bin/perl
#
# Do all the stuff needed to run "match" on a .pht file and 
#   a catalog from the USNOB1.0.  
#
# MWR 3/12/2005
#
# Modified for use with RIT Obs measurements of WZ Sge field
#   and Henden reference photometry file
#
# Modified for use with SN in M74 images from RIT Obs
#
#    usage:  do_astrom.pl  file.pht
#
# MWR 7/27/2013


$debug = 1;

$DEGTORAD = (3.141592653589793/180.0);

if ($#ARGV != 0) {
  printf " usage: do_astrom.pl file.pht ";
  exit(1);
}
$phtfile = $ARGV[0];
if ($debug > 0) {
  printf "phtfile is $phtfile \n";
}

#$catfile = "./gettys_usnob1.dat";
#$cat_ra = 111.59949;
#$cat_dec = 0.77708;
$catfile = "./usno_hand.cat";
$cat_ra = 24.1739458;
$cat_dec = 15.7836619;
$cat_idcol = 0;
$cat_racol = 9;
$cat_deccol = 10;
# mag column is 6 for V
$cat_magcol = 17;

#$phtfile = "./tr08.pht";
#$phtfile = "./sn_i-030.pht";
$pht_idcol = 0;
$pht_xcol = 1;
$pht_ycol = 2;
$pht_magcol = 5;


# run "project_coords" on the catalog file, to turn catalog (RA, Dec)
#  into (xi, eta) around the center of the field.
$proj_file = $catfile;
$proj_file =~ s/.cat/.proj/;
$cmd  = "~/old/tass/match/project_coords ";
$cmd .= " $catfile $cat_racol $cat_deccol $cat_ra $cat_dec ";
$cmd .= " outfile=$proj_file ";
if ($debug > 0) {
  printf "cmd is ..$cmd.. \n";
}
$ret = `$cmd`;
if ($debug > 0) {
  printf "ret is ..$ret.. \n";
}


$mod_pht = $phtfile;
if (0 == 1) {
  # treat the pht file coordinates, which are in (row, col) pixel units,
  #   so that they become closer to the expected (xi, eta) coords:
  #   radians away from the center of the image.
  $ncol = 510;
  $nrow = 340;
  $center_col = $ncol/2;
  $center_row = $nrow/2;
  $field_in_radians = (3.14159/180.0)*(1.0/5.0);
  $mod_pht =~ s/.pht/.mod/;
  open(MODPHT, ">$mod_pht") || die("can't open file $mod_pht");
  open(PHTFILE, "$phtfile") || die("can't open file $phtfile");
  while (<PHTFILE>) {
    $line = " " . $_;
    chomp($line);
    @words = split(/\s+/, $line);
    $row = $words[2];
    $col = $words[3];
    $y = (($row - $center_row)/$nrow)*$field_in_radians;
    $x = (($col - $center_col)/$ncol)*$field_in_radians;
    printf MODPHT " %6d %12.5e %12.5e %6d %7.2f %7.3f %7.3f %6d \n",
             $words[1], $y, $x, $words[4], $words[5],
             $words[6], $words[7], $words[8];
  }
  close(PHTFILE);
  close(MODPHT);
}


# run the "match" program on the projected catalog file and
#  the detected star file
$outfile = $phtfile;
$outfile =~ s/.pht/.mat/;
$cmd  = "~/old/tass/match/match ";
#$cmd .= " $phtfile $pht_xcol $pht_ycol $pht_magcol ";
$cmd .= " $mod_pht $pht_xcol $pht_ycol $pht_magcol ";
$cmd .= " $proj_file $cat_racol $cat_deccol $cat_magcol ";
$cmd .= " matchrad=5e-6 ";
$cmd .= " trirad=0.005 ";
#$cmd .= " id1=$pht_idcol id2=$cat_idcol ";
$cmd .= " outfile=$outfile ";
$cmd .= " nobj=60 ";
$cmd .= " recalc  ";
if ($debug > 0) {
  printf "cmd is ..$cmd.. \n";
}
$ret = `$cmd`;
if ($debug > 0) {
  printf "ret is ..$ret.. \n";
}


# Now apply the transformation to all the objects in the 
#   detected catalog, and restore to RA, Dec.
#
#     first, we get the coeffs of TRANS: a, b, c, d, e, f
$ret_line = $ret;
$ret_line =~ s/=/ /g;
@words = split(/\s+/, $ret_line);
$a = $words[2];
$b = $words[4];
$c = $words[6];
$d = $words[8];
$e = $words[10];
$f = $words[12];
$nm = $words[18];
$sx = $words[20];
$sy = $words[22];
if ($debug > 0) {
  printf "a %12.5e   b %12.5e  c %12.5e \n", $a, $b, $c;
}
# compute the rough scale, using "scale = sqrt(b*b + c*c)"
$scale = sqrt($b*$b + $c*$c);
$scale_recip = 1.0 / $scale;
$arcsec_per_radian = (180.0/3.14159)*3600.0;
$arcsec_per_pixel = $scale * $arcsec_per_radian;
$pixel_per_arcsec = 1.0 / $arcsec_per_pixel;
$sx_arcsec = $sx * $arcsec_per_radian;
$sy_arcsec = $sy * $arcsec_per_radian;
$qual_str =  sprintf " %6.3f arcsec/pix  Nm %2d  sx %6.3f  sy %6.3f  ", 
        $arcsec_per_pixel, $nm, $sx_arcsec, $sy_arcsec;
if ($debug > 0) {
  printf " %s \n", $qual_str;
}

#    now, run the "apply_match" program
$outfile = $phtfile;
$outfile =~ s/.pht/.radec/;
$cmd  = "~/old/tass/match/apply_match ";
$cmd .= " $phtfile $pht_xcol $pht_ycol $cat_ra $cat_dec linear ";
$cmd .= " $a $b $c $d $e $f ";
$cmd .= " outfile=$outfile ";
if ($debug > 0) {
  printf "cmd is ..$cmd.. \n";
}
$ret = `$cmd`;
if ($debug > 0) {
  printf "ret is ..$ret.. \n";
}



# now look for the target, which will be close to this position
$target_ra = 24.2008;
$target_dec = 15.7586;
$cosdec = cos($target_dec*$DEGTORAD);
$match_rad = (1.0 / 3600.0) * 3.0;
$found = 0;
open(OUTFILE, "$outfile") || die("can't open file $outfile");
while (<OUTFILE>) {
  $line = $_;
  @words = split(/\s+/, " " . $line);
  $this_ra = $words[2];
  $this_dec = $words[3];
  $dra = ($this_ra - $target_ra)*$cosdec;
  $ddec = $this_dec - $target_dec;
  $dist = sqrt($dra*$dra + $ddec*$ddec);
  if ($dist < $match_rad) {
    $dist_arcsec = $dist*3600.0;
    $found = 1;
    if ($debug > 0) {
      printf " match dist %6.2f arcsec \n", $dist_arcsec;
      printf "%s", $line;
    }
    $pos_str = sprintf " %9.5f %9.5f  %6.3f ", 
              $this_ra, $this_dec, $dist_arcsec;
    printf "Q %s %s\n", $qual_str, $pos_str;
  }
}
close(OUTFILE);






exit 0;