67 # Jan 25, 2016: - search-radius-doubling heuristic had typo |
67 # Jan 25, 2016: - search-radius-doubling heuristic had typo |
68 # - added %PARAMs |
68 # - added %PARAMs |
69 # Feb 19, 2016: - added support for -l |
69 # Feb 19, 2016: - added support for -l |
70 # - added warning |
70 # - added warning |
71 # Mar 7, 2016: - BUG: editing did not work correctly in all cases |
71 # Mar 7, 2016: - BUG: editing did not work correctly in all cases |
|
72 # Mar 6, 2017: - BUG: assertion in mad_w failed with 2017 P18 DL#206 |
|
73 # Mar 9, 2017: - tightened timelag editing (good_diff: 2->1) |
72 |
74 |
73 # DIFFICULT STATIONS: |
75 # DIFFICULT STATIONS: |
74 # NBP0901#131 this requires the search-radius doubling heuristic |
76 # NBP0901#131 this requires the search-radius doubling heuristic |
75 |
77 |
76 # TODO: |
78 # TODO: |
77 # - better seabed code (from LADCPproc) |
79 # - better seabed code (from LADCPproc) |
78 |
80 |
79 my($TINY) = 1e-6; |
81 my($TINY) = 1e-6; |
80 |
82 |
81 sub mad_w($$$) # mean absolute deviation |
83 sub mad_w($$$) # mean absolute deviation |
82 { |
84 { |
83 my($fe,$le,$so) = @_; # first/last LADCP ens, CTD scan offset |
85 my($fe,$le,$so) = @_; # first/last LADCP ens, CTD scan offset |
84 my($sad) = my($n) = 0; |
|
85 |
86 |
86 my($LADCP_mean_w,$CTD_mean_w,$nsamp) = (0,0,0); |
87 my($LADCP_mean_w,$CTD_mean_w,$nsamp) = (0,0,0); |
87 for (my($e)=$fe; $e<=$le; $e++) { # first, calculate mean w in window |
88 for (my($e)=$fe; $e<=$le; $e++) { # first, calculate mean w in window |
88 my($s) = int(($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[0]) / $CTD{DT} + 0.5); |
89 my($s) = int(($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[0]) / $CTD{DT} + 0.5); |
89 next unless ($s>=0 && $s<=$#{$CTD{ELAPSED}}); |
90 |
|
91 # THE FOLLOWING LINE CAUSES AN ASSERTION FAILURE WITH 2017 P08 DL#206. I AM NOT SURE WHETHER MY |
|
92 # FIX SOLVES THE UNDERLYING PROBLEM OR ONLY THIS SPECIAL CASE. |
|
93 # next unless ($s>=0 && $s<=$#{$CTD{ELAPSED}}); |
|
94 |
|
95 next unless ($s>0 && $s<=$#{$CTD{ELAPSED}}); |
90 die("assertion failed\n" . |
96 die("assertion failed\n" . |
91 "\ttest: abs($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$s]) <= $CTD{DT}/2\n" . |
97 "\ttest: abs($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$s]) <= $CTD{DT}/2\n" . |
92 "\te = $e, s = $s, ensemble = $LADCP{ENSEMBLE}[$e]->{NUMBER}" |
98 "\te = $e, s = $s, ensemble = $LADCP{ENSEMBLE}[$e]->{NUMBER}" |
93 ) unless (abs($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$s]) <= $CTD{DT}/2+$TINY); |
99 ) unless (abs($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[$s]) <= $CTD{DT}/2+$TINY); |
94 next unless numberp($LADCP{ENSEMBLE}[$e]->{REFLR_W}); |
100 next unless numberp($LADCP{ENSEMBLE}[$e]->{REFLR_W}); |
95 my($dw) = $LADCP{ENSEMBLE}[$e]->{REFLR_W}-$LADCP_mean_w - ($CTD{W}[$s+$so]-$CTD_mean_w); |
|
96 next unless (abs($dw) <= $w_max_lim); |
|
97 |
|
98 $LADCP_mean_w += $LADCP{ENSEMBLE}[$e]->{REFLR_W}; |
101 $LADCP_mean_w += $LADCP{ENSEMBLE}[$e]->{REFLR_W}; |
99 $CTD_mean_w += $CTD{W}[$s+$so]; |
102 $CTD_mean_w += $CTD{W}[$s+$so]; |
100 $nsamp++; |
103 $nsamp++; |
101 } |
104 } |
102 return 9e99 unless ($nsamp); |
105 return 9e99 unless ($nsamp); |
103 $LADCP_mean_w /= $nsamp; |
106 $LADCP_mean_w /= $nsamp; |
104 $CTD_mean_w /= $nsamp; |
107 $CTD_mean_w /= $nsamp; |
105 |
108 |
106 for (my($e)=$fe; $e<=$le; $e++) { # now, calculate mad |
109 my($sad) = 0; # now, calculate mad |
|
110 for (my($e)=$fe; $e<=$le; $e++) { |
|
111 next unless numberp($LADCP{ENSEMBLE}[$e]->{REFLR_W}); |
107 my($s) = int(($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[0]) / $CTD{DT} + 0.5); |
112 my($s) = int(($LADCP{ENSEMBLE}[$e]->{ELAPSED} + $CTD{TIME_LAG} - $CTD{ELAPSED}[0]) / $CTD{DT} + 0.5); |
|
113 next unless ($s>=0 && $s<=$#{$CTD{ELAPSED}}); |
|
114 next unless numberp($LADCP{ENSEMBLE}[$e]->{REFLR_W}); |
108 my($dw) = $LADCP{ENSEMBLE}[$e]->{REFLR_W}-$LADCP_mean_w - ($CTD{W}[$s+$so]-$CTD_mean_w); |
115 my($dw) = $LADCP{ENSEMBLE}[$e]->{REFLR_W}-$LADCP_mean_w - ($CTD{W}[$s+$so]-$CTD_mean_w); |
109 next unless numberp($LADCP{ENSEMBLE}[$e]->{REFLR_W}); |
116 # print(STDERR "dw = $dw ($LADCP{ENSEMBLE}[$e]->{REFLR_W}-$LADCP_mean_w - ($CTD{W}[$s+$so]-$CTD_mean_w)\n"); |
110 next unless (abs($dw) <= $w_max_lim); |
117 next unless (abs($dw) <= $w_max_lim); |
111 $sad += abs($dw); |
118 $sad += abs($dw); |
112 $n++; |
119 } |
113 } |
120 return $sad/$nsamp; |
114 return ($n>0) ? $sad/$n : 9e99; # n == 0, e.g. in bottom gap |
|
115 } |
121 } |
116 |
122 |
117 |
123 |
118 sub bestLag($$$$) # find best lag in window |
124 sub bestLag($$$$) # find best lag in window |
119 { |
125 { |
120 my($fe,$le,$ww,$soi) = @_; # first/last LADCP ens, window width, scan-offset increment |
126 my($fe,$le,$ww,$soi) = @_; # first/last LADCP ens, window width, scan-offset increment |
121 my($bestso) = 0; # error at first-guess offset |
127 my($bestso) = 0; # error at first-guess offset |
122 my($bestmad) = mad_w($fe,$le,0); |
128 my($bestmad) = mad_w($fe,$le,0); |
123 |
129 |
|
130 # print(STDERR "bestLag($fe,$le,$ww,$soi)\n"); |
124 for (my($dso) = 1; $dso <= int($ww/2/$CTD{DT} + 0.5); $dso+=$soi) { |
131 for (my($dso) = 1; $dso <= int($ww/2/$CTD{DT} + 0.5); $dso+=$soi) { |
125 my($mad) = mad_w($fe,$le,-$dso); |
132 my($mad) = mad_w($fe,$le,-$dso); |
|
133 # print(STDERR "-$dso $mad\n"); |
126 $bestmad=$mad,$bestso=-$dso if ($mad < $bestmad); |
134 $bestmad=$mad,$bestso=-$dso if ($mad < $bestmad); |
127 $mad = mad_w($fe,$le,$dso); |
135 $mad = mad_w($fe,$le,$dso); |
|
136 # print(STDERR " $dso $mad\n"); |
128 $bestmad=$mad,$bestso=$dso if ($mad < $bestmad); |
137 $bestmad=$mad,$bestso=$dso if ($mad < $bestmad); |
129 } |
138 } |
|
139 # print(STDERR "-> $bestso $bestmad\n"); |
130 return ($bestso,$bestmad); |
140 return ($bestso,$bestmad); |
131 } |
141 } |
132 |
142 |
133 #---------------------------------------------------------------------- |
143 #---------------------------------------------------------------------- |
134 # carry out lag correlations and keep tally of the results |
144 # carry out lag correlations and keep tally of the results |