#!perl 
#
# $Id: ConvQual.pl,v 1.2 2005-11-21 18:54:11 ashted Exp $
#
# Measure the quality of convergence by looking at the variation in lengths
# between struts in a knot (or knots).
#

use strict;
use warnings;
use POSIX;
use Getopt::Std;

our $VERSION = '$Revision: 1.2 $';
$VERSION =~ s/\s*\$//g;

sub HELP_MESSAGE {
  print STDERR <<"EOM";
  Produce a gnuplot input file which displays the variation in strut lengths
so as to show the quality of convergence of the knot.

usage: ConvQual -e <epsilon> <file> ...
  -e <epsilon> How wide the gnuplot boxes should be (default 0.01)
  <file> ...   Files to process.  New file will have the same name with 
               the ".vect" at the end replaced with ".stlen.gnuplot"
EOM
}

print "ConvQual $VERSION\n";

our($opt_e) = 0.01;
$Getopt::Std::STANDARD_HELP_VERSION = 1;
getopts('e:');

foreach my $file (@ARGV) {
  warn("Skipping file which does not end in .vect: $file\n"),next 
    unless $file =~ /.vect$/;
  (my $newfile = $file) =~ s/.vect$/.stlen.gnuplot/;
  open(STRUTS,"struts -e 1e20 --sv $file |")
    or die "Could not run struts on $file, stopped";
  my @lens;
  <STRUTS>;
  <STRUTS>;
  while (<STRUTS>) {
    last if /found/;
    print("skipping :$_:"),next unless /^[-0-9 ]/;
    my @parts = split;
    my $xd = $parts[4] - $parts[0];
    my $yd = $parts[5] - $parts[1];
    my $zd = $parts[6] - $parts[2];
    my $len = sqrt($xd*$xd + $yd*$yd + $zd*$zd);
    push @lens,$len;
  }
  close STRUTS;
  open(GNUPLOT,">$newfile") 
    or die "Could not open $newfile for writing, stopped";
  @lens = sort { $a <=> $b } @lens;
  my($smallest,$largest,$strutcount) = ($lens[0],$lens[-1],scalar(@lens));
  my $prevcnt = $strutcount;
  for (my $len = $smallest; $len < $largest; $len += $opt_e) {
    @lens = grep { $_ > $len + $opt_e } @lens;
    my $found = $prevcnt - scalar(@lens);
    print GNUPLOT $len+($opt_e/2)," $found $opt_e\n" if $found;
    unless ($found) {
      # This slot was empty, skip ahead to a new one
      last unless scalar(@lens);
      my $skip_steps = floor(($lens[0]-$len)/$opt_e);
      $len += $skip_steps*$opt_e;
    }
    $prevcnt = scalar(@lens);
  }
  print GNUPLOT join("\n",@lens),"\n";
  close GNUPLOT;
  my $start = 0;
  print "$file: $strutcount $smallest $largest\n";
}
