#!/usr/bin/perl -w
#      collect_stars.pl    12 March 2002  Tom Droege
#      rev                 20 March 2002  Now ver 1.0   
# Given a photom .cal file, this program brings all the observations
# of the same physical star together in one place.
#   
# This is a modification of Michael Richmond's calc_gammas_cal.pl
# which did other things to the data.  This version just collects
# all the observations of the same star together from a .cal file.  
# 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. The leading 
# number of each line indicates is a star sequence number.
# The second number indicates how many stars are in the group.
#
# usage:    sort -n +1 M*.cal | collect_stars.pl  
#
#   where
#               M*.cal   is a calibrated photometry file created
#                        by the pipeline
#               
# Input format:
#
# star_id ra dec jd "V" vmag vsig vflag "I" imag isig iflag           
#                        
#
# Output format:
# 
# seq_num num_obs avg_ra avg_dec jd_array[] 
# "V" vmag_array[] vsig_array avg_vsig vflag
# "I" imag_array[] isig_array avg_isig iflag
# 
# we need to add an overflow check on num_obs
#use strict;

$debug = 0;

# change these items to add more lists or list length
$number_of_lists = 40; #$ARGV[0];
$list_spacing = 200;
for ($i = 0; $i < $number_of_lists; $i++) {
$num_obs[$i] = 0;
$list_has_been_used[$i] = 0;  
#list_has_been_used is necessary to avoid a problem
#for a star at 0,0 i.e. we cannot test for a 0,0 entry 
#and be sure it is not a real star.
}

# other variables
$active_list_index = 0;
$next_list_to_unload = 0;
$iicl = 0;  # index into current list 
$num_obs = 0;
$line_count = 0;
$match_radius =  40;  $ARGV[0];    # arcseconds  ###
$no_empty_lists = 0;
$match = 0;
$seq_num = 1;

#temp values to make system run

$cur_ra = 9999.9;
$cur_dec = 9999.9;
$num_obs = 0;

# we accumulate values for statistical purposes in these arrays

$star_id[0] = 0;
$cur_ra[0] = 9999.0;    # star RA degrees
$cur_dec[0] = 9999.0;  # star dec degrees
$jd_array[0] = 9999.0; # Julian date of measurement
$vmag_array[0] = 0;     # V filter magnitude
$vsig_array[0] = 0;       # V filter measurement sigma
$vflag_array[0] = 0;      # V measurement error flags
$imag_array[0] = 0;      # I filter magnitude
$isig_array[0] = 0;        # I fliter measurement sigma  
$iflag_array[0] = 0;       # I measurement error flags

#####################################################
#                                                   # 
# Here we start the main loop to get a star         #
# from the list and try to match it with one of     #
#  number_of_lists in memory                        #
#                                                   #
#####################################################


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

##################################################
#                                                # 
#         First time through initialization      #
#         Put the first star on the first list   #
#                                                # 
##################################################

    if ($line_count eq 1) {
        # do it 
        $star_id[0] = $star_id;
        $cur_ra[0] = $ra;
        $cur_dec[0] = $dec;      # degrees
        $jd_array[0] = $jd;
        $vmag_array[0] = $vmag;
        $vsig_array[0] = $vsig;
        $vflag_array[0] = $vflag;
        $imag_array[0] = $imag;
        $isig_array[0] = $isig;
        $iflag_array[0] = $iflag;
        $next_list_to_unload = 0; 
        $num_obs[0] = 1;
        $list_has_been_used[0] = 1;
        next;
    }
###################################################
#                                                 #
#  First, look to see if the current star matches #
#  any list.                                      #
#                                                 #
###################################################
    $match = 0;
  
    for ($active_list_index = 0; $active_list_index < 
        $number_of_lists;            
        $active_list_index++) {
        
        if (coords_match ($ra, $dec, 
        $cur_ra[$active_list_index * $list_spacing], 
        $cur_dec[$active_list_index * $list_spacing],
        $match_radius) == 1) {
           $match = 1;
         
        }

        if ($match == 1) {
            # figure out index into the right list
           
            
          
            $iicl = ($active_list_index * $list_spacing) +                 
            $num_obs[$active_list_index];
      
            # put line on an existing list
            # do it here
            $star_id[$iicl] = $star_id;
            $cur_ra[$iicl] = $ra;
            $cur_dec[$iicl] = $dec;      # degrees
            $jd_array[$iicl] = $jd;
            $vmag_array[$iicl] = $vmag;
            $vsig_array[$iicl] = $vsig;
            $vflag_array[$iicl] = $vflag;
            $imag_array[$iicl] = $imag;
            $isig_array[$iicl] = $isig;
            $iflag_array[$iicl] = $iflag;
            $num_obs[$active_list_index]++;# index after loading
            last;
            #stop the loop.  We have found a match
            #and put the item on a list.  Later we
            #we will write code which catches a match
            #on more than one list       
        }
     }

###################################################
#                                                 #
#  Second, check if there are empty lists         #
#  available which can happen at the start of     #
#  the list filling process                       #
#                                                 #
###################################################                        
                       
                                      
    if (($match == 0) && ($no_empty_lists == 0)) {
    # check if there are empty lists available which can happen at
    #the start of the loading process
        $empty_list_to_use = 0;
        for ($i = 0; $i < $number_of_lists; $i++) {
            if ($list_has_been_used[$i] == 0) {
            $empty_list_to_use = $i;
            $list_has_been_used[$i] = 1;
            last;
            }
        }
        if ($empty_list_to_use == 0) {
        $no_empty_lists = 1;
        } 
       
        if ($empty_list_to_use > 0) {
            # there are empty lists put on next available
            # compute the index to us to first item on list
            $iicl = $empty_list_to_use * $list_spacing;
            # here put on list
            $star_id[$iicl] = $star_id;
            $cur_ra[$iicl] = $ra;
            $cur_dec[$iicl] = $dec;      # degrees
            $jd_array[$iicl] = $jd;
            $vmag_array[$iicl] = $vmag;
            $vsig_array[$iicl] = $vsig;
            $vflag_array[$iicl] = $vflag;
            $imag_array[$iicl] = $imag;
            $isig_array[$iicl] = $isig;
            $iflag_array[$iicl] = $iflag;
            $num_obs[$empty_list_to_use] = 1; #indes after loading
         }             
    } 
   
#######################################################
#                                                     #
# Third, if there are no empty lists, unload the      #
# oldest list and load the new line into the          #
# newly emptied list                                  #
#                                                     #
####################################################### 
               
    if (($match == 0) && ($no_empty_lists == 1)) {    
    # use next_list_to_unload
    # compute index into start of list
    
    $iicl = $next_list_to_unload * $list_spacing;
    
        if (1 == 1) {            
        # unload next_list_to_unload
            for ($i = $iicl; $i < ($iicl + 
                $num_obs[$next_list_to_unload]); $i++) {
                $star_id[$i] = $num_obs[$next_list_to_unload];
            }
            for ($i = $iicl; $i < ($iicl + 
                $num_obs[$next_list_to_unload]); $i++) {
printf  "%4d %3d %8.4f %6.4f %13.5f V %6.3f %4.3f %3d I %6.3f %4.3f %3d \n",
                $seq_num,
                $star_id[$i], 
                $cur_ra[$i], 
                $cur_dec[$i], 
                $jd_array[$i],
                $vmag_array[$i],                                           
                $vsig_array[$i], 
                $vflag_array[$i], 
                $imag_array[$i],
                $isig_array[$i], 
                $iflag_array[$i];
            }
            $seq_num++;
        }
    # load next_list_to_unload with new star
    # we are loading position 0
        $star_id[$iicl] = $star_id;
        $cur_ra[$iicl] = $ra;
        $cur_dec[$iicl] = $dec;      # degrees
        $jd_array[$iicl] = $jd;
        $vmag_array[$iicl] = $vmag;
        $vsig_array[$iicl] = $vsig;
        $vflag_array[$iicl] = $vflag;
        $imag_array[$iicl] = $imag;
        $isig_array[$iicl] = $isig;
        $iflag_array[$iicl] = $iflag;
        $num_obs[$next_list_to_unload] = 1;
        $next_list_to_unload++;
        if ($next_list_to_unload >= $number_of_lists) {
            $next_list_to_unload = 0;
        }
    }
}  # go get next star line 

# we are done but the lists are full, so unload all the lists

for ($j = 0; $j < $number_of_lists; $j++) {
    if ($num_obs[$j] < 1) {
       next; 
    }
    $next_list_to_unload = $j;
    # unload next_list_to_unload
    $iicl = $next_list_to_unload * $list_spacing;
    for ($i = $iicl; $i < ($iicl + 
        $num_obs[$next_list_to_unload]); $i++) {
        $star_id[$i] = $num_obs[$next_list_to_unload];
    }
    for ($i = $iicl; $i < ($iicl + 
        $num_obs[$next_list_to_unload]); $i++) {
printf  "%4d %3d %8.6f %6.6f %13.5f V %6.3f %4.3f %3d I %6.3f %4.3f %3d \n",
        $seq_num,
        $star_id[$i], 
        $cur_ra[$i], 
        $cur_dec[$i], 
        $jd_array[$i],
        $vmag_array[$i],                                           
        $vsig_array[$i], 
        $vflag_array[$i], 
        $imag_array[$i],
        $isig_array[$i], 
        $iflag_array[$i];
    }
    $seq_num++;   
}

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);
#print "in sub", " ", $_[0], " ", $_[1], " ", $_[2], " ", $_[3], "\n";  ###

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