#!/opt/perl/bin/perl
use Smart::Comments;
use GPS::DATES;
use GPS::Geodesy;
use GOA::Stacov;
use File::Basename;
my $progname = basename $0;
my $usage = qq(

$progname stacov_dir comb.stacov

);
my $dir = shift || die $usage;
my $comb = shift || die $usage;

opendir(my $dh,$dir) or die "Couldn't open $dir $!\n";

my @stacovs = map { "$dir/$_" } grep { /\.stacov/ } readdir($dh);
my $sn = GOA::Stacov->new( file => $comb );
my %sn_sites=();
my @sn_sta = $sn->sites;
@sn_sites{@sn_sta}=@sn_sta;
my ($nxyz,$xyz,$cov,$env,$ecov,$sig,$rho,$dy,$cal);
my $header = qq(#SITE  epoch (dec. year)  east (cm)   north (cm) up (cm)  sig_e (cm) sig_n (cm) sig_u (cm) rho_en rho_eu rho_nu 
);
my ( @env , @env_sig );
my ($sta,%rot);
my ($lat,$lon,$h);
foreach my $f ( @stacovs ) { ### working on $cal [===|   ] % done
#foreach my $f ( () ) { ### working on $cal [===|   ] % done
   $cal = unpack( 'a10',basename($f));
   #$cal = get_date($s->date,'cal');
   next if -e "$dir/$cal.envcov";
   my $s = GOA::Stacov->new( file => $f );
   $t = get_date($s->date,'dy');
   open(my $fh,'>',"$dir/$cal.envcov") or die "Couldn't open $dir/$cal.envcov $!\n";
   print $fh $header;
   foreach my $sta ( $s->sites ) {
     next if ! exists $sn_sites{$sta};
     $cov = $s->cov($sta,$sta,'STA');
     $xyz = $s->xyz($sta);
     next if ! $$xyz[0] || ! $$xyz[1] || ! $$xyz[2];
     $nxyz = $sn->xyz($sta);
     if ( ! exists $rot{$sta} ) {
       ($lat,$lon,$h) =   xyz2llh($nxyz);
       $rot{$sta} = env_rot($lat,$lon);
     }
     ($env,$ecov) = xyz2env($nxyz,$xyz-$nxyz,$cov,$rot{$sta});
     @$env = map { $_*1e2 } @$env;
     @$sig = map { sqrt($$ecov[$_][$_])*1e2 } ( 0 .. 2 );
     $$rho[0] = $$ecov[0][1]*1e4/$$sig[0]/$$sig[1];
     $$rho[1] = $$ecov[0][2]*1e4/$$sig[0]/$$sig[2];
     $$rho[2] = $$ecov[1][2]*1e4/$$sig[1]/$$sig[2];
     print $fh "$sta $t @$env @$sig @$rho\n";
   }
   close $fh;
}
my %gds=(); 
foreach my $f ( sort glob("$dir/*.envcov") ) { ### working on envcov $cal [====|   ] % done
   $cal = unpack('a10',$f);
   open(my $fh,'<',$f ) or die "Couldn't open $f $!\n";
   while( <$fh> ) {
     chomp;
     next if /^#/;
     ($sta,$t,@env) = split(' ');
     push @{$gds{$sta}{lon}},"$t $env[0] $env[3] $sta LON $cal\n";
     push @{$gds{$sta}{lat}},"$t $env[1] $env[4] $sta LAT $cal\n";
     push @{$gds{$sta}{rad}},"$t $env[2] $env[5] $sta RAD $cal\n";
    }
}
my $c;
foreach  $c ( qw( lat lon rad ) ) {
  foreach $sta ( sort keys %gds ) { ### working on $sta $c [===|  ] % done
     open( my $fh,'>',"series/$sta.$c") or die "Couldn't open $sta.$c $!\n";
     print {$fh} @{$gds{$sta}{$c}};
     close $fh;
  }
}

