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

($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/libGMT.pl";

require "$ANTS/libSBE.pl";

require "$ANTS/libEOS83.pl";

&antsAddParams('LADCP_w_CTD',"Version $VERSION");

$antsParseHeader = 0;											# usage

&antsUsage('ai:l:orp:s:v:',1,

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

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

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

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

	'[-l)owpass w <cutoff[2s]>]',

	'[profile -i)d <id>]',

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

	'<SBE CNV file>');

&antsFloatOpt(\$opt_l,2);										# defaults

&antsCardOpt(\$opt_s,6);

&antsCardOpt(\$opt_v,0);

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

open(F,&antsFileArg());

&antsActivateOut();												# activate ANTS file

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

# Read Data

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

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

chomp($rec = <F>);

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

	unless defined($rec);

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

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

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

	croak("$CNVfile: unexpected sampling interval $sampint\n")

		unless (abs($sampint-1/24) < 1e-5);

	croak("$CNVfile: no latitude in header\n")

		unless numberp($lat);

	$pressF = fnr('prDM');

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

		$tempF	= fnr('t190C');

		$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	= fnr('t090C');

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

		if (defined($condF)) {

			$condHistRes = 20;

		} else {

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

			$condHistRes = 2;

		}

	}

	$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\n");

		} else {

			printf(STDERR "WARNING: $CNVfile has wrong number of scans in header\n");

		}

		$nscans = @ants_;

	}

} else { 																			# simple ASCII file format:

	$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("$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",$nscans) if ($opt_v > 1);

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

}

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

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

	# 	- nan pressure

	#	- nan conductivity

	#	- conductivity <= 10 mS/cm

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

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

	shift(@ants_),$trimmed++

		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 records (wrong conductivity units?)\n")

		unless (@ants_);			     

	printf(STDERR "\n\t%d initial at-surface records trimmed",$trimmed) if ($opt_v > 1);

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

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

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

	# edit pressure outliers outside contiguous range

	#	- 2dbar resolution

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

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

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

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

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

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

		$hist[$b]++;

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

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

	}

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

		if ($opt_v > 1);

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

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

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

		if ($ants_[$s][$pressF] > $max) { $outliers++; $ants_[$s][$pressF] = nan; }

		if ($ants_[$s][$pressF] < $min) { $outliers++; $ants_[$s][$pressF] = nan; }

	}

	&antsAddParams("pressure_outliers",sprintf("%d",$outliers));

	printf(STDERR "\n\tcontinuous pressure range: %d..%d dbar (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);

	croak("$CNVfile: pressure editing removed too many 'outliers'\n")

		unless ($opt_o || $outliers < 100);

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

	# edit conductivity outliers outside contiguous range

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

	$outliers = $modeSamp = 0;

	undef(@hist);

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

		next unless ($ants_[$s][$condF] > 0);

		my($b) = $ants_[$s][$condF]*$condHistRes;					# 1/10 S/m histogram resolution (1 mS/cm)

		$hist[$b]++;

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

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

	}

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

	$min /= $condHistRes; $max /= $condHistRes;

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

		if ($ants_[$s][$condF] > $max) { $outliers++; $ants_[$s][$condF] = nan; }

		if ($ants_[$s][$condF] < $min) { $outliers++; $ants_[$s][$condF] = nan; }

	}

	&antsAddParams("conductivity_outliers",sprintf("%d",$outliers));

	printf(STDERR "\n\tcontinuous conductivity range: %.1f..%.1f S/m (%d outliers removed)",$min,$max,$outliers) if ($opt_v > 1);

	croak("$CNVfile: conductivity editing removed too many 'outliers'\n")

		unless ($opt_o || $outliers/$nscans < 0.4);

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

	# edit temperature outliers outside contiguous range

	#	- Stan's NBP0901 profiles require resolution of 1deg

	#	- otherwise 0.2deg seems to be fine

	#	- however, on Feb 19, 2016 it was found that the

	#     resolution had been left at 1degC without any

	#     apparent adverse effects

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

	$outliers = $modeSamp = 0;

	undef(@hist);

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

		next unless ($ants_[$s][$tempF] >= -10);

		my($b) = ($ants_[$s][$tempF] + 10) * 0.5;

		$hist[$b]++;

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

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

	}

#	printf(STDERR "\n\ttemperature mode: %.1f degC (%d samples)",$modeBin/0.5-10,$modeSamp)

#		if ($opt_v > 1);

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

	$min = ($min / 0.5) - 10;

	$max = ($max / 0.5) - 10;

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

		if ($ants_[$s][$tempF] > $max) { $outliers++; $ants_[$s][$tempF] = nan; }

		if ($ants_[$s][$tempF] < $min) { $outliers++; $ants_[$s][$tempF] = nan; }

	}

	&antsAddParams("temperature_outliers",sprintf("%d",$outliers));

	printf(STDERR "\n\tcontinuous temperature range: %.1f..%.1f degC (%d outliers removed)",$min,$max,$outliers)

		if ($opt_v > 1);

	croak("$CNVfile: temperature editing removed too many 'outliers'\n")

		unless ($opt_o || $outliers/$nscans < 0.4);

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

	# edit pressure spikes based on gradients

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

	for (my($s)=1; $s<$nscans; $s++) {								# calculate pressure gradients (across gaps)

		if (numberp($ants_[$s][$pressF])) {

			$dp[$s-1] = $ants_[$s][$pressF] - $lvp;

			$lvp = $ants_[$s][$pressF];

		} else {

			$dp[$s-1] = nan;

		}

	}

	my($ns1,$ns2) = (0,0);

	for (my($s)=0; $s<$nscans-2; $s++) {							# consecutive large pressure gradients of opposite sign

		if (($dp[$s]*$dp[$s+1] < 0) &&								# tests return false if either of the dps is not defined

			(abs($dp[$s]) > 10) &&

			(abs($dp[$s+1]) > 10)) {

				$ants_[$s+1][$pressF] = nan;

				$dp[$s] = $dp[$s+1] = undef;

				$ns1++;

		}

	}

	for (my($s)=0; $s<$nscans-3; $s++) {							# 3 consecutive large pressure gradients of opposite sign

		if (($dp[$s]>2	&& $dp[$s+1]<-4 && $dp[$s+2]>2) ||

			($dp[$s]<-2 && $dp[$s+1]>4	&& $dp[$s+2]<-2)) {

				$ants_[$s+1][$pressF] = $ants_[$s+2][$pressF] = nan;

				$dp[$s] = $dp[$s+1] = $dp[$s+2] = undef;

				$ns2+=2;

		}

	}

	&antsAddParams("pressure_spikes_removed",sprintf("%d+%d",$ns1,$ns2));

	printf(STDERR "\n\t%d+%d pressure spikes removed",$ns1,$ns2) if ($opt_v > 1);

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

	# edit conductivity spikes based on large gradients

	#	- $lvc = $ants_[0] is guaranteed numeric here

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

	my($nv) = my($ns) = my($last_dcond) = my($lvcs) = 0;

	my($run_start) = my($run_dcond) = 0;

	for (my($s)=1; $s<$nscans; $s++) {								# calculate conductivity gradients (across gaps)

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

		my($dcond) = $ants_[$s][$condF] - $lvc;						# integrate gradient across runs

#		print(STDERR "ldc/dc/rdc: $last_dcond/$dcond/$run_dcond\n");

		if ($last_dcond*$dcond >= 0) {								# run is continuing (gradient does not change sign)

			$run_dcond += $dcond;

		} else {													# run has ended

#			print(STDERR "new run at $lvcs (run_dcond = $last_dcond)\n");

			$run_start = $lvcs;

			$run_dcond = $dcond;

		}

		if (abs($run_dcond) <= 0.02) {								# small integrated gradient => okay

			$lvc = $ants_[$s][$condF];								# update stored previous values

			$lvcs = $s;

			$last_dcond = $dcond;

			next;													# process next scan

		}

#		print(STDERR "large gradient ($ants_[$run_start][$condF]-$ants_[$s][$condF], $run_dcond) run $run_start-$s\n");

		my($i);

		my($max_spike_length) = 24;

		for ($i=$s; $i<=$run_start+$max_spike_length && $i<$nscans-1; $i++) {	# large gradient => check whether values return within 10s

			next unless (numberp($ants_[$i][$condF]) &&

						 numberp($ants_[$i+1][$condF]));

			last if ((abs($ants_[$i][$condF]-$ants_[$run_start][$condF]) <= 0.005) &&	# 2 vals to avoid large gradients

					 (abs($ants_[$i+1][$condF]-$ants_[$run_start][$condF]) <= 0.005));

		}

		if ($i>$run_start+$max_spike_length || $i==$nscans-1) {		# values don't return => leave data alone

#			print(STDERR "values don't return => new run at $s\n");

			$run_start = $s; $run_dcond = 0;

			$lvc = $ants_[$run_start][$condF];						# start new run

			$lvcs = $run_start;

			$last_dcond = 0;

			next;													# process next scan

		}

		$ns++;

		for (my($j)=$run_start+1; $j<$i; $j++) {					# values return => remove bad spike

			$ants_[$j][$condF] = nan;

			$nv++;

		}

#		print(STDERR "values return at $i ($ants_[$i][$condF]) => deleted $run_start+1-$i; new run\n");

		$run_start = $i; $run_dcond = 0;

		$lvc = $ants_[$run_start][$condF];

		$lvcs = $run_start;

		$last_dcond = 0;

	}

	&antsAddParams("conductivity_spikes_removed",$ns);

	printf(STDERR "\n\t%d conductivity values removed from %d spikes",$nv,$ns) if ($opt_v > 1);

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

} # if $opt_r

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

# Correcting for pressure bias

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

print(STDERR "Correcting for pressure bias...") if ($opt_v);

my($minP) = 9e99;

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

	$minP = $ants_[$s][$pressF]

		if numberp($ants_[$s][$pressF]) && ($ants_[$s][$pressF] < $minP);

}

croak("$CNVfile: no valid CTD pressure data below 25dbar\n")

	unless ($minP < 9e99);

if ($minP < 25) {

	&antsAddParams('pressure_bias',$minP);

	printf(STDERR "\n\tsubtracting %.1f dbar",$minP) if ($opt_v > 1);

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

		$ants_[$s][$pressF] -= $minP

			if numberp($ants_[$s][$pressF]);

	}

} else {

	printf(STDERR "\n\tpartial-depth cast below %.1f dbar (no correction applied)",$minP) if ($opt_v > 1);

}

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

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

# Binning data

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

my($sps) = round(1 / $sampint / $opt_s);

print(STDERR "Creating ${opt_s}Hz time series ($sps samples per bin)...") if ($opt_v);

&antsAddParams('sampling_interval',1/$opt_s);

&antsAddParams('sampling_frequency',$opt_s);

my(@press,@temp,@cond);

my($sp,$np,$st,$nt,$sc,$nc);

$sp = $st = $sc = $np = $nt = $nc = 0;

$slat = $slon = $nGPS = 0;

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

	if ($s*$sampint > $rec/$opt_s) {

		$rec++;

		push(@press,$np>0?$sp/$np:nan);

		push(@temp, $nt>0?$st/$nt:nan);

		push(@cond, $nc>0?$sc/$nc:nan);

		$sp = $st = $sc = $np = $nt = $nc = 0;

		if (defined($latF)) {

			push(@lat,$nGPS>0?$slat/$nGPS:nan);

			push(@lon,$nGPS>0?$slon/$nGPS:nan);

			$slat = $slon = $nGPS = 0;

		}

	}

	$sp+=$ants_[$s][$pressF],$np++ if numberp($ants_[$s][$pressF]);

	$st+=$ants_[$s][$tempF],$nt++ if numberp($ants_[$s][$tempF]);

	$sc+=$ants_[$s][$condF],$nc++ if numberp($ants_[$s][$condF]);

	if (defined($latF) && numberp($ants_[$s][$latF])) {

		$slat += $ants_[$s][$latF];

		$slon += $ants_[$s][$lonF];

		$nGPS++;

	}

}

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

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

# Calculating derived quantities

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

print(STDERR "Calculating vertical package velocity & sound speed...") if ($opt_v);

my($maxP) = -9e99; my($atBtm);

my($min_sspd) = 9e99; my($max_sspd) = -9e99;

for (my($r)=0; $r<@press; $r++) {

	$maxP=$press[$r],$atBtm=$r if ($press[$r] > $maxP);

	$elapsed[$r] = $r/$opt_s;

	$depth[$r] = &depth($press[$r],$lat);

#	print(STDERR "depth[$r] = depth($press[$r],$lat) = $depth[$r]\n");

	$salin[$r] = &salin($cond[$r],$temp[$r],$press[$r]);

	$sspd[$r]  = &sVel($salin[$r],$temp[$r],$press[$r]);

	$min_sspd = $sspd[$r] if ($sspd[$r] < $min_sspd);

	$max_sspd = $sspd[$r] if ($sspd[$r] > $max_sspd);

}

$w[0] = nan;

for (my($r)=1; $r<@depth-1; $r++) {

	$w[$r] = numbersp($depth[$r-1],$depth[$r+1])

		   ? ($depth[$r+1] - $depth[$r-1]) * $opt_s

		   : nan;