#!/usr/bin/perl
#
# Compute color terms which turn ensemble photometry 
#    from RIT Obs into standard Johnson-Cousins magnitudes
#    determined by Landolt.
#
# The ensemble photometry is in a file with a format like this:
#
#       #star      B         V         R         I
#       PG1633   2.084     3.915    4.530     4.519
#
# and the Landolt photometry is in a file with format like this:
#
#       #star      U        B      V       R       I
#       PG1633   13.231  14.205  14.397  14.490  14.609
#
# For output, we want to derive the terms in an equation 
# like this:
#
#       (V - v)  =   k*(b - v)  +  Z
#
# where
#              V      is Landolt mag
#              v,b    are ensemble mags
#              k      is color term
#              Z      is zero-point offset
#
# It would be nice to create a graph showing V vs. (b-v) also.
#
# MWR 8/28/2011
#

use POSIX;

$debug = 0;
$temp_file = "calc_color_terms.tmp";
$cmd_file = "calc_color_terms.cmd";
$gnu_out_file = "calc_color_terms.gnu";
$gnu_out2_file = "calc_color_terms.gnu2";

# Step 0: define the passbands
$passband_array[0] = "U";
$passband_array[1] = "B";
$passband_array[2] = "V";
$passband_array[3] = "R";
$passband_array[4] = "I";

# Step 1: read in the Landolt data
$num_landolt = 0;
$landolt_file = "pg1633_landolt.dat";
open(LANDOLT_FILE, "$landolt_file") || die("can't open file $landolt_file");
while (<LANDOLT_FILE>) {
  $line = $_;
  if ($line =~ /^#/) {
    next;
  }
  @words = split(/\s+/, $line);
  $landolt_name_array[$num_landolt] = $words[0];
  $landolt_mag_array[$num_landolt][0] = $words[1];
  $landolt_mag_array[$num_landolt][1] = $words[2];
  $landolt_mag_array[$num_landolt][2] = $words[3];
  $landolt_mag_array[$num_landolt][3] = $words[4];
  $landolt_mag_array[$num_landolt][4] = $words[5];

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

  $num_landolt++;

}
close(LANDOLT_FILE);


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

  if ($debug > 0) {
    printf " ensemble %10s has mag ", $ensemble_name_array[$num_ensemble];
    for ($i = 0; $i < 5; $i++) {
      printf " %6.3f ", $ensemble_mag_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;
  }

  # We are going to create a temp file with format
  #
  #     ensemble_color    landolt_mag
  #
  # which we will later send to Gnuplot to fit
  open(TEMP_FILE, ">$temp_file") || die("can't open file $temp_file");
  for ($i_star = 0; $i_star < $num_ensemble; $i_star++) {
    
    # make sure that all the necessary magnitudes are valid
    $name = $landolt_name_array[$i_star];
    $cat_mag = $landolt_mag_array[$i_star][$mag_index];
    if ($cat_mag > 90) {
      if ($debug > 0) {
        printf "%10s: cat mag %6.2f is invalid, so skip star \n",
               $name, $cat_mag;
      }
      next;
    }

    $ens_mag = $ensemble_mag_array[$i_star][$mag_index];
    if ($ens_mag > 90) {
      if ($debug > 0) {
        printf "%10s: ens mag %6.2f is invalid, so skip star \n",
               $name, $ens_mag;
      }
      next;
    }
    $ens_color_1 = $ensemble_mag_array[$i_star][$first_index];
    if ($ens_color_1 > 90) {
      if ($debug > 0) {
        printf "%10s: ens color_1 %6.2f is invalid, so skip star \n",
               $name, $ens_color_1;
      }
      next;
    }
    $ens_color_2 = $ensemble_mag_array[$i_star][$second_index];
    if ($ens_color_2 > 90) {
      if ($debug > 0) {
        printf "%10s: ens color_2 %6.2f is invalid, so skip star \n",
               $name, $ens_color_2;
      }
      next;
    }

    # if we got here, this star has valid data, so we can compute a good color
    $color = $ens_color_1 - $ens_color_2;

    # now print a line to the temp file for this star
    printf TEMP_FILE " %6.3f   %6.3f \n",
          $color, $cat_mag - $ens_mag;


  }
  close(TEMP_FILE);


  # next, we prepare a command file for GNUPLOT
  $plot_file = sprintf "calc_%s_%s_%s.ps", $mag_name, $first_name, $second_name;

  open(CMD_FILE, ">$cmd_file") || die("can't open file $cmd_file");
  printf CMD_FILE "f(x) = a + b*x \n";
  printf CMD_FILE "fit f(x) '%s' using 1:2 via a,b \n", $temp_file;
  # and make a graph, too
  printf CMD_FILE "set term postscript color \n";
  printf CMD_FILE "set output '$plot_file' \n";
  printf CMD_FILE "set xlabel '(%s - %s) ensemble color ' \n",
             $first_name, $second_name;
  printf CMD_FILE "set ylabel '(catalog - ensemble) %s mag ' \n",
             $mag_name;
  printf CMD_FILE "set nokey \n";
  printf CMD_FILE "plot '%s' using 1:2 with points pt 5 ps 2.0 t '' ",
             $temp_file;
  printf CMD_FILE " , ";
  printf CMD_FILE "   f(x) with lines lt -1 ";
  printf CMD_FILE " \n";

  printf CMD_FILE "set output \n";
  printf CMD_FILE "quit \n";
  close(CMD_FILE) ;

  
  $cmd = "(gnuplot < $cmd_file > $gnu_out_file) 2> $gnu_out2_file ";
  $ret = exec_cmd($cmd);

  # get the results of the fit: we want to get the values for 
  #   a (and uncert) and b (and uncert)
  $a = -99.0;
  $a_sig = -99.0;
  $b = -99.0;
  $b_sig = -99.0;
  open(GNU2_FILE, "$gnu_out2_file") || die("can't open file $gnu_out2_file");
  while (<GNU2_FILE>) {
    $line = $_;
    if ($line =~ /^a.*=.*\+\/\-*/) {
      printf " line is ..%s.. \n", $line;
      @words = split(/\s+/, $line);
      $a = $words[2];
      $a_sig = $words[4];
    }
    if ($line =~ /^b.*=.*\+\/\-*/) {
      printf " line is ..%s.. \n", $line;
      @words = split(/\s+/, $line);
      $b = $words[2];
      $b_sig = $words[4];
    }
  }
  close(GNU2_FILE);
  if ($debug > 0) {
    printf "  a = %7.4f +/- %7.4f  \n", $a, $a_sig;
    printf "  b = %7.4f +/- %7.4f  \n", $b, $b_sig;
  }
 
  # and convert the Postscript output file to PNG
    $psfile = $plot_file;
    $pngfile = $psfile;
    $pngfile =~ s/.ps/.png/;
    $cmd = "convert -rotate 90 -geometry 720x504 $psfile $pngfile";
    printf "cmd is ..$cmd.. \n";
    $ret = `$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);
}


