.
authorAndreas Thurnherr <ant@ldeo.columbia.edu>
Thu, 07 May 2020 11:43:48 -0400
changeset 41 fa41b3a72c97
parent 40 c1803ae2540f
child 42 22f5d5d35236
.
ANTSlib
ants.pl
libDeines99.pl
libIMP.pl
--- a/ANTSlib
+++ b/ANTSlib
@@ -2,9 +2,9 @@
 #======================================================================
 #                    A N T S L I B 
 #                    doc: Wed May 16 06:19:16 2012
-#                    dlm: Sun Apr  5 16:41:52 2015
+#                    dlm: Mon Apr 13 11:18:17 2020
 #                    (c) 2012 A.M. Thurnherr
-#                    uE-Info: 39 42 NIL 0 0 72 0 2 4 NIL ofnI
+#                    uE-Info: 19 103 NIL 0 0 72 0 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -16,7 +16,7 @@
 #				  - BUG: version consistency check works only on woosher
 
 ($ANTSLIB) = ($0 =~ m{^(.*)/[^/]*$});
-$antsMinLibVersion = 6.0;
+$antsMinLibVersion = 6.0;						# don't change this; change antsLibVersion in [ants.pl]
 
 require "$ANTSLIB/ants.pl";
 #require "$ANTSLIB/libCPT.pl";
--- 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: Tue Nov 27 12:57:36 2018
+#                    dlm: Mon Apr 13 11:17:43 2020
 #                    (c) 1998 A.M. Thurnherr
-#                    uE-Info: 27 60 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 33 21 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -24,12 +24,13 @@
 #  Mar 12, 2017: - updated to V6.8 (for LADCP_w 1.3 release)
 #  Nov 20, 2017: - updated to V6.9 (for DT KVH software)
 #  Dec  8, 2017: - updated to V7.0 (to avoid absolute-path symlink)
-#  Nov 27, 2018: - updaetd to V7.1 (for LADCP_w 1.4 release)
+#  Nov 27, 2018: - updated to V7.1 (for LADCP_w 1.4 release)
+#  Apr 13, 2020: - updated to V7.2 (for MIMP_tools)
 
 exec($ENV{ANTS_PERL},$0,@ARGV),die("$ENV{ANTS_PERL}: $!")
     if (defined($ENV{ANTS_PERL}) && $^X ne $ENV{ANTS_PERL});
 
-$antsLibVersion = 7.1;
+$antsLibVersion = 7.2;
 
 die(sprintf("$0: obsolete library V%.1f; V%.1f required\n",
 	$antsLibVersion,$antsMinLibVersion))
new file mode 100644
--- /dev/null
+++ b/libDeines99.pl
@@ -0,0 +1,22 @@
+#======================================================================
+#                    L I B D E I N E S 9 9 . P L 
+#                    doc: Wed Apr 15 11:57:01 2020
+#                    dlm: Wed Apr 15 11:57:01 2020
+#                    (c) 2020 A.M. Thurnherr
+#                    uE-Info: 14 2 NIL 0 0 70 0 2 4 NIL ofnI
+#======================================================================
+
+sub Sv($$$$$)		  
+{
+	my($temp,$pulse_length,$noise_level,$range,$echo_amplitude) = @_;
+	my($C)		= -143; 				# RDI WHM300 (from Deines)
+	my($Ldbm)	= 10 * log10($pulse_length);
+	my($PdbW)	= 14.0;
+	my($alpha)	= 0.069;
+	my($Kc) 	= 0.45;
+		    
+	return $C + 10*log10(($temp+273)*$range**2) - $Ldbm - $PdbW
+			  + 2*$alpha*$range + $Kc*($echo_amplitude-$noise_level);
+}
+
+1;
--- 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: Sun Apr 12 10:03:38 2020
+#                    dlm: Tue Apr 14 17:23:47 2020
 #                    (c) 2017 A.M. Thurnherr
-#                    uE-Info: 181 61 NIL 0 0 72 0 2 4 NIL ofnI
+#                    uE-Info: 736 58 NIL 0 0 72 0 2 4 NIL ofnI
 #======================================================================
 
 # HISTORY:
@@ -61,6 +61,10 @@
 #				  - BUG: GIMBAL_PITCH had not only been defined for in-water records
 #	Jun  9, 2018: - BUG: some "edge" data (beyond the valid profile) were bogus (does not matter in practice)
 #	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
 
 #----------------------------------------------------------------------
 # gRef() library
@@ -313,6 +317,8 @@
 #--------------------------------------------------------------------------------
 # prep_piro_IMP()
 # 	- calculate pitch/roll offsets & tilt azimuth for IMP
+#	- if first_rec and last_rec are provided, the mean calculation will 
+#	  be restricted to this range (TILT_AZIM, _ANOM calculated for all recs)
 #	- 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,
@@ -324,13 +330,16 @@
 #	  using only the correct time range
 #---------------------------------------------------------------------------------
 
-sub prep_piro_IMP($)
+sub prep_piro_IMP($@)
 {
-	my($verbose) = @_;
+	my($verbose,$first_rec,$last_rec) = @_;
 	my($RDI_max_tilt) = 29; 
 	my($IMP_pitch_mean,$IMP_roll_mean,$nPR) = (0,0,0);
+
+	$first_rec = 0 			unless defined($first_rec);
+	$last_rec  = $#ants_	unless defined($last_rec);
 	
-	for (my($r)=0; $r<@ants_; $r++) {
+	for (my($r)=$first_rec; $r<=$last_rec; $r++) {
 		next unless numbersp($ants_[$r][$pitchF],$ants_[$r][$rollF]);
 		$nPR++;
 		$IMP_pitch_mean += min($ants_[$r][$pitchF],$RDI_max_tilt);
@@ -360,6 +369,9 @@
 {
 	my($verbose) = @_;
 
+	$first_ens = 0 			unless defined($first_ens);
+	$last_ens  = $#ants_	unless defined($last_ens);
+
 	my($LADCP_pitch_mean,$LADCP_roll_mean) = (0,0);
 	for (my($ens)=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
 		$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
@@ -389,6 +401,49 @@
 	return ($LADCP_pitch_mean,$LADCP_roll_mean);
 }
 
+#----------------------------------------------------------------------
+# prep_piro_ADCP()
+#	- calculate pitch/roll offsets & tilt azimuth of moored ADCP
+#	- very similar to prep_piro_LADCP() but using windowing as in
+#	  prep_piro_IMP()
+#	- window ens are indices, not ensemble numbers
+#----------------------------------------------------------------------
+
+sub prep_piro_ADCP($@)
+{
+	my($verbose,$window_first_ens,$window_last_ens) = @_;
+
+	$window_first_ens = 0 					unless defined($window_first_ens);
+	$window_last_ens  = $#{$ADCP{ENSEMBLE}}	unless defined($window_last_ens);
+
+	for (my($ens)=0; $ens<=$#{$ADCP{ENSEMBLE}}; $ens++) {
+		$ADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =
+			gimbal_pitch($ADCP{ENSEMBLE}[$ens]->{PITCH},$ADCP{ENSEMBLE}[$ens]->{ROLL});
+	}
+
+	my($ADCP_pitch_mean,$ADCP_roll_mean) = (0,0);
+	for (my($ens)=$window_first_ens; $ens<=$window_last_ens; $ens++) {
+		$ADCP_pitch_mean += $ADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH};
+		$ADCP_roll_mean  += $ADCP{ENSEMBLE}[$ens]->{ROLL};
+	}
+	$ADCP_pitch_mean /= ($window_last_ens-$window_first_ens+1);
+	$ADCP_roll_mean  /= ($window_last_ens-$window_first_ens+1);
+	printf(STDERR "\n\t\tADCP mean pitch/roll       : %.1f/%.1f deg",$ADCP_pitch_mean,$ADCP_roll_mean)
+			if $verbose;
+	
+	for (my($ens)=0; $ens<=$#{$ADCP{ENSEMBLE}}; $ens++) {
+		$ADCP{ENSEMBLE}[$ens]->{TILT_AZIMUTH} =
+			tilt_azimuth($ADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH}-$ADCP_pitch_mean,
+						 $ADCP{ENSEMBLE}[$ens]->{ROLL}-$ADCP_roll_mean);
+		$ADCP{ENSEMBLE}[$ens]->{TILT_ANOM} =
+			angle_from_vertical($ADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH}-$ADCP_pitch_mean,
+								$ADCP{ENSEMBLE}[$ens]->{ROLL}-$ADCP_roll_mean);
+	}
+	
+	print(STDERR "\n") if $verbose;
+	return ($ADCP_pitch_mean,$ADCP_roll_mean);
+}
+
 #------------------------------------------------------------------
 # sub calc_hdg_offset()
 #	- estimate heading offset from tilt time series
@@ -403,7 +458,7 @@
 	my($plotsize) = '13c';
 	my($xmin,$xmax) = (-180.5,180.5);
 	my($ymin) = 0;
-	my($ymax) = 1.05 * $dhist[$HDG_offset];
+	my($ymax) = 1.05 * $dhist[$HDG_offset/$dhist_binsize];
 	    
 	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");
@@ -415,9 +470,9 @@
 	GMT_unitcoords();
 	GMT_pstext('-F+f14,Helvetica,CornFlowerBlue+jTR -N');
 	if (defined($opt_o)) {
-		printf(GMT "0.99 1.06 -o %g \260 offset (%d%% agreement)\n",angle($HDG_offset),100*$modefrac);
+		printf(GMT "0.99 1.06 -o %g\260 offset (%d%% agreement)\n",angle($HDG_offset),round(100*$modefrac));
 	} else {
-		printf(GMT "0.99 1.06 %g \260 offset (%d%% agreement)\n",angle($HDG_offset),100*$modefrac);
+		printf(GMT "0.99 1.06 %g\260 offset (%d%% agreement)\n",angle($HDG_offset),round(100*$modefrac));
 	}
 	GMT_pstext('-F+f14,Helvetica,blue+jTL -N');
 	printf(GMT "0.01 1.06 $P{profile_id} $opt_a\n");
@@ -487,6 +542,7 @@
 	
 	if (defined($opt_o)) {												# set heading offset, either with -o
 		$HDG_offset = $opt_o;
+		$dhist_binsize = 1;
 	} else {															# or from the data
 		$HDG_offset = 0;
 		for (my($i)=1; $i<@dhist-1; $i++) { 							# make sure mode is not on edge
@@ -495,7 +551,7 @@
 		$HDG_offset *= $dhist_binsize;
 	}
 
-	my($modefrac) = ($dhist[$HDG_offset]+$dhist[$HDG_offset-1]+$dhist[$HDG_offset+1]) / $nhist;
+	my($modefrac) = ($dhist[$HDG_offset/$dhist_binsize]+$dhist[$HDG_offset/$dhist_binsize-1]+$dhist[$HDG_offset/$dhist_binsize+1]) / $nhist;
 	pl_hdg_offset($dhist_binsize,$modefrac,@dhist);
 
 	unless (defined($opt_o)) {	    									# make sure data are consistent (unless -o is used)
@@ -510,10 +566,10 @@
 	
 	printf(STDERR "\n\t") if $verbose;
 	if ($opt_o) {
-		printf(STDERR "IMU heading offset = -o %g deg (%d%% agreement)\n",angle($HDG_offset),100*$modefrac)
+		printf(STDERR "IMU heading offset = -o %g deg (%d%% agreement)\n",angle($HDG_offset),round(100*$modefrac))
 			if $verbose;
 	} else {
-		printf(STDERR "IMU heading offset = %g deg (%d%% agreement)\n",angle($HDG_offset),100*$modefrac)
+		printf(STDERR "IMU heading offset = %g deg (%d%% agreement)\n",angle($HDG_offset),round(100*$modefrac))
 			if $verbose;
 	}
 
@@ -675,9 +731,9 @@
     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);
+    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();