#!/usr/bin/perl

#======================================================================

#                    L A D C P P R O C 

#                    doc: Thu Sep 16 20:36:10 2010

#                    (c) 2010 A.M. Thurnherr & E. Firing

#======================================================================

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

# NOTES:

#	- comments starting with ## are taken from Eric's code

#	- ITS-90 temp field in degC required

#	- salin field prequired

#	- pressure field in dbar required

# 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

($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.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;

	'[use -2)dary CTD sensor pair]',

	'[require -4)-beam LADCP solutions]',

	'[-s)etup <file>]',

	'[enable -p)PI editing]',

	'[-d)iagnostic output]',

	'output: [-t)ime series <file>] [-f)lag <file>] [-b)ottom-track <file>]',

	'<RDI file> <SeaBird file>');

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

&antsFloatOpt($opt_l);

	# 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

&antsCardOpt(\$opt_n,20);

$LADCP_file = &antsFileArg();

$CTD_file 	= &antsFileArg();

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

&antsActivateOut();

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

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

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;

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);

croak("$0: unknown magnetic declination\n")

	unless defined($magdec);

&antsAddParams('magnetic_declination',$magdec);

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

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

printf(STDERR "Pre-processing data...");

printf(STDERR "\n\tCTD...");

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

# clean CTD pressure data

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

my($pSpikes) = 0;

}

print(STDERR "\n\t\t$pSpikes pressure spikes removed")

	if ($pSpikes>0 && $opt_d);

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

# calculate w and find deepest record

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

    }										

}

print(STDERR "\n");

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

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

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) {

			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);

	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]} =

				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");

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

#	- {DEPTH} field is overwritten!

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

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

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

for (my($ens)=$LADCP_start; $ens<=$LADCP_end; $ens++) {

	if ($r < 0 && $ens == $LADCP_start) {

			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;

	}

		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);

	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");

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

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

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");

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

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

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");

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

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

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);

print(STDERR "\n");

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

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

print(STDERR "Getting BT profile...");

getBTprof($LADCP_start,$LADCP_end);

print(STDERR "\n");

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

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

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',

$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);

for (my($gi)=0; $gi<@ush_mu; $gi++) {

	&antsOut(depthOfGI($gi),										# depth in center of bin

			 $dc_ush_mu[$gi],$dc_ush_sig[$gi],

			 $dc_vsh_mu[$gi],$dc_vsh_sig[$gi],

			 $dc_wsh_mu[$gi],$dc_wsh_sig[$gi],

			 $ush_mu[$gi],$ush_sig[$gi],

			 $vsh_mu[$gi],$vsh_sig[$gi],

}

print(STDERR "\n");

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

if (defined($opt_a)) {

	print(STDERR "Writing 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);