#!/usr/bin/perl
#======================================================================
#                    L A D C P 2 C T D M A T C H 
#                    doc: Thu Sep 23 22:01:19 2010
#                    dlm: Tue Sep 28 21:43:29 2010
#                    (c) 2010 A.M. Thurnherr
#                    uE-Info: 312 33 NIL 0 0 70 10 2 4 NIL ofnI
#======================================================================

$antsSummary = 'match LADCP to CTD time series';

# HISTORY:
#	Sep 23, 2010: - incepted

($ANTS) = ($0 =~ m{^(.*)/[^/]*$});
$PERL_TOOLS = '/Data/LADCP/Software/Thurnherr/perl-tools';

require "$ANTS/ants.pl";
require "$PERL_TOOLS/RDI_BB_Read.pl";
require "$PERL_TOOLS/RDI_Coords.pl";
require "$PERL_TOOLS/RDI_Utils.pl";

$antsParseHeader = 0;
&antsUsage('dw:n:',2,
	'[-d)ebug] ',
	'[-w)indow <size[30s]>] [-n) <windows[20]]',
	'<RDI file> <SeaBird file>');

&antsCardOpt(\$opt_w,30);
&antsCardOpt(\$opt_n,20);

$LADCP_file = &antsFileArg();
$CTD_file 	= &antsFileArg();

&antsAddParams('LADCP_file',$LADCP_file,'CTD_file',$CTD_file);
@antsNewLayout = ('elapsed','CTD_w','LADCP_w');
&antsActivateOut();

#----------------------------------------------------------------------
# Step 1: Read Data
#----------------------------------------------------------------------

print(STDERR "Reading LADCP data ($LADCP_file)...");
readData($LADCP_file,\%LADCP);
printf(STDERR "\n\t%d ensembles\n",scalar(@{$LADCP{ENSEMBLE}}));

print(STDERR "Reading CTD data ($CTD_file)...");
open(F,$CTD_file) || croak("$CTD_file: $!\n");
while (1) {											# parse header
	chomp($hdr = <F>);
	$hdr =~ s/\r*$//;
	croak("$0: unexpected EOF (format error)\n") unless defined($hdr);
	last if ($hdr eq '*END*');
	
	$CTD_nfields = $',next if ($hdr =~ /nquan = /);		# Layout
	$CTD_nrecs = $',next if ($hdr =~ /nvalues = /);
	$pressF = $1,next if ($hdr =~ /name (\d+) = prDM:/);

	&antsAddParams('start_time',$1),next			# selected metadata
		if ($hdr =~ /start_time = (.*)/);

	&antsAddParams('station',$1),next
		if ($hdr =~ /Station\s*:\s*(.*)/);
	&antsAddParams('ship',$1),next
		if ($hdr =~ /Ship\s*:\s*(.*)\s*$/);
	&antsAddParams('cruise',$1),next
		if ($hdr =~ /Cruise\s*:\s*(.*)\s*$/);
	&antsAddParams('time',$1),next
		if ($hdr =~ /Time\s*:\s*(.*)/);
	&antsAddParams('date',$1),next
		if ($hdr =~ /Date\s*:\s*(.*)/);

	if ($hdr =~ /Latitude\s*[=:]\s*/) {
		($deg,$min,$NS) = split(/ /,$');
		$lat = $deg + $min/60;
		$lat *= -1 if ($NS eq 'S');
		&antsAddParams('lat',$lat);
		next;
	}
	if ($hdr =~ /Longitude\s*[=:]\s*/) {
		($deg,$min,$EW) = split(/ /,$');
		$lon = $deg + $min/60;
		$lon *= -1 if ($EW eq 'W');
		&antsAddParams('lon',$lon);
		next;
	}
	
	if ($hdr =~ /interval = seconds: /) {
		$CTD_sampint = 1*$';
		&antsAddParams('CTD_interval',1/$CTD_sampint);
		next;
	}
	
	$CTD_badval = $',next
		if ($hdr =~ /bad_flag = /);	
	$CTD_file_type = $',next
		if ($hdr =~ /file_type = /);	
}

croak("$CTD_file: cannot determine CTD file layout\n")
	unless ($CTD_nfields && $CTD_nrecs);
croak("$CTD_file: cannot determine missing value\n")
	unless defined($CTD_badval);
croak("$CTD_file: not a CTD time series file\n")
	unless ($CTD_sampint);
croak("$CTD_file: no pressure field\n")
	unless defined($pressF);

croak("$0: unknown latitude\n") unless defined($lat);

if ($CTD_file_type eq 'ascii') {
	while (1) {
		last unless (@rec = &antsFileIn(F));
		push(@CTD_press,($rec[$pressF] == $CTD_badval) ? nan : $rec[$pressF]);
	}
} elsif ($CTD_file_type eq 'binary') {

	my($fbits) = 8 * length(pack('f',0));
	croak(sprintf("$0: incompatible native CPU float representation (%d instead of 32bits)\n",fbits))
		unless ($fbits == 32);	
		
	croak("$CTD_file: can't read binary data\n")
		unless (read(F,$dta,4*$CTD_nfields*$CTD_nrecs) == 4*$CTD_nfields*$CTD_nrecs);
	print(STDERR "$CTD_file: WARNING: extraneous data at EOF\n") unless eof(F);

	$dta = pack('V*',unpack('N*',$dta))				# big-endian CPU
		if (unpack('h*', pack('s', 1)) =~ /01/);	# c.f. perlport(1)
	
	@dta = unpack("f*",$dta);

	for ($r=0; $r<$CTD_nrecs; $r++) {
		push(@CTD_press,($dta[$r*$CTD_nfields+$pressF] == $CTD_badval) ? nan : $dta[$r*$CTD_nfields+$pressF]);
	}
} else {
	croak("$CTD_file: unknown CTD file type $CTD_file_type\n");
}

printf(STDERR "\n\t%d scans\n",scalar(@CTD_press));

#----------------------------------------------------------------------
# Step 2a: Pre-Process CTD Data
#----------------------------------------------------------------------

printf(STDERR "Pre-processing data...");
printf(STDERR "\n\tCTD...");

#------------------------------------
# calculate w and find deepest record
#------------------------------------
$CTD_maxpress = -9e99;
for (my($r)=1; $r<$CTD_nrecs-1; $r++) {
	print(STDERR '.') if ($r%500==0);
	$CTD_w[$r] = 0.99*($CTD_press[$r+1] - $CTD_press[$r-1]) / (2*$CTD_sampint);
	if ($CTD_press[$r] > $CTD_maxpress) {
		$CTD_maxpress = $CTD_press[$r];
		$CTD_bottom = $r;
    }										
}

#----------------------------------------------------------------------
# Step 2b: Pre-Process LADCP Data
#----------------------------------------------------------------------

print(STDERR "\n\tLADCP...");

#------------------------------------------------------
# construct a depth-vs-time "profile" from the raw data
#------------------------------------------------------

my($LADCP_start,$LADCP_end,$LADCP_bottom,$w_gap_time,$zErr,$maxz) =
    mk_prof(\%LADCP,0,undef,1,6,70,0.1,120);
croak("\n$LADCP_file: no good ensembles found\n")
    unless defined($LADCP_start);
    
if ($opt_d) {
	printf(STDERR "\n\t\tStart of cast    : %s (#%5d) at %6.1fm\n",
	                    $LADCP{ENSEMBLE}[$LADCP_start]->{TIME},
	                    $LADCP{ENSEMBLE}[$LADCP_start]->{NUMBER},
	                    $LADCP{ENSEMBLE}[$LADCP_start]->{DEPTH});
	printf(STDERR "\t\tBottom of cast   : %s (#%5d) at %6.1fm\n",
	                    $LADCP{ENSEMBLE}[$LADCP_bottom]->{TIME},
	                    $LADCP{ENSEMBLE}[$LADCP_bottom]->{NUMBER},
	                    $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
	if (defined($water_depth)) {
	    printf(STDERR "\t\tSeabed           :                      at %6.1fm (+-%dm)\n",$water_depth,$sig_wd);
	} else {
	    print(STDERR "\t\tSeabed           : not found\n");
	}
	printf(STDERR "\t\tEnd of cast      : %s (#%5d) at %6.1fm\n",
	                    $LADCP{ENSEMBLE}[$LADCP_end]->{TIME},
	                    $LADCP{ENSEMBLE}[$LADCP_end]->{NUMBER},
	                    $LADCP{ENSEMBLE}[$LADCP_end]->{DEPTH});
}

#------------------------------------------------------------------------
# find 1st rec & ensemble >=80% down to max depth & make 1st guess at lag
#------------------------------------------------------------------------

my($CTD_80pct_down) = 0;
$CTD_80pct_down++
	while ($CTD_press[$CTD_80pct_down] < 0.8*$CTD_maxpress);

my($LADCP_80pct_down) = 0;
$LADCP_80pct_down++
	while ($LADCP{ENSEMBLE}[$LADCP_80pct_down]->{DEPTH} < 0.8*$LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});


my($first_guess_lag) = $LADCP{ENSEMBLE}[$LADCP_80pct_down]->{ELAPSED_TIME} -
					   $CTD_80pct_down*$CTD_sampint;

printf(STDERR "\t\t1st guess offset = %ds\n",$first_guess_lag*$CTD_sampint)
	if ($opt_d);

#------------------------------------------------------------------------------------
# Linearly interpolate LADCP time series onto a new grid with $CTD_sampint resolution
#	ALSO: apply first_guess_lag to make lags small, which keeps the bestlag data
#		  chunks large
#------------------------------------------------------------------------------------

my($nGaps) = 0;

sub interp_LADCP_w($$)
{
	my($elapsed,$ens) = @_;
	my($sc) = ($elapsed - $LADCP{ENSEMBLE}[$ens-1]->{ELAPSED_TIME}) /
			  ($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} -
					$LADCP{ENSEMBLE}[$ens-1]->{ELAPSED_TIME});
	unless (numberp($LADCP{ENSEMBLE}[$ens-1]->{W})) {
		$nGaps++;
		croak("$0: cannot (yet) handle double gaps --- request feature upgrade\n")
			unless numberp($LADCP{ENSEMBLE}[$ens]->{W});
		return $LADCP{ENSEMBLE}[$ens]->{W};
	}
	return $LADCP{ENSEMBLE}[$ens-1]->{W} +
				$sc * ($LADCP{ENSEMBLE}[$ens]->{W} - $LADCP{ENSEMBLE}[$ens-1]->{W});
}

for (my($ens)=$LADCP_start,my($r)=0; $ens<=$LADCP_end; $ens++) {
	while ($r*$CTD_sampint < $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}) {
		$LADCP_w[$r-$first_guess_lag] = interp_LADCP_w($r*$CTD_sampint,$ens)
			unless ($first_guess_lag > $r);
		$r++;
    }
}

print(STDERR "\t\t$nGaps gaps in w timeseries")
	if ($opt_d);

print(STDERR "\n");

#----------------------------------------------------------------------
# Step 3: Calculate Lags
#----------------------------------------------------------------------

printf(STDERR "Calculating $opt_n lags from %ds-long windows [s]: ",$opt_w);
$opt_w = int($opt_w / $CTD_sampint);

sub bestLag($)
{
	my($ws) = @_;														# window start index

	my($best);
	my($bestmad) = 9e99;												# mean absolute deviation
	for (my($Llag)=-int($opt_w/2); $Llag<int($opt_w/2); $Llag++) {
		my($sad) = my($nad) = 0;
		for (my($Ci)=0; $Ci<$opt_w; $Ci++) {
			my($Li) = $Ci + $Llag;
			next if ($Li<0 || $Li>=$opt_w);
			next unless numberp($CTD_w[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]);
			$sad += abs($CTD_w[$ws+$Ci] - $LADCP_w[$ws+$Li]);
			$nad++;
		}
		next unless ($nad > 0);
		if ($sad/$nad < $bestmad) {
			$best = $Llag;
			$bestmad = $sad/$nad;
		}
	}
	return $best;	
}

#---------------------------------------------------------------
# carry out opt_n lag correlations and keep tally of the results
#---------------------------------------------------------------
my(@nBest);
for (my($window)=0; $window<$opt_n; $window++) {
	my($ws) = $window * int(@LADCP_w/$opt_n);			# window start
	$ws = @LADCP_w-$opt_w if ($ws+$opt_w >= @LADCP_w);

	$bestLag = bestLag($ws);
	if (defined($bestLag)) {
		printf(STDERR "%d ",$bestLag*$CTD_sampint);
		$nBest{$bestLag}++;
	} else {
		printf(STDERR "nan ");
	}
}

#----------------------
# find most popular lag
#----------------------
my($best_lag);
foreach my $i (keys(%nBest)) {
	$best_lag = $i if ($nBest{$i} > $nBest{$best_lag});
}
croak("\n$0: cannot determine a valid lag\n")
	unless ($nBest{$best_lag} > 1);
print(STDERR "\n\n\tWARNING: only $nBest{$best_lag} of the lag estimates agree!\n")
	if ($nBest{$best_lag} < $opt_n/2);

if ($nBest{$best_lag} == $opt_n) {
	printf(STDERR "\n\tunanimous lag = %ds\n",$best_lag*$CTD_sampint);
} else {
	printf(STDERR "\n\tmost popular lag = %ds\n",$best_lag*$CTD_sampint);
}

#----------------------------------------------------------------------
# Step 4: associate CTD data with each LADCP ensemble & write time series
#----------------------------------------------------------------------

print(STDERR "Associating CTD data with LADCP ensembles...\n");

my($overall_dt) = $CTD_sampint * ($best_lag + $first_guess_lag);
for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
	my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} - $overall_dt) / $CTD_sampint);
	$LADCP{ENSEMBLE}[$ens]->{CTD_W} = $CTD_w[$r];
	&antsOut($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME},
			 $LADCP{ENSEMBLE}[$ens]->{CTD_W},
			 $LADCP{ENSEMBLE}[$ens]->{W});
}

&antsExit();