#!/opt/perl/bin/perl
use GOA::Stacov;
use GPS::Geodesy;
use Fortran::F90Format;
use File::Basename;

my $f = shift || die;
my $s = GOA::Stacov->new( file => $f );
my $gmt = Fortran::F90Format->new( fmt => "2(F22.15,1X),5(F9.4,1X),A4");
my ($xyz,$vxyz,$vcov,$lat,$lon,$h,$venv,$vecov,$sig,$rho);
my $out = basename $f;
print STDERR "intput: $out\n";
$out =~ s/\.stacov//;
open(my $hh,'>',"$out.hvel.gmt") or die "Couldn't open $out.hvel.gmt $!\n";
open(my $vh,'>',"$out.vvel.gmt") or die "Couldn't open $out.vvel.gmt $!\n";
foreach my $sta ( $s->sites ) {
  $xyz = $s->xyz($sta,'STA');
  $vxyz = $s->xyz($sta,'VEL');
  $vcov = $s->cov($sta,$sta,'VEL');
  ($lat,$lon,$h) = xyz2llh($xyz);
  ($venv,$vecov) = xyz2env($xyz,$vxyz,$vcov);
  @$sig = map { sqrt($$vecov[$_][$_]) } ( 0 .. 2 );
  $$rho[0] = $$vecov[0][1]/$$sig[0]/$$sig[1];
  $$rho[1] = $$vecov[0][2]/$$sig[0]/$$sig[2];
  $$rho[2] = $$vecov[1][2]/$$sig[1]/$$sig[2];

  print {$hh} $gmt->write($lon,$lat,$$venv[0],$$venv[1],$$sig[0],$$sig[1],$$rho[0],$sta);
  print {$vh} $gmt->write($lon,$lat,0,$$venv[2],0,$$sig[2],0,$sta);
  
  #print "$lat $lon $h @$venv @$sig @$rho\n";
  
}
print STDERR "output: $out.hvel.gmt and $out.vvel.gmt\n";
