#======================================================================

#                    L A D C P P R O C . B E S T L A G 

#                    doc: Tue Sep 28 21:58:48 2010

#                    dlm: Fri Oct 19 11:51:04 2012

#                    (c) 2010 A.M. Thurnherr

#                    uE-Info: 26 46 NIL 0 0 72 2 2 4 NIL ofnI

#======================================================================

# TODO:

#	- first lag is always(?) nan unless CTD is turned on with LADCP in water

# HISTORY:

#	Sep 28, 2010: - created

#	Dec  9, 2010: - adapted to %CTD

#	Dec 10, 2010: - hardened bestlag failure test to require 1/3 agreeing lags

#	Jan  5, 2011: - changed first guess from 80% down to 10% down

#				  - added LADCP time lag to %PARAMs

#				  - added support of -i

#	Jul  7, 2011: - added code to remove window-mean of w before lagging to

#				    make it work in regions of crazy ocean w (IWISE 16007)

#	Jul 15, 2011: - changed screen-output of lag to take first guess lag into

#					account

#	Apr 11, 2011: - removed 1st guess lag consistency check based on large

#					elapsed offsets

#	May 18, 2012: - BUG: window start index was not always calculated correctly

#	Oct 19, 2012: - BUG: opt_i had wrong sign!

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++;

		return $LADCP{ENSEMBLE}[$ens]->{W};

	}

	return $LADCP{ENSEMBLE}[$ens-1]->{W} +

				$sc * ($LADCP{ENSEMBLE}[$ens]->{W} - $LADCP{ENSEMBLE}[$ens-1]->{W});

}

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($mCw,$mLw,$nw) = (0,0,0);									# first calc means

		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]);

			$mCw += $CTD{w}[$ws+$Ci];

			$mLw += $LADCP_w[$ws+$Li];

			$nw++;

		}

		next unless ($nw > 0);

		$mCw /= $nw; $mLw /= $nw;

		my($sad) = my($nad) = 0;										# calc mad with means removed

		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]-$mCw - ($LADCP_w[$ws+$Li]-$mLw));

			$nad++;

		}

		if ($sad/$nad < $bestmad) {

			$best = $Llag;

			$bestmad = $sad/$nad;

		}

	}

	return $best;	

}

sub lagLADCP2CTD()

{

	#------------------------------------------------------------------------

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

	#------------------------------------------------------------------------

	my($first_guess_lag);											# in units of CTD records

	if (defined($opt_i)) {

		$first_guess_lag = $opt_i / $CTD{sampint};

	} else {

		my($CTD_10pct_down) = 0;

		$CTD_10pct_down++											# "until" formulation allows for missing pressures

			until ($CTD{press}[$CTD_10pct_down]-$CTD{press}[0] >= 0.1*($CTD{maxpress}-$CTD{press}[0]));

		my($LADCP_10pct_down) = 0;

		$LADCP_10pct_down++

			while ($LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH} < 0.1*$LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});

		$first_guess_lag = ($LADCP{ENSEMBLE}[$LADCP_10pct_down]->{ELAPSED_TIME} -

							   $CTD_10pct_down*$CTD{sampint}) / $CTD{sampint};

		printf(STDERR "\n\t1st guess offset [CTD pressure, LADCP estimated depth] = %ds [%ddbar, %dm]\n",

				$first_guess_lag*$CTD{sampint},$CTD{press}[$CTD_10pct_down],$LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH})

					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

	#------------------------------------------------------------------------------------

	$nGaps = 0;

	for (my($ens)=$LADCP_start,my($r)=0; $ens<=$LADCP_end; $ens++) {

		while ($r*$CTD{sampint} < $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}) {

			unless ($first_guess_lag > $r) {

				$LADCP_w[$r-$first_guess_lag] = interp_LADCP_w($r*$CTD{sampint},$ens);

				$firstValid = $r-$first_guess_lag unless defined($firstValid);

				$nValid++;

			}

			$r++;

		}

	}

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

		if ($opt_d);

	print(STDERR "\n");

	#----------------------------------------------------------------------

	# Calculate lags

	#----------------------------------------------------------------------

	printf(STDERR "\tcalculating $opt_n lags from %ds-long windows [s]:",$opt_w);

	$opt_w = int($opt_w / $CTD{sampint});

	#---------------------------------------------------------------

	# carry out opt_n lag correlations and keep tally of the results

	#---------------------------------------------------------------

	my(%nBest);

	my($nLags) = 0;

	my($lags) = '';

	my($lastLag) = 9e99; my($nSame) = 1;

	for (my($window)=0; $window<$opt_n; $window++) {

		my($ws) = $firstValid + $window * int($nValid/$opt_n);	# window start

		$ws = @LADCP_w-$opt_w if ($ws+$opt_w >= @LADCP_w);

		$bestLag = bestLag($ws);

		if (defined($bestLag)) {

			if (defined($lastLag) && $bestLag == $lastLag) {

				$nSame++;

			} else {

				printf(STDERR "(x%d)",$nSame)

					if ($nSame > 1);

				printf(STDERR " %d",$bestLag*$CTD{sampint});

				$nSame = 1;

				$lastLag = $bestLag;

			}

			$lags .= sprintf(" %s",$bestLag*$CTD{sampint});

			$nBest{$bestLag}++;

			$nLags++;

		} else {

			if (!defined($lastLag)) {

				$nSame++;

			} else {

				printf(STDERR "(x%d)",$nSame)

					if ($nSame > 1);

				printf(STDERR " nan");

				$nSame = 1;

				$lastLag = $bestLag;

			}

		}

	}

	printf(STDERR "(x%d)",$nSame) if ($nSame > 1);

    &antsAddParams('time_lags',$lags);

	#----------------------

	# 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} > $opt_n/3);

	print(STDERR "\n\n\t\tWARNING: only $nBest{$best_lag} of the lag estimates agree!\n")

		if ($nBest{$best_lag} < $opt_n/2);

	if ($nBest{$best_lag} == $nLags) {

		printf(STDERR "\n\t\tunanimous lag = %ds\n",($first_guess_lag+$best_lag)*$CTD{sampint});

	} else {

		printf(STDERR "\n\t\tmost popular lag = %ds\n",($first_guess_lag+$best_lag)*$CTD{sampint});

	}

	&antsAddParams('LADCP_time_lag',($first_guess_lag + $best_lag) * $CTD{sampint});

	return ($first_guess_lag + $best_lag) * $CTD{sampint};

}

1;