author | A.M. Thurnherr <athurnherr@yahoo.com> |
Wed, 16 Mar 2016 08:25:32 -0400 | |
changeset 40 | 408fc95bcff8 |
parent 36 | be205fe7ad76 |
child 41 | 6bddb82924e3 |
permissions | -rwxr-xr-x |
18 | 1 |
#!/usr/bin/perl |
2 |
#====================================================================== |
|
3 |
# L A D C P _ V K E |
|
4 |
# doc: Tue Oct 14 11:05:16 2014 |
|
36 | 5 |
# dlm: Tue Mar 8 15:55:52 2016 |
18 | 6 |
# (c) 2012 A.M. Thurnherr |
36 | 7 |
# uE-Info: 48 0 NIL 0 0 72 2 2 4 NIL ofnI |
18 | 8 |
#====================================================================== |
9 |
||
10 |
$antsSummary = 'calculate VKE from LADCP-derived vertical-velocity profiles'; |
|
11 |
||
12 |
# HISTORY: |
|
13 |
# Oct 14, 2014: - created from [LADCPfs] |
|
14 |
# Oct 15, 2014: - added parameterization output |
|
15 |
# Oct 17, 2014: - changed parameterization constant $c to 0.021 |
|
16 |
# Nov 6, 2014: - restored from backup and adapted to ANTS V6.0 |
|
17 |
# Nov 7, 2014: - changed parameterization constant $c to 0.0215 |
|
26
b89d4b01fcc5
after implementing new schemes for plotting and active files
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
18
diff
changeset
|
18 |
# Apr 16, 2015: - disabled output activation unless ANTS tools are available |
b89d4b01fcc5
after implementing new schemes for plotting and active files
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
18
diff
changeset
|
19 |
# - removed superfluous $ANTSBIN definition |
29 | 20 |
# May 18, 2015: - added -p)ulse <length> |
21 |
# Jun 11, 2015: - removed w_z code (yfname param requirement) |
|
22 |
# Jun 12, 2015: - adapted to &antsParam() |
|
23 |
# - BUG: %k.0 and %lambda.0 had been required erroneously (.1 are first ones to be used) |
|
24 |
# - made finescale limits optional |
|
25 |
# - renamed -b=>-l, -c=>-o, -s=>-i |
|
26 |
# Jun 14, 2015: - renamed -p=>-a |
|
27 |
# - added -p |
|
28 |
# - removed weird evals |
|
29 |
# Jun 15, 2015: - added plot & other mods |
|
30 |
# Jun 16, 2015: - define default outputs when STDOUT is a tty |
|
31 |
# - adapted to re-added pwrdens.0 in LADCP_wspec output |
|
32 |
# - modified plot label |
|
33 |
# Jun 17, 2015: - added eps.w to plot |
|
31 | 34 |
# Jul 30, 2015: - moved main label and made it blue |
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
35 |
# Sep 28, 2015: - added a p0fit r test to ensure that correlation has the correct sign |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
36 |
# Oct 12, 2015: - BUG: plot was in landscape mode |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
37 |
# - added input run label(s) to plot |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
38 |
# - require ANTSlibs V6.2 for release |
33 | 39 |
# Oct 13, 2015: - adapted to [version.pl] |
40 |
# Nov 30, 2015: - BUG: -a was not allowed in usage |
|
41 |
# - added -f |
|
42 |
# Dec 27, 2015: - reduced minlim for eps to 1e-13 W/kg |
|
34 | 43 |
# Dec 29, 2015: - added 3rd consistency check (p0 limit) |
44 |
# Jan 25, 2016: - added software version %PARAM |
|
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
45 |
# Mar 7, 2016: - removed unused spectral normalization code |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
46 |
# - added latitude constraint for calculation of eps.w |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
47 |
# Mar 8, 2016: - renamed fs_pwr to pwr.fs for consistency |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
48 |
# - renamed eps.w to eps.VKE for consistency |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
49 |
# - changed default output filenames for -d and -u |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
50 |
# - removed ./ from figure label |
36 | 51 |
# - update ANTSlib to 6.3 |
18 | 52 |
|
53 |
($ANTSLIB) = (`which ANTSlib` =~ m{^(.*)/[^/]*$}); |
|
33 | 54 |
($WCALC) = ($0 =~ m{^(.*)/[^/]*$}); |
55 |
$WCALC = '.' if ($WCALC eq ''); |
|
26
b89d4b01fcc5
after implementing new schemes for plotting and active files
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
18
diff
changeset
|
56 |
$ANTS_TOOLS_AVAILABLE = (`which list` ne ''); |
b89d4b01fcc5
after implementing new schemes for plotting and active files
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
18
diff
changeset
|
57 |
|
33 | 58 |
require "$WCALC/version.pl"; |
18 | 59 |
require "$ANTSLIB/ants.pl"; |
60 |
require "$ANTSLIB/libLADCP.pl"; |
|
29 | 61 |
require "$ANTSLIB/libGMT.pl"; |
34 | 62 |
&antsAddParams('LADCP_VKE',"Version $VERSION"); |
29 | 63 |
|
64 |
use FileHandle; |
|
65 |
use IPC::Open2; |
|
18 | 66 |
|
67 |
#---------------------------------------------------------------------- |
|
68 |
# Usage |
|
69 |
#---------------------------------------------------------------------- |
|
70 |
||
29 | 71 |
my($c) = 0.0215; # Thurnherr et al. (GRL 2015) |
72 |
||
33 | 73 |
&antsUsage('a:bc:de:f:g:i:l:mno:p:r:s:tuw:',0, |
29 | 74 |
'[poly-o)rder <n[0]> to de-mean data; -1 to disable>] [apply cosine-t)aper]', |
75 |
'[-d)own/-u)pcast-only] [exclude -b)ottom window]', # LADCP_wspec options |
|
76 |
'[-s)urface <layer depth to exclude[150m]>', |
|
77 |
'[-g)ap <max depth layer to fill with interpolation[40m]>]', |
|
78 |
'[-w)indow <power-of-two input-records[32]>]', |
|
79 |
'[shortwave -c)utoff <kz or lambda[100m]>]', # LADCP_VKE options |
|
80 |
'[o-m)it spectral correction] [spectral-tilt-correction -r)ange <max[0m]>]', |
|
81 |
'[-l) override ADCP bin <length>] [-a) override pulse <length>]', |
|
82 |
"[-e)ps-parameterization <scale[$c]>", |
|
33 | 83 |
'[output -f)iles in <dir>]', |
29 | 84 |
'[output -i)ndividual spectra <basename>]', |
85 |
'[output -p)lot <ps-file>', |
|
86 |
'[file]'); |
|
87 |
||
88 |
#---------------------------------------------------------------------- |
|
89 |
# Calculate VKE spectra with [LADCP_wspec] if input is a w_ocean file |
|
90 |
#---------------------------------------------------------------------- |
|
18 | 91 |
|
29 | 92 |
if (defined(fnrNoErr('dc_w'))) { # pre-process with LADCP_wspec when handed vertical-velocity input |
93 |
my($opts); |
|
94 |
$opts .= ' -d' if ($opt_d); |
|
95 |
$opts .= ' -u' if ($opt_u); |
|
96 |
$opts .= ' -b' if ($opt_b); |
|
97 |
$opts .= ' -t' if ($opt_t); |
|
98 |
$opts .= " -s $opt_s" if defined($opt_s); |
|
99 |
$opts .= " -g $opt_g" if defined($opt_g); |
|
100 |
$opts .= " -w $opt_w" if defined($opt_w); |
|
101 |
$opts .= " -o $opt_o" if defined($opt_o); |
|
102 |
open2(\*FROMCLD,\*TOCLD,"LADCP_wspec $opts") || # spawn sub-process |
|
103 |
croak("LADCP_wspec $opts: $!\n"); |
|
104 |
open(STDIN,"<&FROMCLD") || croak("dup(FROMCLD): $!\n"); |
|
105 |
close(FROMCLD); |
|
106 |
print(TOCLD $antsOldHeaders); undef($antsOldHeaders); # feed already gobbled header & first record to child |
|
107 |
print(TOCLD $antsPeekBuffer); undef($antsPeekBuffer); |
|
108 |
undef(%P); # shouldn't matter, because we'll get the same %PARAMs back |
|
109 |
undef(@antsLayout); # shouldn't matter, because it will get overwritten |
|
110 |
while (<>) { print(TOCLD $_); } # feed remaining data |
|
111 |
close(TOCLD); |
|
112 |
} elsif (defined(fnrNoErr('pwrdens.1'))) { |
|
113 |
croak("$0: -d, -u, -b, -w, -s meaningless when $0 used with spectral input\n") |
|
114 |
if ($opt_d || $opt_u || $opt_b || defined($opt_w) || |
|
115 |
defined($opt_s) || defined($opt_g)); |
|
116 |
} else { |
|
117 |
if ($ARGV[0]) { |
|
118 |
croak("$ARGV[0]: no such file or directory\n"); |
|
119 |
} else { |
|
120 |
croak("$0: empyt input\n"); |
|
121 |
} |
|
122 |
} |
|
18 | 123 |
|
29 | 124 |
#---------------------------------------------------------------------- |
125 |
# Handle LADCP_VKE usage & read data |
|
126 |
#---------------------------------------------------------------------- |
|
127 |
||
128 |
&antsAddParams('ADCP_bin_length',$opt_l) # override bin length |
|
129 |
if defined($opt_l); |
|
130 |
&antsAddParams('ADCP_pulse_length',$opt_a) # override pulse length |
|
131 |
if defined($opt_a); |
|
18 | 132 |
|
133 |
&antsFloatOpt(\$opt_e,$c); # default parameterization |
|
134 |
||
29 | 135 |
if (defined($opt_c)) { # shortwave cutoff supplied |
136 |
$lmin = ($opt_c < 1) ? 2*$PI/$opt_c : $opt_c; |
|
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
137 |
&antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin); # ensure eps.VKE is calculated below |
29 | 138 |
} elsif (defined(antsParam('shortwave_cutoff'))) { # cutoff already applied |
139 |
$lmin = 2*$PI/antsParam('shortwave_cutoff'); |
|
140 |
} else { # use 100m default cutoff |
|
141 |
$lmin = 100; |
|
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
142 |
&antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin); # ensure eps.VKE is calculated below |
29 | 143 |
} |
144 |
$lmax = 9e99; # no longwave cutoff implemented yet |
|
145 |
||
146 |
&antsInstallBufFull(0); # load entire file |
|
147 |
&antsIn(); |
|
148 |
my($Hbuf) = $antsOldHeaders; # save for later |
|
149 |
||
150 |
&antsRequireParam('profile_id'); |
|
151 |
&antsRequireParam('lambda.1'); |
|
152 |
&antsRequireParam('k.1'); |
|
153 |
&antsRequireParam('resolution_bandwidth'); |
|
154 |
&antsRequireParam('input_depth_resolution'); |
|
155 |
&antsRequireParam('output_depth_resolution'); |
|
156 |
&antsRequireParam('ADCP_bin_length'); |
|
157 |
&antsAddParams('ADCP_pulse_length',antsParam('ADCP_bin_length')) |
|
158 |
unless defined(antsParam('ADCP_pulse_length')); |
|
159 |
||
18 | 160 |
$imin = 0; # find frequency bin limits |
29 | 161 |
for ($nfreq=1; defined(antsParam("lambda.$nfreq")); $nfreq++) { |
162 |
$imin = $nfreq if ($imin==0 && antsParam("lambda.$nfreq")<=$lmax); |
|
163 |
$imax = $nfreq if (antsParam("lambda.$nfreq") >= $lmin); |
|
18 | 164 |
} |
165 |
croak("$0: <lambda.min=$lmin> < min(lambda)") |
|
166 |
unless defined($imax); |
|
167 |
||
168 |
$pg_fmin = fnr('pwrdens.0'); # first power field in spectra |
|
169 |
$fs_fmin = $pg_fmin + $imin; # first power field in finescale range |
|
170 |
$fs_fmax = $pg_fmin + $imax; # last power field in finescale range |
|
171 |
||
29 | 172 |
$widxf = fnr('widx'); |
18 | 173 |
$mindf = fnr('depth.min'); |
174 |
$maxdf = fnr('depth.max'); |
|
175 |
||
176 |
#---------------------------------------------------------------------- |
|
29 | 177 |
# Redirect STDOUT & create plot if STDOUT is a tty |
18 | 178 |
#---------------------------------------------------------------------- |
179 |
||
29 | 180 |
if (-t STDOUT) { |
33 | 181 |
$opt_f = '.' unless defined($opt_f); # output directory |
182 |
croak("$opt_f: not a directory\n") unless (-d $opt_f); |
|
183 |
||
29 | 184 |
my($id) = &antsRequireParam('profile_id'); |
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
185 |
|
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
186 |
if ($opt_d) { |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
187 |
$opt_p = sprintf('%s/%03d_dc_VKE.ps',$opt_f,$id) unless defined($opt_p); |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
188 |
$outfile = sprintf('%s/%03d_dc.VKE',$opt_f,$id); |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
189 |
} elsif ($opt_u) { |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
190 |
$opt_p = sprintf('%s/%03d_uc_VKE.ps',$opt_f,$id) unless defined($opt_p); |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
191 |
$outfile = sprintf('%s/%03d_uc.VKE',$opt_f,$id); |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
192 |
} else { |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
193 |
$opt_p = sprintf('%s/%03d_VKE.ps',$opt_f,$id) unless defined($opt_p); |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
194 |
$outfile = sprintf('%s/%03d.VKE',$opt_f,$id); |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
195 |
} |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
196 |
$outfile =~ s@^./@@; |
29 | 197 |
open(STDOUT,">$outfile") || die("$outfile: $!\n"); |
33 | 198 |
} elsif (defined($opt_f)) { |
199 |
croak("-f can only be used without STDOUT redirection\n"); |
|
29 | 200 |
} |
201 |
||
202 |
#---------------------------------------------------------------------- |
|
203 |
# Library |
|
204 |
#---------------------------------------------------------------------- |
|
205 |
||
206 |
sub integrate_fs_power($) # integrate fs spectrum |
|
18 | 207 |
{ |
208 |
my($r) = @_; |
|
209 |
||
210 |
$ants_[$r][$fspwrf] = 0; |
|
211 |
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) { |
|
29 | 212 |
$ants_[$r][$fspwrf] += $ants_[$r][$f]; |
18 | 213 |
} |
29 | 214 |
$ants_[$r][$fspwrf] *= antsParam('resolution_bandwidth'); |
18 | 215 |
} |
216 |
||
217 |
||
29 | 218 |
sub fit_universal_w_spec($) # vertical velocity => p0 |
18 | 219 |
{ |
220 |
my($r) = @_; |
|
221 |
my($nsamp) = $fs_fmax - $fs_fmin + 1; |
|
222 |
||
223 |
#--------------------------------------------------- |
|
224 |
# fit slope-2 line in log-log space (main estimator) |
|
225 |
#--------------------------------------------------- |
|
226 |
||
29 | 227 |
if ($nsamp >= 2) { # require min 2 wavenumber samples |
18 | 228 |
|
229 |
my($sumd,$sumx,$sumy) = (0,0,0); |
|
29 | 230 |
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) { # loop over wavenumbers |
18 | 231 |
my($i) = $f - $pg_fmin; |
29 | 232 |
$sumd += log10($ants_[$r][$f]) + 2*log10(antsRequireParam("k.$i")); |
233 |
$sumx += log10(antsParam("k.$i")); |
|
18 | 234 |
$sumy += log10($ants_[$r][$f]); |
235 |
} |
|
236 |
my($p0) = $sumd/$nsamp; |
|
237 |
$ants_[$r][$p0f] = 10**$p0; |
|
238 |
||
29 | 239 |
my($avgx) = $sumx/$nsamp; # avg for r calc |
18 | 240 |
my($avgy) = $sumy/$nsamp; |
241 |
||
29 | 242 |
my($sumsq,$sxx,$syy,$sxy) = (0,0,0,0); # calc rms misfit & correlation coeff |
18 | 243 |
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) { |
244 |
my($i) = $f - $pg_fmin; |
|
29 | 245 |
my($xt) = log10(antsParam("k.$i")) - $avgx; |
18 | 246 |
my($yt) = log10($ants_[$r][$f]) - $avgy; |
29 | 247 |
$sumsq += ($p0 - 2*log10(antsParam("k.$i")) - log10($ants_[$r][$f]))**2; |
18 | 248 |
$sxx += $xt * $xt; $syy += $yt * $yt; $sxy += $xt * $yt; |
249 |
} |
|
250 |
$ants_[$r][$rmsf] = sqrt($sumsq/$nsamp); |
|
251 |
$ants_[$r][$rf] = $sxy/(sqrt($sxx * $syy) + $SMALL_AMOUNT); |
|
252 |
||
253 |
} else { |
|
254 |
&antsInfo("WARNING: no fit --- need min 2 samples"); |
|
255 |
} |
|
256 |
$ants_[$r][$sampf] = $nsamp; |
|
257 |
} |
|
258 |
||
259 |
#---------------------------------------------------------------------- |
|
260 |
# Process File |
|
261 |
#---------------------------------------------------------------------- |
|
262 |
||
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
263 |
$fspwrf = &antsNewField('pwr.fs'); # derived fields |
18 | 264 |
|
29 | 265 |
$p0f = &antsNewField('p0'); # VKE parameterization |
266 |
$rmsf = &antsNewField('p0fit.rms'); |
|
267 |
$sampf = &antsNewField('p0fit.nsamp'); |
|
268 |
$rf = &antsNewField('p0fit.r'); |
|
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
269 |
$wepsf = &antsNewField('eps.VKE'); |
18 | 270 |
|
271 |
my(@outLayout) = @antsNewLayout; # save for later |
|
272 |
for ($f=0; $f<@outLayout; $f++) { # determine last spectral field in input |
|
273 |
$totf = $f if ($outLayout[$f] eq 'pwr.tot'); |
|
274 |
$tnsf = $f if ($outLayout[$f] eq 'pwr.tot.nsamp'); |
|
275 |
} |
|
276 |
croak("$0: cannot find fields 'pwr.tot' or 'pwr.tot.nsamp' in input\n") |
|
277 |
unless defined($totf) || defined($tnsf); |
|
29 | 278 |
$lsf = defined($tnsf) ? $tnsf : $totf; |
18 | 279 |
|
29 | 280 |
if ($opt_p) { # begin plot |
281 |
$xmin = (&antsParam('output_depth_resolution')>=350) ? 0.012 : 0.018; |
|
282 |
$xmax = 0.2; $ymin = 1; $ymax = 1e4; |
|
283 |
$plotsize = 13; |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
284 |
GMT_begin($opt_p,"-JX${plotsize}l","-R$xmin/$xmax/$ymin/$ymax",'-P -X6 -Y4'); # init plot |
29 | 285 |
GMT_psxy('-G25/255/25 -L'); |
286 |
printf(GMT "%f %f\n",$xmin,$xmin**(-2)*sqrt(2)); |
|
287 |
printf(GMT "%f %f\n",$xmax,$xmax**(-2)*sqrt(2)); |
|
288 |
printf(GMT "%f %f\n",$xmax,$xmax**(-2)/sqrt(2)); |
|
289 |
printf(GMT "%f %f\n",$xmin,$xmin**(-2)/sqrt(2)); |
|
290 |
} |
|
291 |
||
292 |
my(@sumPwr,@nPwr); |
|
293 |
my(@sumGoodPwr,@nGoodPwr); |
|
294 |
my($min_depth) = 9e99; |
|
295 |
my($max_depth) = -9e99; |
|
18 | 296 |
|
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
297 |
my($latM) = abs(&antsRequireParam('lat')); |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
298 |
&antsInfo("WARNING: low latitude-profile no epsilon estimated") |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
299 |
unless ($latM > 5); |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
300 |
|
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
301 |
for (my($r)=0; $r<@ants_; $r++) { # loop over all windows |
18 | 302 |
|
303 |
#-------------------------- |
|
304 |
# apply spectral correction |
|
305 |
#-------------------------- |
|
306 |
||
29 | 307 |
unless ($opt_m) { |
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
308 |
for (my($i)=0; $i<$nfreq; $i++) { # loop over wavenumbers |
29 | 309 |
$ants_[$r][$i+$pg_fmin] *= |
310 |
T_w(antsParam("k.$i"),antsParam('ADCP_bin_length'), |
|
311 |
antsParam('ADCP_pulse_length'),antsParam('input_depth_resolution'), |
|
312 |
$opt_r); |
|
18 | 313 |
} |
314 |
} |
|
315 |
||
316 |
#------------------------ |
|
317 |
# calculate fs quantities |
|
318 |
#------------------------ |
|
319 |
||
29 | 320 |
integrate_fs_power($r); # calculate total finescale power |
321 |
fit_universal_w_spec($r); # fit kz^-2 spectrum |
|
35
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
322 |
|
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
323 |
if ($latM > 5 && # Tests: 0) not equatorial |
54b8bb450e5f
Version 1.2 finished. Out for testing to Jay Hooper.
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
34
diff
changeset
|
324 |
$ants_[$r][$rmsf] <= 0.4 && # 1) rms from kz fit < 0.4 (as in paper) |
33 | 325 |
$ants_[$r][$rf] < 0 && # 2) raw spectra are red spectrum |
326 |
$ants_[$r][$p0f] > 1e-7) { # 3) exclude low-p0 tail apparent at eps<2e-10 W/kg |
|
327 |
$ants_[$r][$wepsf] = ($ants_[$r][$p0f] / $opt_e)**2; |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
328 |
} else { |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
329 |
$ants_[$r][$wepsf] = nan; |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
330 |
} |
29 | 331 |
|
332 |
if (numberp($ants_[$r][$p0f])) { # update min/max depth |
|
333 |
$min_depth = $ants_[$r][$mindf] if ($ants_[$r][$mindf] < $min_depth); |
|
334 |
$max_depth = $ants_[$r][$maxdf] if ($ants_[$r][$maxdf] > $max_depth); |
|
335 |
} |
|
18 | 336 |
|
337 |
#--------------- |
|
338 |
# produce output |
|
339 |
#--------------- |
|
340 |
||
29 | 341 |
if (defined($opt_p)) { # plot current spectrum on -p & calc mean |
342 |
my($R) = 0; # RGB map |
|
343 |
my($G) = int(200*(1-$r/@ants_)); |
|
344 |
my($B) = ($r < @ants_/2) ? 150 : int(100+100*(1-$r/@ants_)); |
|
345 |
GMT_psxy("-W8,$R/$G/$B"); |
|
346 |
for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) { # avg & plot high-quality spectra only |
|
347 |
next unless numberp($ants_[$r][$f]) && numberp($ants_[$r][$p0f]); # - omit zero wavenumber ($pg_fmin) |
|
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
348 |
next unless numberp($ants_[$r][$wepsf]); |
29 | 349 |
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin)); |
350 |
printf(GMT "$k %g\n",$ants_[$r][$f]/$ants_[$r][$p0f]); |
|
351 |
$nGoodPwr[$f]++; |
|
352 |
$sumGoodPwr[$f] += $ants_[$r][$f]/$ants_[$r][$p0f]; |
|
353 |
} |
|
354 |
GMT_psxy("-W4,$R/$G/$B,."); # avg & plot all spectra with dots |
|
355 |
for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) { |
|
356 |
next unless numberp($ants_[$r][$f]) && numberp($ants_[$r][$p0f]); |
|
357 |
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin)); |
|
358 |
printf(GMT "$k %g\n",$ants_[$r][$f]/$ants_[$r][$p0f]); |
|
359 |
$nPwr[$f]++; |
|
360 |
$sumPwr[$f] += $ants_[$r][$f]/$ants_[$r][$p0f]; |
|
361 |
} |
|
362 |
} |
|
363 |
||
364 |
||
365 |
if (defined($opt_i)) { # output current spectrum on -i |
|
18 | 366 |
open(STDOUT_DUP,">&",STDOUT) || croak("$0: cannot dup STDOUT\n"); |
367 |
@antsNewLayout = ('k','lambda','pwrdens','finescale','pwrdens_fit'); |
|
29 | 368 |
if ($ants_[$r][$widxf] > 0) { |
369 |
$ofname = sprintf("${opt_i}_%02d.wspec",$r+1); |
|
370 |
} else { |
|
371 |
$ofname = sprintf("${opt_i}_btm.wspec"); |
|
372 |
} |
|
18 | 373 |
open(STDOUT,">$ofname") || croak("$ofname: $!\n"); |
374 |
undef($antsOldHeaders); |
|
26
b89d4b01fcc5
after implementing new schemes for plotting and active files
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
18
diff
changeset
|
375 |
&antsActivateOut() if ($ANTS_TOOLS_AVAILABLE); |
18 | 376 |
|
377 |
my($saveParams) = $antsCurParams; # add all extra input fields as %PARAMs |
|
378 |
&antsAddParams('binpgrams::depth.min',$ants_[$r][$mindf], |
|
379 |
'binpgrams::depth.max',$ants_[$r][$maxdf]); |
|
380 |
for (my($f)=$lsf+1; $f<@outLayout; $f++) { |
|
381 |
&antsAddParams($outLayout[$f],$ants_[$r][$f]); |
|
382 |
} |
|
383 |
||
29 | 384 |
for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) { |
385 |
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin)); |
|
386 |
my($l) = antsParam(sprintf("lambda.%d",$f-$pg_fmin)); |
|
18 | 387 |
&antsOut($k,$l,$ants_[$r][$f], |
388 |
($f>=$fs_fmin && $f<=$fs_fmax), |
|
389 |
numberp($ants_[$r][$p0f]) ? $ants_[$r][$p0f] * $k**(-2) : nan); |
|
390 |
} |
|
391 |
||
392 |
&antsOut('EOF'); |
|
393 |
open(STDOUT,">&",STDOUT_DUP) || croak("$0: cannot restore STDOUT\n"); |
|
394 |
close(STDOUT_DUP); |
|
395 |
$antsCurParams = $saveParams; |
|
396 |
} |
|
397 |
} |
|
398 |
||
29 | 399 |
if (defined($opt_p)) { # complete plot |
400 |
if (@nPwr) { |
|
401 |
GMT_psxy('-W10,255/100/100,.'); # plot mean spectrum; dotted; entire range |
|
402 |
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) { |
|
403 |
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin)); |
|
404 |
printf(GMT "$k %g\n",$sumPwr[$f]/$nPwr[$f]); |
|
405 |
} |
|
406 |
} |
|
407 |
||
408 |
if (@nGoodPwr) { |
|
409 |
GMT_psxy('-W4/255/100/100'); # plot mean spectrum; entire range |
|
410 |
for (my($f)=$pg_fmin+1; $f<$pg_fmin+$nfreq; $f++) { |
|
411 |
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin)); |
|
412 |
printf(GMT "$k %g\n",$sumGoodPwr[$f]/$nGoodPwr[$f]); |
|
413 |
} |
|
414 |
GMT_psxy('-W16/255/100/100'); # plot mean fs spectrum (heavy) |
|
415 |
for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) { |
|
416 |
my($k) = antsParam(sprintf("k.%d",$f-$pg_fmin)); |
|
417 |
printf(GMT "$k %g\n",$sumGoodPwr[$f]/$nGoodPwr[$f]); |
|
418 |
} |
|
419 |
} |
|
420 |
||
421 |
GMT_psbasemap('-Bf3a2:"Vertical Wavenumber [rad/m]":/' . # annotate axes |
|
422 |
'f3a1:"VKE/p@-0@- [m@+3@+s@+-1/2@+/(rad m@+-1@+)]":WeS'); |
|
423 |
GMT_setJ("-JX-${plotsize}l"); |
|
424 |
GMT_setR(sprintf('-R%f/%f/%f/%f',2*$PI/$xmax,2*$PI/$xmin,$ymin,$ymax)); |
|
425 |
GMT_psbasemap('-B3:"Vertical Wavelength [m]:N"'); |
|
426 |
||
30 | 427 |
GMT_unitcoords_logscale(); # print profile number |
31 | 428 |
GMT_pstext('-Gblue'); |
32
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
429 |
if (defined($outfile)) { printf(GMT "0.98 0.02 14 0 0 BR $outfile [$P{run_label}] %s\n",antsParam('input_data')); } |
6041a20feb39
V1.0 release candidate
A.M. Thurnherr <athurnherr@yahoo.com>
parents:
31
diff
changeset
|
430 |
else { printf(GMT "0.98 0.02 14 0 0 BR %03d [$P{run_label}] %s\n",antsParam('profile_id'),antsParam('input_data')); } |
29 | 431 |
|
432 |
GMT_set('ANNOT_FONT_SIZE_PRIMARY 10','ANNOT_OFFSET_PRIMARY 0.01c', # eps profile inset |
|
433 |
'LABEL_FONT_SIZE 10','LABEL_OFFSET 0'); |
|
434 |
$min_depth = 400 if ($min_depth > 1e4); |
|
435 |
$max_depth = 1500 if ($max_depth < 0); |
|
33 | 436 |
GMT_setR(sprintf("-R5e-13/1e-7/%d/%d",round($min_depth-250,500),round($max_depth+250,500))); |
29 | 437 |
GMT_setJ(sprintf('-JX%fl/-%f',$plotsize/2.2,$plotsize/2.2)); |
438 |
GMT_psbasemap('-X0.3 -Y0.3 -Bg1a1f3p:"@~e@~@-VKE@- [W kg@+-1@+]":/g1000a500f100:"Depth [m]":wEsN'); |
|
439 |
GMT_psxy('-M'); |
|
440 |
for (my($r)=0; $r<@ants_; $r++) { |
|
441 |
next unless numberp($ants_[$r][$wepsf]); |
|
442 |
my($R) = 0; my($G) = int(200*(1-$r/@ants_)); |
|
443 |
my($B) = ($r < @ants_/2) ? 150 : int(100+100*(1-$r/@ants_)); |
|
444 |
print(GMT "> -W12,$R/$G/$B\n"); |
|
445 |
print(GMT "$ants_[$r][$wepsf] $ants_[$r][$mindf]\n"); |
|
446 |
print(GMT "$ants_[$r][$wepsf] $ants_[$r][$maxdf]\n"); |
|
447 |
} |
|
448 |
||
449 |
GMT_end(); # finish plot |
|
450 |
} |
|
451 |
||
18 | 452 |
@antsNewLayout = @outLayout; |
453 |
$antsOldHeaders = $Hbuf; |
|
454 |
$antsHeadersPrinted = 0; |
|
455 |
||
456 |
&antsFlush(); # output record with results |
|
457 |
&antsExit(); |