#!/usr/bin/perl
#Attached is a more useful version of the same script which
#works the same way. I have put in two setting: 1) the minimum
#number of valid pairs before a W-S will be generated as
#output [default: 5], and 2) the baseline photometric sigma
#[default: 0.010]. I think you'll find almost all the stars
#with a W-S statistic greater than 1 will now be variables!

# PW: Added to output: (line printed only if WS > 1)
# - Number of nights observed
# - Significance of Spearman's Rank Correlation Coefficient (%)
#   For each night of I data
# - Significance of Spearman's Rank Correlation Coefficient (%)
#   For each night of V data

#Cheers,
#Doug

#==============================================================
#  Douglas L Welch     | Office/voicemail (905) 525-9140 x23186
#  Physics & Astronomy | FAX              (905) 546-1252
#  McMaster University |
#  Hamilton, Ontario   |
#  Canada L8S 4M1      | E-mail       welch@physics.mcmaster.ca
#==============================================================

#--------------------------------------------
# TASS Welch-Stetson routine by D.L. Welch, May 2002
#    Based on  Welch, D.L & Stetson, P.B. 1993, AJ, 105, 1813-1821
#
#   - requires at least 5 pairs now (min_valid)
#--------------------------------------------
# Format of input file is:
#   0 == seq_num
#   1 == num_obs (number of VI pairs)
#   2 == avg_ra
#   3 == avg_dec
#   4 == jd
#   5 == "V"
#   6 == vmag
#   7 == vsig
#   8 == vflag
#   9 == "I"
#  10 == imag
#  11 == isig
#  12 == iflag
#--------------------------------------------
# Output columns:
#   1 == seq_num
#   2 == num_obs (number of VI pairs)
#   removed: 3 == nok (number of pairs with neither sigma equal to zero)
#   3 == avg_ra
#   4 == avg_dec
#   5 == v_avg
#   6 == vsig_avg
#   7 == i_avg
#   8 == isig_avg
#   9 == total number of nights
#  10 == Welch-Stetson statistic
#  11 == Mean of nightly W-S
#  12 == number of nights used for mean of nightly W-S
#  13 == Largest nightly outlier (expressed in standard deviations from the nightly mean)
#  14 == Difference between maximum and minimum nightly average in V, divided by largest nightly standard deviation in V
#  15 == Difference between maximum and minimum nightly average in I, divided by largest nightly standard deviation in I
#
#  Remarks: 
#    output parameter 11 should be large for short period variables (variation during a night)
#    output parameters 14 and 15 should be large for long period variables (i.e. large differences between nightly means)
#--------------------------------------------
# Settings:
#   min_valid  == minimum number of valid pairs for generating
#                 the Welch-Stetson statistic
#   min_valid_night  == minimum number of valid pairs for generating
#                 the Welch-Stetson statistic for a night
#   base_sigma == photometric error to be added in quadrature
#                 to input errors
#   max_outlier == maximum number of standard deviations an observation
#                 may differ from the mean, for it to be taken into account for the mean W-S, or nightly averages
#--------------------------------------------
#
# Only data with both V and I flags = 0 and Vsig and Isig > 0 are taken into account
#
#--------------------------------------------

$line_num = 0;
$star_num = 0;
$ltvok_num = 0;

$min_valid = 30;
$min_valid_night = 10; # Minimum observations for a night
$base_sigma = 0.01;
$max_outlier = 4.0;

#use VS::Period;

sub Sort
# Quicksort
# Arguments: 
# - Array on which to sort ($SortKey)
# - First index in array
# - Last index in array
# - @SortArray: List of other arrays to be sorted (should not include $SortKey)
{
	my ($SortKey, $l, $r, @SortArray) = @_;
	my $i = $l;
	my $j = $r; 
	my $mid = $$SortKey[int(($l + $r)/2)];
	my $array;
	do
	{
		while ($$SortKey[$i] < $mid) { $i++; }
		while ($mid < $$SortKey[$j]) { $j--; }
		if ($i <= $j)
		{
			for $array ($SortKey, @SortArray)
			{
				my $tmp = $$array[$i];
				$$array[$i] = $$array[$j];
				$$array[$j] = $tmp;
			}
			$i++;
			$j--;
		}
	} until $i > $j;
	if ($l < $j) { Sort($SortKey, $l, $j, @SortArray); }
	if ($i < $r) { Sort($SortKey, $i, $r, @SortArray); }
}


sub treat_night ($$) {
   my ($from, $to) = @_;
   if ($to - $from + 1 >= $min_valid_night) {
      # Make averages for the night
      $v_num_sum = 0.0;
      $v_den_sum = 0.0;
      $v_sqr_sum = 0.0;
      $i_num_sum = 0.0;
      $i_den_sum = 0.0;
      $i_sqr_sum = 0.0;
      $nok = 0;

      for my $i ($from..$to) {
         $vwgt = 1./($v_sig[$i]**2+$base_sigma**2);
         $v_num_sum += $v_obs[$i]*$vwgt;
         $v_den_sum += $vwgt;

         $iwgt = 1./($i_sig[$i]**2+$base_sigma**2);
         $i_num_sum += $i_obs[$i]*$iwgt;
         $i_den_sum += $iwgt;

         $nok++;
      }

      $v_avg = $v_num_sum/$v_den_sum;
      $i_avg = $i_num_sum/$i_den_sum;
      $ws = 0.0;

      # Calculate standard deviations
      for my $i ($from..$to) {
         $vwgt = 1./($v_sig[$i]**2+$base_sigma**2);
         $v_sqr_sum += $vwgt * ($v_obs[$i] - $v_avg)**2;

         $iwgt = 1./($i_sig[$i]**2+$base_sigma**2);
         $i_sqr_sum += $iwgt * ($i_obs[$i] - $i_avg)**2;
      }

      $v_stdev = sqrt($v_sqr_sum/$v_den_sum);
      $i_stdev = sqrt($i_sqr_sum/$i_den_sum);

      # Remove outliers
      for my $i ($from..$to) {
         # Largest deviation ?
         $v_dev = abs($v_obs[$i]-$v_avg)/$v_stdev;
         if ($v_dev > $max_dev) {
            $max_dev = $v_dev;
         }
         $i_dev = abs($i_obs[$i]-$i_avg)/$i_stdev;
         if ($i_dev > $max_dev) {
            $max_dev = $i_dev;
         }
         if ($v_dev > $max_outlier || $i_dev > $max_outlier) {
            $outlier[$i] = 1;
            $vwgt = 1./($v_sig[$i]**2+$base_sigma**2);
            $v_num_sum -= $v_obs[$i]*$vwgt;
            $v_den_sum -= $vwgt;
            $v_sqr_sum -= $vwgt * ($v_obs[$i] - $v_avg)**2;

            $iwgt = 1./($i_sig[$i]**2+$base_sigma**2);
            $i_num_sum -= $i_obs[$i]*$iwgt;
            $i_den_sum -= $iwgt;
            $i_sqr_sum += $iwgt * ($i_obs[$i] - $i_avg)**2;

            $nok--;
         }
      }

      if ($nok >= $min_valid_night) {
         # Recalculate averages
         $v_avg = $v_num_sum/$v_den_sum;
         $i_avg = $i_num_sum/$i_den_sum;
         $v_stdev = sqrt($v_sqr_sum/$v_den_sum);
         $i_stdev = sqrt($i_sqr_sum/$i_den_sum);
         if ($nights_ok) {
            if ($v_avg < $v_max) {
              $v_max = $v_avg;
            }
            if ($v_avg > $v_min) {
              $v_min = $v_avg;
            }
            if ($i_avg < $i_max) {
              $i_max = $i_avg;
            }
            if ($i_avg > $i_min) {
              $i_min = $i_avg;
            }
            if ($v_stdev > $v_max_stdev) {
               $v_max_stdev = $v_stdev;
            }
            if ($i_stdev > $i_max_stdev) {
               $i_max_stdev = $i_stdev;
            }
         } else {
            $v_max = $v_avg;
            $v_min = $v_avg;
            $v_max_stdev = $v_stdev;
            $i_max = $i_avg;
            $i_min = $i_avg;
            $i_max_stdev = $i_stdev;
         }

         # Form residuals and W-S index

         for my $i ($from..$to) {
            if (!$outlier[$i]) {
               $v_sigma = sqrt($v_sig[$i]**2 + $base_sigma**2);
               $i_sigma = sqrt($i_sig[$i]**2 + $base_sigma**2);
               $del_v = ($v_obs[$i]-$v_avg)/$v_sigma;
               $del_i = ($i_obs[$i]-$i_avg)/$i_sigma;
               $ws = $ws + $del_v*$del_i;
            }
         }

         # Normalize W-S index
         $ws2 += $ws/sqrt($nok*($nok-1.0));
         $nights_ok++;
      }

   }
}


# Read line from STDIN

  while (<>) {
    chop;
    ($seq_num, $num_lines, $ra, $dec, $jd,
     $vmag, $vsig, $vflag, $imag, $isig, $iflag)
     = (split(' ',$_))[0..4,6..8,10..12];
    $star_num++;
    $line_num++;

    $sum_ra = $ra;
    $sum_dec = $dec;
    $num_obs = 0;

    if (!$vflag && !$iflag && $vsig && $isig) {
       $jd_obs[++$#jd_obs] = $jd;
       $v_obs[++$#v_obs] = $vmag;
       $v_sig[++$#v_sig] = $vsig;
       $i_obs[++$#i_obs] = $imag;
       $i_sig[++$#i_sig] = $isig;
       $num_obs++;
    }

    for $i (1..$num_lines-1) {
       $tmp = <>;
       chop $tmp;
       ($seq_num, $num_lines, $ra, $dec, $jd,
        $vmag, $vsig, $vflag, $imag, $isig, $iflag)
        = (split(' ',$tmp))[0..4,6..8,10..12];
       $line_num++;

       $sum_ra += $ra;
       $sum_dec += $dec;

       if (!$vflag && !$iflag && $vsig && $isig) {
          $jd_obs[++$#jd_obs] = $jd;
          $v_obs[++$#v_obs] = $vmag;
          $v_sig[++$#v_sig] = $vsig;
          $i_obs[++$#i_obs] = $imag;
          $i_sig[++$#i_sig] = $isig;
          $num_obs++;
       }

    }

    if (!$num_obs) {
       $ltvok_num++;
       next;
    }

    # Calculate WS index per night
    # Sort on JD
    Sort (\@jd_obs, 0, $num_obs-1, \@v_obs, \@v_sig, \@i_obs, \@i_sig);

    $jd = 0;
    $first_j = 0;

    $ws2  = 0.0;
    $nights = 0;
    $nights_ok = 0;
    $max_dev = 0.0;

    for $j (0..$num_obs-1) {
       if ($jd_obs[$j] > $jd + 0.7) {
          # New night starts
          if ($jd) {
             # Treat the previous night
             treat_night($first_j, $j-1);
          }
          $nights++;
          $jd = $jd_obs[$j];
          $first_j = $j;
       }
    }
    # Treat last night
    treat_night($first_j, $num_obs-1);

    # Make the normal W-S index
   
    # Create averages
    $v_num_sum = 0.0;
    $v_den_sum = 0.0;
    $i_num_sum = 0.0;
    $i_den_sum = 0.0;
    for $i (0..$num_obs-1) {
       $vwgt = 1./($v_sig[$i]**2+$base_sigma**2);
       $v_num_sum += $v_obs[$i]*$vwgt;
       $v_den_sum += $vwgt;

       $iwgt = 1./($i_sig[$i]**2+$base_sigma**2);
       $i_num_sum += $i_obs[$i]*$iwgt;
       $i_den_sum += $iwgt;
    }
    $v_avg = $v_num_sum/$v_den_sum;
    $i_avg = $i_num_sum/$i_den_sum;

    # Form residuals and W-S index

    $ws  = 0.0;
    $vsig_sum = 0.0;
    $isig_sum = 0.0;

    for $i (0..$num_obs-1) {
       $v_sigma = sqrt($v_sig[$i]**2 + $base_sigma**2);
       $i_sigma = sqrt($i_sig[$i]**2 + $base_sigma**2);
       $del_v = ($v_obs[$i]-$v_avg)/$v_sigma;
       $del_i = ($i_obs[$i]-$i_avg)/$i_sigma;
       $ws = $ws + $del_v*$del_i;
       $vsig_sum += $v_sig[$i];
       $isig_sum += $i_sig[$i];
    }

    # Normalize W-S index

    printf("%5d %3d %8.4f %8.4f %5.2f %4.2f %5.2f %4.2f %2d",
           $seq_num, $num_obs, $sum_ra/$num_lines, $sum_dec/$num_lines, $v_avg, $vsig_sum/$num_obs, $i_avg, $isig_sum/$num_obs, $nights);
    if ($num_obs >= $min_valid) {
       $ws = $ws/sqrt($num_obs*($num_obs-1.0));
       printf(" %7.3f", $ws);
    } else {
       printf(" %7s", "-");
       $ltvok_num++;
    }
    if ($nights_ok) {
       printf(" %7.3f %2d", $ws2/$nights_ok, $nights_ok);
    } else {
       printf(" %7s %2d", "-", 0);
    }
    if ($max_dev) {
       printf(" %4.1f", $max_dev);
    } else {
       printf(" %4s", "-");
    }
    if ($nights_ok > 1) {
       printf(" %4.1f %4.1f", ($v_min - $v_max)/$v_max_stdev, ($i_min - $i_max)/$i_max_stdev);
    }
    print "\n";
	

    undef @jd_obs;
    undef @outlier;
    undef @v_obs;
    undef @v_sig;
    undef @i_obs;
    undef @i_sig;

  }

  printf(STDERR "Total number of lines processed: %6d\n", $line_num);
  printf(STDERR "Total number of stars processed: %6d\n\n", $star_num);
  printf(STDERR "Number of objects with less than %d valid pairs : %6d\n",
                 $min_valid, $ltvok_num);