|
1 #====================================================================== |
|
2 # L I B L A D C P . P L |
|
3 # doc: Wed Jun 1 20:38:19 2011 |
|
4 # dlm: Mon Apr 1 16:13:46 2019 |
|
5 # (c) 2011 A.M. Thurnherr |
|
6 # uE-Info: 27 65 NIL 0 0 70 2 2 4 NIL ofnI |
|
7 #====================================================================== |
|
8 |
|
9 # HISTORY: |
|
10 # Jun 1, 2011: - created |
|
11 # Jul 29, 2011: - improved |
|
12 # Aug 10, 2011: - made correct |
|
13 # Aug 11, 2011: - added "convenient combinations" |
|
14 # Aug 18, 2011: - made buoyancy frequency non-constant in S() |
|
15 # Jan 4, 2012: - improved T_VI to allow correcting even without superensembles |
|
16 # Jan 5, 2012: - removed S(), which is just pwrdens/N^2 (rather than |
|
17 # pwrdens/N^2/(2pi) as I erroneously thought) |
|
18 # Jan 18, 2012: - added T_VI_alt() to allow assessment of tilt correction extrema |
|
19 # Aug 22, 2012: - added documentation |
|
20 # - added T_w() |
|
21 # Sep 24, 2012: - made "k" argument default in T_w() |
|
22 # Oct 25, 2012: - renamed T_SM() to T_ASM() |
|
23 # Jun 26. 2013: - added T_w_z() |
|
24 # - added parameter checks to processing-specific corrections |
|
25 # May 18, 2015: - added pulse length to T_w() and T_w_z() |
|
26 # Apr 25, 2018: - added eps_VKE() parameterization |
|
27 # Apr 1, 2019: - removed plen parameter from T_w() and T_w_z() |
|
28 |
|
29 require "$ANTS/libvec.pl"; |
|
30 require "$ANTS/libfuns.pl"; |
|
31 |
|
32 #------------------------------------------------------------------------------ |
|
33 # VKE parameterization for epsilon |
|
34 # |
|
35 # NOTES: |
|
36 # - see Thurnherr et al. (GRL 2015) |
|
37 # - calculate eps from p0 |
|
38 # - optional second argument allows free choice of parameterization constant |
|
39 # - default value is from paper, which is slightly lower than the current value |
|
40 # used in [LADCP_VKE], which applies the parameterization only to spectra |
|
41 # passing a few tests |
|
42 #------------------------------------------------------------------------------ |
|
43 |
|
44 sub eps_VKE(@) |
|
45 { |
|
46 my($p0,$c) = |
|
47 &antsFunUsage(-1,'.','<p0[m^2/s^2/(rad/m)]> [c[0.021 [1/sqrt(s)]]]',@_); |
|
48 $c = 0.021 unless defined($c); |
|
49 return numberp($p0) ? ($p0 / $c) ** 2 : nan; |
|
50 } |
|
51 |
|
52 #------------------------------------------------------------------------------ |
|
53 # Spectral corrections for LADCP data |
|
54 # |
|
55 # NOTES: |
|
56 # - see Polzin et al. (JAOT 2002), Thurnherr (JAOT 2012) |
|
57 # - to correct, multiply power densities (or power, I think) with corrections |
|
58 # - apply to down-/up-cast data only |
|
59 #------------------------------------------------------------------------------ |
|
60 #---------------------------------------------------------------------- |
|
61 # 1. Corrections for individual data acquisition and processing steps |
|
62 #---------------------------------------------------------------------- |
|
63 |
|
64 #------------------------------------------------------------------------------ |
|
65 # T_ravg(k,blen[,plen]) |
|
66 # - correct for range averaging due to finite pulse and finite receive window |
|
67 # - when called with 2 arguments, bin-length == pulse-length assumed |
|
68 #------------------------------------------------------------------------------ |
|
69 |
|
70 sub T_ravg(@) |
|
71 { |
|
72 my($k,$blen,$plen) = |
|
73 &antsFunUsage(-2,'ff','<vertical wavenumber[rad/s]> <bin-length[m]> [pulse-length[m]]',@_); |
|
74 $plen = $blen unless defined($plen); |
|
75 return 1 / sinc($k*$blen/2/$PI)**2 / sinc($k*$plen/2/$PI)**2; |
|
76 } |
|
77 |
|
78 #------------------------------------------------------------- |
|
79 # T_fdiff(k,dz) |
|
80 # - correct for first differencing on a grid with dz spacing |
|
81 #------------------------------------------------------------- |
|
82 |
|
83 sub T_fdiff($$) |
|
84 { |
|
85 my($k,$dz) = |
|
86 &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <differencing interval[m]>',@_); |
|
87 return 1 / sinc($k*$dz/2/$PI)**2; |
|
88 } |
|
89 |
|
90 #------------------------------------------------------------ |
|
91 # T_interp(k,blen,dz) |
|
92 # - correct for CODAS gridding-with-interpolation algorithm |
|
93 # - ONLY USED IN UH SOFTWARE |
|
94 #------------------------------------------------------------ |
|
95 |
|
96 sub T_interp($$$) |
|
97 { |
|
98 my($k,$blen,$dz) = |
|
99 &antsFunUsage(3,'fff','<vertical wavenumber[rad/s]> <bin length[m]> <grid resolution[m]>',@_); |
|
100 return 1 / sinc($k*$blen/2/$PI)**4 / sinc($k*$dz/2/$PI)**2; |
|
101 } |
|
102 |
|
103 #------------------------------------------------------------------------- |
|
104 # T_binavg(k,dz) |
|
105 # - correct for simple bin averaging |
|
106 # - Polzin et al. suggest to use blen instead of dz; this must be a typo |
|
107 #------------------------------------------------------------------------- |
|
108 |
|
109 sub T_binavg($$) |
|
110 { |
|
111 my($k,$dz) = |
|
112 &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <grid resolution[m]>',@_); |
|
113 return 1 / sinc($k*$dz/2/$PI)**2; |
|
114 } |
|
115 |
|
116 #-------------------------------------------------------------------------------- |
|
117 # T_tilt(k,d') |
|
118 # - d' is a length scale that depends on tilt stats and range max |
|
119 # - on d' == 0, T_tilt() == 1, i.e. the correction is disabled |
|
120 # - d' = dprime(range_max) |
|
121 # - is from a quadratic fit to three data points given by Polzin et al. |
|
122 # - see Thurnherr (J. Tech. 2012) for notes |
|
123 # - on range_max == 0, d' == 0, i.e. the correction is disabled |
|
124 #-------------------------------------------------------------------------------- |
|
125 |
|
126 sub T_tilt($$) |
|
127 { |
|
128 my($k,$dprime) = |
|
129 &antsFunUsage(2,'ff','<vertical wavenumber[rad/s]> <d-prime[m]>',@_); |
|
130 return $dprime ? 1 / sinc($k*$dprime/2/$PI)**2 : 1; |
|
131 } |
|
132 |
|
133 sub dprime($) |
|
134 { |
|
135 return $_[0] ? -1.2 + 0.0857 * $_[0] - 0.000136 * $_[0]**2 : 0; |
|
136 } |
|
137 |
|
138 #====================================================================== |
|
139 # 2. Processing-Specific Corrections |
|
140 #====================================================================== |
|
141 |
|
142 #---------------------------------------------------------------------- |
|
143 # T_UH(k,blen,dz,range_max) |
|
144 # - UH implementation of the shear method (WOCE standard) |
|
145 # - range_max == 0 disables tilt correction |
|
146 #---------------------------------------------------------------------- |
|
147 |
|
148 sub T_UH($$$$) |
|
149 { |
|
150 my($k,$blen,$dz,$range_max) = |
|
151 &antsFunUsage(4,'ffff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>',@_); |
|
152 croak("T_UH($k,$blen,$dz,$range_max): bad parameters\n") |
|
153 unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0); |
|
154 return T_ravg($k,$blen) * T_fdiff($k,$blen) * T_interp($k,$blen,$dz) * T_tilt($k,dprime($range_max)); |
|
155 } |
|
156 |
|
157 #---------------------------------------------------------------------- |
|
158 # T_ASM(k,blen,dz,range_max) |
|
159 # - re-implemented shear method with simple depth binning |
|
160 # - range_max == 0 disables tilt correction |
|
161 #---------------------------------------------------------------------- |
|
162 |
|
163 sub T_ASM($$$$) |
|
164 { |
|
165 my($k,$blen,$dz,$range_max) = |
|
166 &antsFunUsage(4,'ffff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <shear grid resolution[m]> <range max[m]>',@_); |
|
167 croak("T_ASM($k,$blen,$dz,$range_max): bad parameters\n") |
|
168 unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0); |
|
169 return T_ravg($k,$blen) * T_fdiff($k,$blen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max)); |
|
170 } |
|
171 |
|
172 #------------------------------------------------------------ |
|
173 # T_VI(k,blen,preavg_dz,grid_dz,range_max) |
|
174 # T_VI_alt(k,blen,preavg_dz,grid_dz,dprime) |
|
175 # - velocity inversion method of Visbeck (J. Tech., 2002) |
|
176 # - only valid if pre-averaging into superensembles is used |
|
177 # - range_max == 0 disables tilt correction |
|
178 # - sel == nan disables pre-averaging correction |
|
179 #------------------------------------------------------------ |
|
180 |
|
181 sub T_VI($$$$$) |
|
182 { |
|
183 my($k,$blen,$sel,$dz,$range_max) = |
|
184 &antsFunUsage(5,'ff.ff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <superensemble size[m]|nan> <shear grid resolution[m]> <range max[m]>',@_); |
|
185 return T_VI_alt($k,$blen,$sel,$dz,dprime($range_max)); |
|
186 } |
|
187 |
|
188 sub T_VI_alt($$$$$) |
|
189 { |
|
190 my($k,$blen,$sel,$dz,$dprime) = |
|
191 &antsFunUsage(5,'ff.ff','<vertical wavenumber[rad/s]> <ADCP bin size[m]> <superensemble size[m]|nan> <shear grid resolution[m]> <tilt d-prime[m]>',@_); |
|
192 croak("T_VI_alt($k,$blen,$sel,$dz,$range_max): bad parameters\n") |
|
193 unless ($k>=0 && $blen>0 && $sel ne '' && $dz>0 && $range_max>=0); |
|
194 croak("$0: tilt-dprime outside range [0..$blen]\n") |
|
195 unless ($dprime>=0 && $dprime<=$blen); |
|
196 return ($sel>0) ? T_ravg($k,$blen) * T_binavg($k,$sel) * T_binavg($k,$dz) * T_tilt($k,$dprime) |
|
197 : T_ravg($k,$blen) * T_binavg($k,$dz) * T_tilt($k,$dprime); |
|
198 } |
|
199 |
|
200 #---------------------------------------------------------------------- |
|
201 # T_w(k,blen,dz,range_max) |
|
202 # - vertical-velocity method of Thurnherr (IEEE 2011) |
|
203 # - range_max == 0 disables tilt correction |
|
204 #---------------------------------------------------------------------- |
|
205 |
|
206 { my(@fc); |
|
207 sub T_w(@) |
|
208 { |
|
209 my($k,$blen,$dz,$range_max) = |
|
210 &antsFunUsage(4,'ffff', |
|
211 '[vertical wavenumber[rad/s]] <ADCP bin length[m]> <output grid resolution[m]> <range max[m]>', |
|
212 \@fc,'k',undef,undef,undef,undef,@_); |
|
213 croak("T_w($k,$blen,$dz,$range_max): bad parameters\n") |
|
214 unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0); |
|
215 return T_ravg($k,$blen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max)); |
|
216 } |
|
217 } |
|
218 |
|
219 #---------------------------------------------------------------------- |
|
220 # T_w_z(k,blen,dz,range_max) |
|
221 # - vertical-velocity method of Thurnherr (IEEE 2011) |
|
222 # - first differencing of gridded shear to calculate dw/dz |
|
223 # - NB: grid-scale differentiation assumed |
|
224 # - range_max == 0 disables tilt correction |
|
225 #---------------------------------------------------------------------- |
|
226 |
|
227 { my(@fc); |
|
228 sub T_w_z(@) |
|
229 { |
|
230 my($k,$blen,$dz,$range_max) = |
|
231 &antsFunUsage(4,'ffff', |
|
232 '[vertical wavenumber[rad/s]] <ADCP bin size[m]> <output grid resolution[m]> <range max[m]>', |
|
233 \@fc,'k',undef,undef,undef,undef,@_); |
|
234 croak("T_w_z($k,$blen,$dz,$range_max): bad parameters\n") |
|
235 unless ($k>=0 && $blen>0 && $dz>0 && $range_max>=0); |
|
236 return T_ravg($k,$blen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max)) * T_fdiff($k,$dz); |
|
237 } |
|
238 } |
|
239 |
|
240 1; |