#!/usr/bin/perl -w

use POSIX;

my @x = ();
my @y = ();
my $avgx = 0.0;
my $avgy = 0.0;
my $nstations = 0;
my $unit = 'm';
while(<>){
  if(/^\s*([-0-9.]+)\s+([-0-9.]+)\s+[-0-9.]+/){
    my ($sx, $sy) = ($1, $2);

    push @x, $sx;
    push @y, $sy;

    $nstations = @x;                       # Running means for more accuracy.
    $avgx += ($sx - $avgx) / $nstations;
    $avgy += ($sy - $avgy) / $nstations;
  }
  elsif(/^#.*\s+\(km\)/){
    $unit = 'km';
  }
}

my $medx = value_at_percentile(50.0, @x);
my $medy = value_at_percentile(50.0, @y);

my $avgsep     = 0.0;
my $rmsb     = 0.0;
my $rmsx     = 0.0;
my $rmsy     = 0.0;
my $avgdist    = 0.0;
my $avgdistmav = 0.0;
my $avgdistmed = 0.0;
my $nb = 0;
my @dist    = ();
my @distmav = ();
my @distmed = ();
my @baselines = ();
for(my $i = 0; $i < @x; ++$i){
  my $d    = sqrt($x[$i]**2 + $y[$i]**2);
  my $dmav = sqrt(($x[$i] - $avgx)**2 + ($y[$i] - $avgy)**2);
  my $dmed = sqrt(($x[$i] - $medx)**2 + ($y[$i] - $medy)**2);

  push @dist,    $d;
  push @distmav, $dmav;
  push @distmed, $dmed;
  $avgdist    += ($d    - $avgdist)    / ($i + 1);
  $avgdistmav += ($dmav - $avgdistmav) / ($i + 1);
  $avgdistmed += ($dmed - $avgdistmed) / ($i + 1);

  for(my $j = $i + 1; $j < @x; ++$j){
    ++$nb;
    my $x2 = ($x[$j] - $x[$i])**2;
    my $y2 = ($y[$j] - $y[$i])**2;
    my $b2 = $x2 + $y2;
    my $baseline = sqrt($b2);
    push @baselines, $baseline;
    $avgsep += ($baseline - $avgsep) / $nb;
    $rmsx  += ($x2 - $rmsx) / $nb;
    $rmsy  += ($y2 - $rmsy) / $nb;
    $rmsb  += ($b2 - $rmsb) / $nb;
  }
}
$rmsx = sqrt($rmsx);
$rmsy = sqrt($rmsy);
$rmsb = sqrt($rmsb);

@dist      = sort {$a <=> $b} @dist;
@distmav   = sort {$a <=> $b} @distmav;
@distmed   = sort {$a <=> $b} @distmed;
@baselines = sort {$a <=> $b} @baselines;

if($unit ne 'km' && $avgsep >= 1.0){
  $unit = 'km';
  $avgsep     *= 0.001;
  $rmsb       *= 0.001;
  $rmsx       *= 0.001;
  $rmsy       *= 0.001;
  $avgdist    *= 0.001;
  $avgdistmav *= 0.001;
  $avgdistmed *= 0.001;
  foreach my $item (@dist, @distmav, @distmed, @baselines){
    $item *= 0.001;
  }
}

$sqrtrmsxrmsy = sqrt($rmsx * $rmsy);
printf "Average separation:                     %8.3f $unit\n", $avgsep;
printf "RMS     separation:                     %8.3f $unit\n", $rmsb;
printf "RMS     x:                              %8.3f $unit\n", $rmsx;
printf "RMS     y:                              %8.3f $unit\n", $rmsy;
printf "sqrt(rms_x * rms_y):                    %8.3f $unit\n", $sqrtrmsxrmsy;
printf "Average distance from nominal center:   %8.3f $unit\n", $avgdist;
printf "Average distance from average position: %8.3f $unit\n", $avgdistmav;
printf "Average distance from median position:  %8.3f $unit\n", $avgdistmed;

printf "\nLongest baseline:                       %8.3f $unit\n", $baselines[-1];
printf "Shortest baseline:                      %8.3f $unit\n", $baselines[0];

# Convert baseline in km to resolution in " at 300 GHz.
#               1.22     c        "/deg    deg/pi   300 GHz      pi
my $as_over_b = 1.22 * 2.9978e5 * 3600.0 * 180.0 / (3.00e11 * 3.1415927); 

my $hres = $as_over_b / sqrt($baselines[-1]**2 + $baselines[0]**2);
my $nres = $as_over_b / $avgsep;

# If you calculate the FWHM CLEAN beam axes by fitting a gaussian to the main
# lobe of the dirty beam, i.e. using the 2nd derivatives at 0, you get
# FWHM = sqrt((2Nln2)/sum(b_i**2))/pi.
my $rmsres = 0.374781250258555 * 2.9978e5 * 3600.0 * 180.0 / (3.00e11 *
							      3.1415927 *
							      $sqrtrmsxrmsy);

printf "\nHoldaway resolution: %8.6f\" at 300 GHz\n", $hres;
printf "Nominal resolution acc. to avgsep:  %8.6f\" at 300 GHz\n", $nres;
printf "Nominal resolution acc. to rmssep:  %8.6f\" at 300 GHz\n", $rmsres;

print "\n                    Cumulative distributions ($unit)\n";
print "  %   Dist. from nom. center  Dist. from avg.  Dist. from med.  Baseline\n";
foreach my $p (20, 25, 50, 75, 100){
  printf "%3d       %7.2f               %7.2f          %7.2f        %7.2f\n",
    $p, value_at_percentile($p, @dist), value_at_percentile($p, @distmav),
      value_at_percentile($p, @distmed), value_at_percentile($p, @baselines);
}

# Returns the interpolated value at the percentile % point of @array, which is
# assumed to be sorted.
sub value_at_percentile {
  my ($percentile, @array) = @_;
  my $nelem = @array;
  my $value = 0;

  my $interpolated_index = $percentile * $nelem / 100.0 - 1;
  my $lindex = floor($interpolated_index);
  my $uindex = ceil($interpolated_index);

  if($lindex < 0){
    warn "Extrapolating to unobserved percentile $percentile";
    $value = $array[0] + $interpolated_index * ($array[1] - $array[0]);
  }
  elsif($uindex < $nelem){
    $value = $array[$lindex] + ($interpolated_index - $lindex) *
      ($array[$uindex] - $array[$lindex]);
#    printf "+%7.2f - %7.2f (%d, %7.2f, %d, %d)\n", $array[$uindex],
#      $array[$lindex], $lindex, $interpolated_index, $uindex, $nelem;
  }
  else{
    print "Using total\n";
    $value = $array[-1];
  }
  return $value;
}