#!/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!

#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 == Largest nightly outlier (expressed in standard deviations from the nightly mean)
#  11 == Number of outliers
#  12 == Welch-Stetson statistic
#  13 == Mean of nightly W-S
#  14 == number of nights used for short period variability indicators
#  15 == V short period variability indicator = 1 - mean square successive difference/2/variance
#  16 == I short period variability indicator
#  17 == number of pairs used for the above two indicators
#  18 == Difference between maximum and minimum nightly average in V, divided by largest nightly standard deviation in V (of the minimum or maximum night)
#  19 == Difference between maximum and minimum nightly average in I, divided by largest nightly standard deviation in I (of the minimum or maximum night)
#  20 == number of nights used for long period variability indicators
#
#  Remarks: 
#    output parameters 12, 17 and 18 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_long  == minimum number of valid pairs for a night for taking part
#                 in the long term (internight) variability indicators
#   min_valid_night_short  == minimum number of valid pairs for a night for generating 
#                 short term (intranight) variability indicators
#   It is assumed that $min_valid_night_long <= $min_valid_night_short
#   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_long = 5; 
$min_valid_night_short = 20; 
$base_sigma = 0.01;
$max_outlier = 3.0;

$Nbad = 0;
if ($bad)
{	# A file with a list of bad nights is given (lines with two entries: fromJD toJD)
	((open BAD, $bad) || (die "Bad file $bad does not exist."));
	while (<BAD>)
	{
		($badfrom[++$Nbad], $badto[$Nbad]) = (split(' '))[0..1];
	}
}

sub badnight($)
{
	my $jd = shift;
	for my $i (1..$Nbad)
	{
		if (($jd >= $badfrom[$i]) && ($jd <= $badto[$i]))
		{
			return 1;
		}
	}
	return 0;
}

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_long) {
      # 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;
            $outliers++;
            $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_long) {
         # 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_long) {
            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_max == $v_avg) || ($v_min == $v_avg)) && ($v_stdev > $v_max_stdev)) {
               $v_max_stdev = $v_stdev;
            }
            if ((($v_max == $v_avg) || ($v_min == $v_avg)) && ($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;
         }
         $nights_ok_long++;

         if ($nok >= $min_valid_night_short) {
            # Other short period variability indicators
            # Calculate sum of (squared difference between subsequent observations/variance)

            $vvi_night = 0;
            $ivi_night = 0;

            for my $i0 ($from..$to-1) {
               if (!$outlier[$i0]) {
                  # Get the next measurement
                  $i1 = $i0 + 1;
                  while (($i1 <= $to) && $outlier[$i1]) { $i1++; }
                  if ($i1 > $to) { last; } # No further measurements
                  $vvi_night += ($v_obs[$i1] - $v_obs[$i0])**2;
                  $ivi_night += ($i_obs[$i1] - $i_obs[$i0])**2;
                  $pairs++;
               }
            }

            $vvi += $v_den_sum * $vvi_night/$v_sqr_sum;
            $ivi += $i_den_sum * $ivi_night/$i_sqr_sum;

            # 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_short++;
         }
      }

   }
}

my $tassid = 0; # Current TASS star ID


sub treatstar()
{
    if (!$num_obs) {
       $ltvok_num++;
       return;
    }

    # 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_long = 0;
    $nights_ok_short = 0;
    $max_dev = 0.0;
    $vvi = 0.0;
    $ivi = 0.0;
    $pairs = 0;
    $outliers = 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;
    $ws_obs = 0;
    for $i (0..$num_obs-1) {
      if (!$outlier[$i]) {
        $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;
        $ws_obs++;
      }
    }
    $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) {
      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;
        $vsig_sum += $v_sig[$i];
        $isig_sum += $i_sig[$i];
      }
    }

  if (($nights_ok_long > 1) || ($ws_obs >= $min_valid) || ($nights_ok_short)) {
    printf("%5d %3d %8.4f %8.4f %5.2f %4.2f %5.2f %4.2f %2d",
           $tassid, $num_obs, $sum_ra/$num_lines, $sum_dec/$num_lines, $v_avg, $vsig_sum/$num_obs, $i_avg, $isig_sum/$num_obs, $nights);
    if ($max_dev) {
       printf(" %4.1f %3d", $max_dev, $outliers);
    } else {
       printf(" %4s %3s", "-", "-");
    }
    if ($ws_obs >= $min_valid) {
       # Normalize W-S index
       $ws = $ws/sqrt($ws_obs*($ws_obs-1.0));
       printf(" %6.2f", $ws);
    } else {
       printf(" %6s", "-");
       $ltvok_num++;
    }
    if ($nights_ok_short) {
       printf(" %6.2f %2d", $ws2/$nights_ok_short, $nights_ok_short);
    } else {
       printf(" %6s %2d", "-", 0);
    }
    if ($pairs) {
#       printf(" %4.2f %4.2f %3d", sqrt($pairs/$vvi), sqrt($pairs/$ivi), $pairs);
       printf(" %5.2f %5.2f %3d", 1 - $vvi/$pairs/2, 1 - $ivi/$pairs/2, $pairs);
    } else {
       printf(" %5s %5s %3s", "-", "-", "-");
    }
    if ($nights_ok_long > 1) {
       printf(" %4.1f %4.1f %2d", ($v_min - $v_max)/$v_max_stdev, ($i_min - $i_max)/$i_max_stdev, $nights_ok_long);
    } else {
       printf(" %4s %4s %2d", "-", "-", $nights_ok_long);
    }
    print "\n";
  }

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

}


# Read line from STDIN

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

    if ($seq_num != $tassid)
    {	# Do the calculations for the previous star
	if ($tassid)
	{
		treatstar;
	}

	# New star
	$tassid = $seq_num;
        $star_num++;

        $sum_ra = $ra;
        $sum_dec = $dec;
        $num_obs = 0;
        $num_lines = 1;
    }
    else
    {	# Still the same star
        $sum_ra += $ra;
        $sum_dec += $dec;
        $num_lines += 1;
    }

    $line_num++;

    if (!$vflag && !$iflag && $vsig && $isig && !badnight($jd)) {
       $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++;
    }

  }

  # Do the calculations for the final star
  if ($tassid)
  {
    treatstar;
  }

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