#!/usr/bin/perl

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

#                    L A D C P _ W _ C T D 

#                    doc: Mon Nov  3 17:34:19 2014

#                    (c) 2014 A.M. Thurnherr

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

$antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w';

# HISTORY:

#	Nov  3, 2014: - created

#	Nov  4, 2014: - improved

#	Nov  6, 2014: - BUG: sound speed was not calculated correctly

#				  - added -a

#				  - added conductivity & temperature editing

#	Nov  7, 2014: - loosened outlier editing

#				  - added no-valid-data error message

#				  - modified binning criterion to allow any sampling

#					frequency (not just divisors of 24)

#	Apr 17, 2015: - added in-situ temperature to output

#	Apr 23, 2015: - had to losen unrealistic-soundspeed criterion

#					because of DUK2#45

#	May 13, 2015: - BUG: did not work with casts beginning at depth

#				  - added errors when editing edits too much

#	May 14, 2015: - BUG: depth was wrong for partial-depth casts

#				  - BUG: $badval had not been handled correctly

#				  - BUG: editing errors were too tight

#	Jun 16, 2015: - added STDOUT redirection (default output file) on tty stdout

#				  - added -p)lot option

#	Jun 17, 2015: - added cond to std output

#				  - added conductivity spike filter (different from pressure)

#	Jun 18, 2015: - renamed from SBE_w

#	Jun 19, 2015: - added press & salin to make output suitable for u/v processing

#				  - removed a couple of data assertions

#	Jul 20, 2015: - allow for non-numeric -i

#	Aug 17, 2015: - modified plots

#	Sep 21, 2015: - added lat/lon to output if available in input

#	Sep 22, 2015: - added error on non-numeric station information

#				  - changed plot output name

#	Sep 23, 2015: - added -o

#				  - improved pressure editing

#				  - modified editing algorithm, primarily reducing temperature resolution

#					from 0.1deg to 1deg and finding contiguous range up to single gap in histogram

#					(requring halving of resolution)

#	Oct 12, 2015: - require ANTSlibs V6.2 for release

#	Oct 13, 2015: - added code to deal with CNV files with wrong number of scans in header

#   Oct 13, 2015: - adapted to [version.pl]

#   Jan 25, 2016: - added software version %PARAM

#	Feb 16, 2016: - BUG: %sampling_interval was erroenously called %sampling_frequency

#	Feb 19, 2016: - improved on-surface data detection (conductivity <= 10mS/cm)

#				  - BUG: temperatures < 0 were not allowed

#				  - temperature histogram resolution increased from 1degC to 0.5degC for no good reason

#	Mar 13, 2016: - added simple ASCII file format

#				  - BUG: Layout error when input had no lat/lon

#				  - added lon to simple ASCII format

#	Mar 16, 2016: - adapted to gmt5

#	Mar 19, 2016: - added support for $VERB environment var

#				  - added $libSBE_quiet flag

#	Mar 30, 2016: - added station number of ASCII format

#	Mar 31, 2016: - changed version %PARAM

#	May 26, 2016: - renamed w[.raw] to CTD_w[.raw]

#				  - added winch_w

#	Jul 29, 2016: - change w_CTD plot label to reflect sign convention

#	Oct  3, 2017: - cosmetics

#				  - added -q to suppress plots when generating 1Hz files

#	Dec  9, 2017: - added $antsSuppressCommonOptions = 1;

#	Dec 17, 2017: - added dependencies

#	Mar  8, 2018: - made error messages work with -v

#				  - renamed -l)owpass option to -c)utoff

#				  - added new -l)atitude

#	Mar  9, 2018: - changed -l) to location (include lon as well)

#	May 22, 2018: - added code to remove initial at surface interval

#	Jul 24, 2018: - improved CTD load error message

#	Oct  4, 2018: - minor improvement in log

#				  - BUG: profiles collected with LOTS of wave heave (AAIW 2005, 004)

#						 failed the contiguous pressure test at 2dbar resolution

#					     => increased to 3dbar

#	Oct  5, 2018: - BUG: previous bug "fix" did not solve a real problem, because

#						 problem with AAIW data is lack of pressure resolution

#						 => returned to 2dbar

#				  - added plotting errors

#				  - improved log message

#				  - BUG: initial in-air scans were not handled correctly (nscans not updated)

#	Mar 25, 2019: - changed error message to allow creating 1Hz from 4Hz file (SBE19)

#				  - BUG: ITS was not set. How is this possible?????

#	Apr 21, 2019: - modified code to allow production of 24Hz files (previous code required

#					min 2 samples per bin, allowing for max 12Hz sampling rate)

#	Aug 27, 2019: - began adding correction for dropped CTD scans

#	Aug 28, 2019: - made it work

#				  - added -g

#	Aug 30, 2019: - renamed -g to -m)odulo error correction (not)

#				  - added -f)ill missing data

#	Jun 26, 2020: - added salinity -b)ias correction

# NOTES:

#	w_CTD is positive during the downcast to make the sign of the apparent

#		  water velocity consistent with w_ocean

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

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

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

require "$WCALC/version.pl";

require "$ANTS/ants.pl";

require "$ANTS/fft.pl";

require "$ANTS/libstats.pl";

require "$ANTS/libconv.pl";

require "$ANTS/libGMT.pl";

require "$ANTS/libSBE.pl";

require "$ANTS/libEOS83.pl";

&antsAddParams('LADCP_w_CTD::version',$VERSION);

$antsParseHeader = 0;											# usage

$antsSuppressCommonOptions = 1;

$IS = &antsLoadModel('`','.nminterp','linear');

&antsUsage("ab:c:fgi:l:morp:qs:v:w:$IS_opts",1,

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

	'[use -a)lternate sensor pair]',

	'[correct -S)alinity <bias>]',

	'[-r)etain all data (no editing)] [allow infinite -o)utliers]',

	'[suppress CTD -m)odulo error correction]',

	'[-s)ampling <rate[6Hz]>]',

	'[lowpass w_CTD -c)utoff <limit[2s]>] [-w)inch-speed <granularity[10s]>]',

	'[profile -i)d <id>] [station -l)ocation <lat/lon>]',

	'[-p)lot_basenames <[%03d_w_CTD.ps],[%03d_sspd.ps]>]',

	'[-q)uiet (no plots)]',

	'[-f)ill gaps with linear interpolation]',

	'<SBE CNV file>');

&antsUsageError() unless defined($IS);

&antsFloatOpt(\$opt_c,2);										# default low-pass cutoff for w_CTD

&antsCardOpt(\$opt_s,6);										# default output sampling rate (Hz)

&antsFloatOpt(\$opt_w,10);										# winch velocity granularity

&antsFloatOpt(\$opt_b,0);										# salinity bias

&antsCardOpt(\$opt_v,$ENV{VERB});								# support VERB env variable

$CNVfile = $ARGV[0];											# open CNV file

open(F,&antsFileArg());

&antsAddDeps($CNVfile);

&antsActivateOut();												# activate ANTS file

&antsAddParams('salinity_bias_corrected',$opt_b)

	if ($opt_b != 0);

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

# Read Data

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

sub _croak(@)

{

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

	croak(@_);

}

print(STDERR "Reading $CNVfile...") if ($opt_v);

chomp($rec = <F>);

_croak("$CNVfile: no data\n")

	unless defined($rec);

if ($rec =~ /^\*/) {												# SBE CNV file

	$libSBE_quiet = 1;												# suppress diagnostic messages

	($nfields,$nscans,$sampint,$badval,$ftype,$lat,$lon) =			# decode SBE header 

		SBE_parseHeader(F,0,0); 									# SBE field names, no time check

	_croak("$CNVfile: insufficient time resolution ($sampint s) for ${opt_s}Hz time series\n")

		if (1/$sampint/$opt_s < 1);

	if (defined($opt_l)) {											# set/override station location with -l

		my($slat,$slon) = split('[,/]',$opt_l);

		$lat = GMT2deg($slat);

		$lon = GMT2deg($slon);

		_croak("$0: cannot decode -l $opt_l\n")

			unless numberp($lat) && numberp($lon);

	}

	_croak("$CNVfile: no latitude in header => -l required\n")

		unless numberp($lat);

	&antsAddParams('lat',$lat);

	&antsAddParams('lon',$lon);

	$pressF = fnrNoErr('prdM');										# pressure; SBE19plus

	$pressF = fnr('prDM') unless defined($pressF);					# SBE 9plus

	if ($opt_a) {													# temp/cond alternate sensor pair

		$tempF	= fnr('t190C');

		&antsAddParams('ITS',90);

		$condF	= fnrNoErr('c1S/m');

		if (defined($condF)) {

			$condHistRes = 20;										# 0.2 S/m bins

		} else {

			$condF	= fnr('c1mS/cm');

			$condHistRes = 2;										# 2.0 mS/cm bins

		}

	} else {														# primary sensor pair

		$tempF  = fnrNoErr('tv290C');								# SBE 19plus

		$tempF	= fnr('t090C') unless defined($tempF);				# SBE 9plus

		&antsAddParams('ITS',90);

		$condF	= fnrNoErr('c0S/m');

		if (defined($condF)) {

			$condHistRes = 20;

		} else {

			$condF	= fnr('c0mS/cm');

			$condHistRes = 2;

		}

	}

	unless ($opt_m) {												# set up correction for dropped scans

		$systimeF = &fnrNoErr('timeY');

		$xmerrF  = &fnrNoErr('modError');

		if (defined($systimeF) && defined($xmerrF)) {

			$fill_gaps = 1;

		} else {

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

	    }

	}

	$latF = &fnrNoErr('lat');										# GPS data if available (to make files useful for u/v processing)

	$lonF = &fnrNoErr('lon');

	&antsInstallBufFull(0); 										# read entire CNV file

	&SBEin(F,$ftype,$nfields,$nscans,$badval);

	if (@ants_ != $nscans) {

		if ($opt_v > 1) {

			printf(STDERR "\n\nWARNING: $CNVfile has wrong number of scans in header (%d instead of %d)\n",$nscans,scalar(@ants_));

		} else {

			printf(STDERR "WARNING: $CNVfile has wrong number of scans in header (%d instead of %d)\n",$nscans,scalar(@ants_));

		}

		$nscans = @ants_;

	}

} else { 																			# simple CSV ASCII file format:

	($lat,$lon,$station) = split(',',$rec);

	_croak("$CNVfile: ASCII file format error (1st rec must be lat,lon[,id])\n")		#	header: lat,lon[,station]

		unless numberp($lat) && numberp($lon) &&

				$lat>=-90 && $lat<=90 &&

				$lon>=-360 && $lon<=360;

	&antsAddParams('station',$station) if defined($station);

	&antsAddParams('lat',$lat,'lon',$lon);		

	$sampint = 1/24; $condHistRes = 2;												# 	assumptions: 24Hz

	&antsAddParams('ITS',90,'cond.unit','mS/cm');									#	ITS-90, mS/cm

	$pressF = 1; $tempF = 2; $condF = 3;											# 	Layout: scan press temp cond

	for ($nscans=0; <F>; $nscans++) {

		chomp;

		@{$ants_[$nscans]} = split(',');											# 	CSV format

		_croak(sprintf("$CNVfile: unexpected scan #$ants_[$nscans][0] (%d expected)\n",$nscans+1))

			unless ($ants_[$nscans][0] == $nscans+1);

		$ants_[$nscans][$pressF] = nan unless defined($ants_[$nscans][$pressF]);	# missing values

		$ants_[$nscans][$tempF] = nan unless defined($ants_[$nscans][$tempF]);

		$ants_[$nscans][$condF] = nan unless defined($ants_[$nscans][$condF]);

	}

}

printf(STDERR "\n\t%d scans (%d minutes)",$nscans,round($nscans*$sampint/60))

	if ($opt_v > 1);

printf(STDERR "\n") if ($opt_v);

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

# Fill gaps in CTD time series due to dropped scans (modulo errors)

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

if ($fill_gaps) {

	print(STDERR "Correcting for CTD modulo errors and clock jitter/drift...") if ($opt_v);

	my($scans_filled) = 0;

	my($scans_replaced) = 0;

	my($scans_deleted) = 0;

	for (my($scan)=30; $scan<@ants_; $scan++) {										# start a bit more than 1 second into the cast

		next if ($ants_[$scan][$systimeF] == $ants_[$scan-1][$systimeF]);			# skip forward to next systime second

		if ($ants_[$scan][$systimeF] > $ants_[$scan-1][$systimeF]+1) {				# gap spans at least one full second

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

				$ants_[$scan][$systimeF]-$ants_[$scan-1][$systimeF]-1,

				$ants_[$scan-1][$systimeF]-$ants_[0][$systimeF]);

			my(@splicescan);														# scan to splice in 

			$splicescan[$systimeF] = $ants_[$scan-1][$systimeF] + 1;

			$splicescan[$xmerrF] = $ants_[$scan-1][$xmerrF];

			for (my($i)=0; $i<24; $i++) {

				splice(@ants_,$scan,0,\@splicescan);

				$scans_filled++;

			}

	    }

		my($second_start) = $scan - 1;												# backtrack to beginning of second

		my($scans_this_sec) = 1;

		while ($ants_[$second_start-1][$systimeF] == $ants_[$second_start][$systimeF]) {

			$second_start--; $scans_this_sec++;

		}

		die("\nscan = $scan, second_start = $second_start, st = $ants_[$second_start][$systimeF]")

			unless ($second_start > 0);

		while ($scans_this_sec > 24) {												# CTD clock running fast => remove scans

			splice(@ants_,$scan-1,1);

			$scans_deleted++; $scans_this_sec--; $scan--;

		}

		my($gap_len) = 24 - $scans_this_sec;										# gap length in this second only

		die("gap_len = $gap_len at scan#$scan") unless ($gap_len >= 0);

		my($ngaps) = 0;																# gaps between end of previous and beginning of next section

		for (my($s)=$second_start; $s<=$scan; $s++) {

			$ngaps++ if ($ants_[$s][$xmerrF] > $ants_[$s-1][$xmerrF]);

		}

		my(@splicescan);															# scan to splice in 

		$splicescan[$systimeF] = $ants_[$second_start][$systimeF];

		$splicescan[$xmerrF] = $ants_[$second_start][$xmerrF];

		if ($ngaps==0 && $gap_len>0) {												# clock drift => add scan

#			print(STDERR "adding $gap_len scans to make up for clock drift starting at scan $scan-1 (ss = $second_start) (st = $splicescan[$systimeF])\n");

			for (my($i)=0; $i<$gap_len; $i++) {										# fill gap

				splice(@ants_,$scan-1,0,\@splicescan);

				$scans_filled++; $scan++;

			}

        } elsif ($ngaps == 1) {														# single gap => fill it

			my($gap_location) = $second_start;										# find it

			while ($ants_[$gap_location][$xmerrF] == $ants_[$second_start-1][$xmerrF]) {

				$gap_location++;

			}

			die('INTERNAL ERROR') unless ($gap_location <= $scan);

#			print(STDERR "filling $gap_len gap starting at scan $gap_location (ss = $second_start) (st = $splicescan[$systimeF])\n");

			for (my($i)=0; $i<$gap_len; $i++) {										# fill gap

				splice(@ants_,$gap_location,0,\@splicescan);

				$scans_filled++; $scan++;

			}

		} else {																	# multiple gaps => delete data

			my($last_gap);

			for (my($s)=$second_start; $s<$scan; $s++) {							# replace ambiguous data scans

				next if ($ants_[$s][$xmerrF] == $ants_[$second_start-1][$xmerrF]);

				last if ($ants_[$s][$xmerrF] == $ants_[$scan][$xmerrF]);

				splice(@ants_,$s,1,\@splicescan);									# replace

				$scans_replaced++;

				$last_gap = $s;

            }		

			for (my($i)=0; $i<$gap_len; $i++) {										# fill gaps

				splice(@ants_,$last_gap+1,0,\@splicescan);

				$scans_filled++; $scan++;

			}

		}

		my($scans_this_sec) = 0;

		die("$ants_[$second_start][$systimeF]==$ants_[$second_start-1][$systimeF]+1")

			unless ($ants_[$second_start][$systimeF]==$ants_[$second_start-1][$systimeF]+1);

		for (my($i)=0; $ants_[$second_start+$i][$systimeF]==$ants_[$second_start-1][$systimeF]+1; $i++) {

			$scans_this_sec++;

		}

		die("\n$scans_this_sec scans in second $ants_[$second_start][$systimeF] at scans $second_start .. $scan-1 (gap_len $gap_len, $ngaps gaps)\n")

			unless ($scans_this_sec == 24);

	} # for ($scan

#	my($sts) = 1;

#	for (my($scan)=1; $scan<@ants_; $scan++) {

#		$sts++,next if ($ants_[$scan][$systimeF] == $ants_[$scan-1][$systimeF]);

#		printf(STDERR "%s seconds of data missing after scan %d\n",

#			$ants_[$scan][$systimeF]-$ants_[$scan-1][$systimeF]-1,$scan-1)

#				if ($ants_[$scan][$systimeF]-$ants_[$scan-1][$systimeF]-1 != 0);

#		printf(STDERR "sts = $sts beginning at scan %d (st = $ants_[$scan-1][$systimeF])\n",$scan-$sts)

#			unless ($sts == 24);

#		$sts = 1;

#	}

	if ($opt_v>1 && $scans_filled) {

		printf(STDERR "\n\t%d scans cleared (gap clusters)",$scans_replaced);

	}

	printf(STDERR "\n") if ($opt_v);

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

# Redirect STDOUT to %.6Hz & create %_w_CTD.ps,%_sspd.ps if STDOUT is a tty

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

$id = defined($opt_i) ? $opt_i : &antsParam('station');

_croak("$CNVfile: no station information in header => -i required\n")

	unless defined($id);

_croak("$CNVfile: non-numeric station information <$id> in header => -i required\n")

	unless numberp($id);

if (-t STDOUT) {

	if (numberp($id)) {

		$opt_p = '%03d_w_CTD.ps,%03d_sspd.ps'

			unless defined($opt_p);

		$outfile = sprintf('%03d.%dHz',$id,$opt_s);

	} else {

		$opt_p = '%s_w_CTD.ps,%s_sspd.ps'

			unless defined($opt_p);

		$outfile = sprintf('%s.%dHz',$id,$opt_s);

	}

    open(STDOUT,">$outfile") || die("$outfile: $!\n");

}

undef($opt_p) if $opt_q;											# suppress all plots on -q

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

# Edit Data

#	- pressure outliers & spikes

#	- conductivity outliers & spikes

#	- max_bin+1 is returned to ensure full range

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

sub validRange($)

{

	my($guess_bin) = @_;

	my($min_bin,$max_bin);

	die("assertion failed") unless ($hist[$guess_bin]);

	for ($max_bin=$guess_bin; $hist[$max_bin]; $max_bin++) { }

	for ($min_bin=$guess_bin; $hist[$min_bin]; $min_bin--) { }

	return ($min_bin,$max_bin+1);

}

unless ($opt_r) {

	print(STDERR "Editing Data...") if ($opt_v);

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

	# trim initial scans with

	# 	- nan pressure

	#	- nan conductivity

	#	- conductivity <= 10 mS/cm

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

	my($trimmed) = 0;												# trim leading nan pressures

	shift(@ants_),$trimmed++

#	,printf(STDERR "-> p=$ants_[0][$pressF] c=$ants_[0][$condF]\n")

		until !@ants_ ||

				numberp($ants_[0][$pressF]) &&

			  	numberp($ants_[0][$condF]) &&

			  	(($P{'cond.unit'} eq 'mS/cm' && $ants_[0][$condF] > 10) ||

			     ($P{'cond.unit'} eq 'S/m'   && $ants_[0][$condF] > 1));

	_croak("\n$CNVfile: no valid scans (wrong conductivity units?)\n")

		unless (@ants_);

	$nscans -= $trimmed;

	printf(STDERR "\n\t%d initial in-air scans trimmed",$trimmed) if ($opt_v > 1);

	my($lvp) = $ants_[0][$pressF];

	my($lvc) = $ants_[0][$condF];

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

	# edit pressure outliers outside contiguous range

	#	- 2dbar resolution increased to 3dbar because of 2005 AAIW profile 004

	#	- histogram shifted by 100dbar to allow for negative values

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

	my($press_rez) = 2;

	my($outliers) = my($modeSamp) = 0; my($modeBin,$min,$max); local(@hist);

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

		$ants_[$s][$pressF] = nan unless defined($ants_[$s][$pressF]);

		next unless ($ants_[$s][$pressF]>=-100 && $ants_[$s][$pressF]<=6500);

		my($b) = ($ants_[$s][$pressF]+100) / $press_rez;

		$hist[$b]++;

		next unless ($hist[$b] > $modeSamp);

		$modeSamp = $hist[$b]; $modeBin = $b;

	}

	printf(STDERR "\n\tvalid pressure guess: %d dbar (%d samples)",$press_rez*$modeBin-100,$modeSamp)

		if ($opt_v > 1);

	($min,$max) = validRange($modeBin);

	$min = $press_rez*$min-100; $max = $press_rez*$max-100;

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

		next unless numberp($ants_[$s][$pressF]);