--- a/LADCPintsh
+++ b/LADCPintsh
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P I N T S H
# doc: Thu Oct 14 21:22:50 2010
-# dlm: Fri Dec 10 04:57:50 2010
+# dlm: Wed Feb 16 15:13:46 2011
# (c) 2010 A.M. Thurnherr & E. Firing
-# uE-Info: 207 1 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 401 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'integrate LADCP shear';
@@ -28,6 +28,7 @@
# Oct 24, 2010: - fix spuriously small variances that can occur for BT velocities based on very small
# samples (i.e. primarily when chosing a small -r)
# Dec 9, 2010: - allowed for empty BT file
+# Feb 16, 2011: - BUG: gaps in shear data were not handled correctly in baroclinic solution
($ANTS) = (`which list` =~ m{^(.*)/[^/]*$});
require "$ANTS/ants.pl";
@@ -378,16 +379,26 @@
unless (defined($refU)) {
my($sumU,$sumV,$sumW,$dc_sumU,$dc_sumV,$dc_sumW,$uc_sumU,$uc_sumV,$uc_sumW);
+ my($nSumVel,$dc_nSumVel,$uc_nSumVel);
for (my($r)=0; $r<@depth; $r++) {
- $sumU += $u[$r]; $sumV += $v[$r]; $sumW += $w[$r];
- $dc_sumU += $dc_u[$r]; $dc_sumV += $dc_v[$r]; $dc_sumW += $dc_w[$r];
- $uc_sumU += $uc_u[$r]; $uc_sumV += $uc_v[$r]; $uc_sumW += $uc_w[$r];
+ if (numberp($u[$r])) {
+ $nSumVel++;
+ $sumU += $u[$r]; $sumV += $v[$r]; $sumW += $w[$r];
+ }
+ if (numberp($dc_u[$r])) {
+ $dc_nSumVel++;
+ $dc_sumU += $dc_u[$r]; $dc_sumV += $dc_v[$r]; $dc_sumW += $dc_w[$r];
+ }
+ if (numberp($uc_u[$r])) {
+ $uc_nSumVel++;
+ $uc_sumU += $uc_u[$r]; $uc_sumV += $uc_v[$r]; $uc_sumW += $uc_w[$r];
+ }
}
- $refU = $sumU / @depth; $refV = $sumV / @depth; $refW = $sumW / @depth;
- $dc_refU = $dc_sumU / @depth; $dc_refV = $dc_sumV / @depth; $dc_refW = $dc_sumW / @depth;
- $uc_refU = $uc_sumU / @depth; $uc_refV = $uc_sumV / @depth; $uc_refW = $uc_sumW / @depth;
+ $refU = $sumU / $nSumVel; $refV = $sumV / $nSumVel; $refW = $sumW / $nSumVel;
+ $dc_refU = $dc_sumU / $dc_nSumVel; $dc_refV = $dc_sumV / $dc_nSumVel; $dc_refW = $dc_sumW / $dc_nSumVel;
+ $uc_refU = $uc_sumU / $uc_nSumVel; $uc_refV = $uc_sumV / $uc_nSumVel; $uc_refW = $uc_sumW / $uc_nSumVel;
}
for (my($r)=0; $r<@depth; $r++) { # reference velocities
--- a/LADCPproc
+++ b/LADCPproc
@@ -2,9 +2,9 @@
#======================================================================
# L A D C P P R O C
# doc: Thu Sep 16 20:36:10 2010
-# dlm: Mon Jan 10 13:53:54 2011
+# dlm: Wed Jun 15 15:47:25 2011
# (c) 2010 A.M. Thurnherr & E. Firing
-# uE-Info: 382 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 479 55 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
$antsSummary = 'process LADCP data to get shear, time series';
@@ -32,8 +32,11 @@
# Dec 10, 2010: - change -w) default to 120s
# - changed nshear output to 0 from nan when there are no samples
# Dec 27, 2010: - changed sign of -l to accept lag output from [LADCP_w]
-# Jan 10, 2010: - -o => -k added new -o
+# Jan 10, 2011: - -o => -k added new -o
# - added code to fill CTD sound vel gaps
+# Jan 22, 2011: - added -g) lat,lon
+# - added -c)ompass corr
+# Jun 15, 2011: - added mean backscatter profile to default output
($ANTS) = (`which list` =~ m{^(.*)/[^/]*$});
($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
@@ -52,17 +55,18 @@
require "$PERL_TOOLS/RDI_Utils.pl";
$antsParseHeader = 0;
-&antsUsage('24ab:df:i:kl:n:o:ps:t:w:',2,
+&antsUsage('24ab:c:df:g:i:kl:n:o:ps:t:w:',2,
'[use -2)dary CTD sensor pair]',
'[require -4)-beam LADCP solutions]',
- '[-s)etup <file>]',
+ '[-s)etup <file>] [-g)ps <lat,lon>]',
+ '[-c)ompass-corr <offset,cos-fac,sin-fac>]',
'[enable -p)PI editing]',
'[-o)utput grid <resolution[5m]>]',
'[-i)nitial LADCP time lag <guestimate>]',
'[-l)ag LADCP <by>] [auto-lag -w)indow <size[120s]>] [-n) <auto-lag windows[20]]',
'[-d)iagnostic output]',
'output: [-t)ime series <file>] [-f)lag <file>] [-b)ottom-track <file>]',
- ' [-a)coustic backscatter profiles] [bottom-trac-k) profs]',
+ ' [per-bin -a)coustic backscatter profiles] [bottom-trac-k) profs]',
'<RDI file> <SeaBird file>');
$RDI_Coords::minValidVels = 4 if ($opt_4);
@@ -78,6 +82,19 @@
&antsFloatOpt($opt_i);
&antsCardOpt($opt_o);
+if (defined($opt_g)) {
+ ($CTD{lat},$CTD{lon}) = split(',',$opt_g);
+ croak("$0: cannot decode -g $opt_g\n")
+ unless numberp($CTD{lat}) && numberp($CTD{lon});
+}
+
+if (defined($opt_c)) {
+ ($CC_offset,$CC_cos_fac,$CC_sin_fac) = split(',',$opt_c);
+ croak("$0: cannot decode -c $opt_c\n")
+ unless numberp($CC_offset) && numberp($CC_cos_fac) && numberp($CC_sin_fac);
+}
+
+
$LADCP_file = &antsFileArg();
$CTD_file = &antsFileArg();
@@ -209,10 +226,18 @@
&antsAddParams('3_beam_solutions',"$RDI_Coords::threeBeam_1 $RDI_Coords::threeBeam_2 $RDI_Coords::threeBeam_3 $RDI_Coords::threeBeam_4");
}
} elsif ($LADCP{EARTH_COORDINATES}) {
- printf(STDERR "\n\t\tcorrecting for magnetic declination of %.1f deg...\n",$magdec)
- if ($opt_d);
+ if ($opt_d) {
+ if ($opt_c) {
+ printf(STDERR "\n\t\tcalculating tilt and correcting for compass error and magnetic declination of %.1f deg...\n",$magdec);
+ } else {
+ printf(STDERR "\n\t\tcalculating tilt and correcting for magnetic declination of %.1f deg...\n",$magdec);
+ }
+ }
for (my($ens)=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
$LADCP{ENSEMBLE}[$ens]->{TILT} = &angle_from_vertical($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
+ my($hdg) = rad($LADCP{ENSEMBLE}[$ens]->{HEADING});
+ $LADCP{HEADING_BIAS} = ($CC_offset + $CC_cos_fac*cos($hdg) + $CC_sin_fac*sin($hdg)) - $magdec
+ if ($opt_c);
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]} =
velApplyHdgBias(\%LADCP,$ens,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]});
@@ -421,7 +446,7 @@
@antsNewLayout = ('depth','dc_nshear','dc_u_z','dc_u_z.sig','dc_v_z','dc_v_z.sig','dc_w_z','dc_w_z.sig',
'uc_nshear','uc_u_z','uc_u_z.sig','uc_v_z','uc_v_z.sig','uc_w_z','uc_w_z.sig',
- 'nshear','u_z','u_z.sig','v_z','v_z.sig','w_z','w_z.sig');
+ 'nshear','u_z','u_z.sig','v_z','v_z.sig','w_z','w_z.sig','Sv','Sv.n');
$commonParams = $antsCurParams;
&antsAddParams('ubin_start',$ubin_start,'ubin_end',$ubin_end, # record processing params
@@ -434,7 +459,8 @@
'max_shdev',$max_shdev,'max_shdev_sum',$max_shdev_sum,
'water_depth',round($water_depth),'water_depth.sig',round($sig_water_depth),
'min_hab',round($min_hab),
- 'clip_margin',$clip_margin,'first_clip_bin',$first_clip_bin);
+ 'clip_margin',$clip_margin,'first_clip_bin',$first_clip_bin,
+ 'Svbin_start',$Svbin_start,'Svbin_end',$Svbin_end);
for (my($gi)=0; $gi<@ush_mu; $gi++) {
&antsOut(depthOfGI($gi), # depth in center of bin
@@ -446,10 +472,12 @@
$uc_ush_mu[$gi],$uc_ush_sig[$gi],
$uc_vsh_mu[$gi],$uc_vsh_sig[$gi],
$uc_wsh_mu[$gi],$uc_wsh_sig[$gi],
- $sh_n[$gi],
+ $sh_n[$gi], # combined
$ush_mu[$gi],$ush_sig[$gi],
$vsh_mu[$gi],$vsh_sig[$gi],
$wsh_mu[$gi],$wsh_sig[$gi],
+ $nSv_prof[$gi]?$sSv_prof[$gi]/$nSv_prof[$gi]:nan,
+ $nSv_prof[$gi],
);
}
@@ -458,7 +486,7 @@
#----------------------------------------------------------------------
if (defined($opt_a)) {
- print(STDERR "Writing acoustic backscatter profiles...");
+ print(STDERR "Writing per-bin acoustic backscatter profiles...");
for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
my($fn) = sprintf("bin%02d.Sv",$bin);
--- a/LADCPproc.backscatter
+++ b/LADCPproc.backscatter
@@ -1,15 +1,17 @@
#======================================================================
# L A D C P P R O C . B A C K S C A T T E R
# doc: Wed Oct 20 13:02:27 2010
-# dlm: Fri Dec 10 00:35:45 2010
+# dlm: Wed Jun 15 15:31:43 2011
# (c) 2010 A.M. Thurnherr
-# uE-Info: 12 52 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 67 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Oct 20, 2010: - created
# Dec 10, 2010: - BUG: backscatter above sea surface made code bomb
# when run with uplooker data
+# Jun 15, 2011: - added calculation of backscatter profiles from
+# subset of bins
#----------------------------------------------------------------------
# Volume Scattering Coefficient, following Deines (IEEE 1999)
@@ -67,23 +69,30 @@
my($gi) = int(&depthOfBin($ens,$bin) / $GRID_DZ);
next if ($gi < 0);
my($range_to_bin) = &dzToBin($ens,$bin) / cos(rad($LADCP{BEAM_ANGLE}));
- $sSv[$gi][$bin] += Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
- $LADCP{TRANSMITTED_PULSE_LENGTH},
- $Er[0],$range_to_bin,
- $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][0])/4 +
- Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
- $LADCP{TRANSMITTED_PULSE_LENGTH},
- $Er[1],$range_to_bin,
- $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][1])/4 +
- Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
- $LADCP{TRANSMITTED_PULSE_LENGTH},
- $Er[2],$range_to_bin,
- $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][2])/4 +
- Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
- $LADCP{TRANSMITTED_PULSE_LENGTH},
- $Er[3],$range_to_bin,
- $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][3])/4;
+ my($Sv) = Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[0],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][0])/4 +
+ Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[1],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][1])/4 +
+ Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[2],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][2])/4 +
+ Sv($LADCP{ENSEMBLE}[$ens]->{CTD_TEMP},
+ $LADCP{TRANSMITTED_PULSE_LENGTH},
+ $Er[3],$range_to_bin,
+ $LADCP{ENSEMBLE}[$ens]->{ECHO_AMPLITUDE}[$bin][3])/4;
+
+ $sSv[$gi][$bin] += $Sv;
$nSv[$gi][$bin]++;
+
+ if ($bin>=$Svbin_start && $bin<=$Svbin_end) {
+ $sSv_prof[$gi] += $Sv;
+ $nSv_prof[$gi]++;
+ }
}
}
}
--- a/LADCPproc.defaults
+++ b/LADCPproc.defaults
@@ -1,9 +1,9 @@
#======================================================================
# L A D C P P R O C . D E F A U L T S
# doc: Fri Sep 17 09:44:21 2010
-# dlm: Thu Dec 9 19:40:05 2010
+# dlm: Wed Jun 15 15:15:41 2011
# (c) 2010 A.M. Thurnherr
-# uE-Info: 40 31 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 23 48 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# default parameters for [mkShearProf]
@@ -20,6 +20,7 @@
# HISTORY:
# Sep 17, 2010: - created
# Dec 9, 2010: - added doc for ASCII CTD file support
+# Jun 15, 2011: - added Svbin_start, Svbin_end
#----------------------------------------------------------------------
# ASCII CTD file support
@@ -44,6 +45,23 @@
# $CTD_ASCII_badval = 999;
#----------------------------------------------------------------------
+# Backscatter Profile Parameters
+#----------------------------------------------------------------------
+
+# The output profile of volume-scattering coefficient is calculated
+# from a subset of the bins only. This is based on the observation,
+# from on a single P403 profile, that the volume-scattering coefficient
+# correction of Deines (IEEE 1999) yield similar results only for
+# bins 3-9.
+
+# The seabed-finding routine, on the other hand, uses acoustic
+# backscatter from all bins, as it is possible that the seabed is only
+# seen in the farthest bins.
+
+$Svbin_start = 3;
+$Svbin_end = 9;
+
+#----------------------------------------------------------------------
# Shear Processing Parameters
#----------------------------------------------------------------------
--- a/LADCPproc.loadCTD
+++ b/LADCPproc.loadCTD
@@ -1,16 +1,17 @@
#======================================================================
# L A D C P P R O C . L O A D C T D
# doc: Thu Dec 9 18:39:01 2010
-# dlm: Mon Jan 10 12:13:54 2011
+# dlm: Sat Jan 22 21:40:12 2011
# (c) 2010 A.M. Thurnherr
-# uE-Info: 31 30 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 174 0 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Dec 9, 2010: - exported from LADCPproc
# - added support for ASCII files
# Dec 16, 2010: - BUG cnv read did not work any more
-# Jan 10, 2010: - added code to skip ANTS header
+# Jan 10, 2011: - added code to skip ANTS header
+# Jan 22, 2011: - adapted to new -g
sub readCTD_ASCII($$)
{
@@ -19,8 +20,10 @@
croak("$fn: unknown pressure field\n") unless defined($CTD_ASCII_press_field);
croak("$fn: unknown temperature field\n") unless defined($CTD_ASCII_temp_field);
croak("$fn: unknown salinity field\n") unless defined($CTD_ASCII_salin_field);
- croak("$fn: unknown latitude field\n") unless defined($CTD_ASCII_lat_field);
- croak("$fn: unknown longitude field\n") unless defined($CTD_ASCII_lon_field);
+ unless (numberp($dtaR->{lat})) {
+ croak("$fn: unknown latitude field\n") unless defined($CTD_ASCII_lat_field);
+ croak("$fn: unknown longitude field\n") unless defined($CTD_ASCII_lon_field);
+ }
$CTD_ASCII_badval = 9e99 unless defined($CTD_ASCII_badval);
@@ -33,7 +36,7 @@
push(@{$dtaR->{press}},($rec[$CTD_ASCII_press_field-1] == $CTD_ASCII_badval) ? nan : $rec[$CTD_ASCII_press_field-1]);
push(@{$dtaR->{temp}}, ($rec[$CTD_ASCII_temp_field-1] == $CTD_ASCII_badval) ? nan : $rec[$CTD_ASCII_temp_field-1]);
push(@{$dtaR->{salin}},($rec[$CTD_ASCII_salin_field-1] == $CTD_ASCII_badval) ? nan : $rec[$CTD_ASCII_salin_field-1]);
- unless ($rec[$CTD_ASCII_lat_field-1] == $CTD_ASCII_badval) {
+ unless (!defined($CTD_ASCII_lat_field) || $rec[$CTD_ASCII_lat_field-1] == $CTD_ASCII_badval) {
$nPos++;
$sumLat += $rec[$CTD_ASCII_lat_field-1];
$sumLon += $rec[$CTD_ASCII_lon_field-1];
--- a/TODO
+++ b/TODO
@@ -1,9 +1,9 @@
======================================================================
T O D O
doc: Thu Oct 14 09:15:02 2010
- dlm: Mon Jan 10 13:29:13 2011
+ dlm: Wed Jun 15 15:16:20 2011
(c) 2010 A.M. Thurnherr
- uE-Info: 17 51 NIL 0 0 72 3 2 4 NIL ofnI
+ uE-Info: 40 0 NIL 0 0 72 3 2 4 NIL ofnI
======================================================================
=LADCPintsh=
@@ -36,6 +36,7 @@
- implement TILT_BIT (others?)
- calculate real acoustic backscatter profile
+ - currently, profiles are arbitarily calculated from bins 3-9
- clean up LADCPproc.UHcode:
- remove weirdnesses