mkProfile
changeset 0 229a0d72d2ab
child 3 f3c9dcbbdd68
equal deleted inserted replaced
-1:000000000000 0:229a0d72d2ab
       
     1 #!/usr/bin/perl
       
     2 #======================================================================
       
     3 #                    M K P R O F I L E 
       
     4 #                    doc: Sun Jan 19 18:55:26 2003
       
     5 #                    dlm: Sun May 23 16:34:02 2010
       
     6 #                    (c) 2003 A.M. Thurnherr
       
     7 #                    uE-Info: 788 36 NIL 0 0 72 2 2 4 NIL ofnI
       
     8 #======================================================================
       
     9 
       
    10 # Make an LADCP Profile by Integrating W (similar to Firing's scan*).
       
    11 
       
    12 # HISTORY:
       
    13 #	Jan 19, 2003: - written in order to test the RDI libs
       
    14 #	Jan 20, 2003: - added ensemble number
       
    15 #	Jan 21, 2003: - added horizontal integration
       
    16 #	Jan 22, 2003: - corrected magnetic declination
       
    17 #	Jan 23, 2003: - added -F)ilter
       
    18 #	Jan 24, 2003: - added more %PARAMs; started integration from 1st bin
       
    19 #				  - added -g, -f, battery status
       
    20 #	Jan 25, 2003: - added more %PARAMs
       
    21 #	Feb  1, 2003: - BUG: bottom-track quality checking was bad
       
    22 #	Feb  8, 2003: - allowed for array-indices on -f
       
    23 #	Feb  9, 2003: - added 50% goodvelbin
       
    24 #				  - removed unknown-field err on -f to allow -f W
       
    25 #	Feb 10, 2003: - changed initialization depth to 0m
       
    26 #				  - changed %bottom_depth to %max_depth
       
    27 #	Feb 11, 2003: - changed sign of magnetic declination
       
    28 #	Feb 12, 2003: - corrected BT-range scaling
       
    29 #	Feb 14, 2003: - added %pinging_hours, %min_range
       
    30 #				  - removed magnetic declination from default
       
    31 #	Feb 26, 2004: - added earth coordinates
       
    32 #	Mar  3, 2004: - removed requirement for -M on !-Q
       
    33 #				  - corrected range-stats on earth coordinates
       
    34 #	Mar  4, 2004: - added number of ensebles to output
       
    35 #	Mar 11, 2004: - BUG: rename ACD -> ADC
       
    36 #	Mar 12, 2004: - added %bottom_xmit_{current|voltage}
       
    37 #	Mar 16, 2004: - BUG: on -M u/v/x/y were wrong
       
    38 #	Mar 17, 2004: - added error estimates on u/v/x/y
       
    39 #				  - removed battery stuff (has to be done btw casts)
       
    40 #	Mar 18, 2004: - totally re-did u/v integration
       
    41 #	Mar 19, 2004: - re-designed u/v uncertainty estimation
       
    42 #	Mar 28, 2004: - added MEAN_CORRELATION, MEAN_ECHO_AMPLITUDE
       
    43 #	Sep 15, 2005: - changed BinRead library name
       
    44 #				  - made max gap length variable
       
    45 #	Sep 16, 2005: - re-did u,v,w uncertainties
       
    46 #	Nov  8, 2005: - UNIXTIME => UNIX_TIME
       
    47 #				  - added unix_time, secno, z_BT to default output
       
    48 #	Dec  1, 2005: - moved profile-building code to [RDI_utils.pl]
       
    49 #				  - changed -f syntax to allow name=FIELD
       
    50 #				  - added %bin1_dist, %bin_length
       
    51 #	Dec  8, 2005: - remove spaces from -f argument to allow multiline
       
    52 #				    definitions in Makefiles
       
    53 #	Nov 13, 2006: - BUG: end-of-cast depth had not been reported correctly
       
    54 #				  - cosmetics
       
    55 #	Nov 30, 2007: - adapted to 3-beam solutions
       
    56 #	Dec 11, 2007: - adapted to earlier modifications (Sep 2007) of
       
    57 #					[RDI_BB_Read.pl]
       
    58 #	Dec 14, 2007: - replaced z by depth
       
    59 #	Dec 17, 2007: - BUG: downcast flag was set incorrectly
       
    60 #	Jan 24, 2008: - rotation had been output as degrees/s; to make it more
       
    61 #				    consistent with pitch/roll, I changed it to simple degrees
       
    62 #				  - added net rotations [deployment]/down/up/[recovery]
       
    63 #	Apr  9, 2008: - added profile -B)ottom depth
       
    64 #				  - BUG: depth of first bin was reported as beginning of cast
       
    65 #	Oct 24, 2008: - added RANGE and RANGE_BINS fields
       
    66 #	Mar 18, 2009: - BUG: pitch/roll calculation had typo
       
    67 #				  - calc pitch/roll separately for down-/upcasts
       
    68 #				  - removed approximations in pitch/roll calcs
       
    69 #	Jul 30, 2009: - typo '<' removed from output
       
    70 #				  - NaN => nan
       
    71 
       
    72 # NOTES:
       
    73 #	- the battery values are based on transmission voltages (different
       
    74 #	  from battery voltages) and reported without units (raw 8-bit a2d
       
    75 #	  values)
       
    76 #	- -B with the CTD max depth can be used to linearly scale the depths;
       
    77 #	  even so, the profile can have negative depths, in particular when
       
    78 #	  the CTD is sent to a shallow depth first and then returned to the surface
       
    79 #	  before beginning the cast
       
    80 #	- in one case that I looked at (Anslope ][, cast 82), there are large
       
    81 #	  depth errors, even when -B is used
       
    82 #	- this utility works only approximately for uplookers (profile is
       
    83 #	  roughly ok, but apparently contaminated by surface reflection,
       
    84 #	  but stats are not ok; e.g. NBP0402 037U.prof)
       
    85 
       
    86 $0 =~ m{(.*)/[^/]+}; 
       
    87 require "$1/RDI_BB_Read.pl";
       
    88 require "$1/RDI_Coords.pl";
       
    89 require "$1/RDI_Utils.pl";
       
    90 require "getopts.pl";
       
    91 
       
    92 $USAGE = "$0 @ARGV";
       
    93 die("Usage: $0 " .
       
    94 	"[-A)nts] [-Q)uiet] [-F)ilter <script>] " .
       
    95 	"[-s)uppress checkensemble()] " .
       
    96 	"[require -4)-beam solutions] " .
       
    97 	"[-r)ef-layer <bin|1,bin|6>] [-n) vels <min|2>] " .
       
    98 	"[-e)rr-vel <max|0.1>] [-c)orrelation <min>] " .
       
    99 	"[-m)ax <gap>] " .
       
   100 	"[-d)rift <dx,dy>] [-g)ps <start lat,lon/end lat,lon>] " .
       
   101 	"[output -f)ields <field[,...]> " .
       
   102 	"[-M)agnetic <declination>] [profile -B)ottom <depth>] " .
       
   103 	"<RDI file>\n")
       
   104 		unless (&Getopts("4AB:F:M:Qd:r:n:e:c:g:f:m:s") && @ARGV == 1);
       
   105 
       
   106 die("$0: -Q and -A are mutually exclusive\n")
       
   107 	if ($opt_Q && $opt_A);
       
   108 
       
   109 $RDI_Coords::minValidVels = 4 if ($opt_4);			# no 3-beam solutions
       
   110 
       
   111 require $opt_F if defined($opt_F);					# load filter
       
   112 
       
   113 $opt_r = "1,6" 	unless defined($opt_r);				# defaults
       
   114 $opt_n = 2	   	unless defined($opt_n);
       
   115 $opt_e = 0.1   	unless defined($opt_e);
       
   116 $opt_c = 70	   	unless defined($opt_c);
       
   117 $opt_m = 120	unless defined($opt_m);
       
   118 
       
   119 ($minb,$maxb) = split(',',$opt_r);					# reference layer
       
   120 die("$0: can't decode -r $opt_r\n") unless defined($maxb);
       
   121                                         
       
   122 if ($opt_g) {										# GPS info
       
   123 	($s_lat,$s_lon,$e_lat,$e_lon) = gps_to_deg($opt_g);
       
   124 	$lat = $s_lat/2 + $e_lat/2;
       
   125 	$lon = $s_lon/2 + $e_lon/2;
       
   126 	$ddx = dist($lat,$s_lon,$lat,$e_lon);
       
   127 	$ddy = dist($s_lat,$lon,$e_lat,$lon);
       
   128 }
       
   129 
       
   130 if ($opt_d) {										# ship drift
       
   131 	($ddx,$ddy) = split(',',$opt_d);
       
   132 	die("$0: can't decode -d $opt_d\n") unless defined($ddy);
       
   133 }
       
   134 
       
   135 print(STDERR "Reading $ARGV[0]...");				# read data
       
   136 readData($ARGV[0],\%dta);
       
   137 print(STDERR "done\n");
       
   138 
       
   139 die("$ARGV[0]: not enough bins for choice of -r\n")	# enough bins?
       
   140 	unless ($dta{N_BINS} >= $maxb);
       
   141 if ($dta{BEAM_COORDINATES}) {						# coords used
       
   142 	$beamCoords = 1;
       
   143 } elsif (!$dta{EARTH_COORDINATES}) {
       
   144 	die("$ARGV[0]: only beam and earth coordinates implemented so far\n");
       
   145 }
       
   146 if (defined($opt_M)) {								# magnetic declination
       
   147 	$dta{HEADING_BIAS} = -1*$opt_M;
       
   148 } else {
       
   149 	$dta{HEADING_BIAS} = 0;
       
   150 }
       
   151 
       
   152 ensure_BT_RANGE(\%dta);								# calc if missing
       
   153 
       
   154 if ($opt_f) {										# additional fields
       
   155 	@f = split(',',$opt_f);
       
   156 	foreach $f (@f) {
       
   157 		$f =~ s/\s//g;								# remove spaces
       
   158 		@def = split('=',$f);
       
   159 		if (@def == 2) {							# name=field
       
   160 			$addFields .= " {$def[0]}";
       
   161 			$f = $def[1];
       
   162 		} else {									# field
       
   163 			$addFields .= " {$f}";
       
   164 		}
       
   165 	}
       
   166 #	print(STDERR "addFields = $addFields\n");
       
   167 #	print(STDERR "\@f = @f\n");
       
   168 }
       
   169 
       
   170 #======================================================================
       
   171 # Misc funs used to decode options
       
   172 #======================================================================
       
   173 
       
   174 sub dist($$$$)										# distance
       
   175 {
       
   176 	my($lat1,$lon1,$lat2,$lon2) = @_;
       
   177 	my($a) = 6378139; 								# Earth's radius
       
   178 
       
   179 	$lat1 = rad($lat1); $lon1 = rad($lon1);
       
   180 	$lat2 = rad($lat2); $lon2 = rad($lon2);
       
   181 	$ct1 = cos($lat1); $st1 = sin($lat1); $cp1 = cos($lon1);
       
   182 	$sp1 = sin($lon1); $ct2 = cos($lat2); $st2 = sin($lat2);
       
   183 	$cp2 = cos($lon2); $sp2 = sin($lon2);
       
   184 	$cos = $ct1*$cp1*$ct2*$cp2 + $ct1*$sp1*$ct2*$sp2 + $st1*$st2;
       
   185 	$cos =  1 if ($cos >  1);
       
   186 	$cos = -1 if ($cos < -1);
       
   187 	return $a * acos($cos);
       
   188 }
       
   189 
       
   190 sub deg_to_dec($)									# parse degrees
       
   191 {
       
   192 	my($deg,$min) = split(':',$_[0]);
       
   193 	return $deg + $min/60;
       
   194 }
       
   195 
       
   196 sub gps_to_deg($)									# decode lat/lon
       
   197 {
       
   198 	my($start,$end) = split('/',$_[0]);
       
   199 	my($sa,$so,$ea,$eo);
       
   200 
       
   201 	my($lat,$lon) = split(',',$start);
       
   202 	if ($lat =~ m{N$}) 		{ $sa =  deg_to_dec($`); }
       
   203 	elsif ($lat =~ m{S$}) 	{ $sa = -deg_to_dec($`); }
       
   204 	else 					{ $sa = $lat; }
       
   205 	if ($lon =~ m{E$}) 		{ $so =  deg_to_dec($`); }
       
   206 	elsif ($lon =~ m{W$}) 	{ $so = -deg_to_dec($`); }
       
   207 	else 					{ $so = $lon; }
       
   208 
       
   209 	my($lat,$lon) = split(',',$end);
       
   210 	if ($lat =~ m{N$}) 		{ $ea =  deg_to_dec($`); }
       
   211 	elsif ($lat =~ m{S$}) 	{ $ea = -deg_to_dec($`); }
       
   212 	else 					{ $ea = $lat; }
       
   213 	if ($lon =~ m{E$}) 		{ $eo =  deg_to_dec($`); }
       
   214 	elsif ($lon =~ m{W$}) 	{ $eo = -deg_to_dec($`); }
       
   215 	else 					{ $eo = $lon; }
       
   216 
       
   217 	return ($sa,$so,$ea,$eo);
       
   218 }
       
   219 
       
   220 #======================================================================
       
   221 # Step 1: Integrate w & determine water depth 
       
   222 #======================================================================
       
   223 
       
   224 ($firstgood,$lastgood,$atbottom,$w_gap_time,$zErr,$maxz) =
       
   225 	mk_prof(\%dta,!$opt_s,$opt_F,$minb,$maxb,$opt_c,$opt_e,$opt_m);
       
   226 
       
   227 die("$ARGV[0]: no good ensembles found\n")
       
   228 	unless defined($firstgood);
       
   229 
       
   230 if (defined($opt_B)) {										# scale Z
       
   231 	my($zscale) = $opt_B / ($dta{ENSEMBLE}[$atbottom]->{DEPTH} -# downcast
       
   232 				     	    $dta{ENSEMBLE}[$firstgood]->{DEPTH});
       
   233 #	printf(STDERR "scaling downcast depths by %.2f\n",$zscale);
       
   234 	for (my($e)=$firstgood; $e<$atbottom; $e++) {
       
   235 		next unless defined($dta{ENSEMBLE}[$e]->{DEPTH});
       
   236 		$dta{ENSEMBLE}[$e]->{DEPTH} =
       
   237 			$dta{ENSEMBLE}[$firstgood]->{DEPTH} + $zscale *
       
   238 				($dta{ENSEMBLE}[$e]->{DEPTH}-$dta{ENSEMBLE}[$firstgood]->{DEPTH});
       
   239 	}
       
   240 
       
   241 	$zscale = $opt_B / ($dta{ENSEMBLE}[$atbottom]->{DEPTH} -	# upcast
       
   242 				 		$dta{ENSEMBLE}[$lastgood]->{DEPTH});
       
   243 #	printf(STDERR "scaling upcast depths by %.2f\n",$zscale);
       
   244 	for (my($e)=$atbottom; $e<=$lastgood; $e++) {
       
   245 		next unless defined($dta{ENSEMBLE}[$e]->{DEPTH});
       
   246 		$dta{ENSEMBLE}[$e]->{DEPTH} =
       
   247 			$dta{ENSEMBLE}[$firstgood]->{DEPTH} + $zscale *
       
   248 				($dta{ENSEMBLE}[$e]->{DEPTH}-$dta{ENSEMBLE}[$lastgood]->{DEPTH});
       
   249 	}
       
   250 }
       
   251 
       
   252 ($water_depth,$sig_wd) =									# sea bed
       
   253 	find_seabed(\%dta,$atbottom,$beamCoords);
       
   254 
       
   255 #======================================================================
       
   256 # Step 1a: determine alternate Z by using mean/sigma of w in gaps
       
   257 #======================================================================
       
   258 
       
   259 # This does not make much sense for w, because w is always very close
       
   260 # to zero. It might make sense for u and v, though, and it would
       
   261 # be more consistent with the way the displacement uncertainties are
       
   262 # calculated. However, the way the profiles are calculated at the
       
   263 # moment (using the last valid velocity across the gap) is probably
       
   264 # closer to the truth in most cases.
       
   265 
       
   266 #$dta{ENSEMBLE}[$firstgood]->{ALT_Z} = 0;
       
   267 #$dta{ENSEMBLE}[$firstgood]->{ALT_Z_ERR} = 0;
       
   268 #my($sumVar);
       
   269 #for ($e=$firstgood+1; $e<=$lastgood; $e++) {
       
   270 #	my($dt) = $dta{ENSEMBLE}[$e]->{UNIX_TIME} -
       
   271 #			  $dta{ENSEMBLE}[$e-1]->{UNIX_TIME};
       
   272 #	$dta{ENSEMBLE}[$e]->{ALT_Z} =
       
   273 #		$dta{ENSEMBLE}[$e-1]->{ALT_Z} +
       
   274 #		$dt * (defined($dta{ENSEMBLE}[$e-1]->{W}) ?
       
   275 #				$dta{ENSEMBLE}[$e-1]->{W} : $meanW);
       
   276 #	$sumVar += defined($dta{ENSEMBLE}[$e-1]->{W}) ?
       
   277 #				($dta{ENSEMBLE}[$e-1]->{W_ERR} * $dt)**2 : ($dt**2)*$varW;
       
   278 #	$dta{ENSEMBLE}[$e]->{ALT_Z_ERR} = sqrt($sumVar);
       
   279 #}
       
   280 
       
   281 #======================================================================
       
   282 # Step 2: Integrate u & v
       
   283 #======================================================================
       
   284 
       
   285 sub ref_lr_uv($$$)									# calc ref-level u/v
       
   286 {
       
   287 	my($ens,$z,$water_depth) = @_;
       
   288 	my($i,$n,@v,@goodU,@goodV);
       
   289 
       
   290 	$water_depth = 99999 unless defined($water_depth);
       
   291 
       
   292 	for ($i=$minb; $i<=$maxb; $i++) {
       
   293 		next if ($dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][0] < $opt_c ||
       
   294 				 $dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][1] < $opt_c ||
       
   295 				 $dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][2] < $opt_c ||
       
   296 				 $dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][3] < $opt_c);
       
   297 		if ($beamCoords) {
       
   298 			next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] < 100 ||
       
   299 					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] < 100 ||
       
   300 					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] < 100 ||
       
   301 					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);
       
   302 			@v = velInstrumentToEarth(\%dta,$ens,
       
   303 					velBeamToInstrument(\%dta,
       
   304 						@{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$i]}));
       
   305 		} else {
       
   306 			next if ($dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][0] > 0 ||
       
   307 					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][1] > 0 ||
       
   308 					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][2] > 0 ||
       
   309 					 $dta{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$i][3] < 100);
       
   310 			@v = velApplyHdgBias(\%dta,$ens,
       
   311 					@{$dta{ENSEMBLE}[$ens]->{VELOCITY}[$i]});
       
   312 		}
       
   313 		next if (!defined($v[3]) || abs($v[3]) > $opt_e);
       
   314 
       
   315 #		Martin's BT routines show strong shear just above sea bed
       
   316 #		=> skip lowest 20m.
       
   317 		if (defined($v[0])) {							# valid u,v
       
   318 			if ($dta{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}) {
       
   319 				if ($z - $dta{DISTANCE_TO_BIN1_CENTER}
       
   320 					   - $i*$dta{BIN_LENGTH} > 0) {
       
   321 					push(@goodU,$v[0]); push(@goodV,$v[1]);
       
   322 					$dta{ENSEMBLE}[$ens]->{U} += $v[0];
       
   323 					$dta{ENSEMBLE}[$ens]->{V} += $v[1];
       
   324 					$n++;
       
   325 			   }
       
   326 			} else {
       
   327 				if ($z + $dta{DISTANCE_TO_BIN1_CENTER}
       
   328 					   + $i*$dta{BIN_LENGTH} < $water_depth-20) {
       
   329 					push(@goodU,$v[0]); push(@goodV,$v[1]);
       
   330 					$dta{ENSEMBLE}[$ens]->{U} += $v[0];
       
   331 					$dta{ENSEMBLE}[$ens]->{V} += $v[1];
       
   332 					$n++;
       
   333 			   }
       
   334 			}
       
   335 		}
       
   336 	}
       
   337 
       
   338 	if ($n >= 2) {
       
   339 		my(@sumsq) = (0,0);
       
   340 		$dta{ENSEMBLE}[$ens]->{U} /= $n;
       
   341 		$dta{ENSEMBLE}[$ens]->{V} /= $n;
       
   342 		for ($i=0; $i<$n; $i++) {
       
   343 			$sumsq[0] += ($dta{ENSEMBLE}[$ens]->{U}-$goodU[$i])**2;
       
   344 			$sumsq[1] += ($dta{ENSEMBLE}[$ens]->{V}-$goodV[$i])**2;
       
   345 		}
       
   346 		$dta{ENSEMBLE}[$ens]->{U_ERR} = sqrt($sumsq[0])/($n-1);
       
   347 		$dta{ENSEMBLE}[$ens]->{V_ERR} = sqrt($sumsq[1])/($n-1);
       
   348     } else {
       
   349 		$dta{ENSEMBLE}[$ens]->{U} = undef;
       
   350 		$dta{ENSEMBLE}[$ens]->{V} = undef;
       
   351     }    	
       
   352 }
       
   353 
       
   354 #----------------------------------------------------------------------
       
   355 
       
   356 ($x,$y) = (0,0);											# init
       
   357 
       
   358 $dta{ENSEMBLE}[$firstgood]->{X} = $dta{ENSEMBLE}[$firstgood]->{X_ERR} = 0;
       
   359 $dta{ENSEMBLE}[$firstgood]->{Y} = $dta{ENSEMBLE}[$firstgood]->{Y_ERR} = 0;
       
   360 $prevgood = $firstgood;
       
   361 
       
   362 for ($e=$firstgood+1; defined($opt_M)&&$e<=$lastgood; $e++) {
       
   363 
       
   364 	#--------------------------------------------------
       
   365 	# within profile: both $firstgood and $prevgood set
       
   366 	#--------------------------------------------------
       
   367 
       
   368 	ref_lr_uv($e,$dta{ENSEMBLE}[$e]->{DEPTH},$water_depth)	# instrument vel
       
   369 		if (defined($dta{ENSEMBLE}[$e]->{W}));
       
   370 
       
   371 	if (!defined($dta{ENSEMBLE}[$e]->{U})) {				# gap
       
   372 		$uv_gap_time += $dta{ENSEMBLE}[$e]->{UNIX_TIME} -
       
   373 				   		$dta{ENSEMBLE}[$e-1]->{UNIX_TIME};
       
   374 		next;
       
   375 	}
       
   376 
       
   377 	my($dt) = $dta{ENSEMBLE}[$e]->{UNIX_TIME} -			# time step since
       
   378 			  $dta{ENSEMBLE}[$prevgood]->{UNIX_TIME};		# ...last good ens
       
   379 
       
   380 	#-----------------------------------
       
   381 	# The current ensemble has valid u/v
       
   382 	#-----------------------------------
       
   383 
       
   384 	$x -= $dta{ENSEMBLE}[$prevgood]->{U} * $dt;			# integrate
       
   385 	$xErr += ($dta{ENSEMBLE}[$prevgood]->{U_ERR} * $dt)**2;
       
   386 	$dta{ENSEMBLE}[$e]->{X} = $x;
       
   387 	$dta{ENSEMBLE}[$e]->{X_ERR} = sqrt($xErr);
       
   388 
       
   389 	$y -= $dta{ENSEMBLE}[$prevgood]->{V} * $dt;
       
   390 	$yErr += ($dta{ENSEMBLE}[$prevgood]->{V_ERR} * $dt)**2;
       
   391 	$dta{ENSEMBLE}[$e]->{Y} = $y;
       
   392 	$dta{ENSEMBLE}[$e]->{Y_ERR} = sqrt($yErr);
       
   393 
       
   394 	$prevgood = $e;
       
   395 }
       
   396 
       
   397 unless (defined($dta{ENSEMBLE}[$lastgood]->{X})) {		# last is bad in u/v
       
   398 	my($dt) = $dta{ENSEMBLE}[$lastgood]->{UNIX_TIME} -		# time step since
       
   399 			  $dta{ENSEMBLE}[$prevgood]->{UNIX_TIME};		# ...last good ens
       
   400 
       
   401 	$x -= $dta{ENSEMBLE}[$prevgood]->{U} * $dt;			# integrate
       
   402 	$xErr += ($dta{ENSEMBLE}[$prevgood]->{U_ERR} * $dt)**2;
       
   403 	$dta{ENSEMBLE}[$lastgood]->{X} = $x;
       
   404 	$dta{ENSEMBLE}[$lastgood]->{X_ERR} = sqrt($xErr);
       
   405 
       
   406 	$y -= $dta{ENSEMBLE}[$prevgood]->{V} * $dt;
       
   407 	$yErr += ($dta{ENSEMBLE}[$prevgood]->{V_ERR} * $dt)**2;
       
   408 	$dta{ENSEMBLE}[$lastgood]->{Y} = $y;
       
   409 	$dta{ENSEMBLE}[$lastgood]->{Y_ERR} = sqrt($yErr);
       
   410 }
       
   411 
       
   412 $firstgood++ if ($firstgood == 0);							# centered diff
       
   413 $lastgood-- if ($lastgood == $#{$dta{ENSEMBLE}});			# in step 6
       
   414 
       
   415 #======================================================================
       
   416 # Step 3: Calculate Uncertainties
       
   417 #======================================================================
       
   418 
       
   419 # Time series of W_ERR indicate that errors are very large near the
       
   420 # surface and near the sea bed, perhaps because of reflections.
       
   421 # A reasonable estimate for typical uncertainty is therefore the mode
       
   422 # of the std errors.
       
   423 
       
   424 my(@histUErr,@histVErr,@histWErr);
       
   425 my($histRez) = 1e-4;
       
   426 
       
   427 for ($e=$firstgood; $e<=$lastgood; $e++) {
       
   428 	$histWErr[int($dta{ENSEMBLE}[$e]->{W_ERR}/$histRez+0.5)]++
       
   429 		if defined($dta{ENSEMBLE}[$e]->{W_ERR});
       
   430 	$histUErr[int($dta{ENSEMBLE}[$e]->{U_ERR}/$histRez+0.5)]++
       
   431 		if defined($dta{ENSEMBLE}[$e]->{U_ERR});
       
   432 	$histVErr[int($dta{ENSEMBLE}[$e]->{V_ERR}/$histRez+0.5)]++
       
   433 		if defined($dta{ENSEMBLE}[$e]->{V_ERR});
       
   434 }	
       
   435 
       
   436 my($max) = 0; my($mode);
       
   437 for (my($i)=0; $i<=$#histWErr; $i++) {
       
   438 	next if ($histWErr[$i] < $max);
       
   439 	$max = $histWErr[$i]; $mode = $i;
       
   440 }
       
   441 $wErr = $mode * $histRez if defined($mode);
       
   442 
       
   443 $max = 0; $mode = undef;
       
   444 for (my($i)=0; $i<=$#histUErr; $i++) {
       
   445 	next if ($histUErr[$i] < $max);
       
   446 	$max = $histUErr[$i]; $mode = $i;
       
   447 }
       
   448 $uErr = $mode * $histRez if defined($mode);
       
   449 
       
   450 $max = 0; $mode = undef;
       
   451 for (my($i)=0; $i<=$#histVErr; $i++) {
       
   452 	next if ($histVErr[$i] < $max);
       
   453 	$max = $histVErr[$i]; $mode = $i;
       
   454 }
       
   455 $vErr = $mode * $histRez if defined($mode);
       
   456 
       
   457 #print(STDERR "u: mu = $meanU / sigma = $uErr\n");
       
   458 #print(STDERR "v: mu = $meanV / sigma = $vErr\n");
       
   459 #print(STDERR "w: mu = $meanW / sigma = $wErr\n");
       
   460 
       
   461 if (defined($opt_M)) {									# displacement errors
       
   462 	$x_err = $uErr * $uv_gap_time + $dta{ENSEMBLE}[$lastgood]->{X_ERR};
       
   463 	$y_err = $vErr * $uv_gap_time + $dta{ENSEMBLE}[$lastgood]->{Y_ERR};
       
   464 }
       
   465 $z_err = $wErr * $w_gap_time + $dta{ENSEMBLE}[$lastgood]->{DEPTH_ERR};
       
   466 
       
   467 #printf(STDERR "x_err = $dta{ENSEMBLE}[$lastgood]->{X_ERR} + %g\n",
       
   468 #				$uErr * $uv_gap_time);
       
   469 #printf(STDERR "y_err = $dta{ENSEMBLE}[$lastgood]->{Y_ERR} + %g\n",
       
   470 #				$vErr * $uv_gap_time);
       
   471 #printf(STDERR "z_err = $dta{ENSEMBLE}[$lastgood]->{DEPTH_ERR} + %g\n",
       
   472 #				$wErr * $w_gap_time);
       
   473 
       
   474 #======================================================================
       
   475 # Step 4: Calculate Beam Range Stats
       
   476 #======================================================================
       
   477 
       
   478 my($min_good_bins) = 999;
       
   479 my($worst_beam);
       
   480 
       
   481 sub count_good_vels($)								# count good vels
       
   482 {
       
   483 	my($ens) = @_;
       
   484 	my($good) = -1; my($this_worst_beam);
       
   485 
       
   486 	if ($beamCoords) {
       
   487 		for (my($i)=0; $i<$dta{N_BINS}; $i++) {
       
   488 			for (my($b)=0; $b<4; $b++) {
       
   489 				$good=$i,$this_worst_beam=$b,$nVels[$i][$b]++
       
   490 					if defined($dta{ENSEMBLE}[$ens]->{VELOCITY}[$i][$b]);
       
   491 			}
       
   492 	    }
       
   493 	} else {
       
   494 		for (my($i)=0; $i<$dta{N_BINS}; $i++) {
       
   495 			for (my($b)=0; $b<4; $b++) {
       
   496 				$good=$i,$this_worst_beam=$b,$nVels[$i][$b]++
       
   497 					if ($dta{ENSEMBLE}[$ens]->{CORRELATION}[$i][$b] >=
       
   498 						$dta{MIN_CORRELATION});
       
   499 			}
       
   500 	    }
       
   501 	}
       
   502 	$min_good_ens=$ens, $min_good_bins=$good, $worst_beam=$this_worst_beam
       
   503 		if ((!defined($water_depth) || 
       
   504 			  $dta{ENSEMBLE}[$ens]->{DEPTH} < $water_depth-200)
       
   505 			&& $good >= 0 && $good < $min_good_bins);
       
   506 }
       
   507 
       
   508 #----------------------------------------------------------------------
       
   509 
       
   510 for ($e=$firstgood; $e<=$lastgood; $e++) {					# range
       
   511 	my($i);
       
   512 	for ($i=0; $i<$dta{N_BINS}; $i++) {
       
   513 		last if (defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$i][0]) +
       
   514 				 defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$i][1]) +
       
   515 				 defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$i][2]) +
       
   516 				 defined($dta{ENSEMBLE}[$e]->{VELOCITY}[$i][3]) < 3);
       
   517 	}
       
   518 	$dta{ENSEMBLE}[$e]->{RANGE_BINS} = $i;
       
   519 	$dta{ENSEMBLE}[$e]->{RANGE} =
       
   520 		$dta{DISTANCE_TO_BIN1_CENTER} + $i * $dta{BIN_LENGTH};
       
   521 }
       
   522 
       
   523 for ($e=$firstgood; $e<=$lastgood; $e++) {					# mean corr/amp
       
   524 	$sumcor = $sumamp = $ndata = 0;
       
   525 	for (my($i)=0; $i<$dta{N_BINS}; $i++) {
       
   526 		for (my($b)=0; $b<4; $b++) {
       
   527 			next unless ($dta{ENSEMBLE}[$e]->{CORRELATION}[$i][$b]);
       
   528 			$sumcor += $dta{ENSEMBLE}[$e]->{CORRELATION}[$i][$b];
       
   529 			$sumamp += $dta{ENSEMBLE}[$e]->{ECHO_AMPLITUDE}[$i][$b];
       
   530 			$ndata++;
       
   531 		}
       
   532 	}
       
   533 	$dta{ENSEMBLE}[$e]->{MEAN_CORRELATION} = $sumcor/$ndata;
       
   534 	$dta{ENSEMBLE}[$e]->{MEAN_ECHO_AMPLITUDE} = $sumamp/$ndata;
       
   535 }
       
   536 
       
   537 for ($e=$firstgood+50; $e<=$lastgood-50; $e++) {			# range stats
       
   538 	count_good_vels($e);
       
   539 }
       
   540 for ($i=0; $i<$dta{N_BINS}; $i++) {
       
   541 	for ($b=0; $b<4; $b++) {
       
   542 		$maxVels = $nVels[$i][$b] unless ($maxVels > $nVels[$i][$b]);
       
   543 	}
       
   544 }
       
   545 for ($i=0; $i<$dta{N_BINS}; $i++) {
       
   546 	for ($b=0; $b<4; $b++) {
       
   547 		$gb[$b] = $i if ($nVels[$i][$b] >= 0.8*$maxVels);
       
   548     }
       
   549 }
       
   550 $gb = ($gb[0]+$gb[1]+$gb[2]+$gb[3]) / 4;
       
   551 
       
   552 #======================================================================
       
   553 # Step 5: Remove Ship Drift (probably not useful)
       
   554 #======================================================================
       
   555 
       
   556 if (defined($opt_M) && defined($ddx)) {						# remove barotropic
       
   557 	$du = $ddx / $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};# mean drift vel
       
   558 	$dv = $ddy / $dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};
       
   559 	$iu = $dta{ENSEMBLE}[$lastgood]->{X} /				# mean obs vel
       
   560 			$dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};
       
   561 	$iv = $dta{ENSEMBLE}[$lastgood]->{Y} /
       
   562 			$dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME};
       
   563 	
       
   564 	for ($e=$firstgood; $e<=$lastgood; $e++) {
       
   565 		next unless (defined($dta{ENSEMBLE}[$e]->{X}) &&
       
   566 					 defined($dta{ENSEMBLE}[$e]->{Y}));
       
   567 		$dta{ENSEMBLE}[$e]->{U} -= $du;
       
   568 		$dta{ENSEMBLE}[$e]->{V} -= $dv;
       
   569 		$dta{ENSEMBLE}[$e]->{X} += $dta{ENSEMBLE}[$e]->{ELAPSED_TIME} * ($du-$iu);
       
   570 		$dta{ENSEMBLE}[$e]->{Y} += $dta{ENSEMBLE}[$e]->{ELAPSED_TIME} * ($dv-$iv);
       
   571 	}
       
   572 }
       
   573 
       
   574 #======================================================================
       
   575 # Step 6: Pitch, Roll, Rotation
       
   576 #======================================================================
       
   577 
       
   578 my($prrms,$dnprrms,$upprrms) = (0,0,0);
       
   579 my($rotrms,$prerot,$dnrot,$uprot,$postrot) = (0,0,0,0,0);
       
   580 
       
   581 sub rot($)
       
   582 {
       
   583 	my($e) = @_;
       
   584 	my($rot) = $dta{ENSEMBLE}[$e]->{HEADING} -
       
   585 			   $dta{ENSEMBLE}[$e-1]->{HEADING};
       
   586 	$rot -= 360 if ($rot >  180);
       
   587 	$rot += 360 if ($rot < -180);
       
   588 	return $rot;
       
   589 }
       
   590 
       
   591 for ($e=1; $e<$firstgood; $e++) {				# pre-deployment
       
   592 	$prerot += rot($e);
       
   593 }
       
   594 
       
   595 for (; $e<= $atbottom; $e++) {					# downcast
       
   596 	$dta{ENSEMBLE}[$e]->{PITCHROLL} =
       
   597 		&angle_from_vertical($dta{ENSEMBLE}[$e]->{PITCH},
       
   598 						  	 $dta{ENSEMBLE}[$e]->{ROLL});
       
   599 	$prrms += $dta{ENSEMBLE}[$e]->{PITCHROLL}**2;
       
   600 				 
       
   601 	$dta{ENSEMBLE}[$e]->{ROTATION} = rot($e);
       
   602 	$dnrot += $dta{ENSEMBLE}[$e]->{ROTATION};
       
   603 	$rotrms += $dta{ENSEMBLE}[$e]->{ROTATION}**2;
       
   604 }
       
   605 $dnprrms = $prrms;
       
   606 
       
   607 for (; $e<=$lastgood; $e++) {					# upcast
       
   608 	$dta{ENSEMBLE}[$e]->{PITCHROLL} =
       
   609 		&angle_from_vertical($dta{ENSEMBLE}[$e]->{PITCH},
       
   610 						  	 $dta{ENSEMBLE}[$e]->{ROLL});
       
   611 	$prrms += $dta{ENSEMBLE}[$e]->{PITCHROLL}**2;
       
   612 				 
       
   613 	$dta{ENSEMBLE}[$e]->{ROTATION} = rot($e);
       
   614 	$uprot += $dta{ENSEMBLE}[$e]->{ROTATION};
       
   615 	$rotrms += $dta{ENSEMBLE}[$e]->{ROTATION}**2;
       
   616 }
       
   617 $upprrms = $prrms - $dnprrms;
       
   618 
       
   619 for (; $e<=$#{$dta->{ENSEMBLE}}; $e++) {		# post-recovery
       
   620 	$postrot += rot($e);
       
   621 }
       
   622 
       
   623 $prerot  /= 360;								# rotations, not degrees
       
   624 $dnrot   /= 360;
       
   625 $uprot   /= 360;
       
   626 $postrot /= 360;
       
   627 
       
   628 $prrms 	 = sqrt($prrms/($lastgood-$firstgood));
       
   629 $dnprrms = sqrt($dnprrms/($atbottom-$firstgood));
       
   630 $upprrms = sqrt($upprrms/($lastgood-$atbottom));
       
   631 
       
   632 $rotrms = sqrt($rotrms/($lastgood-$firstgood));
       
   633 
       
   634 #======================================================================
       
   635 # PRODUCE OUTPUT
       
   636 #======================================================================
       
   637 
       
   638 printf(STDERR "# of ensembles   : %d\n",scalar(@{$dta{ENSEMBLE}}));
       
   639 printf(STDERR "Start of cast    : %s (#%5d) at %6.1fm\n",
       
   640 					$dta{ENSEMBLE}[$firstgood]->{TIME},
       
   641 					$dta{ENSEMBLE}[$firstgood]->{NUMBER},
       
   642 					$dta{ENSEMBLE}[$firstgood]->{DEPTH});
       
   643 printf(STDERR "Bottom of cast   : %s (#%5d) at %6.1fm\n",
       
   644 					$dta{ENSEMBLE}[$atbottom]->{TIME},
       
   645 					$dta{ENSEMBLE}[$atbottom]->{NUMBER},
       
   646 					$dta{ENSEMBLE}[$atbottom]->{DEPTH});
       
   647 if (defined($water_depth)) {
       
   648 	printf(STDERR "Seabed           :                      at %6.1fm (+-%dm)\n",$water_depth,$sig_wd);
       
   649 } else {
       
   650 	print(STDERR "Seabed           : not found\n");
       
   651 }
       
   652 printf(STDERR "End of cast      : %s (#%5d) at %6.1fm\n",
       
   653 					$dta{ENSEMBLE}[$lastgood]->{TIME},
       
   654 					$dta{ENSEMBLE}[$lastgood]->{NUMBER},
       
   655 					$dta{ENSEMBLE}[$lastgood]->{DEPTH});
       
   656 
       
   657 printf(STDERR "Rel. Displacement: x = %d(%d)m / y = %d(%d)m\n",
       
   658 					$dta{ENSEMBLE}[$lastgood]->{X}, $x_err, 
       
   659 					$dta{ENSEMBLE}[$lastgood]->{Y}, $y_err, 
       
   660 				) if defined($opt_M);
       
   661 
       
   662 printf(STDERR "Cast Duration    : %.1f hours (pinging for %.1f hours)\n",
       
   663 					$dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME} / 3600,
       
   664 					($dta{ENSEMBLE}[$#{$dta{ENSEMBLE}}]->{UNIX_TIME} -
       
   665 						$dta{ENSEMBLE}[0]->{UNIX_TIME}) / 3600);
       
   666 
       
   667 printf(STDERR "Minimum range    : %dm at ensemble %d, beam %d\n",
       
   668 				$dta{DISTANCE_TO_BIN1_CENTER} +
       
   669 					$min_good_bins*$dta{BIN_LENGTH},
       
   670 				$dta{ENSEMBLE}[$min_good_ens]->{NUMBER},
       
   671 				$worst_beam);
       
   672 printf(STDERR "80%%-valid bins   : %.1f\n",$gb+1);
       
   673 printf(STDERR "80%%-valid range  : %dm\n",
       
   674 				$dta{DISTANCE_TO_BIN1_CENTER} + $gb*$dta{BIN_LENGTH});
       
   675 printf(STDERR "3-beam solutions : $RDI_Coords::threeBeam_1 " .
       
   676 								 "$RDI_Coords::threeBeam_2 " .
       
   677 								 "$RDI_Coords::threeBeam_3 " .
       
   678                                  "$RDI_Coords::threeBeam_4\n")
       
   679 	unless ($opt_4);
       
   680 printf(STDERR "net rotations    : [%d]/%d/%d/[%d]\n",$prerot,$dnrot,$uprot,$postrot);
       
   681 printf(STDERR "rms pitch/roll   : %.1f/%.1f\n",$dnprrms,$upprrms);
       
   682 
       
   683 if ($opt_A) {												# ANTS format
       
   684 	print("#ANTS# [] $USAGE\n");
       
   685 	$uFields = "{u} {u_err} {v} {v_err} {x} {x_err} {y} {y_err}"
       
   686 		if defined($opt_M);
       
   687 	print("#ANTS#FIELDS# {ens} {time} {elapsed} {secno} {downcast} " .
       
   688 						"{w} {w_err} {depth} {depth_err} {depth_BT} " .
       
   689 						"{pitchroll} {rotation} " .
       
   690 		  				"$uFields $addFields\n");
       
   691 
       
   692    	printf("#ANTS#PARAMS# date{$dta{ENSEMBLE}[$firstgood]->{DATE}} " .
       
   693 				   "start_time{$dta{ENSEMBLE}[$firstgood]->{TIME}} " .
       
   694 				  "bottom_time{$dta{ENSEMBLE}[$atbottom]->{TIME}} " .
       
   695 				     "end_time{$dta{ENSEMBLE}[$lastgood]->{TIME}} " .
       
   696 		  "bottom_xmit_voltage{$dta{ENSEMBLE}[$atbottom]->{ADC_XMIT_VOLTAGE}} " .
       
   697 		  "bottom_xmit_current{$dta{ENSEMBLE}[$atbottom]->{ADC_XMIT_CURRENT}} " .
       
   698 			 "pinging_duration{%.1f} " .
       
   699 			    "cast_duration{%.1f} " .
       
   700 	    	   "0.8_valid_bins{%.1f} " .
       
   701 	     	  "0.8_valid_range{%.1f} " .
       
   702 				    "max_depth{%.1f} " .
       
   703 				  "depth_error{%.1f} " .
       
   704 				    "min_range{%d} " .
       
   705 				  "n_ensembles{%d} " .
       
   706 				   "w_gap_time{%d} " .
       
   707 				     "stderr_w{%.4f} " .
       
   708 				"rms_pitchroll{%.1f} " .
       
   709   	   "downcast_rms_pitchroll{%.1f} " .
       
   710 		 "upcast_rms_pitchroll{%.1f} " .
       
   711 				 "rms_rotation{%.2f} " .
       
   712 		 "deployment_rotations{%d} " .
       
   713 		   "downcast_rotations{%d} " .
       
   714 		     "upcast_rotations{%d} " .
       
   715 		   "recovery_rotations{%d} " .
       
   716 				    "bin1_dist{%.1f} " .
       
   717 				   "bin_length{%.1f} " .
       
   718 		  				 "\n",
       
   719 				($dta{ENSEMBLE}[$#{$dta{ENSEMBLE}}]->{UNIX_TIME} -
       
   720 						$dta{ENSEMBLE}[0]->{UNIX_TIME}),
       
   721 				$dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME},
       
   722 				$gb+1,
       
   723 				$dta{DISTANCE_TO_BIN1_CENTER} + $gb*$dta{BIN_LENGTH},
       
   724 				$dta{ENSEMBLE}[$atbottom]->{DEPTH},
       
   725 				$dta{ENSEMBLE}[$lastgood]->{DEPTH} -
       
   726 					$dta{ENSEMBLE}[$firstgood]->{DEPTH},
       
   727 				$dta{DISTANCE_TO_BIN1_CENTER} +
       
   728 					$min_good_bins*$dta{BIN_LENGTH},
       
   729 				scalar(@{$dta{ENSEMBLE}}),
       
   730 				$w_gap_time,$wErr,$prrms,$dnprrms,$upprrms,$rotrms,
       
   731 				$prerot,$dnrot,$uprot,$postrot,
       
   732 				$dta{DISTANCE_TO_BIN1_CENTER},
       
   733 				$dta{BIN_LENGTH},
       
   734 		  );
       
   735 	printf("#ANTS#PARAMS# magnetic_declination{$opt_M} " .
       
   736 							   	  "uv_gap_time{%d} " .
       
   737 						    	       "mean_u{%.4f} " .
       
   738 							         "stderr_u{%.4f} " .
       
   739 							               "dx{%d} " .
       
   740 							           "dx_err{%d} " .
       
   741 							           "mean_v{%.4f} " .
       
   742 							         "stderr_v{%.4f} " .
       
   743 							               "dy{%d} " .
       
   744 							           "dy_err{%d}\n",
       
   745 		$uv_gap_time,
       
   746 		$dta{ENSEMBLE}[$lastgood]->{X} /
       
   747 			$dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME},
       
   748 		$uErr, $dta{ENSEMBLE}[$lastgood]->{X}, $x_err,
       
   749 		$dta{ENSEMBLE}[$lastgood]->{Y} /
       
   750 			$dta{ENSEMBLE}[$lastgood]->{ELAPSED_TIME},
       
   751 		$vErr, $dta{ENSEMBLE}[$lastgood]->{Y}, $y_err,
       
   752 	) if defined ($opt_M);
       
   753 	print("#ANTS#PARAMS# start_lat{$s_lat} start_lon{$s_lon} " .
       
   754 						  "end_lat{$e_lat} end_lon{$e_lon} " .
       
   755 						      "lat{$lat} lon{$lon}\n")
       
   756 		if defined($lat);
       
   757 	if ($dta{TIME_BETWEEN_PINGS} == 0) {
       
   758 		 print("#ANTS#PARAMS# pinging_rate{staggered}\n");
       
   759 	} else {
       
   760 		 printf("#ANTS#PARAMS# pinging_rate{%.2f}\n",
       
   761 			1/$dta{TIME_BETWEEN_PINGS});
       
   762     }		
       
   763     printf("#ANTS#PARAMS# drift_x{%d} drift_y{%d} " .
       
   764 						 "drift_u{%.3f} drift_v{%.3f} " .
       
   765 		   "\n",$ddx,$ddy,$du,$dv) if defined($ddx);
       
   766 	if (defined($water_depth)) {
       
   767 		printf("#ANTS#PARAMS# water_depth{%d} sig-water_depth{%d}\n",
       
   768 					$water_depth,$sig_wd);
       
   769 	} else {
       
   770 		print("#ANTS#PARAMS# water_depth{nan} sig-water_depth{nan}\n");
       
   771 	}
       
   772 }
       
   773 
       
   774 sub p($) { print(defined($_[0])?"$_[0] ":"nan "); }
       
   775 sub pb($) { print($_[0]?"1 ":"0 "); }
       
   776 
       
   777 unless ($opt_Q) {										# write profile
       
   778 	for ($e=$firstgood; $e<=$lastgood; $e++) {
       
   779 		p($dta{ENSEMBLE}[$e]->{NUMBER});
       
   780 		p($dta{ENSEMBLE}[$e]->{UNIX_TIME});
       
   781 		p($dta{ENSEMBLE}[$e]->{ELAPSED_TIME});
       
   782 		p($dta{ENSEMBLE}[$e]->{SECNO});
       
   783 		pb($dta{ENSEMBLE}[$e]->{UNIX_TIME} < $dta{ENSEMBLE}[$atbottom]->{UNIX_TIME});
       
   784 		p($dta{ENSEMBLE}[$e]->{W});
       
   785 		p($dta{ENSEMBLE}[$e]->{W_ERR});
       
   786 		p($dta{ENSEMBLE}[$e]->{DEPTH});
       
   787 		p($dta{ENSEMBLE}[$e]->{DEPTH_ERR});
       
   788 		p($dta{ENSEMBLE}[$e]->{DEPTH_BT});
       
   789 		p($dta{ENSEMBLE}[$e]->{PITCHROLL});
       
   790 		p($dta{ENSEMBLE}[$e]->{ROTATION});
       
   791 		if (defined($opt_M)) {
       
   792 			p($dta{ENSEMBLE}[$e]->{U}); p($dta{ENSEMBLE}[$e]->{U_ERR});
       
   793 			p($dta{ENSEMBLE}[$e]->{V}); p($dta{ENSEMBLE}[$e]->{V_ERR});
       
   794 			p($dta{ENSEMBLE}[$e]->{X}); p($dta{ENSEMBLE}[$e]->{X_ERR});
       
   795 			p($dta{ENSEMBLE}[$e]->{Y}); p($dta{ENSEMBLE}[$e]->{Y_ERR});
       
   796 	    }
       
   797 		if (defined(@f)) {
       
   798 			foreach $f (@f) {
       
   799 				my($fn,$fi) = ($f =~ m{([^[]*)(\[.*)});
       
   800 				$fn = $f unless defined($fn);
       
   801 				p(eval("\$dta{ENSEMBLE}[$e]->{$fn}$fi"));
       
   802 			}
       
   803 		}
       
   804 		print("\n");
       
   805 	}
       
   806 }
       
   807 
       
   808 exit(0);