#!/usr/bin/perl
#======================================================================
#                    L A D C P P R O C 
#                    doc: Thu Sep 16 20:36:10 2010
#                    dlm: Wed Jun 15 15:47:25 2011
#                    (c) 2010 A.M. Thurnherr & E. Firing
#                    uE-Info: 479 55 NIL 0 0 72 2 2 4 NIL ofnI
#======================================================================

$antsSummary = 'process LADCP data to get shear, time series';

# NOTES:
#	- this code is based on merge.c written by Eric Firing
#	- comments starting with ## are taken from Eric's code
#	- for SeaBird files, CTD elapsed time is estimated from recno * CTD{sampint}
#	- CTD{elapsed} is undefined for records before instrument is in the water
#	- ITS-90 temp field in degC required
#	- salin field prequired
#	- pressure field in dbar required
#	- -i should be set to the number that's added to LADCP_elapsed to make the two
#	  time series overplot nicely

# HISTORY:
#	Sep 16, 2010: - incepted
#	Oct 13, 2010: - first working version
#	Oct 14, 2010: - renamed from LADCPshear
#	Oct 19, 2010: - added -a)coustic backscatter profiles
#	Oct 20, 2010: - added -2)dary CTD sensors
#	Oct 23, 2010: - added magnetic-declination correction
#	Oct 26, 2010: - added tilt calculation
#	Dec  9, 2010: - added support for ASCII CTD files
#	Dec 10, 2010: - change -w) default to 120s
#				  - changed nshear output to 0 from nan when there are no samples
#	Dec 27, 2010: - changed sign of -l to accept lag output from [LADCP_w]
#	Jan 10, 2011: - -o => -k added new -o
#				  - added code to fill CTD sound vel gaps
#	Jan 22, 2011: - added -g) lat,lon
#				  - added -c)ompass corr
#	Jun 15, 2011: - added mean backscatter profile to default output

($ANTS) 	  = (`which list` =~ m{^(.*)/[^/]*$});
($PERL_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});
($LADCPPROC)  = ($0 =~ m{^(.*)/[^/]*$});

require "$ANTS/ants.pl";
require "$ANTS/libEOS83.pl";
require "$ANTS/libstats.pl";
require "$LADCPPROC/LADCPproc.loadCTD";
require "$LADCPPROC/LADCPproc.bestLag";
require "$LADCPPROC/LADCPproc.BT";
require "$LADCPPROC/LADCPproc.backscatter";
require "$LADCPPROC/LADCPproc.UHcode";
require "$PERL_TOOLS/RDI_BB_Read.pl";
require "$PERL_TOOLS/RDI_Coords.pl";
require "$PERL_TOOLS/RDI_Utils.pl";

$antsParseHeader = 0;
&antsUsage('24ab:c:df:g:i:kl:n:o:ps:t:w:',2,
	'[use -2)dary CTD sensor pair]',
	'[require -4)-beam LADCP solutions]',
	'[-s)etup <file>] [-g)ps <lat,lon>]',
	'[-c)ompass-corr <offset,cos-fac,sin-fac>]',
	'[enable -p)PI editing]',
	'[-o)utput grid <resolution[5m]>]',
	'[-i)nitial LADCP time lag <guestimate>]',
	'[-l)ag LADCP <by>] [auto-lag -w)indow <size[120s]>] [-n) <auto-lag windows[20]]',
	'[-d)iagnostic output]',
	'output: [-t)ime series <file>] [-f)lag <file>] [-b)ottom-track <file>]',
	'        [per-bin -a)coustic backscatter profiles] [bottom-trac-k) profs]',
	'<RDI file> <SeaBird file>');

$RDI_Coords::minValidVels = 4 if ($opt_4);

&antsFloatOpt($opt_l);
&antsCardOpt(\$opt_w,120);
	# old default of -w 30 does not work if there are significant ambiguity-velocity
	# problems, as is the case, e.g., with 2010_DIMES_US2 station 142
	# old default of -w 60 did not work for DIMES_UK2 station 4 (DL), possibly again
	# related to ambiguity velocity
&antsCardOpt(\$opt_n,20);
&antsFileOpt($opt_s);
&antsFloatOpt($opt_i);
&antsCardOpt($opt_o);

if (defined($opt_g)) {
	($CTD{lat},$CTD{lon}) = split(',',$opt_g);
	croak("$0: cannot decode -g $opt_g\n")
		unless numberp($CTD{lat}) && numberp($CTD{lon});
}

if (defined($opt_c)) {
	($CC_offset,$CC_cos_fac,$CC_sin_fac) = split(',',$opt_c);
	croak("$0: cannot decode -c $opt_c\n")
		unless numberp($CC_offset) && numberp($CC_cos_fac) && numberp($CC_sin_fac);
}
	

$LADCP_file = &antsFileArg();
$CTD_file 	= &antsFileArg();

&antsAddParams('LADCP_file',$LADCP_file,'CTD_file',$CTD_file);
&antsActivateOut();

#----------------------------------------------------------------------
# Step 1: Read LADCP Data
#----------------------------------------------------------------------

print(STDERR "Reading LADCP data ($LADCP_file)...");
readData($LADCP_file,\%LADCP);
printf(STDERR "\n\t%d ensembles\n",scalar(@{$LADCP{ENSEMBLE}}));

#----------------------------------------------------------------------
# Step 2: Set Processing Parameters
#----------------------------------------------------------------------

print(STDERR "Setting processing parameters...\n");

printf(STDERR "\tloading $LADCPPROC/LADCPproc.defaults...\n");
require "$LADCPPROC/LADCPproc.defaults";

if (defined($opt_s)) {
	print(STDERR "\tloading $opt_s...\n");
	require $opt_s;
}

if ($LADCP{BLANKING_DISTANCE} == 0) {
	print(STDERR "\t\tBLANKING_DISTANCE == 0 => excluding all data from bin 1\n")
		if ($opt_d);
	$wbin_start = 2 unless ($wbin_start > 2);
	$ubin_start = 2 unless ($ubin_start > 2);
	$shbin_start = 2 unless ($shbin_start > 2);
}

&antsAddParams('ADCP_orientation',
		$dta->{ENSEMBLE}[0]->{XDUCER_FACING_UP} ? 'uplooker' : 'downlooker');

$SHEAR_PREGRID_DZ = 20;
$GRID_DZ = defined($opt_o) ? $opt_o : 5;

my($year)  = substr($LADCP{ENSEMBLE}[0]->{DATE},6,4);
my($month) = substr($LADCP{ENSEMBLE}[0]->{DATE},0,2);
my($dau  ) = substr($LADCP{ENSEMBLE}[0]->{DATE},3,2);
my($magdec,$maginc,$h_strength,$v_strength) = split('\s+',`magdec $CTD{lon} $CTD{lat} $year $month $day`);

croak("$0: unknown magnetic declination\n")
	unless defined($magdec);

&antsAddParams('magnetic_declination',$magdec);

#----------------------------------------------------------------------
# Step 3: Read CTD data
#----------------------------------------------------------------------

print(STDERR "Reading CTD data ($CTD_file)...");
readCTD($CTD_file,\%CTD);
printf(STDERR "\n\t%d scans\n",scalar(@{$CTD{press}}));

#----------------------------------------------------------------------
# Step 4: Pre-Process CTD & LADCP Data
#----------------------------------------------------------------------

printf(STDERR "Pre-processing data...");
printf(STDERR "\n\tCTD...");

#------------------------
# clean CTD pressure data
#------------------------
my($pSpikes) = 0;
for (my($r)=1; $r<@{$CTD{press}}; $r++) {
	$pSpikes++,$CTD{press}[$r]=nan
		if (abs($CTD{press}[$r]-$CTD{press}[$r-1])/$CTD{sampint} > 2);
}
print(STDERR "\n\t\t$pSpikes pressure spikes removed")
	if ($pSpikes>0 && $opt_d);

#------------------------------------
# calculate w and find deepest record
#------------------------------------
$CTD{maxpress} = -9e99;
for (my($r)=1; $r<@{$CTD{press}}-1; $r++) {
	$CTD{w}[$r] = 0.99*($CTD{press}[$r+1] - $CTD{press}[$r-1]) / (2*$CTD{sampint});
	if ($CTD{press}[$r] > $CTD{maxpress}) {
		$CTD{maxpress} = $CTD{press}[$r];
		$CTD{atbottom} = $r;
    }										
}
printf(STDERR "\n\t\tmax pressure [%ddbar] at scan#%d",$CTD{maxpress},$CTD{atbottom})
	if $opt_d;

print(STDERR "\n");

#----------------------------------------------------------------------
# Step 4b: Pre-Process LADCP Data
#----------------------------------------------------------------------

print(STDERR "\tLADCP...");

#-------------------------------------------
# transform to earth coordinates if required
#-------------------------------------------

$U = 0;		# velocity indices
$V = 1;
$W = 2;
$E = 3;

$LADCP{HEADING_BIAS} = -$magdec;

if ($LADCP{BEAM_COORDINATES}) {
	print(STDERR "\n\t\ttransforming beam to Earth coordinates...")
		if ($opt_d);
	for (my($ens)=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
		$LADCP{ENSEMBLE}[$ens]->{TILT} = &angle_from_vertical($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
		for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
			@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]} =
				velInstrumentToEarth(\%LADCP,$ens,velBeamToInstrument(\%LADCP,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]}));
			@{$LADCP{ENSEMBLE}[$ens]->{PERCENT_GOOD}[$bin]} =					# fake it to fool ref_lr_w
				(0,0,0,defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]) ? 100 : 0);
		}
	}
	$LADCP{BEAM_COORDINATES} = 0;
	$LADCP{EARTH_COORDINATES} = 1;
	unless ($opt_4) {
		print(STDERR "\n\t\t\t3-beam solutions: $RDI_Coords::threeBeam_1 $RDI_Coords::threeBeam_2 $RDI_Coords::threeBeam_3 $RDI_Coords::threeBeam_4\n")
			if ($opt_d);
		&antsAddParams('3_beam_solutions',"$RDI_Coords::threeBeam_1 $RDI_Coords::threeBeam_2 $RDI_Coords::threeBeam_3 $RDI_Coords::threeBeam_4");
	}
} elsif ($LADCP{EARTH_COORDINATES}) {
	if ($opt_d) {
		if ($opt_c) {
			printf(STDERR "\n\t\tcalculating tilt and correcting for compass error and magnetic declination of %.1f deg...\n",$magdec);
		} else {
			printf(STDERR "\n\t\tcalculating tilt and correcting for magnetic declination of %.1f deg...\n",$magdec);
		}
	}
	for (my($ens)=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {
		$LADCP{ENSEMBLE}[$ens]->{TILT} = &angle_from_vertical($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});
		my($hdg) = rad($LADCP{ENSEMBLE}[$ens]->{HEADING});
		$LADCP{HEADING_BIAS} = ($CC_offset + $CC_cos_fac*cos($hdg) + $CC_sin_fac*sin($hdg)) - $magdec
			if ($opt_c);
		for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
			@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]} =
				velApplyHdgBias(\%LADCP,$ens,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]});
		}
	}
} else {
	croak("$0: can only handle beam or earth coordinates\n")
}

#------------------------------------------------------
# construct a depth-vs-time "profile" from the raw data
#------------------------------------------------------

print(STDERR "\t\tconstructing profile time series...")
	if ($opt_d);
	
($LADCP_start,$LADCP_end,$LADCP_bottom,$w_gap_time,$zErr,$maxz) =
    mk_prof(\%LADCP,0,undef,1,6,70,0.1,9999);
croak("\n$LADCP_file: no good ensembles found\n")
    unless defined($LADCP_start);

if ($opt_d) {
	printf(STDERR "\n\t\t\tStart of cast    : %s (#%5d) at %6.1fm\n",
	                    $LADCP{ENSEMBLE}[$LADCP_start]->{TIME},
	                    $LADCP{ENSEMBLE}[$LADCP_start]->{NUMBER},
	                    $LADCP{ENSEMBLE}[$LADCP_start]->{DEPTH});
	printf(STDERR "\t\t\tBottom of cast   : %s (#%5d) at %6.1fm\n",
	                    $LADCP{ENSEMBLE}[$LADCP_bottom]->{TIME},
	                    $LADCP{ENSEMBLE}[$LADCP_bottom]->{NUMBER},
	                    $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});
	printf(STDERR "\t\t\tEnd of cast      : %s (#%5d) at %6.1fm",
	                    $LADCP{ENSEMBLE}[$LADCP_end]->{TIME},
	                    $LADCP{ENSEMBLE}[$LADCP_end]->{NUMBER},
	                    $LADCP{ENSEMBLE}[$LADCP_end]->{DEPTH});
}

print(STDERR "\n");

#----------------------------------------------------------------------
# Step 5: Add CTD to LADCP Data & correct velocities for sound speed
#	- {DEPTH} field is overwritten!
#----------------------------------------------------------------------

print(STDERR "Matching CTD to LADCP time series...");

$opt_l = defined($opt_l) ? -$opt_l : &lagLADCP2CTD();

print(STDERR "Associating CTD data with LADCP ensembles...");

for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
	my($lastSvel);
	my($r) = int(($LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME} - $opt_l) / $CTD{sampint});
	if ($r < 0 && $ens == $LADCP_start) {
		$r = int(($LADCP{ENSEMBLE}[++$ens]->{ELAPSED_TIME} - $opt_l) / $CTD{sampint})
			while ($r < 0);
		printf(STDERR "\n\tCTD data begin with instrument already in water => skipping %ds of LADCP data",
			$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$LADCP_start]->{ELAPSED_TIME});
		$LADCP_start = $ens;
	}
	if ($r > $#{$CTD{press}}) {
		printf(STDERR "\n\tCTD data end while instrument is still in water => truncating %ds of LADCP data",
			$LADCP{ENSEMBLE}[$LADCP_end]->{ELAPSED_TIME}-$LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME});
		$LADCP_end = $ens - 1;
		last;
	}
	my($dr);
	for ($dr=0; !numberp($CTD{press}[$r+$dr]); $dr--) {}
	$LADCP{ENSEMBLE}[$ens]->{DEPTH} = depth($CTD{press}[$r+$dr],$CTD{lat});
	$LADCP{ENSEMBLE}[$ens]->{CTD_W} = $CTD{w}[$r];
	$LADCP{ENSEMBLE}[$ens]->{CTD_TEMP} = $CTD{temp}[$r];
	$LADCP{ENSEMBLE}[$ens]->{CTD_SVEL} = sVel($CTD{salin}[$r],$CTD{temp}[$r],$CTD{press}[$r+$dr]);
	if (numberp($LADCP{ENSEMBLE}[$ens]->{CTD_SVEL})) {
		$lastSvel = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL};
	} else {
		$LADCP{ENSEMBLE}[$ens]->{CTD_SVEL} = $lastSvel;
	}
	my($sscorr) = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL} / $LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND};
	for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
		next unless defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]);
		$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U] *= $sscorr;
		$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V] *= $sscorr;
		$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W] *= $sscorr;
    }
}

&antsAddParams('bottom_depth',round($LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH}),
			   'start_date',$LADCP{ENSEMBLE}[$LADCP_start]->{DATE},
			   'start_time',$LADCP{ENSEMBLE}[$LADCP_start]->{TIME},
			   'bottom_date',$LADCP{ENSEMBLE}[$LADCP_bottom]->{DATE},
			   'bottom_time',$LADCP{ENSEMBLE}[$LADCP_bottom]->{TIME},
			   'end_date',$LADCP{ENSEMBLE}[$LADCP_end]->{DATE},
			   'end_time',$LADCP{ENSEMBLE}[$LADCP_end]->{TIME});

print(STDERR "\n");

#----------------------------------------------------------------------
# Step 6: Calculate Acoustic Backscatter Profile
#----------------------------------------------------------------------

print(STDERR "Finding seabed in acoustic backscatter profiles...");

mk_backscatter_profs($LADCP_start,$LADCP_end);
($water_depth,$sig_water_depth) =
	find_backscatter_seabed($LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH});

$min_hab = $water_depth - $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH};
printf(STDERR "\n\twater depth      = %d(+-%.1f)m",$water_depth,$sig_water_depth);
printf(STDERR "\n\tclosest approach = %dmab",$min_hab);

print(STDERR "\n");

#----------------------------------------------------------------------
# Step 7: Find Seabed
#----------------------------------------------------------------------

print(STDERR "Finding seabed in BT data...");

($BT_water_depth,$sig_BT_water_depth) = 
	find_seabed(\%LADCP,$LADCP_bottom,$LADCP{BEAM_COORDINATES});

if (defined($BT_water_depth)) {
	$min_hab = $BT_water_depth - $LADCP{ENSEMBLE}[$LADCP_bottom]->{DEPTH};
	printf(STDERR "\n\twater depth      = %d(+-%.1f)m",$BT_water_depth,$sig_BT_water_depth);
	printf(STDERR "\n\tclosest approach = %dmab",$min_hab);
#	$water_depth = $BT_water_depth;									# assume BT data are better
#	$sig_water_depth = $sig_BT_water_depth;							# (at least they are higher vertical resolution)
}

unless (defined($water_depth)) {
	print(STDERR "\n\tno seabed found\n");
	print(STDERR "\n\tunknown water depth => PPI editing disabled\n")
		if ($opt_d);
	$clip_margin = 0;
}

print(STDERR "\n");

#----------------------------------------------------------------------
# Step 8: Edit Data
#----------------------------------------------------------------------

print(STDERR "Calculating shear profiles...\n");

$LADCP_start = 1 if ($LADCP_start == 0);							# ensure that there is previous ensemble

print(STDERR "\tdowncast...");
edit_velocity($LADCP_start,$LADCP_bottom);							# downcast
calc_shear($LADCP_start,$LADCP_bottom,$SHEAR_PREGRID_DZ,0);			# pre-grid shear @SHEAR_PREGRID_DZm resolution
calc_shear($LADCP_start,$LADCP_bottom,$GRID_DZ,1);					# calculate final gridded shear profile

@dc_sh_n = @sh_n;													# save downcast results
@dc_ush_mu = @ush_mu; @dc_ush_sig = @ush_sig;
@dc_vsh_mu = @vsh_mu; @dc_vsh_sig = @vsh_sig;
@dc_wsh_mu = @wsh_mu; @dc_wsh_sig = @wsh_sig;

print(STDERR "\n\tupcast...");
edit_velocity($LADCP_end,$LADCP_bottom);							# upcast
calc_shear($LADCP_end,$LADCP_bottom,$SHEAR_PREGRID_DZ,0);
calc_shear($LADCP_end,$LADCP_bottom,$GRID_DZ,1);

@uc_sh_n = @sh_n;													# save upcast results
@uc_ush_mu = @ush_mu; @uc_ush_sig = @ush_sig;
@uc_vsh_mu = @vsh_mu; @uc_vsh_sig = @vsh_sig;
@uc_wsh_mu = @wsh_mu; @uc_wsh_sig = @wsh_sig;

print(STDERR "\n\tcombined...");
for (my($gi)=0; $gi<@dc_ush_mu; $gi++) {
	if ($dc_sh_n[$gi]>0 && $uc_sh_n[$gi]>0) {
		$sh_n[$gi] = $dc_sh_n[$gi] + $uc_sh_n[$gi];
		$ush_mu[$gi] = ($dc_sh_n[$gi]*$dc_ush_mu[$gi] + $uc_sh_n[$gi]*$uc_ush_mu[$gi]) / $sh_n[$gi];
		$vsh_mu[$gi] = ($dc_sh_n[$gi]*$dc_vsh_mu[$gi] + $uc_sh_n[$gi]*$uc_vsh_mu[$gi]) / $sh_n[$gi];
		$wsh_mu[$gi] = ($dc_sh_n[$gi]*$dc_wsh_mu[$gi] + $uc_sh_n[$gi]*$uc_wsh_mu[$gi]) / $sh_n[$gi];
		$ush_sig[$gi] = sqrt(($dc_sh_n[$gi]*$dc_ush_sig[$gi]**2 + $uc_sh_n[$gi]*$uc_ush_sig[$gi]**2) / $sh_n[$gi]);
		$vsh_sig[$gi] = sqrt(($dc_sh_n[$gi]*$dc_vsh_sig[$gi]**2 + $uc_sh_n[$gi]*$uc_vsh_sig[$gi]**2) / $sh_n[$gi]);
		$wsh_sig[$gi] = sqrt(($dc_sh_n[$gi]*$dc_wsh_sig[$gi]**2 + $uc_sh_n[$gi]*$uc_wsh_sig[$gi]**2) / $sh_n[$gi]);
	} elsif ($dc_sh_n[$gi] > 0) {
		$sh_n[$gi] = $dc_sh_n[$gi];
		$ush_mu[$gi]  = $dc_ush_mu[$gi];  $vsh_mu[$gi]  = $dc_vsh_mu[$gi];  $wsh_mu[$gi]  = $dc_wsh_mu[$gi];
		$ush_sig[$gi] = $dc_ush_sig[$gi]; $vsh_sig[$gi] = $dc_vsh_sig[$gi]; $wsh_sig[$gi] = $dc_wsh_sig[$gi];
	} elsif ($uc_sh_n[$gi] > 0) {
		$sh_n[$gi] = $uc_sh_n[$gi];
		$ush_mu[$gi]  = $uc_ush_mu[$gi];  $vsh_mu[$gi]  = $uc_vsh_mu[$gi];  $wsh_mu[$gi]  = $uc_wsh_mu[$gi];
		$ush_sig[$gi] = $uc_ush_sig[$gi]; $vsh_sig[$gi] = $uc_vsh_sig[$gi]; $wsh_sig[$gi] = $uc_wsh_sig[$gi];
	} else {
		$sh_n[$gi] = 0;
		$ush_mu[$gi]  = $vsh_mu[$gi]  = $wsh_mu[$gi]  = nan;
		$ush_sig[$gi] = $vsh_sig[$gi] = $wsh_sig[$gi] = nan;
	}
}

print(STDERR "\n");

#----------------------------------------------------------------------
# Step 9: Get bottom track profile
#----------------------------------------------------------------------

print(STDERR "Getting BT profile...");
getBTprof($LADCP_start,$LADCP_end);
print(STDERR "\n");

#----------------------------------------------------------------------
# Step 10: Write Output
#----------------------------------------------------------------------

print(STDERR "Writing shear profiles...");

@antsNewLayout = ('depth','dc_nshear','dc_u_z','dc_u_z.sig','dc_v_z','dc_v_z.sig','dc_w_z','dc_w_z.sig',
						  'uc_nshear','uc_u_z','uc_u_z.sig','uc_v_z','uc_v_z.sig','uc_w_z','uc_w_z.sig',
						  'nshear','u_z','u_z.sig','v_z','v_z.sig','w_z','w_z.sig','Sv','Sv.n');
						  
$commonParams = $antsCurParams;
&antsAddParams('ubin_start',$ubin_start,'ubin_end',$ubin_end,		# record processing params
			   'wbin_start',$wbin_start,'wbin_end',$wbin_end,
			   'shbin_start',$shbin_start,'shbin_end',$shbin_end,
			   'w_ref_bin',$w_ref_bin,'w_dif',$w_dif,
			   'wake_hd_dif',$wake_hd_dif,'wake_ang_min',$wake_ang_min,
			   'min_wake_w',$min_wake_w,'n_wake_bins',$n_wake_bins,
			   'e_max',$e_max,'min_cor',$min_cor,
			   'max_shdev',$max_shdev,'max_shdev_sum',$max_shdev_sum,
			   'water_depth',round($water_depth),'water_depth.sig',round($sig_water_depth),
			   'min_hab',round($min_hab),
			   'clip_margin',$clip_margin,'first_clip_bin',$first_clip_bin,
			   'Svbin_start',$Svbin_start,'Svbin_end',$Svbin_end);

for (my($gi)=0; $gi<@ush_mu; $gi++) {
	&antsOut(depthOfGI($gi),										# depth in center of bin
			 numberp($dc_sh_n[$gi])?$dc_sh_n[$gi]:0,				# downcast
			 $dc_ush_mu[$gi],$dc_ush_sig[$gi],
			 $dc_vsh_mu[$gi],$dc_vsh_sig[$gi],
			 $dc_wsh_mu[$gi],$dc_wsh_sig[$gi],
			 numberp($uc_sh_n[$gi])?$uc_sh_n[$gi]:0,				# upcast
			 $uc_ush_mu[$gi],$uc_ush_sig[$gi],
			 $uc_vsh_mu[$gi],$uc_vsh_sig[$gi],
			 $uc_wsh_mu[$gi],$uc_wsh_sig[$gi],
			 $sh_n[$gi],											# combined
			 $ush_mu[$gi],$ush_sig[$gi],
			 $vsh_mu[$gi],$vsh_sig[$gi],
			 $wsh_mu[$gi],$wsh_sig[$gi],
			 $nSv_prof[$gi]?$sSv_prof[$gi]/$nSv_prof[$gi]:nan,
			 $nSv_prof[$gi],
	);
}

print(STDERR "\n");

#----------------------------------------------------------------------

if (defined($opt_a)) {
	print(STDERR "Writing per-bin acoustic backscatter profiles...");

	for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
		my($fn) = sprintf("bin%02d.Sv",$bin);
		print(STDERR " $fn");
		
		@antsNewLayout = ('depth','Sv');
		&antsOut('EOF');
		$antsCurParams = $commonParams;
		close(STDOUT);
		open(STDOUT,">$fn") || croak("$fn: $!\n");
	    
		for (my($gi)=0; $gi<@sSv; $gi++) {
		    &antsOut(depthOfGI($gi),
				     $nSv[$gi][$bin] ? $sSv[$gi][$bin]/ $nSv[$gi][$bin] : nan);
	    }
	}
	print(STDERR "\n");
}

#----------------------------------------------------------------------

if (defined($opt_t)) {
	print(STDERR "Writing time series to $opt_t...");
	
	@antsNewLayout = ('ens','elapsed','depth','CTD_w','LADCP_w');
	&antsOut('EOF');
	$antsCurParams = $commonParams;
	close(STDOUT);
	open(STDOUT,">$opt_t") || croak("$opt_t: $!\n");
	
	for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
		&antsOut($LADCP{ENSEMBLE}[$ens]->{NUMBER},
				 $LADCP{ENSEMBLE}[$ens]->{ELAPSED_TIME},
				 $LADCP{ENSEMBLE}[$ens]->{DEPTH},
				 $LADCP{ENSEMBLE}[$ens]->{CTD_W},
				 $LADCP{ENSEMBLE}[$ens]->{W});
	}
	print(STDERR "\n");
}

#----------------------------------------------------------------------

if (defined($opt_b)) {
	print(STDERR "Writing bottom-track data to $opt_b...");
	
	@antsNewLayout = ('depth','u','v','w','u.sig','v.sig','w.sig','ndata');
	&antsOut('EOF');
	$antsCurParams = $commonParams;
	close(STDOUT);
	open(STDOUT,">$opt_b") || croak("$opt_b: $!\n");

	my($skipped);
	for (my($gi)=0; $gi<@BT_nsamp; $gi++) {
		$skipped = 1 if ($BT_nsamp[$gi] > 0);
		next unless ($skipped);
		&antsOut(depthOfGI($gi),$BTu[$gi],$BTv[$gi],$BTw[$gi],$BTu_sig[$gi],$BTv_sig[$gi],$BTw_sig[$gi],$BT_nsamp[$gi]);
	}
	print(STDERR "\n");
}

#----------------------------------------------------------------------

if (defined($opt_f)) {
	print(STDERR "Writing data flags to $opt_f...");
	
	@antsNewLayout = ('ens');
	for (my($i)=1; $i<=$LADCP{N_BINS}; $i++) {
		$antsNewLayout[$i] = "bin$i";
	}
	&antsOut('EOF');
	$antsCurParams = $commonParams;

	close(STDOUT);
	open(STDOUT,">$opt_f") || croak("$opt_f: $!\n");
	
	&antsPrintHeaders(STDOUT,@antsNewLayout);
	for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {
		printf('%4d  ',$LADCP{ENSEMBLE}[$ens]->{NUMBER});
		for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {
			printf("%02x ",$edit_flags[$ens][$bin]);
		}
		print($opt_R);
	}
	
	print(STDERR "\n");
}

&antsExit();