LADCPproc.bestLag
changeset 1 54222c82435f
parent 0 de00d0f32431
child 3 711dd29cb6dd
--- a/LADCPproc.bestLag
+++ b/LADCPproc.bestLag
@@ -1,14 +1,22 @@
 #======================================================================
 #                    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: Tue Oct 26 13:37:19 2010
+#                    dlm: Wed Jan  5 19:44:09 2011
 #                    (c) 2010 A.M. Thurnherr
-#                    uE-Info: 58 94 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 67 28 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
+
 sub interp_LADCP_w($$)
 {
 	my($elapsed,$ens) = @_;
@@ -34,8 +42,8 @@
 		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]);
+			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);
@@ -50,29 +58,34 @@
 sub lagLADCP2CTD()
 {
 	#------------------------------------------------------------------------
-	# find 1st rec & ensemble >=80% down to max depth & make 1st guess at lag
+	# find 1st rec & ensemble >=10% down to max depth & make 1st guess at lag
 	#------------------------------------------------------------------------
 	
-	my($CTD_80pct_down) = 0;
-	$CTD_80pct_down++											# "until" formulation allows for missing pressures
-		until ($CTD_press[$CTD_80pct_down]-$CTD_press[0] >= 0.8*($CTD_maxpress-$CTD_press[0]));
-	
-	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);											# in units of CTD records
 	
-	my($first_guess_lag) = $LADCP{ENSEMBLE}[$LADCP_80pct_down]->{ELAPSED_TIME} -
-						   $CTD_80pct_down*$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_80pct_down],$LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH})
-				if ($opt_d);
-	croak("\n$0: cannot determine valid 1st guess offset ($LADCP_80pct_down,$CTD_80pct_down,$CTD_press[$CTD_80pct_down])\n")
-		if ($first_guess_lag*$CTD_sampint > 1200);
+	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);
+		croak("\n$0: cannot determine valid 1st guess offset ($LADCP_10pct_down,$CTD_10pct_down,$CTD{press}[$CTD_10pct_down])\n")
+	        if ($first_guess_lag*$CTD{sampint} > 1200);
+	}
 	
 	#------------------------------------------------------------------------------------
-	# Linearly interpolate LADCP time series onto a new grid with $CTD_sampint resolution
+	# 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
 	#------------------------------------------------------------------------------------
@@ -80,8 +93,8 @@
 	$nGaps = 0;
 	
 	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)
+		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++;
 		}
@@ -96,8 +109,8 @@
 	# Calculate lags
 	#----------------------------------------------------------------------
 
-	printf(STDERR "\tcalculating $opt_n lags from %ds-long windows [s]: ",$opt_w);
-	$opt_w = int($opt_w / $CTD_sampint);
+	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
@@ -105,20 +118,37 @@
 	my(%nBest);
 	my($nLags) = 0;
 	my($lags) = '';
+	my($lastLag) = 9e99; my($nSame) = 1;
 	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);
-			$lags .= sprintf(" %s",$bestLag*$CTD_sampint);
+			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 {
-			printf(STDERR "nan ");
+			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);
 	
 	#----------------------
@@ -129,17 +159,18 @@
 		$best_lag = $i if ($nBest{$i} > $nBest{$best_lag});
 	}
 	croak("\n$0: cannot determine a valid lag\n")
-		unless ($nBest{$best_lag} > 1);
+		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",$best_lag*$CTD_sampint);
+		printf(STDERR "\n\t\tunanimous lag = %ds\n",$best_lag*$CTD{sampint});
 	} else {
-		printf(STDERR "\n\t\tmost popular lag = %ds\n",$best_lag*$CTD_sampint);
+		printf(STDERR "\n\t\tmost popular lag = %ds\n",$best_lag*$CTD{sampint});
 	}
 
-	return ($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;