#!/opt/perl/bin/perl
use strict;
use File::Basename;
use Pod::Text::Termcap;
use Getopt::Long;
use GOA::StaPos;
use GOA::Stacov;
use GPS::Defaults;
use Matrix;
my $progname = basename $0;
my %units = ( cm => 100, mm => 1000 , m => 1.0 , km => 1/1000.00 );
my ($stacov_in,$sta_pos,$output,$help);
my ($full,$list) = (1,0);
my $t = 0; # treshold cm
my $output_unit = 'm';
my $args = @ARGV;
my $f =''; 

GetOptions( 'i:s' => \$stacov_in,
            'o:s' => \$output,
            'full'   => \$full,
            'list'   => \$list,
            't:f'    => \$t,
            'u:s'    => \$output_unit,
            'pos_file:s' => \$f,
            'help' => \$help,
          ) or die "Error in options\n";

if ( !$args || $help ) {
  usage();
  exit 1;
}
$full = 0 if $list;
my %par = get_defaults( env => \%ENV );
my $u = $units{lc $output_unit};
$f= "$par{STA_INFO_DIR}/sta_pos" if ! $f;
my $dir =dirname $f;
my $file=basename $f;
my $pos = GOA::StaPos->new( file => $f );
my $s   = GOA::Stacov->new( file => $stacov_in );
my ($sta,$delta,$sig,$dxyz,$sx_sig,$o,$px,$sx);
$t*=$u;

format STDOUT =
@<<<< @<<<<<<<<<<<<<<< +- @<<<<<<<<<<<<<< @*
$sta, $delta,             $sig,           $output_unit
.

my %seen=();
foreach $sta ( sort grep { ! $seen{$_}++ } $s->sites  ) { 
  if ( ! $pos->get($sta) ) {
      warn "No sta_pos entry for $sta \n";
      next;
  }
  $px  = $pos->xyz($sta);
  $sx  = $s->xyz($sta);
  $sx_sig = $s->xyz_sig($sta);
  my @sx_size=$sx->dim;
  my @px_size=$px->dim;
  if (  @{$px} != 3 ) {
     warn "No sta_pos entry for $sta\n";
     next;
  }
  if (  @{$sx} != 3 && @{$px} != @{$sx} ) {
     warn "different number of columns for $sta\n".
           "@sx_size != @px_size\n";
     next;
  }
  $dxyz = $px-$sx;
  $delta = sqrt($dxyz.$dxyz); 
  $sig = sqrt($sx_sig.$sx_sig);
  $delta *= $u;
  $sig*=$u;
  $o = $full ? "$sta $delta +- $sig $output_unit " : $sta ;
  $o =  $delta > $t ? $o : next;   
  #print "$o\n";
  write STDOUT;
}


sub usage {
  my $pod=Pod::Text::Termcap->new();
  $pod->parse_from_filehandle(man_page());
}

sub man_page {
my $man =qq(
=head1 NAME

$progname

Computes the distance from the XYZ in the sta_pos file and the XYZ in the input stacov file and shows the ones that are greater than a given treshold (default 10 cm). Use -t 0 to show all.

=head1 OPTIONS

=over 3

=item -i stacov_file

        Input stacov file

=item -t treshold 

        Defines minimal treshold (default 0 )

=item -units [cm | mm | m | km ] 

        Select output unit for distances (default: m)

=item -full  

       Show distances and sigmas with units

=item -list

       Show only the list of stations 

=back

=head1 EXAMPLES

=over 3

=item To show which stations have a difference of more than 1m 
      with respect to their sta_pos entry, do:

=item $progname -i 2008-11-05.comb.stacov -t 1 -unit m


=back

=head1 AUTHOR

Victor Marcelo Santillan <marcelo\@geology.cwu.edu>

=head1 COPYRIGHT

Copyright \(c\) 2009

Central Washington University

=cut
);
open(my $fh,"<",\$man ) or die "$!\n";
return $fh
}

