#!/usr/bin/perl
#
# Go through all the steps required to calibrate the positions of stars
#    detected in RIT Obs images.  The result will be a version of the
#    .pht file in which the (row, col) position of each star is replaced
#    by the (RA, Dec) in decimal degrees.
#
#
# MWR 10/27/2009 

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


# directory with all the "match" software
$match_dir = "~/old/tass/match";
# a temporary file 
$cat_temp_file = "./do_match.cat.temp";
$det_temp_file = "./do_match.det.temp";


# the reference astrometric catalog file, which has (RA, Dec) and some mag
$cat_file = "./mrdel_usnob1.cat";
# zero-indexed columns in the catalog file with RA and Dec and mag
$cat_ra_col = 2;
$cat_dec_col = 3;
$cat_mag_col = 14;
# the central (RA, Dec) around which we project
$center_ra = 307.913;
$center_dec = 5.221;
# name of the projected catalog file, with (xi, eta) coords
$cat_proj_file = sprintf "%s.proj", $cat_file;

# the detected stars in an uncalibrated image
$det_file = "./gidel_2-002V.pht";
# RA and Dec and mag ranges for the detected image
#    use these to discard catalog stars outside the overlap region
$det_min_ra = 307.70; 
$det_max_ra = 308.05;
$det_min_dec = 5.15;
$det_max_dec = 5.38;
$det_max_mag = 16;
#    use these to discard detected stars outside the overlap region
$det_row_min = 0;
$det_row_max = 384;
$det_col_min = 0;
$det_col_max = 510;
# zero-indexed columns in the catalog file with (row, col) and mag
$det_row_col = 1;
$det_col_col = 2;
$det_mag_col = 5;
# some (row, col) position near the center of the image
$det_center_row = 170;
$det_center_col = 250;
# the rough plate scale of the image, in arcsec per pixel
$det_plate_scale = 1.4;


#@det_files = ("./mrdel_2-002R.pht");
@det_files = glob("./mrdel*.pht");
foreach $det_file (@det_files) {
  
  printf " next file is $det_file \n";

  
  ###################
  # Step 1: use "project_coords" to calculate the tangent plane (xi, eta)
  #            coords for stars in the reference catalog
  # first, we remove stars from the catalog which are outside the
  #            detected image's RA and Dec range, or mag range
  $num_cat = 0;
  open(INFILE, "$cat_file") || die("can't open input file $cat_file");
  open(TEMP, ">$cat_temp_file") || die("can't open temp file $cat_temp_file");
  while (<INFILE>) {
    $line = $_;
    chomp($line);
    if ($line =~ /^#/) {
      next;
    }
    $line = " " . $line;
    @words = split(/\s+/, $line);
    $include = 1;
    for ($i = 1; $i <= $#words; $i++) {
      if ((($i - 1) == $cat_ra_col)) {
        $val = $words[$i];
        if (($val < $det_min_ra) || ($val > $det_max_ra)) {
          $include = 0;
        }
      }
      if ((($i - 1) == $cat_dec_col)) {
        $val = $words[$i];
        if (($val < $det_min_dec) || ($val > $det_max_dec)) {
          $include = 0;
        }
      }
      if ((($i - 1) == $cat_mag_col)) {
        $val = $words[$i];
        if ($val > $det_max_mag) {
          $include = 0;
        }
      }
    }
    if ($include == 1) {
      printf TEMP "%s\n", $line;
      $num_cat++;
    }
  }
  close(INFILE);
  close(TEMP);
  if ($debug > 0) {
    printf " we include %d stars in trimmed catalog file \n", $num_cat;
  }
  # now we can run "project_coords" on the trimmed set of stars
  $cmd = "$match_dir/project_coords ";
  $cmd .= " $cat_temp_file $cat_ra_col $cat_dec_col $center_ra $center_dec ";
  $cmd .= " outfile=$cat_proj_file ";
  $ret = exec_cmd($cmd);
  
  
  ###################
  # Step 2: prepare the RIT .pht file for matching, by converting
  #            the (row, col) coords of stars into (x, y), which are
  #                 a) centered on center of image, so (0, 0) near center
  #                 b) scaled so that units are close to radians
  #
  $radians_per_arcsec = (3.14159/180)/3600.0;
  $radians_per_pixel = ($radians_per_arcsec)*($det_plate_scale);
  open(INFILE, "$det_file") || die("can't open input file $det_file");
  open(TEMP, ">$det_temp_file") || die("can't open temp file $det_temp_file");
  $num_det = 0;
  while (<INFILE>) {
    $line = $_;
    chomp($line);
    if ($line =~ /^#/) {
      next;
    }
    $line = " " . $line;
    @words = split(/\s+/, $line);
   
    $include = 1;
    $output_line = "";
    for ($i = 1; $i <= $#words; $i++) {
      if ((($i - 1) == $det_row_col)) {
        if (($words[$i] < $det_row_min) || ($words[$i] > $det_row_max)) {
          $include = 0;
        }
        $val = ($words[$i] - $det_center_row)*$radians_per_pixel;
      }
      elsif ((($i - 1) == $det_col_col)) {
        if (($words[$i] < $det_col_min) || ($words[$i] > $det_col_max)) {
          $include = 0;
        }
        $val = ($words[$i] - $det_center_col)*$radians_per_pixel;
      }
      else {
        $val = $words[$i];
      }
      $output_line = sprintf "%s %s ", $output_line, $val;
    }
    if ($include == 1) {
      printf TEMP "%s\n", $output_line;
      $num_det++;
    }
  }
  close(INFILE);
  close(TEMP);
  if ($debug > 0) {
    printf " num_det is %d \n", $num_det;
  }
  
  
  #########################
  # Step 3: now run "match" on the scaled detected file and the 
  #                projected catalog file
  $cmd = "$match_dir/match ";
  $cmd .= " $det_temp_file $det_row_col $det_col_col $det_mag_col ";
  $cmd .= " $cat_proj_file $cat_ra_col $cat_dec_col $cat_mag_col ";
  $cmd .= " linear recalc ";
  $cmd .= " nobj=60 trirad=0.005 ";
  $cmd .= " matchrad=1.5e-5 ";
  $cmd .= " min_scale=0.5 max_scale=1  ";
  $ret = exec_cmd($cmd);
  @words = split(/\s+/, $ret);
  if ($words[0] ne "TRANS:") {
    printf "error: match fails for file $det_file -- ret follows \n";
    printf "%s", $ret;
    next;
  }
  $ret =~ s/=/ /g;
  @words = split(/\s+/, $ret);
  $a = $words[2];
  $b = $words[4];
  $c = $words[6];
  $d = $words[8];
  $e = $words[10];
  $f = $words[12];
  if ($debug > 0) {
    printf " a %9.4e b %9.4e c %9.4e d %9.4e e %9.4e f %9.4e \n", 
          $a, $b, $c, $d, $e, $f;
    $mag = sqrt($b*$b + $c*$c);
    $rotangle = atan2($c, $b)/$DEGTORAD;
    printf "  mag %6.2f  rot %6.2f \n", $mag, $rotangle;
  }

  #########################
  # Step 4: apply the TRANS to the detected stars, so that we end up
  #                with (RA, Dec) for each star
  
  # we have to recreate the modified detected file,
  #    this time we keep all the stars (instead of discarding those
  #    outside the overlap region)
  open(INFILE, "$det_file") || die("can't open input file $det_file");
  open(TEMP, ">$det_temp_file") || die("can't open temp file $det_temp_file");
  while (<INFILE>) {
    $line = $_;
    chomp($line);
    if ($line =~ /^#/) {
      next;
    }
    $line = " " . $line;
    @words = split(/\s+/, $line);
   
    $include = 1;
    $output_line = "";
    for ($i = 1; $i <= $#words; $i++) {
      if ((($i - 1) == $det_row_col)) {
        $val = ($words[$i] - $det_center_row)*$radians_per_pixel;
      }
      elsif ((($i - 1) == $det_col_col)) {
        $val = ($words[$i] - $det_center_col)*$radians_per_pixel;
      }
      else {
        $val = $words[$i];
      }
      $output_line = sprintf "%s %s ", $output_line, $val;
    }
    if ($include == 1) {
      printf TEMP "%s\n", $output_line;
    }
  }
  close(INFILE);
  close(TEMP);
  #
  # okay, now we can use the "apply_match" program to turn these modified
  #    coordinates into (RA, Dec)
  $radec_file = $det_file . ".radec";
  $cmd = "$match_dir/apply_match $det_temp_file ";
  $cmd .= " $det_row_col $det_col_col ";
  $cmd .= " $center_ra $center_dec linear ";
  $cmd .= " $a $b $c $d $e $f ";
  $cmd .= " outfile=$radec_file ";
  $ret = exec_cmd($cmd);

  
  
  
  
}

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