#!/usr/bin/perl

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

#                    L A D C P _ W 

#                    doc: Fri Dec 17 18:11:13 2010

#                    (c) 2010 A.M. Thurnherr

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

# TODO:

#	  detection so that editing stats make sense

#	- own seabed detection (P403)

#	- make upcast-flag valid for yoyo casts

#	- make diagnostic output 3-beam field work for Earth coordinates

#	- remove water-depth from BT code, which is not really used and a bit of an outlier

#	  because mean and stddev are used instead of median/mad

#	- read assumed ADCP soundspeed from data file, instead of assuming 1500m/s

#	? use instrument tilt in sidelobe editing?

$antsSummary = 'calculate vertical velocities from LADCP & CTD time series';

$antsMinLibVersion = 6.0;

# HISTORY:

#	Dec 17, 2010: - created from [mergeCTD+LADCP]

#	Dec 18, 2010: - made to work

#	Dec 19, 2010: - improved considerably

#	Dec 20, 2010: - onward

#				  - BUG: depth-binning was off by 1 bin?!

#				  - added binning correction for instrument tilt

#	Dec 21, 2010: - added -h (seafloor depth)

#	Dec 22, 2010: - BUG: had not applied soundspeed-correction to w

#				  - debugged opt_d

#	Dec 23, 2010: - continued implementation of soundspeed corrections

#	Dec 24, 2010: - added winch_w, wave_w

#				  - removed beampair velocities from code

#	Dec 25, 2010: - adapted for surface-wave correction in terms of acceleration (CTD_w_t)

#				  - removed elapsed_mismatch

#				  - removed winch_w, wave_w

#	Dec 26, 2010: - made -p output layout independent of -x to avoid Makefile problems

#	Dec 30, 2010: - cleaned up some

#				  - folded-in backscatter calculation from shear method

#				  - folded-in BT calculation from shear method

#	Dec 31, 2010: - added weighted mean w profile to output

#	Jan  2, 2010: - BUG: BT_w.mad could bomb with division by 0

#				  - BUG: division by zero if no valid data

#	Jan  5, 2010: - adapted to allow for gaps in CTD time series

#	Feb 16, 2011: - cosmetics

#	Jun 22, 2011: - cosmetics

#	Jun 23, 2011: - disabled error on large rms reflr w

#				  - added -l

#				  - removed CTD headers from output

#	Jun 26, 2011: - added -u

#				  - changed package correction from acceleration to velocity, because of

#					Stan's Antarctic data set where accelerations are zero but package effects are

#					there

#	Jul  2, 2011: - increased tilt default to 15 degrees

#	Jul  3, 2011: - replaced old package-velocity correction -x code by new beamvel correction

#				  - removed -p

#				  - replaced -d by residual (diagnostics) output

#	Jul  4, 2011: - saved a snapshot in [Jul_04_2011]

#				  - removed output of binned profile (but not calculation code, which is required for residual)

#				  - BUG: firstGoodEns or lastGoodEns could end up in a reflr_w gap when valid LADCP data begin

#						 or end with CTD in water

#				  - moved rarely used option -s to -k

#				  - added -s)kip <ens> option

#				  - had to very very slightly relax an assertion (by 1e-10 seconds...)

#	Aug  3, 2011: - added -z)ap to truncate range, based on Stan's 2009 data set (017-019 most clearly)

#				    obs that final 2 velocities are outliers; could be due to color bar, should check with

#					residuals

#	Sep 22, 2011: - removed (commented out) CTD acceleration

#				  - added basename to Parameter File Matching

#	Sep 23, 2011: - cosmetics (lots)

#				  - added -d)

#	Oct  6, 2011: - removed old commented-out code (CTD acceleration & old -x)

#				  - renamed time-series field names

#				  - added uncorrected reflr w to time-series output

#	Oct 10, 2011: - BUG: LADCP_w in .tds output had not been soundspeed corrected

#				  - removed -x

#				  - added "false-positives" editing filter

#	Oct 11, 2011: - adapted to the new parameter file, removing a few options

#	Oct 12, 2011: - added surface-layer editing

#				  - made %min/max_depth/elapsed more accurate

#	Oct 13, 2011: - fiddled

#	Oct 14, 2011: - renamed .prof fields from .N to .nsamp

#				  - added %output_basename

#				  - removed all ">" from open to allow use in pipelines

#				  - made <run-label> argument optional

#				  - replaced stdout output by $w_out

#				  - renamed _out to out_; output_basename to out_basename

#	Oct 15, 2011: - added editWOutliers()

#				  - added step to remove single-ping noise

#	Oct 16, 2011: - BUG: ensemble and bin numbers in output files were off by 1

#	Oct 17, 2011: - added default run label

#				  - added out_BR

#				  - BUG: closed STDOUT caused problems with tee in plotting scripts

#				  - added %run_label

#				  - added PostProcess.sh

#	Oct 18, 2011: - BUG: %{min,max}_ens used ensemble indices, rather than numbers

#				  - disabled false-positives editing (removes too much good data on 2011_IWISE)

#	Oct 19, 2011: - BUG: time-series output included ensembles without valid w

#				  - increased post-process sleep time to 30s

#				  - added corr, Sv and echo_amplitude as new fields to w output

#	Oct 20, 2011: - added editFarBins()

#				  - added downcast flag to time-series output

#				  - added ensemble number to LADCP-time-series output

#	Oct 21, 2011: - BUG: Sv code was not enabled for uplooker data

#	Oct 24, 2011: - added code to copy %lat/%lon from CTD data

#				  - added %ADCP_freuquency, %ADCP_blanking_distance

#				  - %LADCP_bin_length => %ADCP_bin_length

#	Oct 26, 2011: - added $first_guess_timelag

#	Oct 27, 2011: - added w12, w34, ww (pitch/roll-weighted average) to full output

#				  - removed ww because it apparently does not improve the solutions

#				  - moved editTilt() to after time-series calculation to allow

#				    time matching to work even with strict tilt limit

#				  - added correctAttitude()

#	May 22, 2012: - adapted to ANTS V5

#	Oct 15, 2012: - added $edit_data_hook

#				  - HG COMMIT

#				  - separated dc/uc time lagging

#				  - removed support for TLhist

#	Oct 16, 2012: - added support for dc/uc only solutions

#	Oct 17, 2012: - renamed $edit_data_hook to $post_merge_hook

#				  - HG COMMIT

#	Feb 13, 2013: - BUG: CTD_neg_press_offset did not work for CTD time series with -ve starting depth

#	Mar 23, 2013: - cosmetics

#				  - HG COMMIT

#	Apr 22, 2013: - removed opt_? variable aliases

#	May  8, 2013: - replaced default run label (default to profiles) to make more readable directory structure

#	May 14, 2013: - opt_m => w_max_lim

#				  - BUG: time-series output had two messed-up fields

#	May 15, 2013: - added w12 and w34 (2-beam solutions) to profile output

#	May 16, 2013: - added CTD_w_tt to time-series & all-sample outputs

#			      - -a => -d (CTD depth offset)

#				  - implemented pressure-sensor acceleration correction (-a)

#				  - added re-gridding of full profile after ping-coherent error removal

#	Jun  5, 2013: - renamed $discard_velocities_from_beam to $bad_beam

#				  - BUG: $bad_beam did not discard BT_VELOCITY data

#	Jun  6, 2013: - BUG: error message had -a instead of -d

#	Sep  5, 2013: - BUG: w12/w34 do not work for earth-coordinate data, of course

#	Apr 17, 2014: - BUG: edit_tilt was never called when all recorded bins are valid

#	Apr 21, 2014: - updated comments

#	May 19, 2014: - began adding support for PPI filtering

#	May 20, 2014: - changed volume_scattering_coeff to Sv in output

#				  - added editPPI()

#	Jul  6, 2014: - BUG: nan water depth had been interpreted as known

#				  - BUG: sVelProf[] was not allowed to have any gaps

#	Jul  9, 2014: - BUG: Jul 6 bug fixes had been applied to older

#				  		 version

#				  - BUG: code meant to ensure gap-freee svel profiles did not work correctly

#	Jul 12, 2014: - finally made output files executable

#	Apr  5, 2015: - added check for required software

#				  - BUG: removed dc/uc mean w fields from .prof again

#	Apr  7, 2015: - made LADCP_w callable from installation directory

#				  - BUG: -v default was wrong in usage message

#				  - replaced 'ens' in output files by 'ensemble'

# CTD REQUIREMENTS

#	- elapsed		elapsed seconds; see note below

#	- depth

#	- ss			sound speed

#	- w				ddepth/dt

#	- temp			OPTIONAL; used for backscatter calculation (i.e. not very important)

#	- %lat/%lon		OPTIONAL

# 2-BEAM SOLUTIONS

#	- for beam-coordinate data, two separate two-beam solutions (w12 & w34) are calculated

#	- discrepancies can be indicative of package wake effects (e.g. MIXET data)

#	- for earth-coordinate data, w12 and w34 are simply copies of w

# NUMERICAL OPTIONS

#	- the first option in the list cannot be numerical!

#	- if need be, use -v 1 as a dummy option

# ELAPSED TIMES

#	- there are 2 different elapsed times used in this program:

#		1) elapsed based on firstgoodens in the LADCP time series

#		2) CTD elapsed time

#	- CTD elapsed time does not have to start with zero!

#	- do not use the Seabird elapsed field, which is only reported to

#	  3 significant digits, causing significant jitter in dt; however,

#	  at least up to 2010 Seabird simply calculates elapsed time by 

#	  assuming a 24Hz sampling rate and no record drop; therefore,

#	  it is best to calculate elapsed time as %RECNO/24

#	- the elapsed field of the output is the elapsed time from the CTD

#	  file; this is required in order to be able to compare the times

#	  from the uplooker and downlooker-derived vertical velocity

#	  profiles

#	- the .TL files use the LADCP elapsed time in output; as a result,

#	  the DL/UL .TL files have different time axes!

#	- as a result, a profile only starts with elapsed==0 if the CTD

#	  is turned on when the LADCP is already in the water

# BINNED-PROFILE OUTPUT NOTES

#	- residuals are calculated with respect to down-/upcast medians

#	- w12 and w34 2-beam solutions are reported without considering

#	  -k (min samples)

#	- equivalence, assuming -o 10:

#		1) 004DL.prof dc_w depth

#		2) bindata -Sdowncast:1 -Fw.median,depth -n 20 depth 10 004DL.w

# VELOCITY AMBIGUITY ERRORS

#	- quite extensive tests with DIMES US2 station 146, which has a lot of

#	  ambiguity velocity errors, reveal that -m 1 catches those errors

#	  quite nicely

#	- even when the errors are not filtered with -m 1, they do not

#	  affect the w profiles, as long as the median bin values are used

($ANTS) 	  = (`which ANTSlib`   =~ m{^(.*)/[^/]*$});

($ADCP_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$});

($WCALC) 	  = ($0                =~ m{^(.*)/[^/]*$});

$WCALC = '.' if ($WCALC eq '');

die("$0: Korn shell (/bin/ksh) required but not found\n")

	unless (-x '/bin/ksh');

die("$0: Generic Mapping Tools (GMT) required but not found (bad \$PATH?)\n")

	unless (`which psxy` ne '');

die("$0: ANTSlib required but not found (bad \$PATH?)\n")

	unless ($ANTS ne '');

die("$0: ADCP Tools required but not found (bad \$PATH?)\n")

	unless ($ADCP_TOOLS ne '');

require "$ANTS/ants.pl";

require "$ANTS/libstats.pl";

require "$WCALC/edit_data.pl";

require "$WCALC/time_series.pl";

require "$WCALC/time_lag.pl";

require "$WCALC/find_seabed.pl";

require "$WCALC/svel_corrections.pl";

require "$WCALC/acoustic_backscatter.pl";

require "$WCALC/bottom_tracking.pl";

require "$ADCP_TOOLS/RDI_BB_Read.pl";

require "$ADCP_TOOLS/RDI_Coords.pl";

use IO::Handle;

my($TINY) = 1e-16;

#------

# Usage

#------

@ARGS = @ARGV;														# save opts & arguments

$antsParseHeader = 0;

&antsUsage('3:4a:b:c:e:g:h:i:k:m:n:o:p:qs:t:uv:w:x:',1,

	'[-v)erbosity <level[2]>]',

	'[-q)uick (no single-ping denoising)]',

    '[require -4)-beam solutions] [apply beamvel-m)ask <file> if it exists]',

	'[valid LADCP -b)ins <bin[2],bin[*]>',

	'[-c)orrelation <min[70]>] [-t)ilt <max[10deg]> [-e)rr-vel <max[0.1m/s]>]',

	'[-h water <depth>]',

	'[max LADCP time-series -g)ap <length[60s]>]',

	'[offset CTD -d)epth <by[0m]>]',

	'[-i)nitial CTD time offset <guestimate> [-u)se as final]]',

	'[calculate -n) <lags,lags[10,100]>] [lag -w)indow <sz,sz[240s,20s]>] [lag-p)iece <CTD_elapsed_min|+[,...]>]',

	'[require top-3) lags to account for <frac[0.6]> of all]',

	'[pressure-sensor -a)cceleration correction <residual/CTD_w_tt>]',

	'[-o)utput bin <resolution[10m]>] [-k) require <min[20]> samples]',

	'[e-x)ecute <perl-expr>]',

	'<station-id> [run-label]');

&antsUsageError() if ($opt_u && !defined($opt_i));

&antsCardOpt(\$opt_s,0);											# skip # initial ensembles

$opt_p = '+' unless defined($opt_p);								# separate uc/dc time lagging

&antsFloatOpt(\$opt_a,1);											# pressure acceleration correction

$STN = &antsCardArg();												# station id

if (@ARGV) {														# run label

	$RUN = $ARGV[0];

	shift;

} else {

	$RUN = 'profiles';

}

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

# Handle Processing Parameters

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

require "$WCALC/defaults.pl";											# load default/option parameters

require "$processing_param_file";										# load processing parameters

eval($opt_x) if defined($opt_x);										# eval cmd-line expression to override anything

$RDI_Coords::minValidVels = 4 if ($opt_4);								# decode cmd-line arguments

($LADCP_firstBin,$LADCP_lastBin) = split(',',$opt_b);

	croak("$0: cannot decode -b $opt_b\n")

    	unless (numberp($LADCP_firstBin) &&

        	    ($LADCP_lastBin eq '*' || numberp($LADCP_lastBin)));

@number_of_timelag_windows = split(',',$opt_n);

	croak("$0: cannot decode -n $opt_n\n")

		unless numberp($number_of_timelag_windows[0]) && numberp($number_of_timelag_windows[1]);

@length_of_timelag_windows = split(',',$opt_w);

	croak("$0: cannot decode -w $opt_w\n")

		unless numberp($length_of_timelag_windows[0]) && numberp($length_of_timelag_windows[1]);

croak("$0: \$out_basename undefined\n")									# all plotting routines use this

	unless defined($out_basename);

&antsAddParams('out_basename',$out_basename);

&antsAddParams('run_label',$RUN);

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

# Screen Logging

#	- warning levels:

#		0 probably unimportant, e.g. nonsensical parameters that probably won't affect solution

#		1 may be somewhat important

#		2 important

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

open(LOGF,">$out_log") || croak("$out_log: $!\n")

	if defined($out_log);

LOGF->autoflush(1);

sub progress(@)

{

	printf(LOGF @_) if defined($out_log);

	printf(STDERR @_) if ($opt_v > 1);

}

sub info(@)

{

	printf(LOGF "\t",@_) if defined($out_log);

	print(STDERR ($opt_v > 1) ? "\t" : "$LADCP_file: ");

	printf(STDERR @_)  if ($opt_v > 0);

}

sub warning(@)

{

	my($lvl,@msg) = @_;

	if (defined($out_log)) {

		print(LOGF "\nWARNING (level $lvl): ");

		printf(LOGF @msg);

		print(LOGF "\n");

	}

	return if ($opt_v == 0);

	print(STDERR "\n") if ($opt_v > 1);

	print(STDERR "WARNING (level $lvl): ");

	printf(STDERR @msg);

	print(STDERR "\n") if ($opt_v > 1);

}

sub debugmsg(@)

{ printf(STDERR @_) if ($opt_v > 2); }

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

# Read LADCP data

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

progress("Reading LADCP data from <$LADCP_file>...\n");

readData($LADCP_file,\%LADCP);

progress("\t%d ensembles\n",scalar(@{$LADCP{ENSEMBLE}}));

croak("$LADCP_file: not enough LADCP bins ($LADCP{N_BINS}) for choice of -r\n")

	unless ($LADCP{N_BINS} >= $refLr_lastBin);

croak("$0: first reference-layer bin outside valid range\n")

	unless ($refLr_firstBin>=1 && $refLr_firstBin<=$LADCP{N_BINS});

croak("$0: last reference-layer bin outside valid range\n")

	unless ($refLr_lastBin>=1 && $refLr_lastBin<=$LADCP{N_BINS});

croak("$0: first reference-layer bin > last reference-layer bin\n")

	unless ($refLr_firstBin <= $refLr_lastBin);

$LADCP_lastBin = $LADCP{N_BINS}-1

	if ($LADCP_lastBin eq '*');

croak("$0: first valid LADCP bin outside valid range\n")

	unless ($LADCP_firstBin>=1 && $LADCP_firstBin<=$LADCP{N_BINS});

croak("$0: last valid LADCP bin outside valid range\n")

	unless ($LADCP_lastBin>=1 && $LADCP_lastBin<=$LADCP{N_BINS});

croak("$0: first valid LADCP bin > last valid LADCP bin\n")

	unless ($LADCP_firstBin <= $LADCP_lastBin);

warning(0,"first reference-layer bin < first valid LADCP bin\n")

	unless ($refLr_firstBin >= $LADCP_firstBin);

warning(0,"last reference-layer bin > last valid LADCP bin\n")

	unless ($refLr_lastBin <= $LADCP_lastBin);

warning(1,"if at all, bin 1 should not be used for short blank-after-transmit values\n")

	if ($LADCP{BLANKING_DISTANCE}<$LADCP{BIN_LENGTH} && $refLr_firstBin==1);

&antsAddParams('ADCP_bin_length',$LADCP{BIN_LENGTH},

			   'ADCP_frequency',$LADCP{BEAM_FREQUENCY},

			   'ADCP_blanking_distance',$LADCP{BLANKING_DISTANCE});

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

# Edit beam-velocity data

#	0) beam-vel mask on -m

#	1) correlation threshold

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

if ($LADCP{BEAM_COORDINATES}) {

	progress("Editing beam-velocity data...\n");

	if (defined($opt_m) && -r $opt_m) {

		progress("\tmasking beam velocities according to $opt_m:");

		my($nee) = 0;

		open(BVM,$opt_m) || die("$opt_m: $!\n");

		while (<BVM>) {

			s/#.*//;

			s/^\s*$//;

			next if ($_ eq '');

			my($fe,$te,$db) = split;

			die("$opt_m: cannot decode $_\n")

				unless (numberp($fe) && numberp($te) && $te>=$fe && $db>=1 && $db<=4);

			die("$0: assertion failed")

				unless ($LADCP{ENSEMBLE}[$fe-1]->{NUMBER} == $fe &&

						$LADCP{ENSEMBLE}[$te-1]->{NUMBER} == $te);

			for (my($ens)=$fe-1; $ens<=$te-1; $ens++) {

				$nee++;

				for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {

					undef($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$db-1]);

				}

			}

		}

		print(STDERR " $nee ensembles edited\n");

		close(BVM);

	}

	$nvv = $cte = 0;

	for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {

		correctAttitude($ens,$pitch_bias,$roll_bias,$heading_bias);

		$nvv += countValidBeamVels($ens);

		$cte += editCorr($ens,$opt_c);

#		$pte += editTilt($ens,$opt_t);

	}

	croak("$LADCP_file: no valid data\n") unless ($nvv > 0);

	progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$opt_c,$cte,round(100*$cte/$nvv));

#	progress("\tattitude threshold (-t %d deg): %d velocites removed (%d%% of total)\n",$opt_t,$pte,round(100*$pte/$nvv));

} else {

	progress("Editing velocity data...\n");

	croak("$LADCP_file: cannot apply beamvel-mask $opt_m to earth-coordinate data\n")

		if defined($opt_m);

	$nvv = $cte = 0;

	for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {

		$nvv += countValidBeamVels($ens);

		$cte += editCorr_Earthcoords($ens,$opt_c);

#		$pte += editTilt($ens,$opt_t);

	}

	croak("$LADCP_file: no valid data\n") unless ($nvv > 0);

	progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$opt_c,$cte,round(100*$cte/$nvv));

#	progress("\tattitude threshold (-t %d deg): %d velocites removed (%d%% of total)\n",$opt_t,$pte,round(100*$pte/$nvv));

}

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

# Calculate earth velocities

#	- this is done for all bins (not just valid ones), to allow

#	  useless possibility that invalid bins are used for reflr calcs

#	- also calculate separate beam-pair velocities

#	- the UNEDITED velocities are saved for the BT calculations

#	  (W is required, U & V are only used for stats that have not

#	  been very useful so far)

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

if ($LADCP{BEAM_COORDINATES}) {

	my($dummy);

	progress("Calculating earth-coordinate velocities...\n");

	progress("\tdiscarding velocities from beam $bad_beam\n")

		if ($bad_beam);

	$nvw = 0;

	for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) {

		for (my($bin)=$LADCP_firstBin-1; $bin<=$LADCP_lastBin-1; $bin++) {

			if ($bad_beam) {

				undef($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$bad_beam-1]);

				undef($LADCP{ENSEMBLE}[$ens]->{BT_VELOCITY}[$bin][$bad_beam-1]);

			}

			($LADCP{ENSEMBLE}[$ens]->{U}[$bin],

			 $LADCP{ENSEMBLE}[$ens]->{V}[$bin],

			 $LADCP{ENSEMBLE}[$ens]->{W}[$bin],

			 $LADCP{ENSEMBLE}[$ens]->{ERRVEL}[$bin]) =

				velInstrumentToEarth(\%LADCP,$ens,

					velBeamToInstrument(\%LADCP,

										@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]}));

			if (defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin])) {

				$per_bin_nsamp[$bin]++;

				$nvw++;

			}

			$LADCP{ENSEMBLE}[$ens]->{GIMBAL_PITCH} =

			 	gimbal_pitch($LADCP{ENSEMBLE}[$ens]->{PITCH},$LADCP{ENSEMBLE}[$ens]->{ROLL});

			$LADCP{ENSEMBLE}[$ens]->{U_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{U}[$bin];

			$LADCP{ENSEMBLE}[$ens]->{V_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{V}[$bin];

			$LADCP{ENSEMBLE}[$ens]->{W_UNEDITED}[$bin] = $LADCP{ENSEMBLE}[$ens]->{W}[$bin];

			($dummy,$LADCP{ENSEMBLE}[$ens]->{W12}[$bin],

			 $dummy,$LADCP{ENSEMBLE}[$ens]->{W34}[$bin]) =

			 	velBeamToBPEarth(\%LADCP,$ens,@{$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin]});

		}

	}

	progress("\t$nvw valid velocities in bins $LADCP_firstBin-$LADCP_lastBin\n");

	progress("\t3-beam solutions : $RDI_Coords::threeBeam_1 " .

								  "$RDI_Coords::threeBeam_2 " .

								  "$RDI_Coords::threeBeam_3 " .

								  "$RDI_Coords::threeBeam_4\n")