1 #====================================================================== |
1 #====================================================================== |
2 # L A D C P P R O C . B E S T L A G |
2 # L A D C P P R O C . B E S T L A G |
3 # doc: Tue Sep 28 21:58:48 2010 |
3 # doc: Tue Sep 28 21:58:48 2010 |
4 # dlm: Wed Jan 5 19:44:09 2011 |
4 # dlm: Fri Jul 8 02:54:29 2011 |
5 # (c) 2010 A.M. Thurnherr |
5 # (c) 2010 A.M. Thurnherr |
6 # uE-Info: 67 28 NIL 0 0 72 2 2 4 NIL ofnI |
6 # uE-Info: 55 101 NIL 0 0 72 2 2 4 NIL ofnI |
7 #====================================================================== |
7 #====================================================================== |
8 |
8 |
9 # TODO: |
9 # TODO: |
10 # - first lag is always(?) nan unless CTD is turned on with LADCP in water |
10 # - first lag is always(?) nan unless CTD is turned on with LADCP in water |
11 |
11 |
14 # Dec 9, 2010: - adapted to %CTD |
14 # Dec 9, 2010: - adapted to %CTD |
15 # Dec 10, 2010: - hardened bestlag failure test to require 1/3 agreeing lags |
15 # Dec 10, 2010: - hardened bestlag failure test to require 1/3 agreeing lags |
16 # Jan 5, 2011: - changed first guess from 80% down to 10% down |
16 # Jan 5, 2011: - changed first guess from 80% down to 10% down |
17 # - added LADCP time lag to %PARAMs |
17 # - added LADCP time lag to %PARAMs |
18 # - added support of -i |
18 # - added support of -i |
|
19 # Jul 7, 2011: - added code to remove window-mean of w before lagging to |
|
20 # make it work in regions of crazy ocean w (IWISE 16007) |
19 |
21 |
20 sub interp_LADCP_w($$) |
22 sub interp_LADCP_w($$) |
21 { |
23 { |
22 my($elapsed,$ens) = @_; |
24 my($elapsed,$ens) = @_; |
23 my($sc) = ($elapsed - $LADCP{ENSEMBLE}[$ens-1]->{ELAPSED_TIME}) / |
25 my($sc) = ($elapsed - $LADCP{ENSEMBLE}[$ens-1]->{ELAPSED_TIME}) / |
36 my($ws) = @_; # window start index |
38 my($ws) = @_; # window start index |
37 |
39 |
38 my($best); |
40 my($best); |
39 my($bestmad) = 9e99; # mean absolute deviation |
41 my($bestmad) = 9e99; # mean absolute deviation |
40 for (my($Llag)=-int($opt_w/2); $Llag<int($opt_w/2); $Llag++) { |
42 for (my($Llag)=-int($opt_w/2); $Llag<int($opt_w/2); $Llag++) { |
41 my($sad) = my($nad) = 0; |
43 my($mCw,$mLw,$nw) = (0,0,0); # first calc means |
42 for (my($Ci)=0; $Ci<$opt_w; $Ci++) { |
44 for (my($Ci)=0; $Ci<$opt_w; $Ci++) { |
43 my($Li) = $Ci + $Llag; |
45 my($Li) = $Ci + $Llag; |
44 next if ($Li<0 || $Li>=$opt_w); |
46 next if ($Li<0 || $Li>=$opt_w); |
45 next unless numberp($CTD{w}[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]); |
47 next unless numberp($CTD{w}[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]); |
46 $sad += abs($CTD{w}[$ws+$Ci] - $LADCP_w[$ws+$Li]); |
48 $mCw += $CTD{w}[$ws+$Ci]; |
|
49 $mLw += $LADCP_w[$ws+$Li]; |
|
50 $nw++; |
|
51 } |
|
52 next unless ($nw > 0); |
|
53 $mCw /= $nw; $mLw /= $nw; |
|
54 |
|
55 my($sad) = my($nad) = 0; # calc mad with means removed |
|
56 for (my($Ci)=0; $Ci<$opt_w; $Ci++) { |
|
57 my($Li) = $Ci + $Llag; |
|
58 next if ($Li<0 || $Li>=$opt_w); |
|
59 next unless numberp($CTD{w}[$ws+$Ci]) && numberp($LADCP_w[$ws+$Li]); |
|
60 $sad += abs($CTD{w}[$ws+$Ci]-$mCw - ($LADCP_w[$ws+$Li]-$mLw)); |
47 $nad++; |
61 $nad++; |
48 } |
62 } |
49 next unless ($nad > 0); |
|
50 if ($sad/$nad < $bestmad) { |
63 if ($sad/$nad < $bestmad) { |
51 $best = $Llag; |
64 $best = $Llag; |
52 $bestmad = $sad/$nad; |
65 $bestmad = $sad/$nad; |
53 } |
66 } |
54 } |
67 } |