#!/usr/bin/perl
# 
# Given a photom .cal file, this program brings all the observations
#   of a same physical star together in one place, and analyzes the
#   properties of those observations as a group.
#
# The user must SORT INPUT FILE BY RA in order for this program
#   to work properly.
#
# It prints out an ASCII text file with one line per observation,
#   grouping the lines for the same star together.
#
# usage:    sort -n +1 M*.cal | calc_gammas_cal.pl  radecfile  
#
#   where
#               M*.cal         is a calibrated photometry file created
#                                    by the pipeline
#               radecfile      is a file with imagename, JD, RA, Dec,
#                                    created by grab_pos.pl
#

$debug = 0;

$radecfile = $ARGV[0];
if (-e $radecfile == 0) {
  printf "calc_gammas_cal.pl: radecfile ..$radecfile.. doesn't exist?\n";
  exit 1;
}

$line_count = 0;
$cur_ra = 9999.0;       # degrees
$cur_dec = 9999.0;      # degrees
$cur_jd = 9999.0;
$cur_calv = 99.0;
$match_radius = 5.0;    # arcseconds

# we accumulate values for statistical purposes in these arrays
$num_obs = 0;
$vmag_array[0] = 0;
$vsig_array[0] = 0;
$imag_array[0] = 0;
$isig_array[0] = 0;
$jd_array[0] = 0;

# we set this flag to 1 if the next star we read is a new one
$is_new_star = 0;

while (<STDIN>) {
  $line_count++;
  $line = $_;
  @words = split(/\s+/, $line);
  $star_id = $words[1];
  $ra = $words[2];
  $dec = $words[3];
  $jd = $words[4];
  $vmag = $words[6];
  $vsig = $words[7];
  $vflag = $words[8];
  $imag = $words[10];
  $isig = $words[11];
  $iflag = $words[12];

  if (0 == 1) {
    printf "%9.4f %9.4f %13.5f  %6.3f %6.3f  %6.3f %6.3f \n", 
       $ra, $dec, $jd, $vmag, $vsig, $imag, $isig;
  }

  # ignore stars if flags indicate possible problems
  if (0 == 1) {
    if (($vflag != 0) || ($iflag != 0)) {
      next;
    }
  }


  # ignore stars which are outside of the specified Dec limits
  #   (which can be used to isolate a certain portion of the field of view)
  if (0 == 1) {
    $center_dec = -0.16;
    $dec_width = 0.5;
    if (($dec > $center_dec + $dec_width) || ($dec < $center_dec - $dec_width)) {
      next;
    }
  }


 if (1 == 1) {
  # check -- is this a new star?
  if (coords_match($ra, $dec, $cur_ra, $cur_dec, $match_radius) == 0) {
    $is_new_star = 1;

    # we should perform analysis on old star's value here, if there
    #   _is_ an old star to analyze 
    if ($line_count > 1) {
  

      # printf "last star had %d obs \n", $num_obs;
      $avg_v = calc_mean($num_obs, @vmag_array);
      $sig_v = calc_sigma($num_obs, @vmag_array);
      $avg_vsig = calc_mean($num_obs, @vsig_array);
      $avg_i = calc_mean($num_obs, @imag_array);
      $sig_i = calc_sigma($num_obs, @imag_array);
      $avg_isig = calc_mean($num_obs, @isig_array);

      for ($i = 0; $i < $num_obs; $i++) {

        # find the central (RA, Dec) values for this V-band and I-band image
        ($central_rav, $central_decv, $central_rai, $central_deci) = 
              find_central_radec($jd_array[$i], $radecfile);
        if ($debug > 0) {
          printf "central  %9.4f %9.4f   %9.4f %9.4f \n", 
                   $central_rav, $central_decv,
                   $central_rai, $central_deci;
        }
        ($delta_rav, $delta_decv, $delta_rai, $delta_deci) = 
             calc_delta_pos($cur_ra, $cur_dec, $central_rav, $central_decv,
  	                                     $central_rai, $central_deci);

        $gamma_v = $vmag_array[$i] - $avg_v;
        $gamma_i = $imag_array[$i] - $avg_i;
        printf "%-2d %6.3f %6.3f %6.3f %6.3f    %6.3f %6.3f %6.3f %6.3f %9.4f %9.4f %9.4f %9.4f \n", 
                              $num_obs, $avg_v, $gamma_v, $avg_vsig, $sig_v, 
                                        $avg_i, $gamma_i, $avg_isig, $sig_i,
                                        $delta_rav, $delta_decv,
  	  		                $delta_rai, $delta_deci;
      }

    }

    # and reset, fill zero'th elements with values for new star
    $num_obs = 0;
    $vmag_array[$num_obs] = $vmag;
    $vsig_array[$num_obs] = $vsig;
    $imag_array[$num_obs] = $imag;
    $isig_array[$num_obs] = $isig;
    $jd_array[$num_obs] = $jd;
    $num_obs++;
    $cur_ra = $ra;
    $cur_dec = $dec;
  }
  else {
    # nope, this is still same old star, just add elements to arrays
    $vmag_array[$num_obs] = $vmag;
    $vsig_array[$num_obs] = $vsig;
    $imag_array[$num_obs] = $imag;
    $isig_array[$num_obs] = $isig;
    $jd_array[$num_obs] = $jd;
    $num_obs++;
  }
 }


}

  # we should perform analysis on the last star's value here...

    # we should perform analysis on old star's value here...
    # printf "last star had %d obs \n", $num_obs;
    $avg_v = calc_mean($num_obs, @vmag_array);
    $sig_v = calc_sigma($num_obs, @vmag_array);
    $avg_vsig = calc_mean($num_obs, @vsig_array);
    $avg_i = calc_mean($num_obs, @imag_array);
    $sig_i = calc_sigma($num_obs, @imag_array);
    $avg_isig = calc_mean($num_obs, @isig_array);
    for ($i = 0; $i < $num_obs; $i++) {

      # find the central (RA, Dec) values for the V-band and I-band images
      ($central_rav, $central_decv, $central_rai, $central_deci) = 
            find_central_radec($jd_array[$i], $radecfile);
      if ($debug > 0) {
        printf "central  %9.4f %9.4f   %9.4f %9.4f \n", 
               $central_rav, $central_decv,
               $central_rai, $central_deci;
      }
      if ($central_rav == -99) {
        next;
      }
      
      ($delta_rav, $delta_decv, $delta_rai, $delta_deci) = 
           calc_delta_pos($cur_ra, $cur_dec, $central_rav, $central_decv,
	                                     $central_rai, $central_deci);

      $gamma_v = $vmag_array[$i] - $avg_v;
      $gamma_i = $imag_array[$i] - $avg_i;
      printf "%-2d %6.3f %6.3f %6.3f %6.3f    %6.3f %6.3f %6.3f %6.3f %9.4f %9.4f %9.4f %9.4f \n", 
                            $num_obs, $avg_v, $gamma_v, $sig_v, $avg_vsig,
                                      $avg_i, $gamma_i, $sig_i, $avg_isig,
                                      $delta_rav, $delta_decv,
				      $delta_rai, $delta_deci;
    }



exit 0;


#####################################################################
# PROCEDURE: coords_match
#
# DESCRIPTION: given two sets of coords (both in decimal degrees),
#              calculate the distance between them.  If the distance
#              is 
#
#                   <= match_radius         then we match: return 1
#                   >  match_radius         then no match: return 0
#
# Note that the "match_radius" argument is in arcseconds, 
#   not degrees.
#
# USAGE:
#         coords_match   this_ra this_dec  old_ra old_dec  match_radius
#
sub coords_match {
 
  my($this_ra, $this_dec,
     $old_ra,  $old_dec,
     $match_radius,
     $delta_ra, $delta_dec,
     $avg_dec, $factor,
     $DEGTORAD);

  $this_ra = $_[0];
  $this_dec = $_[1];
  $old_ra = $_[2];
  $old_dec = $_[3];
  $match_radius = $_[4];

  # convert matching radius to degrees
  $match_radius /= 3600.0;

  $delta_ra = abs($this_ra - $old_ra);
  $delta_dec = abs($this_dec - $old_dec);

  # printf "match_radius %9.5lf  delta_ra %9.5lf  delta_dec %9.5lf \n",
  #        $match_radius, $delta_ra, $delta_dec;

  # do a quick check first 
  if ($delta_dec > $match_radius) {
    return 0;
  }

  # we need to modify the RA distance due to cos(Dec) factor.
  #   Use the avg dec to calc this factor
  $avg_dec = ($this_dec + $old_dec) / 2.0;
  $DEGTORAD = 0.017453293;
  $factor = cos($avg_dec*$DEGTORAD);
  $delta_ra *= $factor;
  if ($delta_ra > $match_radius) {
    return 0;
  }

  # well, if we get this far, the two coords do match within
  #   a square box of the given matching radius.
  #   Call it a match
  return 1;

}


########################################################################
# PROCEDURE: calc_mean
#
# DESCRIPTION: given the number of elements in an array, and the
#              array itself, calculate the mean value of the array.
#
# RETURNS:
#              the mean value
#
sub calc_mean {

  my($num, @array);
  my($i, $sum);
  my($mean);

  $num = $_[0];
  for ($i = 0; $i < $num; $i++) {
    $array[$i] = $_[$i+1];
  }

  if ($num == 0) {
    $mean = 0.0;
    return($mean);
  }
 
  # printf "%2d ", $num;
  for ($i = 0; $i < $num; $i++) {
  #   printf "%6.3f ", $array[$i];
  }

  $sum = 0.0;
  for ($i = 0; $i < $num; $i++) {
    $sum += $array[$i];
  }
  $mean = $sum/$num;
  # printf(" %7.4f \n", $mean);

  return($mean);
}


########################################################################
# PROCEDURE: calc_sigma
#
# DESCRIPTION: given the number of elements in an array, and the
#              array itself, calculate the standard deviation from the
#              mean of an array
#
# RETURNS:
#              the stdev
#
sub calc_sigma {

  my($num, @array);
  my($i, $sum);
  my($sumsq);
  my($mean, $stdev);

  $num = $_[0];
  for ($i = 0; $i < $num; $i++) {
    $array[$i] = $_[$i+1];
  }

  if ($num < 2) {
    $stdev = 0.0;
    return($stdev);
  }
 
  # printf "%2d ", $num;
  for ($i = 0; $i < $num; $i++) {
  #   printf "%6.3f ", $array[$i];
  }

  $sum = 0.0;
  $sumsq = 0.0;
  for ($i = 0; $i < $num; $i++) {
    $sum += $array[$i];
    $sumsq += $array[$i]*$array[$i];
  }
  $mean = $sum/$num;
  $stdev = sqrt(($sumsq - $num*$mean*$mean)/($num - 1));
  # printf(" %7.4f \n", $stdev);

  return($stdev);
}


########################################################################
# PROCEDURE: calc_delta_pos
#
# DESCRIPTION: given the (RA, Dec) of a star, and the (RA, Dec) centers
#              of the V-band and I-band images in which it appears,
#              calculate the difference between its RA and the
#              central RA for its images, and between its Dec and 
#              the central Dec for its images.
#
# RETURNS:
#              We return a list with 4 values:
#
#                    this ra  - central RA   in V-band
#                    this dec - central Dec  in V-band
#                    this ra  - central RA   in I-band
#                    this dec - central Dec  in I-band
#
#
sub calc_delta_pos {

  my($ra, $dec);
  my($rav_center, $decv_center);
  my($rai_center, $deci_center);
  my($delta_rav, $delta_decv);
  my($delta_rai, $delta_deci);

  $ra = $_[0];
  $dec = $_[1];
  $rav_center = $_[2];
  $decv_center = $_[3];
  $rai_center = $_[4];
  $deci_center = $_[5];

  # calc deltas in each band
  $delta_rav = $ra - $rav_center;
  $delta_decv = $dec - $decv_center;
  $delta_rai = $ra - $rai_center;
  $delta_deci = $dec - $deci_center;

  return($delta_rav, $delta_decv, $delta_rai, $delta_deci);
}


#############################################################################
# Given a JD and a file name with JD/RA/Dec cross reference
#   find the central (RA, Dec) for the V-band and I-band
#   images taken at that time.  Return those values as a 4-element list.
#
# If we can't find them, return values of -99 for all offsets.
#
sub find_central_radec {
  my($jd);
  my($file);
  my($rav, $decv, $rai, $deci);


  $jd = $_[0];
  $file = $_[1];

  if ($debug > 0) {
    printf "find_central_radec: JD %13.5f  file %s  \n", $jd, $file;
  }

  # initialize values to nonsense
  $rav = -99;
  $decv = -99;
  $rai = -99;
  $deci = -99;

  # now, look through the given file for the given JD
  open(JDFILE, $file) || die ("find_central_radec: can't open file $file");
  while (<JDFILE>) {
    @words = split(/\s+/);
    $this_jd = $words[2];
    if (abs($jd - $this_jd) < 0.0002) {
      if ($words[1] =~ /hira/) {
        $rai = $words[3];
	$deci = $words[4];
      }
      if ($words[1] =~ /hvra/) {
        $rav = $words[3];
	$decv = $words[4];
      }
    }
  }
  close(JDFILE);

  # check to make sure that we did find the given JD
  if (($rav == -99) || ($rai == -99)) {
    printf STDERR "find_central_radec: couldn't find jd $jd\n";
    $rav = -99; $decv = -99; $rai = -99; $deci = -99;
  }

  @retlist = ($rav, $decv, $rai, $deci);

  return(@retlist);
}