#!/usr/bin/perl
#
# Given a "good" data file, with photometric data like so:
#
#       55     272.24 336.11    8 12115.57264   0.000 -0.129 0.121
#
# where 
#          col 5             is JD - 2,440,000
#          col 7             is differential mag
#
# create a file which has three different times and one differential
#   magnitude, like so:
#
#     col 1   fractional UT date                      Jul_25.07264
#
#     col 2   Julian Date - 2,450,000                 2115.57264
#
#     col 3   Heliocentric Julian Date - 2,450,000    2115.57234
#
#     col 4   differential mag                          -0.129

#     col 5   uncertainty in diff mag                    0.012
#
# -----
# Subtract the magnitude of a comparison star from the magnitude in the
#   "good" data file before printing the result.  This allows us to
#   set the zeropoint as we wish.
# MWR 7/4/2003
#
# Modified to allow me to set the zeropoint by hand at a value
#    estimated from nearby star.  I had trouble with scattered
#    light and a rotating field :-(
# MWR 6/22/2004
#
# Modified to add a fixed "uncertainty" column, the same for all
#    magnitudes.
# MWR 6/23/2005
#
# Fixed some fundamental errors in calculation of dates
#    - got rid of all shortened versions of JD and HJD -- now
#          we print the whole thing, darnit
#    - fixed the UT_day by adding -0.5 day offset, which is necessary
#          for RIT Obs
#
# MWR 6/23/2009
#
# Modified for use with Kitt Peak WIYN data from 2016



# need this for "floor"
use POSIX;
$target = "GSC_05140-03255";
$ra = "19:05:48.04";
$dec = "-06:05:49";
$utprefix = "Sep";
$utday = 19;
$filter = "v";
$filter_lc = lc($filter);

$JDBASE = 2458000.0;
$hjd_dir = "/w/spiff2012/home/richmond/bait";

# must change the following 3 items for each night
#  (get uncert_mag from datafile itself, words[7]
# $uncert_mag = 0.035;
# or, in this case, this is the amount by which to shift the differential
#   magnitudes to get a very very rough V-band mag
$aper = "4";
$arcsec_per_pix = (1.32*(2.0/2.0));
$aper_arcsec_nice = sprintf "%.1f", $aper*$arcsec_per_pix;
$add_to_mag = 0.0;

# these are 
#       a) the calibrated magnitude of a comparison star on standard scale
#            (star = 000-BNG-662 at 19:06:00.95   -06:01:14.9 
#                   has V=13.350 according to AAVSO
#       b) the mean magnitude of the comparison star in the 
#             ensemble solution
$comparison_name = "AAVSO 000-BNG-662";
$true_comparison_mag = 13.350;
$tcm = $true_comparison_mag;
$ensemble_comparison_mag = 3.392;



# print the "header" information
print "# Measurements of $target made at RIT Obs, UT 2019 $utprefix $utday, \n";
print "#    in good conditions, \n";
print "#    by Michael Richmond, \n";
print "#    using Meade 12-inch LX200 and ATIK 11000. \n";
print "# Exposures 30 seconds long, $filter filter. \n";
print "# Tabulated times are midexposure (FITS header time - half exposure length) \n";
print "#    and accurate only to +/- 1 second (??). \n";
print "# 'mag' is a differential magnitude based on ensemble photometry \n";
print "#    using a circular aperture of radius $aper pix = $aper_arcsec_nice arcseconds.  \n";
print "#    which has been shifted so $comparison_name has mag=$tcm \n"; 
print "#    which is its V-band magnitude according to AAVSO.  \n";
print "# \n";
print "# UT_day             JD            HJD        mag    uncert\n";
#print "# UT day      JD-2,450,000  HJD-2,450,000   mag    \n";


$file = sprintf "star_gsc_%s.dat%d", $filter_lc, $aper;
open(INFILE, "$file") || die ("can't open file $file");
while (<INFILE>) {
  @words = split(/\s+/);
  $jd = $words[4];
  $mag = $words[6] + ($true_comparison_mag - $ensemble_comparison_mag);
  #$mag = $words[6] + $add_to_mag;
  $uncert_mag = $words[7];

  # calculate the fractional UT day
  #   this section very different for SDSS data!
  $fraction = $jd - floor($jd);
  $fraction -= 0.5;
  $fraction += $utday;
  $utdate = sprintf("%s%08.5f", $utprefix, $fraction);

  # don't mess around with shortened values of the dates
  if (0 == 1) {
    # calculate the Julian Date - 2,450,000
    $full_julian_date = $jd + 2450000.0;
    #$full_julian_date = $jd;
    $short_jd = $full_julian_date - 2450000.0;
  }

  # calculate the heliocentric julian date
  $full_hjd = `$hjd_dir/hjd jd=$jd ra=$ra dec=$dec`;
  if (0 == 1) {
    $short_hjd = $full_hjd - 2450000.0;
  }

  # print out the results
  #   this version if I have a real uncertainty estimate
  # printf "%-15s %11.5f  %11.5f  %6.3f %6.3f \n", $utdate, $short_jd, $short_hjd, 
   #                  $mag, $uncert_mag;
  #   this version if I don't
  printf "%-15s %12.5f  %12.5f  %6.3f %6.3f \n", 
           $utdate, $jd, $full_hjd, $mag, $uncert_mag;

}

exit 0;


