#!/usr/local/bin/perl 

# snedit.pl  Version 1.0  - Jim Lovell, ATNF, 22-Apr-1999
# snedit.pl  Version 1.2  - Glen Langston, NRAO, 15-July-1999
# snedit.pl  Version 1.5  - Glen Langston, NRAO, 22-July-1999
#

# An SN table editing program because I could never get SNSMO to do
# what I wanted and SNEDT is difficult to use. Requires a SN table
# written by AIPS using TBOUT with the DOCRT adverb set large enough
# so that all data for a single (IF, Antenna, Time) are on the same
# line. DOCRT = 1000 is a safe setting.

#
# The SN table may then be flagged in a similar Manner to flagging
# in Difmap. Different sources have different plotting symbols
#
# usage
#      snedit.pl  <infile> <outfile>
#
print "\nsnedit v1.5: measure and edit an AIPS SN table \n";

use PGPLOT;

my $progname = "snedit.pl";
my ($day, $time, $tsys, $dummy, $i);
my @subs;
my %color = (black => 0,white => 1,red => 2,green => 3,blue => 4, 
             cyan => 5, magenta => 6, yellow => 7, orange => 8);

my $nargs = @ARGV;

if ($nargs < 2){
    print "\nUsage: $progname [input SN file] [output SN file]\n";
    die;
}

my $infile = $ARGV[0];
my $outfile = $ARGV[1];

open(INLOG,"$infile") || die "cannot open $infile: $!\n";
open(OLOG,">$outfile") || die "cannot open $outfile: $!\n";

my $nl = 0;

$line = <INLOG>;
# just write out the header stuff without changing it
while ( $line !~ /BEGIN\*PASS/ ) {
 
  print OLOG $line;
  $line = <INLOG>;
  next;
}

print OLOG $line;

  $line = <INLOG>;

# now we're into the data
while ($line !~ /END\*PASS/ ) {

  # we asume two IF's
  ($if1ro[$nl], $if1td[$nl], $if1ti[$nl], $if1so[$nl], $if1an[$nl], 
   $if1su[$nl], $if1fi[$nl], $if1fr[$nl], $if1no[$nl], $if1mb[$nl],
   $if1re[$nl], $if1im[$nl], $if1de[$nl], $if1ra[$nl], $if1wt[$nl],
   $if1ref[$nl]) = split ' ', $line;

  $line = <INLOG>;
  ($if2ro[$nl], $if2td[$nl], $if2ti[$nl], $if2so[$nl], $if2an[$nl], 
   $if2su[$nl], $if2fi[$nl], $if2fr[$nl], $if2no[$nl], $if2mb[$nl],
   $if2re[$nl], $if2im[$nl], $if2de[$nl], $if2ra[$nl], $if2wt[$nl],
   $if2ref[$nl]) = split ' ', $line;

  $if1td[$nl] =~ s/D/E/;
  # convert to hours
  $if1td[$nl] = $if1td[$nl]*24.;
  $if2td[$nl] = $if1td[$nl];
  $if2ti[$nl] = $if1ti[$nl];
  $if2so[$nl] = $if1so[$nl];
  $if2an[$nl] = $if1an[$nl];
  $if2su[$nl] = $if1su[$nl];
  $if2fi[$nl] = $if1fi[$nl];
  $if2fr[$nl] = $if1fr[$nl];
  $if2no[$nl] = $if1no[$nl];
  $if2mb[$nl] = $if1mb[$nl];

  # third IF is difference of two
  $if3su[$nl] =   $if1su[$nl] - $if2su[$nl];
  $if3fi[$nl] =   $if1fi[$nl] - $if2fi[$nl];
  $if3fr[$nl] =   $if1fr[$nl] - $if2fr[$nl];
  $if3no[$nl] =   $if1no[$nl] - $if2no[$nl];
  $if3mb[$nl] =   $if1mb[$nl] - $if2mb[$nl];
  $if3de[$nl] =   $if1de[$nl] - $if2de[$nl];
  $if3ra[$nl] =   $if1ra[$nl] - $if2ra[$nl];
  if ($if2wt[$nl] < $if1wt[$nl]) {
    $if3wt[$nl] =  $if2wt[$nl];
  } else {
    $if3wt[$nl] =  $if1wt[$nl];
  }
  
  $if3td[$nl] = $if1td[$nl];
  $if3so[$nl] = $if1so[$nl];
  $if3an[$nl] = $if1an[$nl];
  $if3su[$nl] = $if1su[$nl];
  $if3fi[$nl] = $if1fi[$nl];
  $if3fr[$nl] = $if1fr[$nl];
  $if3no[$nl] = $if1no[$nl];
  $if3mb[$nl] = $if1mb[$nl];
  # calculate the phase from the angle of the real and imaginary parts
  $if1pha[$nl] = atan2( $if1im[$nl], $if1re[$nl]);
  $if2pha[$nl] = atan2( $if2im[$nl], $if2re[$nl]);
  $if3pha[$nl] = $if1pha[$nl] - $if2pha[$nl];
  $if3re[$nl]  = cos( $if3pha[$nl]);
  $if3im[$nl]  = sin( $if3pha[$nl]);
  $if3pha[$nl] = atan2( $if3im[$nl], $if3re[$nl]);
  # convert phases to degrees
  $if1pha[$nl] = $if1pha[$nl] * 57.21685979;
  $if2pha[$nl] = $if2pha[$nl] * 57.21685979;
  $if3pha[$nl] = $if3pha[$nl] * 57.21685979;

  $nl++;
  $line = <INLOG>;
} # end for all input lines

$nl--;    #assume lastline is ***END*PASS***, ignore
 
close(INLOG);

# find out how many antennas
($minant, $maxant) = minmax(\@if1an, 0.0);

# begin plotting.
#
# Divide screen into 3 with the top being SNR, middle Delay, bottom Rate
# pgbegin(0,"/x",1,3);
pgbegin(0,"/xserve",1,3);
pgask(0);

# start with IF 1, antenna 1
($o_tmin, $o_tmax) = minmax(\@if1td, 0.1);
#($smin, $smax) = minmax(\@if1wt, 0.1);
#($dmin, $dmax) = minmax(\@if1de, 0.1);
#($rmin, $rmax) = minmax(\@if1ra, 0.1);
$tmin = $o_tmin;
$tmax = $o_tmax;

# main loop
my $ch = "";
my $x1 = 0;
my $x2 = 0;
my $y1 = 0;
my $y2 = 0;
my $redraw = 1;
my $ant = 1;
my $if = 1;
my $igflag = 0;
my $bothif = 0;

while ($ch ne "x") {

  if ($redraw) {
    if ($if == 1) {
      plot_sdr($ant, $tmin, $tmax, $nl, \@if1td, \@if1wt, 
	       \@if1pha, \@if1de, \@if1ra, \@if1an, \@if1so);
    } elsif ($if == 3) {
      plot_sdr($ant, $tmin, $tmax, $nl, \@if3td, \@if3wt, 
	       \@if3pha, \@if3de, \@if3ra, \@if3an, \@if3so);
    } else {
      plot_sdr($ant, $tmin, $tmax, $nl, \@if2td, \@if2wt, 
	       \@if2pha, \@if2de, \@if2ra, \@if2an, \@if2so);
    }
    $redraw = 0;
  } # end if time to redraw

  pgband(6,0,0.0,0.0, $x1, $y1, $ch); $ch = uc2lc($ch);

  if (($ch eq "h") || ($ch eq "?")) {
    print "$progname help\n";
    print "--------------\n\n";
    print "X : (right mouse button) Exit\n";
    print "N : Display the next antenna.\n";
    print "P : Display the previous antenna.\n";
    print "I : Display the next IF.\n";
    print "A : (Left mouse button) flag/restore the nearest point in time.\n";
    print "M : Display the mean values of a time interval.\n";
    print "C : Flag over a time range\n";
    print "R : Restore solutions over a time range.\n";
    print "B : Toggle flagging both or a single IF at once\n";
    print "U : Zoom in on a narrower time range.\n";
    print "L : Redraw the screen.\n";
    print "S : Toggle to ignore flagged data when scaling.\n";
    print "\n";

  } elsif ($ch eq "u") {
    # timerange
    print "Display a time range of solutions.\n";
    if (($tmin > $o_tmin) || ($tmax < $o_tmax)) {
      $tmin = $o_tmin;
      $tmax = $o_tmax;
    } else {
      print "\tUse left mouse button to select minimum and maximum times\n";
      print "\tor press X (or right mouse to cancel)\n";
      print "Press U again to display the full range\n";
      pgband(6, 0, $x1, $y1, $x1, $y1, $ch); $ch = uc2lc($ch);
      if ($ch eq "a") {
        pgband(4, 0, $x1, $y1, $x2, $y2, $ch); $ch = uc2lc($ch);
        if ($ch eq "a") {
          ($tmin, $tmax) = sort numerically ($x1, $x2);
        }
      }
    }
    $redraw = 1;
    $ch = "";
  } elsif ($ch eq "b") {
    $bothif = !$bothif;
    if ($bothif) {
      print "Both IFs will be flagged/restored simultaneously.\n";
    } else {
      print "IF's will be flagged/restored separately.\n";
    }
  } elsif ($ch eq "c") {
    # timerange
    print "Flagging a range of solutions.\n";
    print "\tUse the left mouse button to select the minimum and maximum times\n";
    print "\tor press X (or right mouse to cancel)\n";
    pgsci(2);
    pgband(6, 0, $x1, $y1, $x1, $y1, $ch); $ch = uc2lc($ch);
    if ($ch ne "x") {
      pgband(4, 0, $x1, $y1, $x2, $y2, $ch); $ch = uc2lc($ch);
    } else {
      print "Cancelled\n";
      $ch = "";
    }
    if ($ch ne "x") {
      if ($bothif) {
	  flag_range($x1, $x2, \@if1td, \@if1wt, \@if1an, $ant, 1);
	  flag_range($x1, $x2, \@if2td, \@if2wt, \@if2an, $ant, 1);
	  flag_range($x1, $x2, \@if3td, \@if3wt, \@if3an, $ant, 1);
      } else {
	if ($if == 1) {
	  flag_range($x1, $x2, \@if1td, \@if1wt, \@if1an, $ant, 1);
	} elsif ($if == 2) {
	  flag_range($x1, $x2, \@if2td, \@if2wt, \@if2an, $ant, 1);
	} elsif ($if == 3) {
	  flag_range($x1, $x2, \@if3td, \@if3wt, \@if3an, $ant, 1);
	}
      }
      $redraw = 1;
    } else { 
      print "Cancelled\n";
      $ch = "";
    }
    pgsci(1);
  } elsif ($ch eq "m") {
    # timerange
    print "Taking the mean of a range of solutions.\n";
    print "\tUse the left mouse button to select the minimum and maximum times\n";
    print "\tor press X (or right mouse to cancel)\n";
    pgsci(2);
    pgband(6, 0, $x1, $y1, $x1, $y1, $ch); $ch = uc2lc($ch);
    if ($ch ne "x") {
      pgband(4, 0, $x1, $y1, $x2, $y2, $ch); $ch = uc2lc($ch);
    } else {
      print "Cancelled\n";
      $ch = "";
    }
    if ($ch ne "x") {
  mean_range($x1, $x2, \@if1pha, \@if1de, \@if1ra, \@if1td, \@if1wt, \@if1an, $ant, 1);
  mean_range($x1, $x2, \@if2pha, \@if2de, \@if2ra, \@if2td, \@if2wt, \@if2an, $ant, 2);
  mean_range($x1, $x2, \@if3pha, \@if3de, \@if3ra, \@if3td, \@if3wt, \@if3an, $ant, 3);
      $redraw = 1;
    } else { 
      print "Cancelled\n";
      $ch = "";
    }
    pgsci(1);

  } elsif ($ch eq "r") {
    # timerange
    print "Restoring a range of solutions.\n";
    print "\tUse the left mouse button to select the minimum and maximum times\n";
    print "\tor press X (or right mouse to cancel)\n";
    pgsci(3);
    pgband(6, 0, $x1, $y1, $x1, $y1, $ch); $ch = uc2lc($ch);
    if ($ch ne "x") {
      pgband(4, 0, $x1, $y1, $x2, $y2, $ch); $ch = uc2lc($ch);
    } else {
      print "Cancelled\n";
      $ch = "";
    }

    if ($ch ne "x") {
      if ($bothif) {
	  flag_range($x1, $x2, \@if1td, \@if1wt, \@if1an, $ant, 0);
	  flag_range($x1, $x2, \@if2td, \@if2wt, \@if2an, $ant, 0);
	  flag_range($x1, $x2, \@if3td, \@if3wt, \@if3an, $ant, 0);
      } else {
	if ($if == 1) {
	  flag_range($x1, $x2, \@if1td, \@if1wt, \@if1an, $ant, 0);
	} elsif ($if == 2) {
	  flag_range($x1, $x2, \@if2td, \@if2wt, \@if2an, $ant, 0);
        } else {
	  # flag both ifs if difference is bad
	  flag_range($x1, $x2, \@if1td, \@if1wt, \@if1an, $ant, 0);
	  flag_range($x1, $x2, \@if2td, \@if2wt, \@if2an, $ant, 0);
	  flag_range($x1, $x2, \@if3td, \@if3wt, \@if3an, $ant, 0);
	}
      }
      $redraw = 1;
    }  else {
      print "Cancelled\n";
      $ch = "";
    }

    pgsci(1);
  } elsif ($ch eq "a") {
    # flag a point
    if ($bothif) {
      flag_nearest($x1, \@if1td, \@if1wt, \@if1an, $ant);
      flag_nearest($x1, \@if2td, \@if2wt, \@if2an, $ant);
      flag_nearest($x1, \@if3td, \@if3wt, \@if3an, $ant);
    } else {
      if ($if == 1) {
	flag_nearest($x1, \@if1td, \@if1wt, \@if1an, $ant);
      } elsif ($if == 2) {
	flag_nearest($x1, \@if2td, \@if2wt, \@if2an, $ant);
      } elsif ($if == 3) {
	# flag both ifs if difference is bad
	flag_nearest($x1, \@if1td, \@if1wt, \@if1an, $ant);
	flag_nearest($x1, \@if2td, \@if2wt, \@if2an, $ant);
	flag_nearest($x1, \@if3td, \@if3wt, \@if3an, $ant);
      }
    }
    $redraw = 1;
  } elsif ($ch eq "l") {
    #redraw
    $redraw = 1;
  } elsif ($ch eq "n") {
    $ant++;
    if ($ant > $maxant) {$ant = 1;}
    print "Displaying solutions for Antenna $ant.\n";
    #redraw
    $redraw = 1;
  } elsif ($ch eq "p") {
    $ant--;
    if ($ant < 1) {$ant = $maxant;}
    print "Displaying solutions for Antenna $ant.\n";
    #redraw
    $redraw = 1;
  } elsif ($ch eq "i") {
    $if++;
    if ($if > 3.0) {$if = 1;}
    if ($if > 2.0) {
	print "Displaying IF 2 - IF 1\n";
    } else {
	print "Changing to IF $if.\n";
    }
    $redraw = 1;
  } elsif ($ch eq "s") {
    $igflag = !$igflag;
    if ($igflag) {
      print "Flagged solutions will be ignored when scaling plots.\n";
    } else {
      print "Flagged solutions will be used in scaling plots.\n";
    }
    $redraw = 1;
  }
} # end while not exiting (not "X")

# stop pgplot service
pgend;

# now write out the data in TBIN format
for($i = 0; $i <= $nl; $i++) {
  #convert back to days;
  $if1td[$i] = $if1td[$i]/24.;
  $if2td[$i] = $if2td[$i]/24.;
  $if1td[$i] =~ s/E/D/;
  $if2td[$i] =~ s/E/D/;
  printf OLOG "%8d   %s   %s %10i %10i %10i %10i   %s %10i   %s  %s   %s". 
    "   %s   %s%15.6e %10i\n", $if1ro[$i], $if1td[$i], $if1ti[$i], $if1so[$i], 
    $if1an[$i], $if1su[$i], $if1fi[$i],  $if1fr[$i], $if1no[$i], $if1mb[$i],
    $if1re[$i], $if1im[$i], $if1de[$i], $if1ra[$i],   $if1wt[$i],  $if1ref[$i];
  printf OLOG "%8d                      ''             ''         ''         ''".
    "         ''         ''             ''         ''             ''    %s   %s". 
    "   %s   %s%15.6e %10i\n", $if2ro[$i], 
    $if2re[$i], $if2im[$i], $if2de[$i], $if2ra[$i],   $if2wt[$i],  $if2ref[$i];
}
  print OLOG "***END*PASS***\n";
close(OLOG);

######################## end of main function

sub plot_sdr {
  # plot_sdr() produces 3 plots for one IF of one Antenna
  # the plots are phase, delay and rate, using auto-scaling
  my ($ant, $x1, $x2, $nl, 
      $xref, $sref, $sPhase, $dref, $rref, $antref, $srcref) = @_;
  my ($xtxt, $ytxt, $tiqtle, @time, @snr, @phase, @delay, @rate, @srcok, 
      @srcbad, $n); 

  my @symb_ok  = (2, 13, 16, 17);
  my @symb_bad = (5,  7, 6,  4 );

  # form an array of time, snr, delay and rate for $ant
  $n = 0;
  for ($i = 0; $i <= $nl; $i++) {
    if ($$antref[$i] == $ant && $tmin < $$xref[$i] && $tmax > $$xref[$i]) {
      $time[$n]  = $$xref[$i];
      $snr[$n]   = $$sref[$i];
      $phase[$n] = $$sPhase[$i];
      $delay[$n] = 1.0e+9 * $$dref[$i];
      $rate[$n]  = 1.0e+12 * $$rref[$i];
      $srcok[$n] = $symb_ok[($$srcref[$i] % 4) - 1];
      $srcbad[$n] = $symb_bad[($$srcref[$i] % 4) - 1];
      $n++;
    }
  } # end for all input lines
  $n--;

  if ($n < 0) {
    print "No valid solutions found for antenna $ant.\n";
    print "Flagged pints will be plotted next time.\n";
    $igflag = 0;
    return;
  }

  # find the maximum time range in all 3 arrays
  ($xmin, $xmax) = minmax2(\@time, \@snr, 0.1);
  
  pgsch(2.0);
  $xtxt = "";
  $ytxt = "Phase";
  $title = "Fringe Fitting Solutions. Antenna $ant, IF $if";
  if ($if > 2) {
    $title = "Fringe Fitting Solutions. Antenna $ant, IF 1 - 2";
  }

  pgsci(7);
  ($ymin, $ymax) = minmax2(\@phase, \@snr, 0.1);
  if ($ymin == $ymax) {
      $ymax = $ymax + 10;
      $ymin = $ymin - 10;
  }
  pgenv($xmin, $xmax, $ymin, $ymax, 0, 0);
  pgsci(1);
  pglabel($xtxt, $ytxt, $title);
  for ($i = 0; $i <= $n; $i++) {
    if ($snr[$i] > 0.0) {
      pgsci(3);
      pgpoint(1,$time[$i],$phase[$i],$srcok[$i]);
    } else {
      pgsci(2);
      pgpoint(1,$time[$i],$phase[$i],$srcbad[$i]);
    }
  }

  $xtxt = "";
  $ytxt = "Delay (ns)";
  $title = "";
#  pgpage;
  pgsci(7);
  ($ymin, $ymax) = minmax2(\@delay, \@snr, 0.1);
  if ($ymin == $ymax) {
    if ($ymin == 0.0) {
      $ymin = -0.1; $ymax = 0.1;
    } else {
      $ymax = $ymax + 0.1*$ymax;
      $ymin = $ymin - 0.1*$ymin;
    }
  }

  pgenv($xmin, $xmax, $ymin, $ymax, 0, 0);
  pgsci(1);
  pglabel($xtxt, $ytxt, $title);
  for ($i = 0; $i <= $n; $i++) {
    if ($snr[$i] > 0.0) {
      pgsci(3);
      pgpoint(1,$time[$i],$delay[$i],$srcok[$i]);
    } else {
      pgsci(2);
      pgpoint(1,$time[$i],$delay[$i],$srcbad[$i]);
    }
  }

  $xtxt = "Time (hours)";
  $ytxt = "Rate (ps/s)";
  $title = "";
  pgsci(7);
  ($ymin, $ymax) = minmax2(\@rate, \@snr, 0.1);
  if ($ymin == $ymax) {
    $ymax = $ymax + (0.1*$ymax) + 1;
    $ymin = $ymin - (0.1*$ymin) - 1;
  }
  pgenv($xmin, $xmax, $ymin, $ymax, 0, 0);
  pgsci(1);
  pglabel($xtxt, $ytxt, $title);
  for ($i = 0; $i <= $n; $i++) {
    if ($snr[$i] > 0.0) {
      pgsci(3);
      pgpoint(1,$time[$i],$rate[$i],$srcok[$i]);
    } else {
      pgsci(2);
      pgpoint(1,$time[$i],$rate[$i],$srcbad[$i]);
    }
  }

  pgsci(1);
} # end of plot_sdr()


sub minmax {
  # takes an array and provides the minimum and maximum. If
  # $expand_fac is > 0, decrease $min and increase $max by 
  # $expand_fac * ($max - $min).
  my ($arr_ref, $expand_fac) = @_;
  my @sorted = sort numerically (@$arr_ref);
  my $min = $sorted[0];
  my $max = $sorted[$#sorted];

  if ($expand_fac > 0.0) {
    my $fact = $expand_fac * ($max - $min);
    $max = $max + $fact;
    $min = $min - $fact;
  }
  return ($min, $max);
} # end of minmax()

sub minmax2 {
  # does same as minmax but accompanying array is used to determine if
  # point is used to determine the range

  my ($arr_ref, $flg_ref, $expand_fac) = @_;
  my ($min, $max, $i);
  my @arr = @$arr_ref;
  my @newarr;
  my @flag = @$flg_ref;
  my $n = 0;
  for ($i = 0; $i <= $#arr; $i++) {
    if (($flag[$i] > 0.0) || (($flag[$i] <= 0.0) && (!$igflag))) {
      $newarr[$n] = $arr[$i];
      $n++;
    }
  }
  $n--;
  if ($n < 0) {
    ($min, $max) = minmax(\@arr, $expand_fac);
  } else {
    my @sorted = sort numerically (@newarr);
    $min = $sorted[0];
    $max = $sorted[$#sorted];
    if ($expand_fac > 0.0) {
      my $fact = $expand_fac * ($max - $min);
      $max = $max + $fact;
      $min = $min - $fact;
    }
  } # end if several values 

  return ($min, $max);
} #end of minmax2()

sub numerically { $a <=> $b; }

sub uc2lc {
  # converts upper case to lower case
  # call $lc_string = uc2lc($string);
    local($lc_string);
    $lc_string = "\L$_[0]\E";
}

sub flag_nearest {
  # flag_nearest() finds the closest point an flags the values
  # flagging is done by making the weights negative
  my ($x, $xref, $okref, $antref, $ant) = @_;
  my @xa = @$xref;
  my @ants = @$antref;
  my $dist = 99999.0;
  my $newdist = $dist;
  my $nearest = 0;

  for ($i = 0; $i <= $#xa; $i++) {
    if ($ants[$i] == $ant) {
      $newdist = ($xa[$i] - $x)**2;

      if ($newdist < $dist) {
	$nearest = $i;
	$dist = $newdist;
      }
    }
  } #end for all points

  $$okref[$nearest] = -abs($$okref[$nearest]);
} # end of flag_nearest

sub flag_range {
  # flag range will toggle the flagged values in a range
  # flagging is done by making the weights negative
  my ($x1, $x2, $xref, $okref, $antref, $ant, $flag) = @_;
  ($x1, $x2) = sort numerically ($x1, $x2);

  my @xa = @$xref;
  my @ants = @$antref;
  my $i;
  for ($i = 0; $i <= $#xa; $i++) {
    # if time (hours) is within range 
    if (($xa[$i] >= $x1) && ($xa[$i] <= $x2))  {
      if ($ants[$i] == $ant) {
	if ($flag) {
	  # Flag it
	  $$okref[$i] = -abs($$okref[$i]);
	} else {
	  # Restore it
	  $$okref[$i] = abs($$okref[$i]);
	} # end else if un-flagging
      } # end if required antenna
    } # end if in required time range
  } # end for all data values
} # end of flag range

sub mean_range {
  # mean_range() gets the mean value of parameters as in a time range
  # when processing IF 3, it is assumed that IFs 1 and 2 have already
  # been processed by this function, so that saved average values may be used
  my ($x1, $x2, $inPhase, $inDelay, $inRate, $xref, $okref, $antref, $ant, $ifNumber) = @_;
  ($x1, $x2) = sort numerically ($x1, $x2);

  my @xa = @$xref;
  my @ants = @$antref;
  my @Ws = @$okref;
  my @phase = @$inPhase;
  my @rate  = @$inRate;
  my @delay = @$inDelay;
  my $i = 0, $j = 0;
  my $npoints = 0, $nAcc = 0;
  my $nphases = 0;
  my $aveDelay = 0, $aveRate = 0, $avePhase = 0, $aveAcc = 0;
  my $aveTDel  = 0, $aveTRat = 0, $aveTPha = 0, $aveTAcc = 0;
  my $refFreq  = 4.93E9, $minTime = -1;
  my $a1 = 0, $a2 = 0, $a3 = 0, $a4 = 0, $dt = 0;

  for ($i = 0; $i < $#xa-1; $i++) {
    for ($j = $i + 1; $j < $#xa; $j++) {
      # if sn entry in time range and also antenna required
      if (($xa[$i] >= $x1) && ($xa[$i] <= $x2) && ($ants[$i] == $ant)) {
        if (($xa[$j] >= $x1) && ($xa[$j] <= $x2) && ($ants[$j] == $ant)) {
	  # if not flagged
	  if ($Ws[$i] > 3 && $Ws[$j] > 3) {
	      # convert from hours to seconds time difference
	      $dt = ($xa[$j] - $xa[$i])*3600.;
	      if ($dt > 0) {
		  $acc[$nAcc] = ($rate[$j] - $rate[$i])/$dt;
		  $accT[$nAcc] = ($xa[$j] + $xa[$i])/2.;
		  $nAcc++;
	      }
	  } # end if not flagged 
        } # end if second measurement is OK
      } # end if in time range and antenna of interest
    } # end for all second measurments
  } # end for all values in SN table
  print "N acc =", $nAcc, "\n";

  ($aveTRat,$aveRate) =aveRangeAnt($x1,$x2, \@xa, \@rate,  \@Ws, \@ants, $ant);
  ($aveTDel,$aveDelay)=aveRangeAnt($x1,$x2, \@xa, \@delay, \@Ws, \@ants, $ant);

  # save values for averaging, when replacing with a model
  $ifRate[$ifNumber] = $aveRate;
  $ifTRat[$ifNumber] = $aveTRat;
  $ifDelay[$ifNumber]= $aveDelay;
  $ifTDel[$ifNumber] = $aveTRat;

  print "Ant ", $ant, " IF ", $ifNumber;
  printf " Delay %10.3e %s", $aveDelay, " (s)" ;
  printf " Rate  %10.3e %s", $aveRate, " (s/s)\n"; 

  # if displaying difference of IFs, then calculate phase average
  if ($ifNumber == 3) {
    ($aveTPha,$avePhase) =aveRangeAnt($x1,$x2, \@xa, \@phase, \@Ws, \@ants, $ant);
    printf " Phase %8.3f %s at %8.3f (h\n)", $avePhase, " (d)", $aveTPha; 
  } # end if calculating phase value

  if ($ifNumber == 3) {
      print "\n Type Y to replace with model values; else Ignore \n";

      pgband(6, 0, $a1, $a2, $a3, $a4, $ch); $ch = uc2lc($ch);  

      if ($ch eq "y") {

	  $aveRate  = ($ifRate[1]  + $ifRate[2])/2;
	  $aveTRat  = ($ifTRat[1]  + $ifTRat[2])/2;
	  $aveDelay = ($ifDelay[1] + $ifDelay[2])/2;
	  $aveTDel  = ($ifTDel[1]  + $ifTDel[2])/2;

	  for ($i = 0; $i <= $#xa; $i++) {
	      # if sn entry in time range and also antenna required
	      if (($xa[$i] >= $x1) && ($xa[$i] <= $x2) && ($ants[$i] == $ant)){

		  $if1ra[$i] = $aveRate +(($xa[$i]-$aveTAcc)*3600.*$aveAcc);
		  $if1de[$i] = $aveDelay+(($xa[$i]-$aveTRat)*3600.*$if1ra[$i]);
		  $if2ra[$i] = $aveRate +(($xa[$i]-$aveTAcc)*3600.*$aveAcc);
		  $if2de[$i] = $aveDelay+(($xa[$i]-$aveTRat)*3600.*$if2ra[$i]);
		  $if3de[$i] =   $if1de[$i] - $if2de[$i];
		  $if3ra[$i] =   $if1ra[$i] - $if2ra[$i];
	      } # end if in time range, and correct antenna
	  } # end for all values;
          print "Updated model:\n";
          printf " Delay %10.3e %s", $aveDelay, " (s)" ;
          printf " Rate  %10.3e %s", $aveRate, " (s/s)\n"; 
      } else { 
          print "Data not changed\n";
      }
  } # end if replacing with model 
} # end of mean_range()

sub aveRangeAnt {
  # aveRangeAnt() returns the average of some set of values for a given
  # input range ($x1..$x2) and antenna $ant.
  # output is a pair of values, the X of the average and the average value Y
  # the two extreama Y values in the range are discarded before averaging
  my ($x1, $x2, $inXs, $inYs, $okref, $antref, $ant) = @_;
  ($x1, $x2) = sort numerically ($x1, $x2);

  my @Xs = @$inXs;
  my @Ys = @$inYs;
  my @ants = @$antref;
  my @Ws = @$okref;
  my $i = 0, $npoints = 0;
  my $sumY = 0, $maxY = 0, $minY = 0;
  my $sumX = 0, $maxX = 0, $minX = 0;

  # for all input values
  for ($i = 0; $i <= $#Xs; $i++) {
      # if sn entry in X range and also antenna required
      if (($Xs[$i] >= $x1) && ($Xs[$i] <= $x2) && ($ants[$i] == $ant)) {
	  # if weight ok (not flagged)
	  if ($Ws[$i] > .1) {
	      $sumX = $sumX + $Xs[$i];
	      $sumY = $sumY + $Ys[$i];
	      # on first point, init extreama values
	      if ($npoints == 0) {
		  $minY = $Ys[$i];
		  $maxY = $Ys[$i];
		  $minX = $Xs[$i];
		  $maxX = $Xs[$i];
	      } else {
		  # else if an extreama, save
		  if ($Ys[$i] < $minY) {
		      $minX = $Xs[$i];
		      $minY = $Ys[$i];
		  } elsif ($Ys[$i] > $maxY) {
		      $maxX = $Xs[$i];
		      $maxY = $Ys[$i];
		  }
	      } # end else not first point 
	      $npoints++;
	  } # end if not flagged 
      } # end if in X range and antenna of interest
  } # end for all Ys

  # if several points, exclude min and max
  if ($npoints > 3) {
      $sumX = ($sumX - $maxX - $minX) / ($npoints - 2);
      $sumY = ($sumY - $maxY - $minY) / ($npoints - 2);
  } elsif ($npoints > 0) {
      $sumX = $sumX / $npoints;
      $sumY = $sumY / $npoints;
  }
  print "sum X, Y: ", $sumX, ",", $sumY, " (", $npoints,")\n";
  return ($sumX, $sumY);
} # end of aveRangeAnt() 

sub aveRange {
  # aveRange() returns the average of some set of values for a given
  # for an input set of coordianate array pair
  # output is a pair of values, the X of the average and the average value Y
  # the two extreama Y values in the range are discarded before averaging
  my ($xref, $values) = @_;

  my @Xs = @$xref;
  my @Ys = @$values;
  my $i = 0, $npoints = 0;
  my $sumY = 0, $maxY = 0, $minY = 0;
  my $sumX = 0, $maxX = 0, $minX = 0;

  # for all input values
  for ($i = 0; $i <= $#Xs; $i++) {
      $sumX = $sumX + $Xs[$i];
      $sumY = $sumY + $Ys[$i];
      # on first point, init extreama values
      if ($npoints == 0) {
	  $minY = $Ys[$i];
	  $maxY = $Ys[$i];
	  $minX = $Xs[$i];
	  $maxX = $Xs[$i];
      } else {
          # else if an extreama, save
	  if ($Ys[$i] < $minY) {
	      $minX = $Xs[$i];
	      $minY = $Ys[$i];
	  } elsif ($Ys[$i] > $maxY) {
	      $maxX = $Xs[$i];
	      $maxY = $Ys[$i];
	  }
      } # end else not first point 
      $npoints++;
  } # end for all Ys

  # if several points, exclude min and max
  if ($npoints > 3) {
      $sumX = ($sumX - $maxX - $minX) / ($npoints - 2);
      $sumY = ($sumY - $maxY - $minY) / ($npoints - 2);
  } elsif ($npoints > 0) {
      $sumX = $sumX / $npoints;
      $sumY = $sumY / $npoints;
  }
  print "sum X, Y: ", $sumX, ",", $sumY, " (", $npoints, ")\n";
  return ($sumX, $sumY);
} # end of aveRange()