#!/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

#	- make timelagging work for short casts (make sure 10% is not more than half window size)

#	- 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

#	? use instrument tilt in sidelobe editing?

#	? remove profile-binning code but add option to output BT-referenced data

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

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

# 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

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

# 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

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

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

# OUTPUT

#	- default is "long" output with all the data; required to generate Visbeck-type plots

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

#	- -p requests old profile output

#	- BT-referenced profiles are only found in -p output

#	- equivalence, assuming -o 10:

#		1) 004DL.prof dc_w depth

#		2) bindata -Sdowncast:1 -Fw.median,depth -n 20 depth 0 1000 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

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

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

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

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 "$PERL_TOOLS/RDI_BB_Read.pl";

require "$PERL_TOOLS/RDI_Coords.pl";

@ARGS = @ARGV;														# save opts & arguments

$antsParseHeader = 0;

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

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

	'[-m)ax vertical <velocity[1m/s]>',

	'[-a)djust CTD depth <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]>]',

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

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

&antsFloatOpt($opt_i);												# externally supplied lag

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

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

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

# Screen Logging

#	- warning levels:

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

#		1 may be somewhat important

#		2 important

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

sub progress(@)

sub info(@)

{

}

sub warning(@)

{

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

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

	printf(STDERR @msg);

}

sub debugmsg(@)

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

# Read LADCP data

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

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

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

# Edit beam-velocity data

#	1) correlation threshold

#	2) tilt threshold (also sets TILT field in all ensembles)

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

if ($LADCP{BEAM_COORDINATES}) {

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

	$nvv = $cte = 0;

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

		$nvv += countValidBeamVels($ens);

	}

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

} else {

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

	$nvv = $cte = 0;

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

		$nvv += countValidBeamVels($ens);

	}

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

}

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

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

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

	$nvw = 0;

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

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

			($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]}));

			$nvw += defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);

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

		}

	}

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

	    unless ($opt_4);

} else {

	progress("Counting valid vertical velocities...\n");

	$nvw = 0;

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

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

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

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

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

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

			$nvw += defined($LADCP{ENSEMBLE}[$ens]->{W}[$bin]);

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

		}

	}

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

}

croak("$LADCP_file: insufficient valid velocities\n") unless ($nvw > 1000);

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

# S1 STEP: Edit earth-coordinate -velocity data

#	1) error-velocity threshold

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

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

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

	$fpte += editFalsePositives($ens);

}

progress("\terror-velocity threshold (-e %.1f m/s): %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",

progress("\t\"false-positives\" filter: %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",

	$fpte,round(100*$fpte/$nvw));

progress("\trange truncation (-z %d bins): %d velocites removed (%d%% of total in bins $LADCP_firstBin-$LADCP_lastBin)\n",

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

# Calculate LADCP time series

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

progress("Calculating LADCP time-series...\n");

($firstGoodEns,$lastGoodEns,$LADCP_atbottom,$LADCP_w_gap_time) =

croak("$LADCP_file: no good ensembles\n")

	unless defined($firstGoodEns) && ($lastGoodEns-$firstGoodEns > 0);

my($cast_duration) = $LADCP{ENSEMBLE}[$lastGoodEns]->{ELAPSED} -

				     $LADCP{ENSEMBLE}[$firstGoodEns]->{ELAPSED};

croak("$0: implausibly short cast ($cast_duration seconds)\n")

	unless ($cast_duration > 600);

$LADCP{MEAN_DT} = $cast_duration / ($lastGoodEns-$firstGoodEns-1);

progress("\tStart of cast     : %s (#%5d)\n",

					$LADCP{ENSEMBLE}[$firstGoodEns]->{TIME},

					$LADCP{ENSEMBLE}[$firstGoodEns]->{NUMBER});

progress("\tBottom of cast    : %s (#%5d) @ dz=%6dm\n",

					$LADCP{ENSEMBLE}[$LADCP_atbottom]->{TIME},

					$LADCP{ENSEMBLE}[$LADCP_atbottom]->{NUMBER},

					$LADCP{ENSEMBLE}[$LADCP_atbottom]->{DEPTH});

progress("\tEnd of cast       : %s (#%5d) @ dz=%6dm\n",

					$LADCP{ENSEMBLE}[$lastGoodEns]->{TIME},

					$LADCP{ENSEMBLE}[$lastGoodEns]->{NUMBER},

					$LADCP{ENSEMBLE}[$lastGoodEns]->{DEPTH});

progress("\tCast duration     : %.1f hours (pinging for %.1f hours)\n",

					$cast_duration / 3600,

					($LADCP{ENSEMBLE}[$#{$LADCP{ENSEMBLE}}]->{UNIX_TIME} -

						$LADCP{ENSEMBLE}[0]->{UNIX_TIME}) / 3600);

progress("\tMean ping interval: %.1f seconds\n",$LADCP{MEAN_DT});

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

# Read CTD data

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

undef($antsOldHeaders);

($CTD_elapsed,$CTD_depth,$CTD_svel,$CTD_w) = &fnr('elapsed','depth','ss','w');

$CTD_temp = &fnrNoErr('temp');

$CTD_maxdepth = -1;

do {

	croak("$0: cannot deal with non-numeric CTD elapsed time\n")

		unless &antsNumbers($CTD_elapsed);

	push(@{$CTD{ELAPSED}},$ants_[0][$CTD_elapsed]);

	push(@{$CTD{SVEL}},   $ants_[0][$CTD_svel]);

	push(@{$CTD{W}},      $ants_[0][$CTD_w]);

	push(@{$CTD{TEMP}},	  $ants_[0][$CTD_temp]) if defined($CTD_temp);

		$CTD_atbottom = $#{$CTD{DEPTH}};

	}

} while (&antsIn());

$CTD{DT} = ($CTD{ELAPSED}[$#{$CTD{ELAPSED}}] - $CTD{ELAPSED}[0]) / $#{$CTD{ELAPSED}};

progress("\t%d scans at %.1fHz\n",scalar(@{$CTD{DEPTH}}),1/$CTD{DT});

progress("\tstart depth = %.1fm\n",$CTD{DEPTH}[0]);

croak("$0: CTD start depth must be numeric\n")

	unless numberp($CTD{DEPTH}[0]);

progress("\tmax depth   = %dm (# $CTD_atbottom)\n",$CTD_maxdepth);

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

# Construct sound-speed correction profile from CTD 1Hz downcast data

#	very simple algorithm that stores the last value found

#	in each 1m bin

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

progress("Constructing sound-speed correction profile\n");

my($scans_per_sec) = int(1/$CTD{DT}+0.5);

for (my($s)=0; $s<=$CTD_atbottom; $s+=$scans_per_sec) {

	next unless ($CTD{DEPTH}[$s] >= 0 && numberp($CTD{SVEL}[$s]));

	$sVelProf[int($CTD{DEPTH}[$s])] = $CTD{SVEL}[$s];

}

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

# Determine time lag

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

if (defined($opt_i)) {

	progress("Setting initial time lag...\n");

	$CTD{TIME_LAG} = $opt_i;

	progress("\t-i => elapsed(CTD) ~ elapsed(LADCP) + %.1fs\n",$CTD{TIME_LAG});

} else {

	progress("Guestimating time lag...\n");

	my($CTD_10pct_down) = 0;

	my($LADCP_10pct_down) = $firstGoodEns;

	$CTD_10pct_down++

		until ($CTD{DEPTH}[$CTD_10pct_down]-$CTD{DEPTH}[0] >= 0.1*($CTD_maxdepth-$CTD{DEPTH}[0]));

	$LADCP_10pct_down++

		until ($LADCP{ENSEMBLE}[$LADCP_10pct_down]->{DEPTH} >= 0.1*$LADCP{ENSEMBLE}[$LADCP_atbottom]->{DEPTH});

	$CTD{TIME_LAG} =

		$CTD{ELAPSED}[$CTD_10pct_down] - $LADCP{ENSEMBLE}[$LADCP_10pct_down]->{ELAPSED};

	progress("\telapsed(dz(CTD)=%.1fm) ~ elapsed(dz(LADCP)=%.1fm) + %.1fs\n",