#!/usr/bin/perl
use Getopt::Std;
use Math::Trig;
# option -f for ETERNA output format added  2016-10-19 / HWz
# sign for southern hemisphere corrected    2018-01-16 / HWz

$blq=0;
# $$ COLUMN ORDER:  M2  S2  N2  K2  K1  O1  P1  Q1  MF  MM SSA
@phf  = ( 2, 2, 2, 2, 1, 1, 1, 1, 0, 0,  0 );
@indx = ( 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11 ); 

# wave group selection for ETERNA:
# @eternawaves = ("Q1", "O1", "K1", "P1", "N2", "M2", "S2", "K2", "MF", "MM", "SSA");
@tsoftwaves  = ("SSA","MM","MF","Q1","O1","P1","K1","N2","M2","S2","K2");
@eternawaves = ("Q1", "O1", "K1", "P1", "N2", "M2", "S2", "K2");
$amf=1.e-6;
$ampl=6;
$phal=7; $phald=9;
$endl=7;
$amfdm=1.e-5;
$show="GRAVITY [m/s^2]";
$tilt=0; $tiltamf=1.e-9;
$sign=1.;
$phaa=0.;
$st="";
$model="";

# equatorial band to reverse long period phases (based on GRS80)
$e2=0.00669438002290;        # 1st excentricity
$th=pi/2-acos(sqrt(1/3));   # equatorial band (geozentric)
$gam=atan(tan($th)/(1-$e2)); # above, but ellipsoidal

getopts('Raehnrsf');
if ( $opt_h )
{ 
    print "\nUSAGE:\n";
    print "          olgte.pl [options] blqfile\n\n";
    print "PURPOSE:\n";
    print "          Converts ocean loading gravity parameters\n";
    print "          from Automatic ocean loading provider format\n";
    print "          (a.k.a. BLQ) to TSOFT or, optionally, to ETERNA format.\n\n";
    print "          Only one file can be processed at one time,\n";
    print "          and the file must be headered, i.e. the complete\n";
    print "          email that http://holt.oso.chalmers.se/loading\n";
    print "          upon a user request must be input.\n\n";
    print "          This procedure copes with the following BLQ variants:\n";
    print "          One-component    (Gravity only)\n";
    print "          Three-components (Gravity and Tilt)\n";
    print "          nm/s**2 um/s**2 or mGal units\n\n";
    print "\n          Revised January 2018.\n";
    print "OPTIONS:\n";
    print " -e  -n   - Extract the tilt information, -e for east, -n for north.\n";
    print "            Default is gravity. -e and -n are mutually exclusive.\n\n";
    print " -r       - Reverse sign of amplitudes\n";
    print " -R       - Reverse sign by adding 180 degrees\n";
    print "            You can specify both -r -R but to no avail\n\n";
    print " -s       - Add the \"s\" before the tide symbol (old tradition)\n\n";
    print " -f       - Format ETERNA LOAD-Info in ini-file (ET34-Vxx) \n\n";
    print "OUTPUT:   Header record 1 contains the component and its units,\n";
    print "          the tide symbols and the name of the ocean model, \n";
    print "          taken from the BLQ input header.\n\n";
    print "          Header record 2 is issued per site segment of the\n";
    print "           BLQ file and contains\n";
    print "          station name, position, and a CompCode\n\n";
    print "          The CompCode gives you information about the ocean\n";
    print "          model, the date of the computation and whether\n";
    print "          the post-processor was used\n";
    print "          CompCode examples:\n";
    print "          HAMTIDE_PP = the Hamburg ocean model, post-processed\n";
    print "          SCHW_1S    = the Schwiderski model in one step\n\n";
    print "          Post processing is applied when stations have been\n";
    print "          located near a coast, else the parameters come from\n";
    print "          the 1-step solution\n\n";
    exit 
}
# PROCESS OPTIONS 
#
if ( $opt_e ){ $ampl=7; $phald=10; $show="EAST TILT [nrad]" }
if ( $opt_n ){ $ampl=8; $phald=11; $show="NORTH TILT [nrad]" }
if ( $opt_r ){ $sign=-1; }
if ( $opt_R ){ $phaa=180. }
if ( $opt_s ){ $st="s" }
#___________________________________________________________________________
#
# PROCESS FILE 
#
while (<>){
    chomp;
    if ( ! $blq ) # _____________________________________ PROCESS BLQ HEADER
    {
	if ( m|COLUMN ORDER:(.*)| )
	{
	    $tt=$1;
	    $tt=~s/\W(\S+)/$1 /g;
	    $tt=~s/MF/Mf/; 
	    $tt=~s/MM/Mm/;
	    $tt=~s/SSA/Ssa/;
	    $tt=~s/^\s+//; # remove leading blanks
	    $tt=~s/\s+$//; # remove trailing blanks
	    @t=( split(/ +/,$tt) )
	}
#
	if ( m|\$\$   TILT E| )
	{ 
	    $tilt=1; 
	    $phal=$phald; $endl=11; 
	    if ( $opt_e || $opt_n ){ $amf=$tiltamf; $amfdm=$tiltamf }
	}
#
	if ( m|\$\$ Ocean tide model: (.*)| )
	{ $model=$1 }
#
	if ( m|\$\$ (.*)_comp:| )
	{ $model=$1 }
#
	if ( m|\(mGal| )
	{ $amf=$amfdm; 
#          print "# BLQ Amplitudes were in mGal, using factor $amf\n" 
        }
#
	if ( m|\(nm/s| )
	{ $amf=1e-9; 
#          print "# BLQ Amplitudes were in nm/s^2, using factor $amf\n" 
        }
#
	if ( m|0\.1 urad| )
	{ $tiltamf = 1.e-6 }
#
	if ( m|END HEADER| )
	{ 
	    if ( ( $opt_n || $opt_e ) && ! $tilt )
	    { die "No tilt information in this file\n" }

	    $blq=1; $nr=-1; 
	    print "# Loading Tide in $show: @t from $model\n"; 
	}
    }
#
    if ( $blq )  #_________________________________________ CONVERT THE DATA 

    {   
	$nr++; 
        if ( $nr > $endl ) # _____________________________ PRINT THE SEGMENT
	{ 
	    print "# Station $s long.= $x lat.= $y height= $h  CompCode=$m\n";
            foreach $j ( @indx )
	    {
		$i=$j-1;  ## array starts at zero
		$name = uc $st.$t[$i];

		# normalize longitude
		$x+= 360. if ( $x < -180 );
		$x-= 360. if ( $x >  180 );

		# check phase
		$phab=0;
		# long
		$phab = 180 if ($phf[$i]==0 & abs($y) > $gam);
		# diurnal
		$phab = 180 if ($phf[$i]==1 & $y < 0);
		
		$pv = -$p[$i] - $phf[$i]*$x + $phaa - $phab;
		#printf "pv: $pv\n";
		
		# normalize phases
		$pv+= 360. if ( $pv < -180 );
		$pv-= 360. if ( $pv >  180 );

		if ( $opt_f )
		{
		    if ( grep { $t[$i] } @eternawaves )
		    {
			$model=~s/.*: //;
			$model=~s/,.*//;
			
			printf "LOADVECTOR=%10.3f%10.2f  %-8s  %-4s\n", 
			$a[$i]*$amf*$sign*1e9, $pv, $model, $name;
		    }
		}
		else
		{   
 		    if ( grep { $t[$i] } @tsoftwaves )
		    {
   		       printf "%-7s: %16.4e%15.2f\n",
		       $st.$t[$i], $a[$i]*$amf*$sign, $pv
		     }
		}
	    }
	    $nr=1; 
	    @a = ( );  @p = ( )
	}
#
	if ( m|lon/lat:\s+(\S+)\s+(\S+)\s+(\S+)| )  
	{ $x=$1; $y=$2; $h=$3 }
#
        if ( $nr == 3 )
	{
	    ## newer blq-files don't have descriptor 'Ocean tide model' 
	    ## in header
            s|(_PP)||; if ( $1 ){ $pp=$1 }else{$pp="_1S"}; 
#
	    m|^\$\$\s*(\S+).*ID:[ ]*(\S+) |; $model=$1.$pp; $m=$model."_on_".$2;
#
	    $model =~ s/_.*// if ( $model == "" );
#
	    $m =~ s|Complete(.*)|$1|;
	    $m =~ s|ID:||;
            $m =~ s|_PP|$pp|;
	}
#
	if ( $_ =~ /^ / )
	{   
	    
	    s/^\s+//; # remove leading blanks

	    if ( $nr == 2 )
	    { 
		m|^(.*)|; $s=$1
	    }
	    elsif ( $nr == $ampl )
	    {
		@a = ( @a, split(/ +/) );
	    }
	    elsif ( $nr == $phal )
	    {
		@p = ( @p, split(/ +/) );
	    }
	}
    }
}
exit