#!/usr/bin/perl
#
# Carry out the matching operations on two images of NGC 4526,
#   then subtract the template from the SN image.
#
# This isn't a really general-purpose script, since it makes the 
#   assumption in the convolution stage that the template image
#   has better seeing.  See the spot marked "ASSUMPTION" below.
#
# MWR 6/6/2014

use POSIX;

$debug = 1;

$xvista_dir = "/home/richmond/xvista";

# the image WITH the supernova
$sn_img = "./May42ggb.fts";
# the template image
$templ_img = "./94d_tempr.fts";

# subtract sky from both images
$cmd = "$xvista_dir/sky $sn_img bin=3 ";
$ret = exec_cmd($cmd);
$ret =~ s/=/ /g;
@words = split(/\s+/, $ret);
$sn_sky = $words[1];
if ($debug > 0) {
  printf " SN sky is %s \n", $sn_sky;
}
$cmd = "$xvista_dir/sky $templ_img bin=3 ";
$ret = exec_cmd($cmd);
$ret =~ s/=/ /g;
@words = split(/\s+/, $ret);
$templ_sky = $words[1];
if ($debug > 0) {
  printf " template sky is %s \n", $templ_sky;
}

$cmd = "cp $sn_img sn_skysub.fts";
$ret = exec_cmd($cmd);
$cmd = "$xvista_dir/sub sn_skysub.fts const=$sn_sky";
$ret = exec_cmd($cmd);

$cmd = "cp $templ_img templ_skysub.fts";
$ret = exec_cmd($cmd);
$cmd = "$xvista_dir/sub templ_skysub.fts const=$templ_sky";
$ret = exec_cmd($cmd);


# measure the light from the bright star "A" at right side of image
#   first, in the SN image
$cmd = "$xvista_dir/box 1 cr=266 cc=421 nr=60 nc=60 ";
$ret = exec_cmd($cmd);
$cmd = "$xvista_dir/stars sn_skysub sky=0 skysig=5 minsig=10 box=1 ";
$cmd .= "   outfile=sn_skysub.coo ";
$ret = exec_cmd($cmd);
$cmd = "$xvista_dir/phot sn_skysub infile=sn_skysub.coo ";
$cmd .= "   outfile=sn_skysub.pht aper=10 skyinner=15 skyouter=20 nomag ";
$ret = exec_cmd($cmd);
#
#   now, in the template image
$cmd = "$xvista_dir/stars templ_skysub sky=0 skysig=5 minsig=10 box=1 ";
$cmd .= "   outfile=templ_skysub.coo ";
$ret = exec_cmd($cmd);
$cmd = "$xvista_dir/phot templ_skysub infile=templ_skysub.coo ";
$cmd .= "   outfile=templ_skysub.pht aper=10 skyinner=15 skyouter=20 nomag ";
$ret = exec_cmd($cmd);

# determine the ratio of brightness of star A
#
#  ratio  =  (flux of star A in SN image) / (flux of star A in template image)
#
$cmd = "cat sn_skysub.pht";
$ret = exec_cmd($cmd);
@words = split(/\s+/, $ret);
$target_flux = $words[6];
$cmd = "cat templ_skysub.pht";
$ret = exec_cmd($cmd);
@words = split(/\s+/, $ret);
$templ_flux = $words[6];
$ratio = sprintf "%.5f", $target_flux / $templ_flux;
if ($debug > 0) {
  printf " target flux %12.1f  templ flux %12.1f   ratio %8.4f \n", 
      $target_flux, $templ_flux, $ratio;
}

# tweak the ratio to get better results (if desired)
$ratio = sprintf "%.5f", $ratio*1.00;


# multiply the template image by the ratio so it matches SN image in flux
$scale_factor = 10000;
$mul_factor = sprintf "%.0f", $ratio*$scale_factor;
$cmd = "cp templ_skysub.fts templ_flux.fts";
$ret = exec_cmd($cmd);
$cmd = "$xvista_dir/mul templ_flux.fts const=$mul_factor div=$scale_factor ";
$ret = exec_cmd($cmd);


# figure out the FWHM of each image, and then determine the factor
#    by which we should convolve the template to match the SN image
$cmd = "cat sn_skysub.coo";
$ret = exec_cmd($cmd);
@words = split(/\s+/, $ret);
$sn_fwhm = $words[5];
$cmd = "cat templ_skysub.coo";
$ret = exec_cmd($cmd);
@words = split(/\s+/, $ret);
$templ_fwhm = $words[5];
# ASSUMPTION
if ($templ_fwhm >= $sn_fwhm) {
  printf " we can't handle template FWHM %7.2f > SN FWHM %7.2f \n",
        $templ_fwhm, $sn_fwhm;
  printf " halting now \n";
  exit(1);
}
$conv_fwhm = sprintf "%.2f", sqrt($sn_fwhm*$sn_fwhm - $templ_fwhm*$templ_fwhm);
if ($debug > 0) {
  printf " SN FWHM %7.2f    templ FWHM %7.2f  conv FWHM %7.2f \n",
      $sn_fwhm, $templ_fwhm, $conv_fwhm;
}

# now, convolve the template image by the appropriate amount
$cmd = "cp templ_flux.fts templ_fwhm.fts ";
$ret = exec_cmd($cmd);
$cmd = "$xvista_dir/smooth templ_fwhm.fts fw=$conv_fwhm ";
$ret = exec_cmd($cmd);



# use the position of star A in each image to determine the offset
#    by which template image should be shifted
#
#     offset =  (SN image pos) - (templ image pos)
#
$cmd = "cat sn_skysub.pht";
$ret = exec_cmd($cmd);
@words = split(/\s+/, $ret);
$sn_row = $words[2];
$sn_col = $words[3];
$cmd = "cat templ_skysub.pht";
$ret = exec_cmd($cmd);
@words = split(/\s+/, $ret);
$templ_row = $words[2];
$templ_col = $words[3];
$d_row = sprintf "%.2f", $sn_row - $templ_row;
$d_col = sprintf "%.2f", $sn_col - $templ_col;
if ($debug > 0) {
  printf " SN pos %7.2f %7.2f   templ pos %7.2f %7.2f   offset %7.2f %7.2f \n", 
      $sn_row, $sn_col, $templ_row, $templ_col, $d_row, $d_col;
}

# now, shift the template image to match the SN image
$cmd = "cp templ_fwhm.fts templ_pos.fts ";
$ret = exec_cmd($cmd);
$cmd = "$xvista_dir/imshift templ_pos.fts dr=$d_row dc=$d_col";
$ret = exec_cmd($cmd);


# finally, subtract the template from the SN image
$cmd = "cp sn_skysub.fts residual.fts ";
$ret = exec_cmd($cmd);
$cmd = "$xvista_dir/sub residual.fts templ_pos.fts";
$ret = exec_cmd($cmd);


# now, measure the light of the SN in the residual image 
#   in exactly the same manner that we measured the light of star "A"
#   in the original target image above
$cmd = "$xvista_dir/box 1 cr=263 cc=276 nr=60 nc=60 ";
$ret = exec_cmd($cmd);
$cmd = "$xvista_dir/stars residual.fts sky=0 skysig=5 minsig=10 box=1 ";
$cmd .= "   outfile=residual.coo ";
$ret = exec_cmd($cmd);
$cmd = "$xvista_dir/phot residual.fts infile=residual.coo ";
$cmd .= "   outfile=residual.pht aper=10 skyinner=15 skyouter=20 nomag ";
$ret = exec_cmd($cmd);

$cmd = "cat residual.pht";
$ret = exec_cmd($cmd);
@words = split(/\s+/, $ret);
$sn_flux = $words[6];

# use the ratio of flux (SN / A) to compute a difference in magnitudes
$ratio = $sn_flux / $target_flux;
$delta_mag = -2.5*log10($ratio);
if ($debug > 0) {
  printf "  SN flux %9.1f   A flux %9.1f     delta_mag %6.3f \n",
       $sn_flux, $target_flux, $delta_mag;
}






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


