#!/usr/bin/perl
#
# Create a file with photometry of the star "D", which is apparently
#    a variable star.  The format must be
#
#      col 1         HJD
#      col 2         (D - C) mag
#      col 3         (D - C) uncert
#      col 4         (D - B) mag
#      col 5         (D - B) uncert
#
# This is the format requested by Enrique de Miguel Agustino
#    demiguel@uhu.es
# who is doing the analysis.
#
# I'll simply add the estimated uncertainty in mag for each
#   star in quadrature to get estimated uncertainty in (D-C) and (D-B).
#
# MWR 7/9/2010

$debug = 1;

$ra = "16:25:20.3";
$dec = "+12:03:09";
$hjd_dir = "~/old/bait/";


$d_file = "starD.dat";
$b_file = "starB.dat";
$c_file = "starC.dat";


# first, we read in all the data for star D
$num_d_obs = 0;
$d_hjd_array[0] = 0.0;
$d_mag_array[0] = 0.0;
$d_err_array[0] = 0.0;
open(INFILE, $d_file) || die("can't open file $d_file");
while (<INFILE>) {
  $line = $_;
  chomp($line);
  $line = " " . $line;
  @words = split(/\s+/, $line);
  $jd = $words[5];
  $mag = $words[7];
  $err = $words[8];
  
  # correct the JD to HJD
  $hjd = `$hjd_dir/hjd jd=$jd ra=$ra dec=$dec`;

  # save the values for later reference
  $d_hjd_array[$num_d_obs] = $hjd;
  $d_mag_array[$num_d_obs] = $mag;
  $d_err_array[$num_d_obs] = $err;

  $num_d_obs++;
}
close(INFILE);
if ($debug > 0) {
  printf " read %5d measures of star D \n", $num_d_obs;
}



# next, we read in all the data for star B
$num_b_obs = 0;
$b_hjd_array[0] = 0.0;
$b_mag_array[0] = 0.0;
$b_err_array[0] = 0.0;
open(INFILE, $b_file) || die("can't open file $b_file");
while (<INFILE>) {
  $line = $_;
  chomp($line);
  $line = " " . $line;
  @words = split(/\s+/, $line);
  $jd = $words[5];
  $mag = $words[7];
  $err = $words[8];
  
  # correct the JD to HJD
  $hjd = `$hjd_dir/hjd jd=$jd ra=$ra dec=$dec`;

  # save the values for later reference
  $b_hjd_array[$num_b_obs] = $hjd;
  $b_mag_array[$num_b_obs] = $mag;
  $b_err_array[$num_b_obs] = $err;

  $num_b_obs++;
}
close(INFILE);
if ($debug > 0) {
  printf " read %5d measures of star B \n", $num_b_obs;
}

# next, we read in all the data for star C
$num_c_obs = 0;
$c_hjd_array[0] = 0.0;
$c_mag_array[0] = 0.0;
$c_err_array[0] = 0.0;
open(INFILE, $c_file) || die("can't open file $c_file");
while (<INFILE>) {
  $line = $_;
  chomp($line);
  $line = " " . $line;
  @words = split(/\s+/, $line);
  $jd = $words[5];
  $mag = $words[7];
  $err = $words[8];
  
  # correct the JD to HJD
  $hjd = `$hjd_dir/hjd jd=$jd ra=$ra dec=$dec`;

  # save the values for later reference
  $c_hjd_array[$num_c_obs] = $hjd;
  $c_mag_array[$num_c_obs] = $mag;
  $c_err_array[$num_c_obs] = $err;

  $num_c_obs++;
}
close(INFILE);
if ($debug > 0) {
  printf " read %5d measures of star C \n", $num_c_obs;
}



# now, we walk through all the measurements of star "D".
#   For each one, we look for a matching-in-time measurement of C and B.
#   If either C or B is missing, skip this time
#   If both C and B have measurements, then compute the differences
#      and print a line with the format
#         HJD     (D - C)_mag  (D - C)_uncert  (D - B)_mag (D - B)_uncert
#
for ($i_d = 0; $i_d < $num_d_obs; $i_d++) {
  $this_hjd = $d_hjd_array[$i_d];

  $found_b_flag = 0;
  for ($i_b = 0; $i_b < $num_b_obs; $i_b++) {
    if ($b_hjd_array[$i_b] == $this_hjd) {
      $found_b_flag = 1;
      last;
    }
  }

  $found_c_flag = 0;
  for ($i_c = 0; $i_c < $num_c_obs; $i_c++) {
    if ($c_hjd_array[$i_c] == $this_hjd) {
      $found_c_flag = 1;
      last;
    }
  }

  if (($found_b_flag == 0) || ($found_c_flag == 0)) {
    if ($debug > 0) {
      printf " no match for HJD %12.5f \n", $this_hjd;
    }
    next;
  }

  # compute the differences and uncertainties
  $d_c_mag = $d_mag_array[$i_d] - $c_mag_array[$i_c];
  $d_c_err = sqrt(  $d_err_array[$i_d]*$d_err_array[$i_d] +
                    $c_err_array[$i_c]*$c_err_array[$i_c] );
  $d_b_mag = $d_mag_array[$i_d] - $b_mag_array[$i_b];
  $d_b_err = sqrt(  $d_err_array[$i_d]*$d_err_array[$i_d] +
                    $b_err_array[$i_b]*$b_err_array[$i_b] );

  # print out the result
  printf " %12.5f  %6.3f %6.3f   %6.3f %6.3f \n",
         $this_hjd, $d_c_mag, $d_c_err, $d_b_mag, $d_b_err;

}




exit 0;