#!/usr/bin/perl
#
# Compute magnitude values for SN 2011fe,
#    using raw ensemble mags, catalog mags for stars in the field,
#    and a set of color terms.
#
# MWR 8/28/2011
#
# Modified for SN 2013ej use.
#
# MWR 8/1/2013


use POSIX;

$debug = 0;
$color_terms_file = "color_terms.dat";
$temp_file = "calc_sn_mags.tmp";

# Step 0: define the passbands
$num_passband = 5;
$passband_array[0] = "U";
$passband_array[1] = "B";
$passband_array[2] = "V";
$passband_array[3] = "R";
$passband_array[4] = "I";
# and define the ensemble, instrumental mags
for ($i = 0; $i < $num_passband; $i++) {
  $raw_passband_array[$i] = lc($passband_array[$i]);
  if ($debug > 0) { 
    printf "  raw passband %3d is %s \n", $i, $raw_passband_array[$i];
  }
}

# Step 0.1: compute average JD for each passband
#     based on the times associated with each image in .img file
if (calc_avg_jd() != 0) {
  printf " calc_avg_jd fails \n";
  exit(1);
}


# Step 0.5: read in color term info
$num_colorterm = 0;
#$colorterm_file = "color_terms.dat";
$colorterm_file = "color_terms_aug05_2013.dat";
open(COLORTERM_FILE, "$colorterm_file") || die("can't open file $colorterm_file");
while (<COLORTERM_FILE>) {
  $line = $_;
  if ($line =~ /^#/) {
    next;
  }
  $line = " " . $line;
  @words = split(/\s+/, $line);
  $colorterm_cat_array[$num_colorterm] = $words[1];
  $colorterm_ens_array[$num_colorterm] = $words[2];
  $colorterm_first_array[$num_colorterm] = $words[3];
  $colorterm_second_array[$num_colorterm] = $words[4];
  $colorterm_k_array[$num_colorterm] = $words[5];
  $colorterm_ksig_array[$num_colorterm] = $words[6];

  if ($debug > 0) {
    printf " colorterm %2d  has ", $num_colorterm;
    printf "   %4s  %4s  %4s %4s   %7.3f %7.3f ", 
       $colorterm_cat_array[$num_colorterm],
       $colorterm_ens_array[$num_colorterm],
       $colorterm_first_array[$num_colorterm],
       $colorterm_second_array[$num_colorterm],
       $colorterm_k_array[$num_colorterm],
       $colorterm_ksig_array[$num_colorterm];
    printf "\n";
  }

  $num_colorterm++;

}
close(COLORTERM_FILE);




# Step 1: read in the catalog data
$num_catalog = 0;
$catalog_file = "aavso_only_hand.dat";
open(CATALOG_FILE, "$catalog_file") || die("can't open file $catalog_file");
while (<CATALOG_FILE>) {
  $line = $_;
  if ($line =~ /^#/) {
    next;
  }
  @words = split(/\s+/, " " . $line);
  $catalog_name_array[$num_catalog] = $words[2];
  $catalog_mag_array[$num_catalog][0] = 99.0;
  $catalog_mag_array[$num_catalog][1] = $words[5];
  $catalog_mag_array[$num_catalog][2] = $words[7];
  $catalog_mag_array[$num_catalog][3] = $words[9];
  $catalog_mag_array[$num_catalog][4] = $words[11];

  if ($debug > 0) {
    printf " star %10s has mag ", $catalog_name_array[$num_catalog];
    for ($i = 0; $i < 5; $i++) {
      printf " %6.3f ", $catalog_mag_array[$num_catalog][$i];
    }
    printf "\n";
  }

  $num_catalog++;

}
close(CATALOG_FILE);




# STep 2: read in the ensemble data
$num_ensemble = 0;
$ensemble_file = "m74_ensemble.dat";
open(ENSEMBLE_FILE, "$ensemble_file") || die("can't open file $ensemble_file");
while (<ENSEMBLE_FILE>) {
  $line = $_;
  if ($line =~ /^#/) {
    next;
  }
  if ($line =~ /^$/) {
    next;
  }
  @words = split(/\s+/, " " . $line);
  $ensemble_name_array[$num_ensemble] = $words[1];
  $ensemble_mag_array[$num_ensemble][1] = $words[2];
  $ensemble_magerr_array[$num_ensemble][1] = $words[3];
  $ensemble_mag_array[$num_ensemble][2] = $words[4];
  $ensemble_magerr_array[$num_ensemble][2] = $words[5];
  $ensemble_mag_array[$num_ensemble][3] = $words[6];
  $ensemble_magerr_array[$num_ensemble][3] = $words[7];
  $ensemble_mag_array[$num_ensemble][4] = $words[8];
  $ensemble_magerr_array[$num_ensemble][4] = $words[9];

  if ($debug > 0) {
    printf " ensemble %10s has mag ", $ensemble_name_array[$num_ensemble];
    for ($i = 0; $i < 5; $i++) {
      printf " %6.3f %5.3f   ", 
         $ensemble_mag_array[$num_ensemble][$i],
         $ensemble_magerr_array[$num_ensemble][$i];
    }
    printf "\n";
  }

  $num_ensemble++;

}
close(ENSEMBLE_FILE);




# Step 3: loop over the desired magnitudes and colors,
#             for each one, compute best-fit coefficients
#             (and make a graph)

for ($mag_index = 1; $mag_index <= 4; $mag_index++) {

  # choose the color
  #    B  -> (B-V)    V (V-R)    R (R-I)    I (R-I)  
  if (($mag_index >= 1) && ($mag_index <= 3)) {
    $first_index = $mag_index;
    $second_index = $mag_index + 1;
  }
  else {  
    $first_index = $mag_index - 1;
    $second_index = $mag_index;
  }
  $mag_name = $passband_array[$mag_index];
  $first_name = $passband_array[$first_index];
  $second_name = $passband_array[$second_index];
  if ($debug > 0) {
    printf "  working on filter %s   as func of (%s - %s) \n",
        $mag_name, $first_name, $second_name;
  }


  # prepare to compute the "A" term for this filter
  $sum = 0.0;
  $sumsq = 0.0;
  $num = 0;

  open(TEMP_FILE, ">$temp_file") || die("can't open file $temp_file");
  for ($i_star = 0; $i_star < $num_ensemble; $i_star++) {
    
    $ens_name = $ensemble_name_array[$i_star];
    $ens_mag = $ensemble_mag_array[$i_star][$mag_index];
    $ens_jd = $ensemble_jd_array[$i_star][$mag_index];
    if ($ens_mag > 90) {
      if ($debug > 0) {
        printf "%10s: ens mag %6.2f is invalid, so skip star \n",
               $ens_name, $ens_mag;
      }
      next;
    }
    $first_mag = $ensemble_mag_array[$i_star][$first_index];
    if ($first_mag > 90) {
      if ($debug > 0) {
        printf "%10s: ens first mag %6.2f is invalid, so skip star \n",
               $ens_name, $first_mag;
      }
      next;
    }
    $second_mag = $ensemble_mag_array[$i_star][$second_index];
    if ($second_mag > 90) {
      if ($debug > 0) {
        printf "%10s: ens second mag %6.2f is invalid, so skip star \n",
               $ens_name, $second_mag;
      }
      next;
    }



    # find the corresponding catalog star (if it exists)
    $found = 0;
    for ($j_star = 0; $j_star < $num_catalog; $j_star++) {
      $cat_name = $catalog_name_array[$j_star];
      if ($cat_name eq $ens_name) {
        $found = 1;
        $cat_mag = $catalog_mag_array[$j_star][$mag_index];
        last;
      }
    }
    if ($found == 0) {
      next;
    }
    if ($cat_mag > 90) {
      if ($debug > 0) {
        printf "%10s: cat mag %6.2f is invalid, so skip star \n",
               $cat_name, $cat_mag;
      }
      next;
    }

    if ($debug > 0) {
      printf " star %10s  cat_mag %7.3f  ens_mag %7.3f  %7.3f %7.3f \n",
             $cat_name, $cat_mag, $ens_mag, $first_mag, $second_mag;
    }


    # figure out which k-term goes with this magnitude
    $found = 0;
    for ($i_term = 0; $i_term < $num_colorterm; $i_term++) {
      if (  ($passband_array[$mag_index] eq $colorterm_cat_array[$i_term]) &&
            ($raw_passband_array[$mag_index] eq $colorterm_ens_array[$i_term]) &&
            ($raw_passband_array[$first_index] eq $colorterm_first_array[$i_term]) &&
            ($raw_passband_array[$second_index] eq $colorterm_second_array[$i_term])) {
        $found = 1;
        $k_term = $colorterm_k_array[$i_term];
        last;
      }
    }
    if ($found == 0) {
      printf " error: can't find colorterm \n";
      exit(1);
    }
    if ($debug > 0) {
      printf " colorterm for %s as func of (%s - %s) is %8.4f \n",
               $passband_array[$mag_index],
               $raw_passband_array[$first_index],
               $raw_passband_array[$second_index],
               $k_term;
    }
       
    # compute the contribution of this star to the average value of
    #
    #     A  =  (cat_mag - ens_mag) - k*(first_mag - second_mag)
    #
    $diff = ($cat_mag - $ens_mag) - $k_term*($first_mag - $second_mag);
    $sum += $diff;
    $sumsq += ($diff*$diff);
    if ($debug > 0) {
      printf "  diff is %8.4f   sum %8.4f  sumsq %8.4f \n",
             $diff, $sum, $sumsq;
    }
    $num++;

  }
  
  # now we compute the average value of "A", and its stdev
  if ($num == 0) {
    printf " error, no mags to average? \n";
    $sum = 99.0;
    $num = 1;
    # exit(1);
  }
  $mean = $sum / $num;
  if ($num == 1) {
    $stdev = 0.0;
  }
  else {
    $stdev = sqrt(($sumsq - $num*$mean*$mean)/($num - 1.0));
  }
  printf " %s = %s + %7.4f * (%s - %s) + %7.4f +/- %7.4f   N = %d \n",
               $passband_array[$mag_index],
               $raw_passband_array[$mag_index],
               $k_term,
               $raw_passband_array[$first_index],
               $raw_passband_array[$second_index],
               $mean, 
               $stdev, 
               $num;
 

  # now compute the magnitude of the SN in this passband
  $found = 0;
  for ($i_star = 0; $i_star < $num_ensemble; $i_star++) {
    
    $ens_name = $ensemble_name_array[$i_star];
    if ($ens_name eq "SN") {
      $found = 1;
   
      $ens_mag = $ensemble_mag_array[$i_star][$mag_index];
      $ens_jd = $ensemble_jd_array[$i_star][$mag_index];
      if ($ens_mag > 90) {
        if ($debug > 0) {
          printf "%10s: ens mag %6.2f is invalid, so skip star \n",
                 $ens_name, $ens_mag;
        }
        next;
      }
      $first_mag = $ensemble_mag_array[$i_star][$first_index];
      if ($first_mag > 90) {
        if ($debug > 0) {
          printf "%10s: ens first mag %6.2f is invalid, so skip star \n",
                 $ens_name, $first_mag;
        }
        next;
      }
      $second_mag = $ensemble_mag_array[$i_star][$second_index];
      if ($first_mag > 90) {
        if ($debug > 0) {
          printf "%10s: ens second mag %6.2f is invalid, so skip star \n",
                 $ens_name, $second_mag;
        }
        next;
      }

      # Compute the overall uncertainty in mag
      $ens_magerr = $ensemble_magerr_array[$i_star][$mag_index];
      $total_magerr = sqrt($stdev*$stdev + $ens_magerr*$ens_magerr);


      # Okay, we can compute corrected SN mag now
      $sn_mag = $ens_mag + $k_term*($first_mag - $second_mag) + $mean;
      printf "   SN  %s = %8.3f   +/-  %6.3f  (ens %6.3f zp %6.3f)    %12.5f \n\n",
               $passband_array[$mag_index], $sn_mag, 
               $total_magerr, $ens_magerr, $stdev,
               $jd_array[$mag_index];

      last;

    }

  }
  if ($found == 0) {
    printf " error, couldn't find SN? \n";
    exit(1);
  }

  
  
  if ($mag_index == 5) {
    exit(0);
  }


}







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




#####################################################################
# PROCEDURE: calc_avg_jd() 
#
# DESCRIPTION: Figure out the average JD of the SN measurement
#              in each passband, based on the times of 
#              mid-exposure of each image used in the ensemble
#              solution.
#
#              Look in the files like
#
#                       solve_sn_4_b.img
#
#              and compute the average value of the JDs. 
#              Yes, this will be wrong if some image(s)
#              were not included in the solution ... big deal.
#
#              This procedure sets values in the global array
#
#                    jd_array[]
#
# RETURNS:     0           if all goes well
#              1           if an error occurs
#
sub calc_avg_jd {

  my($filter);
  my($sum, $num, $avg);
  my($line, @words);
  my($jd);
  my($img_base);

  $img_base = "./solve_sn_4_";
  for ($i = 0; $i < $num_passband; $i++) {
    $filter_lc = lc($passband_array[$i]);
    $file = sprintf "%s%s.img", $img_base, $filter_lc;
    if (!-f $file) {
      $jd_array[$i] = -1;
      if ($debug > 0) {
        printf " calc_avg_jd: can't find file %s, so JD %s = %12.5f \n",
            $file, $filter_lc, $jd_array[$i];
      }
      next;
    }

    $sum = 0.0;
    $num = 0;
    open(FILE, "$file") || die("can't open file $file");
    while (<FILE>) {
      $line = $_;
      if ($line =~ /^#/) {
        next;
      }
      @words = split(/\s+/, " " . $line);
      $jd = $words[4];
      $sum += $jd;
      $num++;
      if ($debug > 0) {
        printf " filter %s   add JD %12.5f  so num %3d \n",
               $filter_lc, $jd, $num;
      }
    }
    if ($num < 1) {
      printf " calc_avg_jd: zero images in file %s ?! \n", $file;
      $avg = -1;
    }
    else {
      $avg = $sum / $num;
    }
    if ($debug > 0) {
      printf " calc_avg_jd: filter %s  num %3d  avg %12.5f \n",
            $filter_lc, $num, $avg;
    }

    $jd_array[$i] = $avg;
  }

  return(0);
}