#!/usr/bin/perl
#
# Stack a bunch of FITS images after shifting to register.
# Use this script in the following manner:
#
#    1. run multipht to generate a list of images 
#               and offsets.  Save output into 
#               a text file
#
#    2. run this script, which will parse the
#               multipht file, locate the FITS
#               images, make copies, shift the copies
#               to register, and co-add them
#               to create two output files
#                 e1) sum all images, normalizing by number of images 
#                 e2) median of all images
# 
# MWR 4/14/2006

$debug = 1;

# name of multipht output file
$multipht_output = "./do_multi_b.out";

# the images are located in a directory 
#    which has this prefix
$image_dir = "../";

# ordinary temp files have names starting with this prefix
$output_base = "QQ";
# group temp files have names starting with this prefix
$group_base = "RR";

# maximum number of files in a group
$max_group_number = 50;

# remove temporary files when we are done?
$cleanup = 1;


#####################
# First, we read the multipht output file
#
$image_index = 0;
$group_index = 0;
$group_start[$group_index] = 0;

open(MULTIPHT, "$multipht_output") || die("can't open file $multipht_output");
while (<MULTIPHT>) {
  $line = " " . $_;
  chomp($line);
  @words = split(/\s+/, $line);

  # this is the template image, which needs no shift
  if ($image_index == 0) {
    $template_fits = $words[2];
    $template_fits =~ s/.pht/.fit/;
    $template_fits = sprintf "%s%s", $image_dir, $template_fits;
    if ($debug > 0) { 
      printf "template image is ..%s.. \n", $template_fits;
    }
    if (!(-e $template_fits)) {
      printf "cannot find file ..$template_fits..  will skip it \n";
      $template_fits = "";
    }
    else {
      if ($debug > 0) {
        printf "found image file ..$template_fits.. \n";
      }
    }
  }


  $pht_file = $words[3];
  $fits_file = $pht_file;
  $fits_file =~ s/.pht/.fit/;
  $fits_file = sprintf "%s%s", $image_dir, $fits_file;
  if ($debug > 0) { 
    printf "next image is ..%s.. \n", $fits_file;
  }

  if (!(-e $fits_file)) {
    printf "cannot find file ..$fits_file..  will skip it \n";
    next;
  }
  else {
    if ($debug > 0) {
      printf "found image file ..$fits_file.. \n";
    }
  }

  # figure out the amount by which to shift the image
  $dr = $words[6];
  $dc = $words[9];
  if ($debug > 0) {
    printf "   dr %7.3f  dc %7.3f \n", $dr, $dc;
  }


  # in the special case of the first time through the
  #   loop, we deal with the template image here.
  if ($image_index == 0) {
    if ($template_fits ne "") {
      $temp = sprintf "%s_%03d.fit", $output_base, $image_index;
      $temp_array[$image_index] = $temp;
      $cmd = "cp $template_fits $temp";
      $ret = exec_cmd($cmd);
      $cmd = "chmod 600 $temp";
      $ret = exec_cmd($cmd);
      $image_index++;
    }
  }


  # here, we make a shifted copy of the image
  $temp = sprintf "%s_%03d.fit", $output_base, $image_index;
  $temp_array[$image_index] = $temp;
  if ($debug > 0) {
    printf "temp_array[$image_index] is ..$temp_array[$image_index].. \n";
  }
  # we make a copy of the image, then shift it to match template
  $cmd = "cp $fits_file $temp";
  $ret = exec_cmd($cmd);
  $cmd = "chmod 600 $temp";
  $ret = exec_cmd($cmd);
  $dr = sprintf "%-.2f", $dr;
  $dc = sprintf "%-.2f", $dc;
  $cmd = "imshift $temp dr=$dr dc=$dc";
  $ret = exec_cmd($cmd);
  
  $image_index++;


  # do we need to start a new group at this point?
  if (($image_index - $group_start[$group_index]) >= $max_group_number) {
    $group_end[$group_index] = $image_index - 1;
    if ($debug > 0) {
      printf "  end of group %2d  which runs %5d to %5d \n",
          $group_index, $group_start[$group_index], $group_end[$group_index];
    }
    $group_index++;
    $group_start[$group_index] = $image_index;
    $group_end[$group_index] = $image_index;
  }
  

  # if ($image_index > 3) {
  #  last;
  #  }

}

# end the last group
$group_end[$group_index] = $image_index - 1;
if ($debug > 0) {
  printf "  end of group %2d  which runs %5d to %5d \n",
      $group_index, $group_start[$group_index], $group_end[$group_index];
}






###############################################################3
#  Now we do the stacking, in groups.
#     First, stack the files within each group to form
#          group output stacks.
#     Next, stack the group output stacks to form 
#          the final stacks.
#
#  We first add together the images to form sums.
#  Afterwards, before we combine to make medians,
#     we determine sky value in each image and subtract it
#     from the shifted image before combining.

for ($group_i = 0; $group_i <= $group_index; $group_i++) {

  $group_sum[$group_i] = sprintf "%s_sum_%03d.fit", $group_base, $group_i;
  $group_median[$group_i] = sprintf "%s_med_%03d.fit", $group_base, $group_i;
  if ($debug > 0) {
    printf "about to stack group %d into %s %s \n", $group_i, 
                $group_sum[$group_i], $group_median[$group_i];
  }

  # create the stack via sums
  $cmd = "add ";
  for ($i = $group_start[$group_i]; $i <= $group_end[$group_i]; $i++) {
    $cmd .= " $temp_array[$i] " ;
  }
  $cmd .= " norm outfile=$group_sum[$group_i] ";
  if ($debug > 0) {
    printf "cmd is ..$cmd.. \n";
  }
  $ret = `$cmd`;
  if ($debug > 0) {
    printf "ret is ..$ret.. \n";
  }
  
  # create the stack via median
  #
  # first, subtract sky from each temp shifted file
  for ($i = $group_start[$group_i]; $i <= $group_end[$group_i]; $i++) {
    $cmd = "sky $temp_array[$i] bin=6 ";
    $ret = exec_cmd($cmd);
    $ret =~ s/=/ /g;
    @words = split(/\s+/, $ret);
    $skyval = $words[1];
    $skysig = $words[3];
    if ($debug > 0) {
      printf "sky is %10.2f  skysig %8.2f \n", $skyval, $skysig;
    }
    $cmd = "sub $temp_array[$i] const=$skyval ";
    $ret = exec_cmd($cmd);
  }
  # now we can combine the sky-subtracted images
  $cmd = "median ";
  for ($i = $group_start[$group_i]; $i <= $group_end[$group_i]; $i++) {
    $cmd .= " $temp_array[$i] " ;
  }
  $cmd .= " nomean verbose outfile=$group_median[$group_i] ";
  if ($debug > 0) {
    printf "cmd is ..$cmd.. \n";
  }
  $ret = `$cmd`;
  if ($debug > 0) {
    printf "ret is ..$ret.. \n";
  }

}

##############################################
# now, deal with the group results.
#   if there was only one group, 
#         simply rename its output files as final output files
#   if there was more than one group,
#         run another round of "add" and "median"
  
if ($group_index == 0) {
  $cmd = "mv $group_sum[0] shift_sum.fit ";
  exec_cmd($cmd);
  $cmd = "mv $group_median[0] shift_median.fit "; 
  exec_cmd($cmd);
}
else {

  # add together all the group sums, again normalizing
  $cmd = "add ";
  for ($i = 0; $i <= $group_index; $i++) {
    $cmd .= " $group_sum[$i] " ;
  }
  $cmd .= " norm outfile=shift_sum.fit ";
  $ret = exec_cmd($cmd);
  if ($debug > 0) {
    printf "ret is ..$ret.. \n";
  }
  
  # combine together all the group medians
  $cmd = "median ";
  for ($i = 0; $i <= $group_index; $i++) {
    $cmd .= " $group_median[$i] " ;
  }
  $cmd .= " nomean verbose outfile=shift_median.fit ";
  $ret = exec_cmd($cmd);
  if ($debug > 0) {
    printf "ret is ..$ret.. \n";
  }


}
  
  

#####################
# if desired, remove all temporary files
if ($cleanup == 1) {
  @remove_files = glob("$output_base*");
  foreach $remove_file (@remove_files) {
    unlink($remove_file);
  }
  @remove_files = glob("$group_base*");
  foreach $remove_file (@remove_files) {
    unlink($remove_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);
}