#!/opt/perl/bin/perl  
#
# cat together individual station stacov files
#
# Marcelo Santillan 2006-03-06
# Based on catcov, to allow parameter names used by posmrg
# AAS : amplitude annual sin
# AAC : amplitude annual cos 
# ASS : amplitude semi-annual sin
# ASC : amplitude semi-annual cos 
# JMP : Antenna jumps
# OFF : Offset events 
#
################################################
if ($#ARGV < 1) {
  print "\n  catcov useage:  catstacov file1 file2 ...\n";
  print "        where the input files are in stacov format \n";
  print "        output is to standard out.\n";
  exit 2;
}
$i = 0;
$j = 0;
$q = 0;
$prevref = "no_frame";
$ant_ht_src = "no_ant_ht_src";
open(SVEC,"> /tmp/svec.$$") || die "Can't open input file ",$_,"\n";
################################################
for(@ARGV){
  open(INPUT,$ARGV[$i]) || die "Can't open input file ",$_,"\n";
  $infile=$_;
  chop($line = <INPUT>);
  ($nparm[$i], $junk, $junk, $ttag[$i]) = split(' ',$line);
  for ($k = 0; $k < $nparm[$i]; $k++){
    chop($line = <INPUT>);
    ($iparam[$j],$sta,$type,$coord,$value[$j],$junk,$sigma[$j]) = split(' ',$line);
    $pname[$j] = $sta." ".$type." ".$coord."      ";
    if ($junk ne '+-') {
      ($iparam[$j],$sta,$type,$coord,$event,$value[$j],$junk,$sigma[$j]) = split(' ',$line);
      $pname[$j] = $sta." ".$type." ".$coord." ".$event;
    }
    if( $type eq "STA"){
     $psta[$s] = $j;
     $s++;
    }
    elsif ( $type eq "VEL"){
     $pvel[$v] = $j;
     $v++;
    }
    elsif ( $type eq "AAS"){
     $paas[$aas] = $j;
     $aas++;
    }
    elsif ( $type eq "AAC"){
     $paac[$aac] = $j;
     $aac++;
    }
    elsif ( $type eq "ASS"){
     $pass[$ass] = $j;
     $ass++;
    }
    elsif ( $type eq "ASC"){
     $pasc[$asc] = $j;
     $asc++;
    }
    elsif ( $type eq "JMP"){
     $pjmp[$jmp] = $j;
     $jmp++;
    }
    elsif ( $type eq "OFF"){
     $poff[$off] = $j;
     $off++;
    }
    $j++;
  }
#############################################################
  while(<INPUT>){
    chop;
    (@dline) = split(' ');
    if ( $#dline == 2) {
     $id1[$q] =  $dline[0] + $j - $nparm[$i];
     $id2[$q] =  $dline[1] + $j - $nparm[$i];
     $corr[$q] = $dline[2];
     $q++;
    }
    elsif ( $dline[1] eq "ANTENNA") {
      print SVEC $_,"\n";
    }
    else {
     $t++;
    }
    if (/Reference/) {
      if ($prevref ne "no_frame" && $prevref ne $dline[$#dline]) {
	 print STDERR "WARNING: different reference frames. $prevref != $dline[$#dline] in $infile\n";
      }
      $prevref = $dline[$#dline];
    }
    if (/ant hght/) {
      if ($_ ne $ant_ht_src && $ant_ht_src ne "no_ant_ht_src") {
	 print STDERR "NOTE: different antenna height sources in $infile\n";
	 $ant_ht_src = " !possibly inconsistent antenna height info source\n";
      } else {
        $ant_ht_src=$_
      }
    }
  }
  $i++;
}
close (SVEC);
################################################
if ($s+$v > 9999) {
  die "stacov format cant deal with more than 9999 parameters\n";
}
printf(" %4d PARAMETERS ON %s\n",$j, $ttag[0]);
# first positions
for ( $l = 0; $l < $s; $l++) {
   $p = $psta[$l];
   printf(" %4d  %s  %22.15e  +-  %21.15e\n", 
   $l+1, $pname[$p], $value[$p],$sigma[$p]);
   $ni[$p] = $l + 1;
#   printf("lpn %3d %3d %3d\n", $l , $p, $ni[$p]);
}
# next velocities
for ( $l = 0; $l < $v; $l++) {
   $p = $pvel[$l];
   printf(" %4d  %s  %22.15e  +-  %21.15e\n",
   $l + $s + 1, $pname[$p], $value[$p],$sigma[$p]);
   $ni[$p] = $l + $s + 1;
#   printf("lpn %3d %3d %3d\n", $l , $p, $ni[$p]);
}
for ( $l = 0; $l < $aas; $l++) {
   $p = $paas[$l];
   printf(" %4d  %s  %22.15e  +-  %21.15e\n",
   $l + $s + $v + 1, $pname[$p], $value[$p],$sigma[$p]);
   $ni[$p] = $l + $s + $v + 1;
#   printf("lpn %3d %3d %3d\n", $l , $p, $ni[$p]);
}
for ( $l = 0; $l < $aac; $l++) {
   $p = $paac[$l];
   printf(" %4d  %s  %22.15e  +-  %21.15e\n",
   $l + $s + $v + $aas + 1, $pname[$p], $value[$p],$sigma[$p]);
   $ni[$p] = $l + $s + $v + $aas + 1;
#   printf("lpn %3d %3d %3d\n", $l , $p, $ni[$p]);
}
for ( $l = 0; $l < $ass; $l++) {
   $p = $pass[$l];
   printf(" %4d  %s  %22.15e  +-  %21.15e\n",
   $l + $s + $v + $aas + $aac + 1, $pname[$p], $value[$p],$sigma[$p]);
   $ni[$p] = $l + $s + $v + $aas + $aac + 1;
#   printf("lpn %3d %3d %3d\n", $l , $p, $ni[$p]);
}
for ( $l = 0; $l < $asc; $l++) {
   $p = $pasc[$l];
   printf(" %4d  %s  %22.15e  +-  %21.15e\n",
   $l + $s + $v + $aas + $aac + $ass + 1, $pname[$p], $value[$p],$sigma[$p]);
   $ni[$p] = $l + $s + $v + $aas + $aac + $ass + 1;
#   printf("lpn %3d %3d %3d\n", $l , $p, $ni[$p]);
}
for ( $l = 0; $l < $jmp; $l++) {
   $p = $pjmp[$l];
   printf(" %4d  %s  %22.15e  +-  %21.15e\n",
   $l + $s + $v + $aas + $aac + $ass + $asc + 1, $pname[$p], $value[$p],$sigma[$p]);
   $ni[$p] = $l + $s + $v + $aas + $aac + $ass + $asc + 1;
#   printf("lpn %3d %3d %3d\n", $l , $p, $ni[$p]);
}
for ( $l = 0; $l < $off; $l++) {
   $p = $poff[$l];
   printf(" %4d  %s  %22.15e  +-  %21.15e\n",
   $l + $s + $v + $aas + $aac + $ass + $asc + $jmp +1, $pname[$p], $value[$p],$sigma[$p]);
   $ni[$p] = $l + $s + $v + $aas + $aac + $ass + $asc + $jmp + 1;
#   printf("lpn %3d %3d %3d\n", $l , $p, $ni[$p]);
}
################################################
for ( $l = 0; $l < $q; $l++){
   $ii = $ni[$id1[$l] - 1];
   $jj = $ni[$id2[$l] - 1];
   $k = $j*($jj - 1) + $ii;
   $corr2[$k] =  $corr[$l]; 
#   printf("%3d %3d   %3d %3d %3d  %10.4f\n", 
#      $id1[$l], $id2[$l], $ii, $jj, $k,  $corr[$l]);
}
$dim = $j;
for ( $ii = 1; $ii < $dim; $ii++) {
  for ( $jj = 0; $jj < $ii; $jj++) {
     $k = $dim*($jj ) + $ii + 1;
     if ( abs($corr2[$k]) > 1.0e-20){ 
       printf(" %4d  %4d  %22.14e\n", $ii + 1, $jj + 1, $corr2[$k]) ;
     }
  }
}
################################################
open(SVEC,"/tmp/svec.$$") || die "Can't open input file ",$_,"\n";
while (<SVEC>){
 print $_;
}
close(SVEC);
`rm /tmp/svec.$$`;
################################################
if ($prevref ne "no_frame") {
  print " ! Reference frame: $prevref\n";
}
################################################
if ($ant_ht_src) {
  print "$ant_ht_src\n";
}
################################################
exit;
################################################

