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