diff --git a/2011_07_10/LADCPproc.bestLag b/2011_07_10/LADCPproc.bestLag deleted file mode 100644 --- a/2011_07_10/LADCPproc.bestLag +++ /dev/null @@ -1,189 +0,0 @@ -#====================================================================== -# 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 Jul 8 02:54:29 2011 -# (c) 2010 A.M. Thurnherr -# uE-Info: 55 101 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) - -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=$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); - 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 - # 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}) { - $LADCP_w[$r-$first_guess_lag] = interp_LADCP_w($r*$CTD{sampint},$ens) - unless ($first_guess_lag > $r); - $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) = $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)) { - 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",$best_lag*$CTD{sampint}); - } else { - printf(STDERR "\n\t\tmost popular lag = %ds\n",$best_lag*$CTD{sampint}); - } - - &antsAddParams('LADCP_time_lag',($first_guess_lag + $best_lag) * $CTD{sampint}); - return ($first_guess_lag + $best_lag) * $CTD{sampint}; -} - -1;