--- a/HISTORY
+++ b/HISTORY
@@ -1,9 +1,9 @@
#======================================================================
# H I S T O R Y
# doc: Thu May 7 13:12:05 2015
-# dlm: Tue Nov 27 13:14:17 2018
+# dlm: Sat Jul 24 09:37:46 2021
# (c) 2015 A.M. Thurnherr
-# uE-Info: 188 70 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 257 70 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
May 7, 2015:
@@ -186,3 +186,72 @@
- published V7.1
----------------------------------------------------------------------
+V7.2 (published Jun 29, 2021)
+ - various bug fixes and minor additions
+ - new features
+ - TEOS10 support
+ - GMT6 support
+ - significant changes
+ - dependency warnings disablecd by default
+----------------------------------------------------------------------
+
+[antsio.pl]:
+ Apr 10, 2019: - disabled dependency warnings
+
+[antsutils.pl]
+# Apr 5, 2019: - disabled weird line of code in antsFunUsage() (see comment)
+# - improved error messages in antsFunUsage()
+# - BUG: antsFunUsage did not work with -ve argc (variable argument funs)
+# Aug 30, 2019: - BUG: antsLoadModel() did not respect $ANTS
+
+[libGM76.pl]
+# Feb 2, 2019: - renamed from libGM.pl to libGM76.pl
+# - replaced beta with m
+# - replaced old code based on Gregg + Kunze formalism with
+# expressions from Thurnherr et al. (GRL 2015)
+# - added GM_strain
+# - BUG: Sw usage message had wrong parameter order
+# - renamed Sw => GM_VKE; Su_z => GM_shear
+# Mar 31, 2019: - updated doc for shear and strain spectra
+# Apr 1, 2019: - added GM_vdiv
+# Apr 5, 2019: - BUG: GM_VKE was erroneous (argument shifting was wrong)
+# - adapted to improved antsFunUsage()
+# Apr 6, 2019: - fiddled during debugging of fs_finestructure issue
+
+[libGMT.pl]
+# Apr 10, 2021: - adapted to GMT6 (suppress warnings)
+# Apr 11, 2021: - added gmt set GMT_AUTO_DOWNLOAD off
+
+[libIMP.pl]
+# Apr 12, 2020: - modified rot_vecs to pass all records with non-numerical elapsed
+# Apr 13, 2020: - added prep_piro_ADCP
+# - BUG: dhist_binsize != 1 did not work
+# - BUG: dhist agreement % was not rounded
+# Apr 14, 2020: - cosmetics
+
+[libRDI_Coords.pl]
+# Jun 5, 2020: - added sscorr_w & sscorr_w_mooring
+# Jun 29, 2020: - added comments for sscorr_w, which conflicts with LADCP_w_ocean
+
+[libSBE.pl]
+# Jan 3, 2019: - BUG: SBE_parseHeader() did not correctly detect missing lat/lon
+
+[libTEOS10.pl]
+# Jan 18, 2019: - created
+
+[libconv.pl]
+# Jan 17, 2019: - added ISO_Datetime()
+# Feb 15, 3019: - added deg2lat, deg2lon
+
+[libstats.pl]
+# Mar 26, 2019: - added regress()
+
+[libvec.pl]
+# Mar 1, 2021: - adapted rotation_ts and angle_ts to deal with nans
+
+
+----------------------------------------------------------------------
+V7.3 (published Jul 24, 2021)
+ - bug fixes [libGHP.pl] [libGMT.pl]
+ - minor improvement to [libIMP.pl]
+----------------------------------------------------------------------
deleted file mode 100644
--- a/OLD
+++ /dev/null
@@ -1,721 +0,0 @@
-#======================================================================
-# O L D
-# doc: Tue Nov 26 21:59:40 2013
-# dlm: Mon Jul 1 15:15:19 2019
-# (c) 2017 A.M. Thurnherr
-# uE-Info: 698 42 NIL 0 0 70 2 2 4 NIL ofnI
-#======================================================================
-
-# HISTORY:
-# Nov 26, 2013: - created
-# Dec 1, 2013: - removed HDG_ROT
-# - added support for IMP data gaps
-# Mar 4, 2014: - added support for DATA_SOURCE_ID
-# May 22, 2014: - use +/-2deg to assess quality of heading offset
-# May 23, 2014: - added $dhist_binsize, $dhist_min_pirom
-# Jul 27, 2016: - updated variable names for consistency
-# Jul 28, 2016: - major re-write if merging routines
-# Jul 29, 2016: - cosmetics
-# - increased heading-offset resolution from 2 to 1 degrees
-# - BUG: inconsistent heading definition used (from old IMP with
-# confused coordinates)
-# - modified initial timelag guess (there was a bug and it is
-# likely more robust based on end time rather than start time)
-# Aug 5, 2016: - BUG: weird statement accessing LADCP_begin-1
-# - BUG: DSID of first ensemble was not left original
-# Aug 22, 2016: - major changes to timelagging
-# Aug 23, 2016: - changed semantics for removing ensembles with bad attitudes:
-# instead of setting attitude to undef (or large pitch/roll),
-# clearEns() is used
-# Aug 24, 2016: - overhauled time-lagging
-# Aug 25, 2016: - significant code cleanup
-# Aug 26, 2016: - added _hdg_err.ps output plot
-# Oct 13, 2016: - made hdg nan for invalid records (BUG with current versions of IMP+LADCP, IMPatchPD0)
-# Nov 22, 2016: - added heading-offset plot
-# - added sensor info to plots
-# Nov 29, 2016: - added stats to compass error plot
-# Dec 29, 2016: - improved histogram plot
-# Nov 16, 2017: - adapted rot_vecs() to KVM coordinates
-# - made sensor information optional in
-# Nov 20, 2017: - major code cleanup
-# Nov 22, 2017: - replaced "IMP" output in routines used by KVH by "IMU"
-# Dec 8, 2017: - replaced remaing "IMP" output (e.g. in plot) by "IMU"
-# May 22, 2018: - added data trimming to rot_vec
-# May 23, 2018: - added horizontal field strength to mag calib plot
-# - added support for S/N 8 board inside UL WH300 (neg_piro == 2)
-# May 24, 2018: - continued working (coord_trans == 2)
-# May 25, 2018: - continued working
-# May 30, 2018: - BUG: trimming did not re-calculate elapsed time
-# Jun 5, 2018: - BUG: on -c main magnetometer and accelerometer vectors were
-# modified twice
-# Jun 7, 2018: - relaxed definition of -e
-# Jun 9, 2018: - BUG: some "edge" data (beyond the valid profile) were bogus (does not matter in practice)
-
-#----------------------------------------------------------------------
-# gRef() library
-#----------------------------------------------------------------------
-
-sub pl_mag_calib_begin($$$) # initialize mag_calib plot
-{
- my($pfn,$plotsize,$axlim) = @_;
- my($xmin,$xmax) = (-$axlim,$axlim);
- my($ymin,$ymax) = (-$axlim,$axlim);
- GMT_begin($pfn,"-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
- GMT_psxy('-Sc0.05');
-}
-
-sub pl_mag_calib_plot($$$) # plot data point
-{
- my($valid,$magX,$magY) = @_;
- if ($valid) { print(GMT "> -Wred -Gred\n$magX $magY\n"); }
- else { print(GMT "> -Wgreen -Ggreen\n$magX $magY\n"); }
-}
-
-sub pl_mag_calib_end($$) # finish mag_calib plot
-{
- my($axlim,$HF_mag,$sensor_info) = @_;
-
- GMT_psxy('-Sc0.1 -Gblue'); # calibration circle
- for (my($a)=0; $a<2*$pi; $a+=0.075) {
- printf(GMT "%f %f\n",$HF_mag*sin($a),$HF_mag*cos($a));
- }
-
- if ($axlim <= 0.1) { # axes labels
- GMT_psbasemap('-Bg1a.04f.001:"X Magnetic Field [Gauss]":/g1a0.02f0.001:"Y Magnetic Field [Gauss]":WeSn');
- } else {
- GMT_psbasemap('-Bg1a.1f.01:"X Magnetic Field [Gauss]":/g1a0.1f0.01:"Y Magnetic Field [Gauss]":WeSn');
- }
- GMT_unitcoords(); # horizontal field strength
- GMT_pstext('-F+f12,Helvetica,blue+jTR -N');
- printf(GMT "0.98 0.98 HF = %.2f Gauss\n",$HF_mag);
-
- printf(GMT "0.98 0.94 $sensor_info\n") # sensor info
- if ($sensor_info ne '');
-
- GMT_pstext('-F+f14,Helvetica,blue+jTL -N'); # profile id
- printf(GMT "0.01 1.06 $P{profile_id}\n");
- GMT_end();
-}
-
-sub rot_vecs(@) # rotate & output IMU vector data
-{
- my($coord_trans,$min_elapsed,$max_elapsed,$plot_milapsed,$plot_malapsed) = @_; # negate KVH pitch/roll data if first arg set to 1
- $min_elapsed = 0 unless defined($min_elapsed);
- $max_elapsed = 9e99 unless defined($max_elapsed);
- $plot_milapsed = $min_elapsed unless defined($plot_milapsed);
- $plot_malapsed = $max_elapsed unless defined($plot_malapsed);
-
- while (&antsIn()) {
- next if ($ants_[0][$elapsedF] < $min_elapsed); # trim data
- last if ($ants_[0][$elapsedF] > $max_elapsed);
-
- my($cpiro) = -1; # current pitch/roll accelerometer
- my(@R); # rotation matrix
- for (my($i)=0; $i<@vecs; $i++) { # rotate vector data
- if ($piro[$i][0] != $cpiro) { # next sensor chip
- $cpiro = $piro[$i][0];
-
- my($accX) = $ants_[0][$vecs[$cpiro][0]];
- my($accY) = $ants_[0][$vecs[$cpiro][1]];
- my($accZ) = $ants_[0][$vecs[$cpiro][2]];
-
- if ($coord_trans == 2) { # S/N 8 inside UL WH300
- $accY *= -1; $accZ *= -1;
- }
-
- my($roll) = atan2($accY,$accZ); # eqn 25 from Freescale AN3461
- my($pitch) = atan2($accX,sqrt($accY**2+$accZ**2)); # eqn 26 from <Freescale AN3461
- if ($coord_trans == 1) { # KVH
- $pitch *= -1;
- $roll *= -1;
- }
- $ants_[0][$piro[$i][1]] = deg($pitch); # add pitch/roll to data
- $ants_[0][$piro[$i][2]] = deg($roll);
-
- my($sp) = sin($pitch); my($cp) = cos($pitch); # define rotation matrix
- my($sr) = sin($roll); my($cr) = cos($roll);
- @R = ([ $cp, 0, -$sp ],
- [-$sp*$sr, $cr, -$cp*$sr],
- [ $sp*$cr, $sr, $cp*$cr]);
- }
- my($xval) = $ants_[0][$vecs[$i][0]];
- my($yval) = $ants_[0][$vecs[$i][1]];
- my($zval) = $ants_[0][$vecs[$i][2]];
-
- if ($coord_trans == 2) { # S/N 8 inside UL WH300
- $yval *= -1; $zval *= -1;
- }
-
- $ants_[0][$vecs[$i][0]] = ($xval-$bias[$i][0]) * $R[0][0] +
- ($yval-$bias[$i][1]) * $R[0][1] +
- ($zval-$bias[$i][2]) * $R[0][2];
- $ants_[0][$vecs[$i][1]] = ($xval-$bias[$i][0]) * $R[1][0] +
- ($yval-$bias[$i][1]) * $R[1][1] +
- ($zval-$bias[$i][2]) * $R[1][2];
- $ants_[0][$vecs[$i][2]] = ($xval-$bias[$i][0]) * $R[2][0] +
- ($yval-$bias[$i][1]) * $R[2][1] +
- ($zval-$bias[$i][2]) * $R[2][2];
- }
-
- my($magX) = $ants_[0][$magXF];
- my($magY) = $ants_[0][$magYF];
- my($magZ) = $ants_[0][$magZF];
- my($accX) = $ants_[0][$accXF];
- my($accY) = $ants_[0][$accYF];
- my($accZ) = $ants_[0][$accZF];
-
- my($HF) = sqrt($magX**2+$magY**2);
- my($valid)= ($HF >= $minfac*$HF_mag) && ($HF <= $maxfac*$HF_mag);
- my($hdg) = $valid ? mag_heading($magX,$magY) : nan;
-
- &antsOut($ants_[0][$elapsedF]-$min_elapsed,$ants_[0][$tempF],
- RDI_pitch($ants_[0][$pitchF],$ants_[0][$rollF]),$ants_[0][$rollF],
- $hdg,$accX,$accY,$accZ,$magX,$magY,$magZ,
- vel_u($HF,$hdg),vel_v($HF,$hdg),$valid);
-
- pl_mag_calib_plot($valid,$magX,$magY)
- if defined($P{profile_id}) &&
- ($ants_[0][$elapsedF] >= $plot_milapsed) &&
- ($ants_[0][$elapsedF] <= $plot_malapsed);
- }
-}
-
-#----------------------------------------------------------------------
-# LADCP merging library
-#----------------------------------------------------------------------
-
-#-----------------------------------------------------
-# Instrument Offset Estimation
-#
-# 1: resolution of histogram in deg
-# 1 deg okay for good sensors
-# 2 deg sometimes required (2016 P18 003 2nd sensor)
-# 2: min tilt anom to consider for offset estimation
-# 0.3 deg works even for calm casts
-# increased values improve histogram
-# 2.0 deg is too high for quiet casts (2016 P18 003)
-# 3: minimum fraction of hist mode required
-# 10% default
-# decreasing histogram resolution is better than
-# decreasing this value, I think
-#-----------------------------------------------------
-
-#----------------------------------------------------------------------
-# trim_out_of_water()
-# - attempt to remove out-of-water records from time-series of
-# horizontal acceleration
-#----------------------------------------------------------------------
-
-sub trim_out_of_water($)
-{
- my($verbose) = @_;
-
- #--------------------------------------------------------------------------
- # first-difference horizontal acceleration at full resolution to pass-filter
- # dangling motion
- #--------------------------------------------------------------------------
-
- $IMP{Ah}[0] = sqrt($ants_[0][$accXF]**2+$ants_[0][$accYF]**2);
- $IMP{dAhdt}[0] = nan;
- for (my($r)=1; $r<@ants_; $r++) {
- $IMP{Ah}[$r] = sqrt($ants_[$r][$accXF]**2+$ants_[$r][$accYF]**2);
- $IMP{dAhdt}[$r] = ($IMP{Ah}[$r]-$IMP{Ah}[$r-1]) / ($ants_[$r][$elapsedF]-$ants_[$r-1][$elapsedF]);
- }
-
- #--------------------------------------------------------------------------------------
- # create 10-s binned time series to calculate rms values of this quantity (dAhdt), and
- # scale this with cos(sqrt($$pitch**2+$$roll**2)) to dampen underwater peaks (when the
- # instrument has a large tilt because it is being dragged)
- # NB: dAhdt, pitch and roll are only set up to last full bin (not so, sum and n)
- #--------------------------------------------------------------------------------------
-
- my(@dAhdt,@pitch,@roll,@dAhdt_rms);
- my(@sum) = my(@sume) = my(@sump) = my(@sumr) = my(@n) = (0);
-
- my($bin_start) = $ants_[0][$elapsedF];
- for (my($r)=1; $r<@ants_; $r++) {
-
- if ($ants_[$r][$elapsedF] - $bin_start <= 10) { # within 10-s bin
- $sum[$#sum] += $IMP{dAhdt}[$r]; # sums
- $sume[$#sume] += $ants_[$r][$elapsedF];
- $sump[$#sump] += $ants_[$r][$pitchF];
- $sumr[$#sumr] += $ants_[$r][$rollF];
- $n[$#n]++;
- next;
- }
-
- $dAhdt[$#sum] = $sum[$#sum] / $n[$#n]; # bin done => means
- $elapsed[$#sum] = $sume[$#sume] / $n[$#n];
- $pitch[$#sum] = $sump[$#sump] / $n[$#n];
- $roll[$#sum] = $sumr[$#sumr] / $n[$#n];
-
- my($sumsq) = 0; # sum of squares for rms(accel)
- for (my($rr)=$r-$n[$#n]; $rr<$r; $rr++) {
- $sumsq += ($IMP{dAhdt}[$rr] - $dAhdt[$#sum])**2;
- }
- $dAhdt_rms[$#sum] = sqrt($sumsq / $n[$#n]);
-
- push(@sum,$IMP{dAhdt}[$r]); # begin next bin
- push(@sume,$ants_[$r][$elapsedF]);
- push(@sump,$ants_[$r][$pitchF]);
- push(@sumr,$ants_[$r][$rollF]);
- push(@n,1);
- $bin_start = $ants_[$r][$elapsedF];
- }
-
- #--------------------------------------------
- # trim beginning/end when IMP is out of water
- #--------------------------------------------
-
- my($i,$si);
- for ($i=int(@dAhdt_rms/2); $i>0; $i--) {
- last if ($dAhdt_rms[$i] * cos(rad(sqrt($pitch[$i]**2+$roll[$i]**2))) > 1.0);
- }
- if ($dAhdt_rms[$i] * cos(rad(sqrt($pitch[$i]**2+$roll[$i]**2))) > 1.0) {
- for ($si=0; $ants_[$si][$elapsedF]<=$elapsed[$i]; $si++) {}
- splice(@ants_,0,$si);
- printf(STDERR "\n\t\t%5d leading out-of-water IMU records removed",$si)
- if ($si>0 && $verbose);
- } else {
- print(STDERR "\n\t\tWARNING: no leading out-of-water IMU records detected/removed") if $verbose;
- }
-
- for ($i=int(@dAhdt_rms/2); $i<@dAhdt_rms; $i++) {
- last if ($dAhdt_rms[$i] * cos(rad(sqrt($pitch[$i]**2+$roll[$i]**2))) > 1.0);
- }
- if ($dAhdt_rms[$i] * cos(rad(sqrt($pitch[$i]**2+$roll[$i]**2))) > 1.0) {
- for ($si=$#ants_; $ants_[$si][$elapsedF]>=$elapsed[$i]; $si--) {}
- my($rem) = @ants_ - $si;
- splice(@ants_,$si);
- printf(STDERR "\n\t\t%5d trailing out-of-water IMU records removed",$rem)
- if ($rem>0 && $verbose);
- } else {
- print(STDERR "\n\t\tWARNING: no trailing out-of-water IMU records detected/removed") if $verbose;
- }
-
- printf(STDERR "\n\t\tcast duration : %.1f min",
- ($ants_[$#ants_][$elapsedF] - $ants_[0][$elapsedF]) / 60)
- if $verbose;
-}
-
-#--------------------------------------------------------------------------------
-# prep_piro_IMP()
-# - calculate pitch/roll offsets & tilt azimuth for IMP
-# - during an attempt to improve time lagging for the 2015 IOPAS data set,
-# it was noticed that one particular instrument, WHM300#12973 maxes out
-# pitch at 27.36 degrees, whereas the roll may not be maxed out at 28.76 deg,
-# the max observed during the cruise.
-# - therefore, IMP{TILT_AZIM} and IMP{TILT_ANOM} are calculated here, first,
-# with a pitch/roll cutoff value of 29 degrees
-# - after the time lagging, when the LADCP start and end times are known,
-# the TILT values are re-calculated without the pitch/roll limit, and
-# using only the correct time range
-#---------------------------------------------------------------------------------
-
-sub prep_piro_IMP($)
-{
- my($verbose) = @_;
- my($RDI_max_tilt) = 29;
- my($IMP_pitch_mean,$IMP_roll_mean,$nPR) = (0,0,0);
-
- for (my($r)=0; $r<@ants_; $r++) {
- next unless numbersp($ants_[$r][$pitchF],$ants_[$r][$rollF]);
- $nPR++;
- $IMP_pitch_mean += min($ants_[$r][$pitchF],$RDI_max_tilt);
- $IMP_roll_mean += min($ants_[$r][$rollF],$RDI_max_tilt);
- }
- $IMP_pitch_mean /= $nPR;
- $IMP_roll_mean /= $nPR;
- printf(STDERR "\n\t\tIMU mean pitch/roll : %.1f/%.1f deg",$IMP_pitch_mean,$IMP_roll_mean)
- if $verbose;
-
- for (my($r)=0; $r<@ants_; $r++) {
- next unless numbersp($ants_[$r][$pitchF],$ants_[$r][$rollF]);
- $IMP{TILT_AZIMUTH}[$r] = tilt_azimuth(min($ants_[$r][$pitchF],$RDI_max_tilt)-$IMP_pitch_mean,
- min($ants_[$r][$rollF],$RDI_max_tilt) -$IMP_roll_mean);
- $IMP{TILT_ANOM}[$r] = angle_from_vertical(min($ants_[$r][$pitchF],$RDI_max_tilt)-$IMP_pitch_mean,
- min($ants_[$r][$rollF],$RDI_max_tilt) -$IMP_roll_mean);
- }
- return ($IMP_pitch_mean,$IMP_roll_mean);
-}
-
-#----------------------------------------------------------------------
-# prep_piro_LADCP()
-# - calculate pitch/roll offsets & tilt azimuth of LADCP
-#----------------------------------------------------------------------
-
-sub prep_piro_LADCP($)
-{
- my($verbose) = @_;
-
- my($LADCP_pitch_mean,$LADCP_roll_mean) = (0,0);
- for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
- $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
- gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
- $LADCP_pitch_mean += $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH};
- $LADCP_roll_mean += $LADCP{ENSEMBLE}[$ens]->{ROLL};
- }
- $LADCP_pitch_mean /= ($LADCP_end-$LADCP_begin+1);
- $LADCP_roll_mean /= ($LADCP_end-$LADCP_begin+1);
- printf(STDERR "\n\t\tLADCP mean pitch/roll : %.1f/%.1f deg",$LADCP_pitch_mean,$LADCP_roll_mean)
- if $verbose;
-
- for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
- $LADCP{ENSEMBLE}[$ens]->{TILT_AZIMUTH} =
- tilt_azimuth($LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH}-$LADCP_pitch_mean,
- $LADCP{ENSEMBLE}[$ens]->{ROLL}-$LADCP_roll_mean);
- $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} =
- angle_from_vertical($LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH}-$LADCP_pitch_mean,
- $LADCP{ENSEMBLE}[$ens]->{ROLL}-$LADCP_roll_mean);
- }
-
- print(STDERR "\n") if $verbose;
- return ($LADCP_pitch_mean,$LADCP_roll_mean);
-}
-
-#------------------------------------------------------------------
-# sub calc_hdg_offset()
-# - estimate heading offset from tilt time series
-# - returns heading offset and updated IMP mean tilts
-# - also creates diagnostic plot with pl_hdg_offset()
-#------------------------------------------------------------------
-
-sub pl_hdg_offset($@)
-{
- my($dhist_binsize,$modefrac,@dhist) = @_;
-
- my($plotsize) = '13c';
- my($xmin,$xmax) = (-180.5,180.5);
- my($ymin) = 0;
- my($ymax) = 1.05 * $dhist[$HDG_offset];
-
- GMT_begin("$P{profile_id}${opt_a}_hdg_offset.ps","-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
- GMT_psxy("-Sb${dhist_binsize}u -GCornFlowerBlue");
- for (my($i)=0; $i<@dhist; $i++) {
- next unless $dhist[$i];
- printf(GMT "%f $dhist[$i]\n",$i*$dhist_binsize>180 ? $i*$dhist_binsize-360 : $i*$dhist_binsize);
- }
- GMT_psbasemap('-Bg45a90f15:"IMU Heading Offset [\260]":/ga100f10:"Frequency":WeSn');
- GMT_unitcoords();
- GMT_pstext('-F+f14,Helvetica,CornFlowerBlue+jTR -N');
- printf(GMT "0.99 1.06 %g \260 offset (%d%% agreement)\n",angle($HDG_offset),100*$modefrac);
- GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
- printf(GMT "0.01 1.06 $P{profile_id} $opt_a\n");
- GMT_end();
-}
-
-sub calc_hdg_offset($)
-{
- my($verbose) = @_;
-
- print(STDERR "\n\tRe-calculating IMU pitch/roll anomalies") if $verbose;
- ($IMP_pitch_mean,$IMP_roll_mean,$nPR) = (0,0,0);
- for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
- my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
- if ($r < 0 && $ens == $LADCP_begin) {
- $r = int(($LADCP{ENSEMBLE}[++$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT})
- while ($r < 0);
- printf(STDERR "\n\tIMU data begin with instrument already in water => skipping %ds of LADCP data",
- $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME})
- if ($verbose);
- $LADCP_begin = $ens;
- }
- if ($r > $#ants_) {
- printf(STDERR "\n\tIMU data end while instrument is still in water => truncating %ds of LADCP data",
- $LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME})
- if ($verbose);
- $LADCP_end = $ens - 1;
- last;
- }
- next unless numberp($IMP{TILT_AZIMUTH}[$r]);
- $nPR++;
- $IMP_pitch_mean += $ants_[$r][$pitchF];
- $IMP_roll_mean += $ants_[$r][$rollF];
- }
- $IMP_pitch_mean /= $nPR;
- $IMP_roll_mean /= $nPR;
-
- for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
- my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
- next unless numberp($IMP{TILT_AZIMUTH}[$r]);
- $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_AZIMUTH} =
- $IMP{TILT_AZIMUTH}[$r] = tilt_azimuth($ants_[$r][$pitchF]-$IMP_pitch_mean,
- $ants_[$r][$rollF] -$IMP_roll_mean);
- $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM} =
- $IMP{TILT_ANOM}[$r] = angle_from_vertical($ants_[$r][$pitchF]-$IMP_pitch_mean,
- $ants_[$r][$rollF] -$IMP_roll_mean);
- }
-
- printf(STDERR "\n\t\tIMU mean pitch/roll: %.1f/%.1f deg",$IMP_pitch_mean,$IMP_roll_mean)
- if $verbose;
-
- if (defined($opt_s)) {
- $HDG_offset = $opt_s;
- printf(STDERR "\n\tHEADING_OFFSET = %.1f (set by -s)\n",$HDG_offset) if $verbose;
- } else {
- my($dhist_binsize,$dhist_min_pirom,$dhist_min_mfrac) = split(/,/,$opt_e);
- croak("$0: cannot decode -e $opt_e\n")
- unless ($dhist_binsize > 0 && $dhist_min_pirom > 0 && $dhist_min_mfrac >= 0);
-
- my(@dhist); my($nhist) = my($modeFreq) = 0;
- for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
- my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
- next unless numberp($IMP{TILT_AZIMUTH}[$r]);
- next unless (abs($ants_[$r][$pitchF]-$IMP_pitch_mean) >= $dhist_min_pirom &&
- abs($ants_[$r][$rollF] -$IMP_roll_mean) >= $dhist_min_pirom);
- $dhist[int(angle_pos($LADCP{ENSEMBLE}[$ens]->{TILT_AZIMUTH}-$IMP{TILT_AZIMUTH}[$r])/$dhist_binsize+0.5)]++;
- $nhist++;
- }
- croak("$0: empty histogram\n")
- unless ($nhist);
-
- $HDG_offset = 0;
- for (my($i)=1; $i<@dhist-1; $i++) { # make sure mode is not on edge
- $HDG_offset = $i if ($dhist[$i] >= $dhist[$HDG_offset]);
- }
- my($modefrac) = ($dhist[$HDG_offset]+$dhist[$HDG_offset-1]+$dhist[$HDG_offset+1]) / $nhist;
-
- $HDG_offset *= $dhist_binsize;
- pl_hdg_offset($dhist_binsize,$modefrac,@dhist);
-
- if ($opt_f) {
- printf(STDERR "\n\nIGNORED WARNING (-f): Cannot determine reliable heading offset; $HDG_offset+/-$dhist_binsize deg accounts for only %f%% of total\n",$modefrac*100)
- if ($modefrac < $dhist_min_mfrac);
- } else {
- croak(sprintf("\n$0: Cannot determine reliable heading offset; $HDG_offset+/-$dhist_binsize deg accounts for only %f%% of total\n",$modefrac*100))
- if ($modefrac < $dhist_min_mfrac);
- }
-
- printf(STDERR "\n\t") if $verbose;
- printf(STDERR "IMU heading offset = %g deg (%d%% agreement)\n",angle($HDG_offset),100*$modefrac) if $verbose;
- }
- return ($HDG_offset,$IMP_pitch_mean,$IMP_roll_mean);
-}
-
-#-----------------------------------------------------------
-# rot_IMP()
-# - rotate IMP Data Into LADCP Instrument Coords
-# - also replaced pitch/roll by corresponding anomalies!!!
-#-----------------------------------------------------------
-
-sub rot_IMP($)
-{
- my($verbose) = @_;
- my($crho) = cos(rad($HDG_offset));
- my($srho) = sin(rad($HDG_offset));
-
- for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
- my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
-
- if (numbersp($ants_[$r][$pitchF],$ants_[$r][$rollF])) { # pitch/roll
- my($rot_p) = (($ants_[$r][$pitchF]-$IMP_pitch_mean)*$crho +
- ($ants_[$r][$rollF]-$IMP_roll_mean)*$srho);
- my($rot_r) = (-($ants_[$r][$pitchF]-$IMP_pitch_mean)*$srho +
- ($ants_[$r][$rollF]-$IMP_roll_mean)*$crho);
- $ants_[$r][$pitchF] = $rot_p;
- $ants_[$r][$rollF] = $rot_r;
- }
-
- $ants_[$r][$hdgF] = angle_pos($ants_[$r][$hdgF] - $HDG_offset)
- if numberp($ants_[$r][$hdgF]);
- }
-
- my($rot_p) = $IMP_pitch_mean * $crho + $IMP_roll_mean * $srho; # mean pitch roll
- my($rot_r) = -$IMP_pitch_mean * $srho + $IMP_roll_mean * $crho;
- $IMP_pitch_mean = $rot_p;
- $IMP_roll_mean = $rot_r;
-
- print(STDERR "\n") if $verbose;
- return ($IMP_pitch_mean,$IMP_roll_mean);
-}
-
-#----------------------------------------------------------------------
-# create_merge_plots()
-# - tilt time series (*_time_lag.ps)
-#----------------------------------------------------------------------
-
-sub create_merge_plots($$$)
-{
- my($basename,$plotsize,$verbose) = @_;
-
- #---------------------------------
- # Tilt Time Series (*_time_lag.ps)
- #---------------------------------
-
- my($mint,$maxt) = (99,-99);
- for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
- if (($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME} <= 180) ||
- ($LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} <= 180)) {
- $mint = $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM} if ($LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM} < $mint);
- $maxt = $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM} if ($LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM} > $maxt);
- $mint = $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} if ($LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} < $mint);
- $maxt = $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} if ($LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} > $maxt);
- }
- }
-
- my($xmin,$xmax) = (-90,90);
- my($ymin) = round($mint-0.5);
- my($ymax) = round($maxt+0.5);
-
- GMT_begin("${basename}_time_lag.ps","-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
-
- GMT_psxy('-W2,coral');
- for (my($ens) = $LADCP_begin + 5;
- $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME} <= 90;
- $ens++) {
- printf(GMT "%f %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME},
- $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM});
- }
- GMT_psxy('-W1');
- for (my($ens) = $LADCP_begin + 5;
- $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME} <= 90;
- $ens++) {
- printf(GMT "%f %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_begin]->{ELAPSED_TIME},
- $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM});
- }
- GMT_psxy('-W2,SeaGreen');
- for (my($ens) = $LADCP_end - 5;
- $LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} <= 90;
- $ens--) {
- printf(GMT "%f %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME},
- $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM});
- }
- GMT_psxy('-W1');
- for (my($ens) = $LADCP_end - 5;
- $LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} <= 90;
- $ens--) {
- printf(GMT "%f %f\n",$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME},
- $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM});
- }
-
- GMT_psbasemap('-Bg30a30f10:"Elapsed Time [sec]":/g5a5f1:"Tilt Magnitude [\260]":WeSn');
- GMT_unitcoords();
- GMT_pstext('-F+f14,Helvetica,Coral+jTL');
- printf(GMT "0.52 0.98 downcast\n");
- GMT_pstext('-F+f14,Helvetica,SeaGreen+jTR');
- printf(GMT "0.48 0.98 upcast\n");
- GMT_psxy('-W4,LightSkyBlue');
- printf(GMT "0.5 0\n0.5 1\n");
- GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
- printf(GMT "0.01 1.06 $P{profile_id} $opt_a\n");
- GMT_end();
-
- #------------------------------
- # Heading Errors (*_hdg_err.ps)
- #------------------------------
-
- my(@err_binned,@err_nsamp);
- my($sumErr) = 0; my($nErr) = $LADCP_end - $LADCP_begin + 1;
- for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
- next unless ($LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} < 10);
- my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
- my($bi) = int($ants_[$r][$hdgF]/5);
- my($err) = angle_diff($ants_[$r][$hdgF],$LADCP{ENSEMBLE}[$ens]->{HEADING});
- next unless numberp($err);
- $err_binned[$bi] += $err; $sumErr += $err;
- $err_nsamp[$bi]++;
- }
- for (my($bi)=0; $bi<@err_nsamp; $bi++) {
- $err_binned[$bi] = ($err_nsamp[$bi] >= 5)
- ? $err_binned[$bi]/$err_nsamp[$bi]
- : undef;
- }
- my(@err_dssq);
- for (my($ens)=$LADCP_begin; $ens<=$LADCP_end; $ens++) {
- next unless ($LADCP{ENSEMBLE}[$ens]->{TILT_ANOM} < 10);
- my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
- my($bi) = int($ants_[$r][$hdgF]/5);
- my($err) = angle_diff($ants_[$r][$hdgF],$LADCP{ENSEMBLE}[$ens]->{HEADING});
- next unless numberp($err);
- $err_dssq[$bi] += ($err-$err_binned[$bi])**2;
- }
-
- my($xmin,$xmax) = (0,360);
- my($ymin,$ymax) = (-45,45);
- GMT_begin("${basename}_hdg_err.ps","-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
- GMT_psxy('-Ey/3,CornFlowerBlue');
- my($sumSq,$sumBe) = my($nSq,$nBe) = (0,0);
- for (my($bi)=0; $bi<@err_binned; $bi++) {
- next unless ($err_nsamp[$bi] >= 2);
- next unless numberp($err_binned[$bi]);
- $sumSq += $err_binned[$bi]**2; $nSq++;
- $sumBe += $err_binned[$bi]; $nBe++;
-# printf(GMT "%f %f\n",2.5+5*$bi,$err_binned[$bi]);
- printf(GMT "%f %f %f\n",2.5+5*$bi,$err_binned[$bi],sqrt($err_dssq[$bi]/($err_nsamp[$bi]-1)));
- }
- GMT_psbasemap('-Bg90a45f5:"ADCP Heading [\260]":/g15a15f5:"ADCP Compass Error [\260]":WeSn');
- GMT_unitcoords();
- GMT_pstext('-F+f12,Helvetica,CornFlowerBlue+jTR -Gwhite -C25%');
- printf(GMT "0.98 0.98 rms error = %7.1f \260\n",sqrt($sumSq/$nSq));
- printf(GMT "0.98 0.94 time-averaged error = %7.1f \260\n",$sumErr/$nErr);
- printf(GMT "0.98 0.90 heading-averaged error = %7.1f \260\n",$sumBe/$nBe);
- GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
- printf(GMT "0.01 1.06 $P{profile_id} $opt_a\n");
- GMT_end();
-
- print(STDERR "\n") if $verbose;
-}
-
-#----------------------------------------------------------------------
-# output_merged()
-# - output merged data
-#----------------------------------------------------------------------
-
-sub output_merged($)
-{
- my($verbose) = @_;
-
- my($tazimF) = &antsNewField('tilt_azimuth');
- my($tanomF) = &antsNewField('tilt_magnitude');
- my($L_tazimF) = &antsNewField('LADCP_tilt_azimuth');
- my($L_tanomF) = &antsNewField('LADCP_tilt_magnitude');
- my($L_elapsedF) = &antsNewField('LADCP_elapsed');
- my($L_ensF) = &antsNewField('LADCP_ens');
- my($L_depthF) = &antsNewField('LADCP_depth_estimate');
- my($L_pitchF) = &antsNewField('LADCP_pitch');
- my($L_rollF) = &antsNewField('LADCP_roll');
- my($L_hdgF) = &antsNewField('LADCP_hdg');
- my($dcF) = &antsNewField('downcast');
-
- for (my($ens)=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
- my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});
-# print(STDERR "$r\[$ens\,$LADCP{ENSEMBLE}[$ens]->{NUMBER}] = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG} - $ants_[0][$elapsedF]) / $IMP{DT});\n");
- if ($r<0 || $r>$#ants_) { # ensemble beyond limits of IMU data
- my(@out);
- if ($ens >= $LADCP_begin && $ens <= $LADCP_end) {
- $out[$elapsedF] = $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} + $IMP{TIME_LAG};
- $out[$L_elapsedF] = $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME};
- }
- $out[$L_ensF] = $LADCP{ENSEMBLE}[$ens]->{NUMBER};
- $out[$dcF] = ($ens <= $LADCP_bottom);
- &antsOut(@out);
- } elsif ($ens < $LADCP_begin || $ens > $LADCP_end) { # pre deplyment or post recovery
- $ants_[$r][$L_elapsedF] = undef;
- $ants_[$r][$L_ensF] = $LADCP{ENSEMBLE}[$ens]->{NUMBER};
- $ants_[$r][$L_pitchF] = $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} - $LADCP_pitch_mean;
- $ants_[$r][$L_rollF] = $LADCP{ENSEMBLE}[$ens]->{ROLL} - $LADCP_roll_mean;
- $ants_[$r][$L_hdgF] = $LADCP{ENSEMBLE}[$ens]->{HEADING};
- $ants_[$r][$dcF] = nan;
- &antsOut(@{$ants_[$r]});
- } else {
- $ants_[$r][$tazimF] = angle($LADCP{ENSEMBLE}[$ens]->{IMP_TILT_AZIMUTH} + $HDG_offset);
- $ants_[$r][$tanomF] = $LADCP{ENSEMBLE}[$ens]->{IMP_TILT_ANOM};
- $ants_[$r][$L_tazimF] = $LADCP{ENSEMBLE}[$ens]->{TILT_AZIMUTH};
- $ants_[$r][$L_tanomF] = $LADCP{ENSEMBLE}[$ens]->{TILT_ANOM};
- $ants_[$r][$L_elapsedF] = $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME};
- $ants_[$r][$L_ensF] = $LADCP{ENSEMBLE}[$ens]->{NUMBER};
- $ants_[$r][$L_depthF] = $LADCP{ENSEMBLE}[$ens]->{DEPTH};
- $ants_[$r][$L_pitchF] = $LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} - $LADCP_pitch_mean;
- $ants_[$r][$L_rollF] = $LADCP{ENSEMBLE}[$ens]->{ROLL} - $LADCP_roll_mean;
- $ants_[$r][$L_hdgF] = $LADCP{ENSEMBLE}[$ens]->{HEADING};
- $ants_[$r][$dcF] = ($ens <= $LADCP_bottom);
- &antsOut(@{$ants_[$r]});
- }
- }
-
- print(STDERR "\n") if $verbose;
-}
-
-#----------------------------------------------------------------------
-
-1; # return true for all the world to see
--- a/ants.pl
+++ b/ants.pl
@@ -2,9 +2,9 @@
#======================================================================
# A N T S . P L
# doc: Fri Jun 19 14:01:06 1998
-# dlm: Mon Apr 13 11:17:43 2020
+# dlm: Thu Jul 1 18:46:26 2021
# (c) 1998 A.M. Thurnherr
-# uE-Info: 33 21 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 29 51 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -26,11 +26,12 @@
# Dec 8, 2017: - updated to V7.0 (to avoid absolute-path symlink)
# Nov 27, 2018: - updated to V7.1 (for LADCP_w 1.4 release)
# Apr 13, 2020: - updated to V7.2 (for MIMP_tools)
+# Jul 1, 2021: - updated to V7.3 because of a bug
exec($ENV{ANTS_PERL},$0,@ARGV),die("$ENV{ANTS_PERL}: $!")
if (defined($ENV{ANTS_PERL}) && $^X ne $ENV{ANTS_PERL});
-$antsLibVersion = 7.2;
+$antsLibVersion = 7.3;
die(sprintf("$0: obsolete library V%.1f; V%.1f required\n",
$antsLibVersion,$antsMinLibVersion))
--- a/antsio.pl
+++ b/antsio.pl
@@ -2,9 +2,9 @@
#======================================================================
# A N T S I O . P L
# doc: Fri Jun 19 19:22:51 1998
-# dlm: Wed Apr 10 16:57:59 2019
+# dlm: Tue Nov 30 12:28:03 2021
# (c) 1998 A.M. Thurnherr
-# uE-Info: 217 48 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 220 50 NIL 0 0 70 10 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -215,6 +215,10 @@
# Apr 5, 2017: - BUG: stale file mtime dependency info was not printed correctly
# Apr 23, 2018: - BUG: @antsLayout was not kept up-to-date when layout-changes are allowed
# Apr 10, 2019: - disabled dependency warnings
+# Oct 14, 2021: - changed nan to NaN, to make gnuplot work nicely again
+# Nov 30, 2021: - finally made -F not croak on division-by-zero, etc.
+# - BUG: some NaNs were still nans
+# HISTORY END
# GENERAL NOTES:
# - %P was named without an ants-prefix because associative arrays
@@ -491,8 +495,14 @@
# Handle Layout changes:
# - only allow when $antsAllowEmbeddedLayoutChange is set
# - ensure that antsLayout always contains up-to-date Layout
- croak("$0: embedded layout change when reading file $ARGV <@antsLayout> -> <@Layout>")
- if (!$antsAllowEmbeddedLayoutChange && @Layout && @antsLayout && ("@Layout" ne "@antsLayout"));
+ if (!$antsAllowEmbeddedLayoutChange && @Layout && @antsLayout && ("@Layout" ne "@antsLayout")) {
+ my(@OL) = @antsLayout;
+ my(@NL) = @Layout;
+ while ($OL[0] eq $NL[0]) {
+ shift(@OL); shift(@NL);
+ }
+ croak("$0: embedded layout change when reading file $ARGV <@OL[0..3] ...> -> <@NL[0..3] ...>")
+ }
@antsLayout = @Layout if (@Layout);
$P{RECNO} = -1 unless defined($P{RECNO}); # set pseudo %PARAMs
@@ -525,17 +535,17 @@
$P{RECNO}++; # update pseudo %PARAMs
$P{LINENO} = ($ARGV eq $lastFile) ? $P{LINENO}+1 : 0;
- s/[Nn][Aa][Nn]/nan/g; # make all nans lower case
+ s/[Nn][Aa][Nn]/NaN/g; # turn all nan into NaN
local(@in) = split($opt_I); # needs to be local for -S
if (@in > $antsBufNFields) { # increase # of fields to expect
$antsBufNFields = @in;
for ($i=0; $i<$#ants_; $i++) { # update recs already in buffer
- push(@{$ants_[$i]},nan)
+ push(@{$ants_[$i]},NaN)
while ($#{$ants_[$i]}+1 < $antsBufNFields);
}
}
- push(@in,nan) # pad current record
+ push(@in,NaN) # pad current record
while (@in<$antsBufNFields);
if (defined($opt_S)) { # -S)elect
@@ -570,12 +580,12 @@
### $antsBufNFields = $#{$ants_[$#ants_]} + 1;
#### print("antsBufNFields := $antsBufNFields --- $_");
### for ($i=0; $i<$#ants_; $i++) {
-### push(@{$ants_[$i]},nan)
+### push(@{$ants_[$i]},NaN)
### while ($#{$ants_[$i]}+1 < $antsBufNFields);
### }
### }
### }
-### push(@{$ants_[$#ants_]},nan) # pad this
+### push(@{$ants_[$#ants_]},NaN) # pad this
### while ($#{$ants_[$#ants_]}+1 < $antsBufNFields);
}
@@ -670,7 +680,8 @@
} elsif (defined($OEfield[$f])) {
$out[$f] = $out_buf[$OEfield[$f]];
} else {
- $out[$f] = &{$OEexpr[$f]};
+ eval(\'$out[$f] = &{$OEexpr[$f]};\'); # print errors and continue
+ print(STDERR $@) unless ($@ eq "");
}
}
}
@@ -692,14 +703,14 @@
# STEP 5: PRINT DATA
$antsPadOut = @ofn if ($antsPadOut >= 0 && @ofn);
- push(@out,nan) while (@out < $antsPadOut);
+ push(@out,NaN) while (@out < $antsPadOut);
my($outStr);
for (my($fnr)=0; $fnr<=$#out; $fnr++) {
$out[$fnr] =
fmtNum($out[$fnr],
@antsNewLayout ? $antsNewLayout[$fnr] : $antsLayout[$fnr]);
- $outStr .= (defined($out[$fnr]) && $out[$fnr] ne "" ? $out[$fnr] : nan)
+ $outStr .= (defined($out[$fnr]) && $out[$fnr] ne "" ? $out[$fnr] : NaN)
. ($fnr == $#out ? $opt_R : $opt_O);
}
print($outStr);
@@ -755,7 +766,7 @@
$antsBufNFields = $f+1 # auto extension
if ($antsBufNFields-1 < $f);
while ($#{$ants_[$r]} < $f-1) {
- push(@{$ants_[$r]},nan);
+ push(@{$ants_[$r]},NaN);
}
$ants_[$r][$f] = $v;
}
--- a/antsnc.pl
+++ b/antsnc.pl
@@ -23,13 +23,18 @@
# Mar 20, 2008: - added progress output to NC_stringify
# Jul 21, 2009: - allowed for suppression of %PARAMs
# Jan 15, 2016: - BUG: %DEPS pseudo-%PARAM was encoded
+# Apr 5, 2022: - added NC_writeMDataMulti()
+# - implemented disappeared NetCDF::DOUBLE etc.
+# - enabled fill value handling (library bug has gone away)
+# Apr 22, 2022: - added "_coordinate" to abscissa dimension to
+# allow reading with python xr
+# HISTORY END
# NOTES:
# - multi-valued attribs are not loaded by getInfo()
# - spaces in NC strings are replaced by underscores
-# - data filling is disabled, because of a bug in the NetCDF library
-# NetCDF Library Bug:
+# NetCDF Library Bug: NO LONGER ACTIVE AS OF APR 2022
# The library appears to have incorrect default _FillValue types for
# integer data types. The error appears if the "setfill" line is commented
# out and the following command is run:
@@ -40,6 +45,48 @@
use NetCDF;
+#----------------------------------------------------------------------
+# NetCDF Constants
+# - as of April, 2022 these are no longer part of NetCDF???
+# - manually encoded from .h file
+#----------------------------------------------------------------------
+
+sub NetCDF::NAT { return 0; }
+sub NetCDF::BYTE { return 1; }
+sub NetCDF::CHAR { return 2; }
+sub NetCDF::SHORT { return 3; }
+sub NetCDF::INT { return 4; }
+sub NetCDF::LONG { return 4; }
+sub NetCDF::FLOAT { return 5; }
+sub NetCDF::DOUBLE { return 6; }
+sub NetCDF::UBYTE { return 7; }
+sub NetCDF::USHORT { return 8; }
+sub NetCDF::UINT { return 9; }
+sub NetCDF::INT64 { return 10; }
+sub NetCDF::UINT64 { return 11; }
+sub NetCDF::STRING { return 12; }
+
+sub NetCDF::FILL_BYTE { return -127; }
+sub NetCDF::FILL_CHAR { return "\0"; }
+sub NetCDF::FILL_SHORT { return -32767; }
+sub NetCDF::FILL_INT { return -2147483647; }
+sub NetCDF::FILL_LONG { return -2147483647; }
+sub NetCDF::FILL_FLOAT { return 9.9692099683868690e+36; }
+sub NetCDF::FILL_DOUBLE { return 9.9692099683868690e+36; }
+sub NetCDF::FILL_UBYTE { return 255; }
+sub NetCDF::FILL_USHORT { return 65535; }
+sub NetCDF::FILL_UINT { return 4294967295; }
+sub NetCDF::FILL_INT64 { return -9223372036854775806; }
+sub NetCDF::FILL_UINT64 { return 18446744073709551614; }
+sub NetCDF::FILL_STRING { return ''; }
+
+sub NetCDF::NOWRITE { return 0x0000; }
+sub NetCDF::CLOBBER { return 0x0000; }
+sub NetCDF::NOFILL { return 0x100; }
+sub NetCDF::GLOBAL { return -1; }
+sub NetCDF::UNLIMITED { return 0; }
+
+
#----------------------------------
# string representation of NC types
#----------------------------------
@@ -265,7 +312,7 @@
}
croak("$0: varid != fnr (implementation restriction)")
unless ($vid == $f);
- foreach my $anm (keys(%P)) { # variable attributes
+ foreach my $anm (keys(%P)) { # VARIABLE ATTRIBUTES
next unless defined($P{$anm});
my($var,$attr) = ($anm =~ m{([^:]+):(.*)});
next unless ($var eq $antsLayout[$f]);
--- a/antsusage.pl
+++ b/antsusage.pl
@@ -2,9 +2,9 @@
#======================================================================
# A N T S U S A G E . P L
# doc: Fri Jun 19 13:43:05 1998
-# dlm: Wed Dec 13 09:09:58 2017
+# dlm: Tue May 17 15:35:36 2022
# (c) 1998 A.M. Thurnherr
-# uE-Info: 167 49 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 168 89 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -165,6 +165,8 @@
# Dec 9, 2017: - added -E, $antsSuppressCommonOptions
# - common options cosmetics
# Dec 13, 2017: - BUG: common options cosmetics
+# May 17, 2022: - BUG: antsFileError() did not report permission errors in a useful way
+# HISTORY END
# NOTES:
# - ksh expands {}-arguments with commas in them!!! Use + instead
@@ -542,6 +544,8 @@
sub antsNoFileErr($$)
{
croak("$0: $_[0] $_[1] is not a valid file\n")
+ unless (-f $_[1]);
+ croak("$0: $_[0] $_[1] is not readable\n")
unless (-r $_[1]);
&antsAddDeps($_[1]);
}
--- a/antsutils.pl
+++ b/antsutils.pl
@@ -2,9 +2,9 @@
#======================================================================
# A N T S U T I L S . P L
# doc: Fri Jun 19 23:25:50 1998
-# dlm: Fri Apr 5 16:21:54 2019
+# dlm: Tue Apr 5 21:20:29 2022
# (c) 1998 A.M. Thurnherr
-# uE-Info: 106 55 NIL 0 0 70 10 2 4 NIL ofnI
+# uE-Info: 157 0 NIL 0 0 70 10 2 4 NIL ofnI
#======================================================================
# Miscellaneous auxillary functions
@@ -106,6 +106,8 @@
# - improved error messages in antsFunUsage()
# - BUG: antsFunUsage did not work with -ve argc (variable argument funs)
# Aug 30, 2019: - BUG: antsLoadModel() did not respect $ANTS
+# Nov 29, 2021: - made fmtNum() return NaN on undefined input
+# HISTORY END
# fnr notes:
# - matches field names starting with the string given, i.e. "sig" is
@@ -240,6 +242,7 @@
{
my($num,$fname) = @_;
+ $num = NaN unless defined($num);
$num = 0 if ($num eq '-0'); # perl 5.8.8: 0*-0.1 = -0, which is
# not handled correctly by all progs
$num = str2num($num) if ($opt_C);
--- a/libGHP.pl
+++ b/libGHP.pl
@@ -1,17 +1,19 @@
#======================================================================
# L I B G H P . P L
# doc: Fri Sep 7 09:56:08 2012
-# dlm: Mon Oct 22 13:10:47 2012
+# dlm: Wed Aug 11 13:12:52 2021
# (c) 2012 A.M. Thurnherr
-# uE-Info: 11 0 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 38 0 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
# Sep 7, 2012: - created
# Oct 22, 2012: - cosmetics
+# Jul 14, 2021: - BUG: adapted to new libGM name
+# Aug 11, 2021: - modified j() to handle N<f
require "$ANTS/libfuns.pl"; # arccosh
-require "$ANTS/libGM.pl"; # GM_N0
+require "$ANTS/libGM76.pl"; # GM_N0
#----------------------------------------------------------------------------
# h(R_omega) correction factor for different shear/strain (R_omega) ratios
@@ -29,17 +31,30 @@
return 3*($R_omega+1) / (2*sqrt(2)*$R_omega*sqrt($R_omega-1));
}
+sub h2($)
+{
+ my($R_omega) = @_;
+ return 1/(6*sqrt(2)) * $R_omega*($R_omega+1) / sqrt($R_omega-1);
+}
+
#----------------------------------------------------------------------------
# j(f,N) correction factor for latitude
# - this version from Kunze et al. (2006)
+# - if N<f, N/f = 1 assumed
#----------------------------------------------------------------------------
sub j(@)
{
my($f,$N) = @_;
- return 0 if ($f == 0);
+
+ $f = abs($f);
+ return 0 if ($f < 1e-6);
+
+ my($f30) = &f(30);
+
$N = $GM_N0 unless defined($N);
- my($f30) = &f(30);
+ $N = $f unless ($N > $f);
+
return $f*acosh($N/$f) / ($f30*acosh($GM_N0/$f30));
}
--- a/libGMT.pl
+++ b/libGMT.pl
@@ -1,9 +1,9 @@
#======================================================================
# L I B G M T . P L
# doc: Sun Jun 14 13:45:47 2015
-# dlm: Sun Apr 11 09:55:22 2021
+# dlm: Thu Jul 1 18:45:15 2021
# (c) 2015 A.M. Thurnherr
-# uE-Info: 47 34 NIL 0 0 72 2 2 4 NIL ofnI
+# uE-Info: 48 66 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================
# perl implementation of /Data/Makefiles/Makefile.GMT
@@ -45,6 +45,7 @@
# Mar 17, 2016: - added check for gmt5 on load
# Apr 10, 2021: - adapted to GMT6 (suppress warnings)
# Apr 11, 2021: - added gmt set GMT_AUTO_DOWNLOAD off
+# Jul 1, 2021: - BUG: gmt check was based on psxy, not gmt psxy
$DEBUG = 0;
@@ -53,7 +54,7 @@
#----------------------------------------------------------------------
if (`which gmt` eq '') {
- if (`which psxy` eq '') {
+ if (`which gmt` eq '') {
croak("$0: [libGMT.pl] GMT version 6 required\n");
} else {
croak("$0: [libGMT.pl] GMT version 6 required (gmt4 installed)\n");
--- a/libIMP.pl
+++ b/libIMP.pl
@@ -1,9 +1,9 @@
#======================================================================
# L I B I M P . P L
# doc: Tue Nov 26 21:59:40 2013
-# dlm: Tue Apr 14 17:23:47 2020
+# dlm: Sat Jul 24 09:37:27 2021
# (c) 2017 A.M. Thurnherr
-# uE-Info: 736 58 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 70 13 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -65,6 +65,9 @@
# - BUG: dhist_binsize != 1 did not work
# - BUG: dhist agreement % was not rounded
# Apr 14, 2020: - cosmetics
+# Jul 12, 2021: - improvements to histogram plot
+# Jul 24, 2021: - updated docu
+# HISTORY END
#----------------------------------------------------------------------
# gRef() library
@@ -458,9 +461,18 @@
my($plotsize) = '13c';
my($xmin,$xmax) = (-180.5,180.5);
my($ymin) = 0;
- my($ymax) = 1.05 * $dhist[$HDG_offset/$dhist_binsize];
+ my($ymax) = 0;
+
+ for (my($i)=0; $i<@dhist; $i++) {
+ $ymax = $dhist[$i] if ($dhist[$i] > $ymax);
+ }
+ $ymax = 1.05 * $ymax;
GMT_begin("$P{profile_id}${opt_a}_hdg_offset.ps","-JX${plotsize}","-R$xmin/$xmax/$ymin/$ymax",'-X6 -Y4 -P');
+ if (defined($opt_o)) {
+ GMT_psxy("-W2,red");
+ printf(GMT "%f %f\n%f %f\n",$opt_o,0,$opt_o,$ymax);
+ }
GMT_psxy("-Sb${dhist_binsize}u -GCornFlowerBlue");
for (my($i)=0; $i<@dhist; $i++) {
next unless $dhist[$i];
--- a/libLADCP.pl
+++ b/libLADCP.pl
@@ -1,9 +1,9 @@
#======================================================================
-# L I B L A D C P . P L
+# . . / L I B / L I B L A D C P . P L
# doc: Wed Jun 1 20:38:19 2011
-# dlm: Sat Apr 6 19:42:07 2019
+# dlm: Tue Sep 14 08:17:30 2021
# (c) 2011 A.M. Thurnherr
-# uE-Info: 211 34 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 201 0 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -20,12 +20,15 @@
# - added T_w()
# Sep 24, 2012: - made "k" argument default in T_w()
# Oct 25, 2012: - renamed T_SM() to T_ASM()
-# Jun 26. 2013: - added T_w_z()
+# Jun 26, 2013: - added T_w_z()
# - added parameter checks to processing-specific corrections
# May 18, 2015: - added pulse length to T_w() and T_w_z()
# Apr 25, 2018: - added eps_VKE() parameterization
# Apr 5, 2018: - adapted to improved antsFunUsage()
# - BUG eps_VKE() had erroneous string
+# Sep 9, 2021: - removed T_w_z(); consistent with Thurnherr (JAOT 2012),
+# vertical divergence spectra should be corrected with
+# T_w()
require "$ANTS/libvec.pl";
require "$ANTS/libfuns.pl";
@@ -202,6 +205,14 @@
# T_w(k,blen,plen,dz,range_max)
# - vertical-velocity method of Thurnherr (IEEE 2011)
# - range_max == 0 disables tilt correction
+# - this is the expression that should be used to correct spectra
+# of VKE and vertical divergence (w_z)
+# - I am not sure why the finite differencing of the processed
+# velocities does not have to be corrected for, but this is
+# consistent with:
+# - Thurnherr (JAOT 2012) where the shear is corrected with T_VI
+# - the strain correction of Kunze et al. (2006), which corrects
+# for bin averaging, but not for finite differencing(?)
#----------------------------------------------------------------------
{ my(@fc);
@@ -217,25 +228,25 @@
}
}
-#----------------------------------------------------------------------
-# T_w_z(k,blen,plen,dz,range_max)
-# - vertical-velocity method of Thurnherr (IEEE 2011)
-# - first differencing of gridded shear to calculate dw/dz
-# - NB: grid-scale differentiation assumed
-# - range_max == 0 disables tilt correction
-#----------------------------------------------------------------------
-
-{ my(@fc);
- sub T_w_z(@)
- {
- my($k,$blen,$plen,$dz,$range_max) =
- &antsFunUsage(-4,'ffff',
- 'T_w_z([vertical wavenumber[rad/s]] <ADCP bin size[m]> <pulse length[m]> <output grid resolution[m]> <range max[m]>)',
- \@fc,'k',undef,undef,undef,@_);
- croak("T_w_z($k,$blen,$plen,$dz,$range_max): bad parameters\n")
- unless ($k>=0 && $blen>0 && $plen>0 && $dz>0 && $range_max>=0);
- return T_ravg($k,$blen,$plen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max)) * T_fdiff($k,$dz);
- }
-}
+##----------------------------------------------------------------------
+## T_w_z(k,blen,plen,dz,range_max)
+## - vertical-velocity method of Thurnherr (IEEE 2011)
+## - first differencing of gridded shear to calculate dw/dz
+## - NB: grid-scale differentiation assumed
+## - range_max == 0 disables tilt correction
+##----------------------------------------------------------------------
+#
+#{ my(@fc);
+# sub T_w_z(@)
+# {
+# my($k,$blen,$plen,$dz,$range_max) =
+# &antsFunUsage(-4,'ffff',
+# 'T_w_z([vertical wavenumber[rad/s]] <ADCP bin size[m]> <pulse length[m]> <output grid resolution[m]> <range max[m]>)',
+# \@fc,'k',undef,undef,undef,@_);
+# croak("T_w_z($k,$blen,$plen,$dz,$range_max): bad parameters\n")
+# unless ($k>=0 && $blen>0 && $plen>0 && $dz>0 && $range_max>=0);
+# return T_ravg($k,$blen,$plen) * T_binavg($k,$dz) * T_tilt($k,dprime($range_max)) * T_fdiff($k,$dz);
+# }
+#}
1;
new file mode 100644
--- /dev/null
+++ b/libStokes.pl
@@ -0,0 +1,30 @@
+#======================================================================
+# L I B S T O K E S . P L
+# doc: Tue Feb 15 12:41:27 2022
+# dlm: Tue Feb 15 12:55:44 2022
+# (c) 2022 A.M. Thurnherr
+# uE-Info: 30 2 NIL 0 0 70 2 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Feb 15, 2022: - created
+
+#----------------------------------------------------------------------
+# Stokes Settling Velocity
+#----------------------------------------------------------------------
+
+sub Stv(@)
+{
+ my($rho_p,$r_p,$g,$rho_f,$mu_f) = @_;
+
+ $g = 9.81 unless defined($g); # acceleration due to gravity
+ $rho_f = 1000 unless defined($rho_f); # fluid density
+ $mu_f = 1e-3 unless defined($mu_f); # dynamic viscosity (order of magnitude)
+
+ croak("Usage: Stv(rho_p,r_p[,g[,rho_f[,mu_f]]]\n")
+ unless numbersp($rho_p,$r_p,$g,$rho_f,$mu_f);
+
+ return 2/9 * ($rho_p - $rho_f)/$mu_f * $g * $r_p**2;
+}
+
+1;
--- a/libconv.pl
+++ b/libconv.pl
@@ -1,9 +1,9 @@
#======================================================================
# L I B C O N V . P L
# doc: Sat Dec 4 13:03:49 1999
-# dlm: Fri Feb 15 13:21:08 2019
+# dlm: Wed Apr 6 18:41:02 2022
# (c) 1999 A.M. Thurnherr
-# uE-Info: 529 24 NIL 0 0 70 2 2 4 NIL ofnI
+# uE-Info: 73 59 NIL 0 0 70 2 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -67,7 +67,11 @@
# Dec 18, 2017: - removed ambiguous-date warning
# May 22, 2018: - added NMEA2dec_time()
# Jan 17, 2019: - added ISO_Datetime()
-# Feb 15, 3019: - added deg2lat, deg2lon
+# Feb 15, 2019: - added deg2lat, deg2lon
+# Aug 27, 2021: - improved Date to take dn from %PARAMs
+# Nov 15, 2021: - modified O2 coversion routines to return nan on missing press, temp, salin
+# Apr 6, 2022: - made =Date work with %dn params as well
+# HISTORY END
require "$ANTS/libEOS83.pl"; # &sigma()
require "$ANTS/libPOSIX.pl"; # &floor()
@@ -329,6 +333,24 @@
}
}
+ unless (defined($dnf)) {
+ foreach my $p (keys(%P)) {
+ next unless ($p =~ /^dn(\d\d)$/);
+ $year = ($1 > 70) ? 1900 + $1 : 2000 + $1;
+ $day = int($P{$p});
+ while ($day > 365+&leapYearP($year)) { # adjust year
+ $day -= 365 + &leapYearP($year);
+ $year++;
+ }
+ my($month) = 1;
+ while ($day > &monthLength($year,$month)) {
+ $day -= &monthLength($year,$month);
+ $month++;
+ }
+ return sprintf('%04d/%02d/%02d',$year,$month,$day);
+ }
+ }
+
my($year,$day) = &antsFunUsage(2,"cf","epoch, dayNo",\@fc,undef,$dnf,@_);
$year += ($year < 50) ? 2000 : 1900 # Y2K
@@ -361,6 +383,26 @@
last;
}
+ unless (defined($dnf)) {
+ foreach my $p (keys(%P)) {
+ next unless ($p =~ /^dn(\d\d)$/);
+ my($fday) = $P{$p};
+ my($day) = int($fday);
+ $fday -= $day;
+
+ my($hour) = int(24*$fday);
+ $fday -= $hour/24;
+ my($min) = int(24*60*$fday);
+ $fday -= $min/24/60;
+ my($sec) = round(24*3600*$fday);
+ $min++,$sec=0 if ($sec == 60);
+ $hour++,$min=0 if ($min == 60);
+ $day++,$hour=0 if ($hour == 24);
+
+ return sprintf('%02d:%02d:%02d',$hour,$min,$sec);
+ }
+ }
+
my($fday) = &antsFunUsage(1,"f","dayNo",\@fc,$dnf,@_);
my($day) = int($fday);
$fday -= $day;
@@ -594,7 +636,9 @@
{ my(@fc);
sub O2mlpl2umpkg(@)
{
- return nan if isnan($_[3]);
+# return nan if isnan($_[3]);
+ return nan
+ if (@_ == 4 && !numberp($_[0]+$_[1]+$_[2]+$_[3]));
my($S,$T,$P,$mlpl) =
&antsFunUsage(4,'ffff','[S, T, P [dbar], O2 [ml/l]]',
\@fc,'salin','temp','press','O2',@_);
@@ -605,7 +649,9 @@
{ my(@fc);
sub O2umpkg2mlpl(@)
{
- return nan if isnan($_[3]);
+# return nan if isnan($_[3]);
+ return nan
+ if (@_ == 4 && !numberp($_[0]+$_[1]+$_[2]+$_[3]));
my($S,$T,$P,$umpkg) =
&antsFunUsage(4,'ffff','[S, T, P [dbar], O2 [ml/l]]',
\@fc,'salin','temp','press','O2',@_);
--- a/libstats.pl
+++ b/libstats.pl
@@ -1,9 +1,9 @@
#======================================================================
-# . . / L I B / L I B S T A T S . P L
+# L I B S T A T S . P L
# doc: Wed Mar 24 13:59:27 1999
-# dlm: Tue Mar 26 12:59:32 2019
+# dlm: Sat Aug 14 13:08:30 2021
# (c) 1999 A.M. Thurnherr
-# uE-Info: 50 31 NIL 0 0 72 0 2 4 NIL ofnI
+# uE-Info: 86 1 NIL 0 0 72 0 2 4 NIL ofnI
#======================================================================
# HISTORY:
@@ -79,6 +79,12 @@
return $sig / sqrt($dof);
}
+sub dof($)
+{
+ my($nsamp) = &antsFunUsage(1,"c","nsamp",@_);
+ return ($nsamp > 1) ? $nsamp-1 : 1;
+}
+
#----------------------------------------------------------------------
# calc standard stats from vector of vals
# - used, e.g., in [rfilt]
new file mode 100644
--- /dev/null
+++ b/libwaterdepth.pl
@@ -0,0 +1,24 @@
+#======================================================================
+# L I B W A T E R D E P T H . P L
+# doc: Mon Feb 7 16:06:56 2022
+# dlm: Mon Feb 7 16:20:09 2022
+# (c) 2022 A.M. Thurnherr
+# uE-Info: 10 27 NIL 0 0 70 0 2 4 NIL ofnI
+#======================================================================
+
+# HISTORY:
+# Feb 7, 2022: - created
+
+sub waterdepth(@)
+{
+ my($lon,$lat) = &antsFunUsage(2,"ff","lon, lat",@_);
+ open(my $F, "-|", "waterdepth $lon $lat") || croak("waterdepth: $!\n");
+ my($wd);
+ chomp($wd = <$F>);
+ croak("Cannot decode waterdepth output ($wd)\n")
+ unless numberp($wd);
+ close($F);
+ return $wd;
+}
+
+1;