listW
changeset 0 229a0d72d2ab
child 14 8c79b38a7086
new file mode 100755
--- /dev/null
+++ b/listW
@@ -0,0 +1,221 @@
+#!/usr/bin/perl
+#======================================================================
+#                    L I S T W 
+#                    doc: Wed Mar 24 06:45:09 2004
+#                    dlm: Thu Jul 30 17:42:33 2009
+#                    (c) 2004 A.M. Thurnherr
+#                    uE-Info: 205 53 NIL 0 0 72 2 2 4 NIL ofnI
+#======================================================================
+
+# dump vertical velocities
+
+# NB: currently broken
+
+# HISTORY:
+#	Mar 24, 2004: - created from [listens] and [mkprofile]
+#	Mar 27, 2004: - added elapsed field
+#				  - floatized time
+#	Apr  3, 2004: - cosmetics
+#	Nov  8, 2005: - UNIXTIME => UNIX_TIME
+#	Sep 19, 2007: - adapted to new [RDI_BB_Read.pl] (not tested)
+#	Jul 30, 2009: - NaN => nan
+
+$0 =~ m{(.*)/[^/]+}; 
+require "$1/WorkhorseBinRead.pl";
+require "$1/WorkhorseCoords.pl";
+require "$1/WorkhorseUtils.pl";
+
+use Getopt::Std;
+
+$USAGE = "$0 @ARGV";
+die("Usage: $0 " .
+		"[-A)nts] " .
+		"[-F)ilter <script>] " .
+		"[bin -r)ange  <bin|0,bin|*>] " .
+		"[-e)rr-vel <max|0.1>] [-c)orrelation <min|70>] " .
+		"[-S)alin <val|35>] [-t)emp <bias>] " .
+		"[output -f)ields <field[,...]> " .
+		"<RDI file>\n")
+	unless (&getopts("Ac:e:F:f:r:S:t:") && @ARGV == 1);
+
+$opt_e = 0.1 unless defined($opt_e);				# defaults
+$opt_c = 70	 unless defined($opt_c);
+$opt_S = 35  unless defined($opt_S);
+print(STDERR "WARNING: Using uncalibrated ADCP temperature!\n"),$opt_t = 0
+	 unless defined($opt_t);
+
+require $opt_F if defined($opt_F);					# load filter
+
+if ($opt_f) {										# additional fields
+	@f = split(',',$opt_f);
+	foreach $f (@f) {
+		my($f) = $f;								# copy it
+		$f =~ s{\[.*$}{};							# remove indices
+		$addFields .= " {$f}";
+	}
+}
+
+#----------------------------------------------------------------------
+
+print(STDERR "Reading $ARGV[0]...");				# read data
+readData($ARGV[0],\%dta);
+print(STDERR "done\n");
+
+if (defined($opt_r)) {								# bin range
+	($minb,$maxb) = split(',',$opt_r);
+	die("$0: can't decode -r $opt_r\n") unless defined($maxb);
+} else {
+	$minb = 0;
+	$maxb = $dta{N_BINS} - 1;
+}
+
+die("$ARGV[0]: not enough bins for choice of -r\n")	# enough bins?
+	unless ($dta{N_BINS} >= $maxb);
+
+if ($dta{BEAM_COORDINATES}) {						# coords used
+	$beamCoords = 1;
+} elsif (!$dta{EARTH_COORDINATES}) {
+	die("$ARGV[0]: only beam and earth coordinates implemented so far\n");
+}
+
+#----------------------------------------------------------------------
+# Reference-Layer w (from [mkprofile])
+#	- also sets W field when valid
+#----------------------------------------------------------------------
+
+sub ref_lr_w($)
+{
+	my($ens) = @_;
+	my($i,$n,$w);
+
+	for ($i=$minb; $i<=$maxb; $i++) {
+		next if ($dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][0] < $opt_c ||
+				 $dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][1] < $opt_c ||
+				 $dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][2] < $opt_c ||
+				 $dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][3] < $opt_c);
+		if ($beamCoords) {
+			next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] < 100 ||
+					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] < 100 ||
+					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] < 100 ||
+					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);
+			@v = velInstrumentToEarth(\%dta,$ens,
+					velBeamToInstrument(\%dta,
+						@{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$i]}));
+		} else {
+			next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] > 0 ||
+					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] > 0 ||
+					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] > 0 ||
+					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);
+			@v = @{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$i]};
+		}
+		next if (!defined($v[3]) || abs($v[3]) > $opt_e);
+
+		if (defined($v[2])) {							# valid w
+			$dta{ENSEMBLE}[$ens]->{W}[$i] = $v[2];
+			$w += $v[2]; $n++;
+		}
+	}
+	return $n ? $w/$n : undef;
+}
+
+#----------------------------------------------------------------------
+
+print(STDERR "Generating profile by integrating w...");
+
+for ($e=0; $e<=$#{$dta{ENSEMBLE}}; $e++) {
+	checkEnsemble(\%dta,$e);								# sanity checks
+	filterEnsemble(\%dta,$e)								# filter ensemble 
+		if (defined($opt_F) &&
+			$dta{ENSEMBLE}[$e]->{PERCENT_GOOD}[0][0] > 0);
+
+	$dta{ENSEMBLE}[$e]->{REFW} = ref_lr_w($e);
+	next unless defined($dta{ENSEMBLE}[$e]->{REFW});
+
+	unless (defined($firstgood)) {							# init profile
+		$firstgood = $lastgood = $e;			
+		$depth = 0;
+	}
+
+	my($dt) = $dta{ENSEMBLE}[$e]->{UNIX_TIME} -			# time step since
+			  $dta{ENSEMBLE}[$lastgood]->{UNIX_TIME};		# ... last good ens
+	if ($dt > 120) {
+		printf(STDERR "\nWARNING: %d-s gap too long, profile restarted at ensemble $e...",$dt);
+		$firstgood = $lastgood = $e;			
+		$dt = $depth = 0;
+	}
+
+	$depth += $dta{ENSEMBLE}[$lastgood]->{REFW} * $dt		# integrate depth
+		if ($dt > 0);
+	$dta{ENSEMBLE}[$e]->{DEPTH} = $depth;
+	$atbottom = $e, $maxdepth = $depth if ($depth > $maxdepth);	
+
+	my($ss) = soundSpeed($opt_S,							# sound-speed corr
+				  	     $dta{ENSEMBLE}[$e]->{TEMPERATURE}-$opt_t,
+				   	     $depth);
+	$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} =
+		$ss / $dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND};
+	$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND} = $ss;
+	$dta{ENSEMBLE}[$e]->{REFW} *=
+		$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION};
+
+	$lastgood = $e;
+}
+
+printf(STDERR "done (max depth = %.1fm, depth at end of cast = %.1fm)\n",
+			$maxdepth,$depth);
+
+filterEnsembleStats() if defined($opt_F);
+
+#----------------------------------------------------------------------
+
+print(STDERR "Writing output...");
+
+if ($opt_A) {
+	print("#ANTS# [] $USAGE\n");
+	print("#ANTS#FIELDS# {ens} {unix_time} {time} {bin} {depth} {dz} {w} {ref_w} {dw} $addFields\n");
+	printf("#ANTS#PARAMS# maxdepth{$max_depth} bottom_time{%d}\n",
+	    $dta{ENSEMBLE}[$atbottom]->{UNIX_TIME} - 
+		 	$dta{ENSEMBLE}[$firstgood]->{UNIX_TIME});
+} else {
+	print("# ens-no time elapsed bin-no depth dz w ref-w dw $addFields\n");
+	print("#----------------------------------------------------------------------\n");
+}
+	
+for ($e=$firstgood; $e<=$lastgood; $e++) {
+
+	for ($i=$minb; $i<=$maxb; $i++) {						# dump valid
+		next unless defined($dta{ENSEMBLE}[$e]->{W}[$i]);
+
+		printf("%d %f %f %d %.1f %.1f %g %g %g ",
+			$e,$dta{ENSEMBLE}[$e]->{UNIX_TIME},
+			$dta{ENSEMBLE}[$e]->{UNIX_TIME} -
+				$dta{ENSEMBLE}[$firstgood]->{UNIX_TIME},$i,
+			$dta{ENSEMBLE}[$e]->{DEPTH} +
+				$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *
+					($dta{DISTANCE_TO_BIN1_CENTER} + $i*$dta{BIN_LENGTH}),
+			$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *
+				($dta{DISTANCE_TO_BIN1_CENTER} + $i*$dta{BIN_LENGTH}),
+			$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *
+				$dta{ENSEMBLE}[$e]->{W}[$i],
+			$dta{ENSEMBLE}[$e]->{REFW},
+			$dta{ENSEMBLE}[$e]->{SPEED_OF_SOUND_CORRECTION} *
+				$dta{ENSEMBLE}[$e]->{W}[$i] - $dta{ENSEMBLE}[$e]->{REFW},
+		);
+
+		sub p($) { print(defined($_[0])?"$_[0] ":"nan "); }
+		
+		if (defined(@f)) {
+			foreach $f (@f) {
+				my($fn,$fi) = ($f =~ m{([^[]*)(\[.*)});
+				$fn = $f unless defined($fn);
+				p(eval("\$dta{ENSEMBLE}[$e]->{$fn}$fi"));
+			}
+		}
+		print("\n");
+	}
+}
+
+print(STDERR "done\n");
+
+exit(0);	
+