#!/usr/bin/perl

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

#                    S B E _ W 

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

#                    dlm: Fri Nov  7 15:52:27 2014

#                    (c) 2014 A.M. Thurnherr

#                    uE-Info: 22 0 NIL 0 0 72 2 2 4 NIL ofnI

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

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

$antsMinLibVersion = 6.0;

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

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

require "$ANTS/ants.pl";

require "$ANTS/fft.pl";

require "$ANTS/libSBE.pl";

require "$ANTS/libEOS83.pl";

$antsParseHeader = 0;											# usage

&antsUsage('al:rs:v:',1,

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

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

	'[-r)etain all data (no editing)]',

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

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

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

($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: unknown latitude\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 = 10;		# 0.1 bins

	} else {

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

		$condHistRes = 1;		# 1.0 bins

	}

} else {														# primary sensor pair

	$tempF	= fnr('t090C');

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

	if (defined($condF)) {

		$condHistRes = 10;

	} else {

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

		$condHistRes = 1;

	}

}

&antsInstallBufFull(0);											# read entire CNV file

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

printf(STDERR "\n\t%d scans",$nscans) if ($opt_v > 1);

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

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

# Edit Data

#	- pressure outliers & spikes

#	- conductivity outliers & spikes

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

unless ($opt_r) {

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

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

	# trim initial records with nan pressures

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

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

	shift(@ants_),$trimmed++

		until numberp($ants_[0][$pressF]);

	printf(STDERR "\n\t%d initial records with nan pressure and/or conductivity trimmed",$trimmed) if ($opt_v > 1);

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

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

	# edit pressure outliers outside contiguous range

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

	my($outliers) = 0; my($min,$max);

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

		$phist[$ants_[$s][$pressF]+1000]++

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

	}

	for ($max=25; 	$phist[$max+1000]||$phist[$max+1001]; $max++) {}	# outliers after 2 gaps

	for ($min=$max; $phist[$min+1000]||$phist[$min+ 999]; $min--) {}

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

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

	# edit conductivity outliers outside contiguous range

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

	$outliers = 0;

	my($modeSamp)=0;

	undef(@phist);

	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)

		$phist[$b]++;

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

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

	}

	for ($max=$modeBin; $phist[$max]||$phist[$max+1]; $max++) {}	# outliers after 2 gaps

	for ($min=$modeBin; $phist[$min]||$phist[$min-1]; $min--) {}

	$max /= $condHistRes; $min /= $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);

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

	# edit temperature outliers outside contiguous range

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

	$outliers = 0;

	my($modeSamp)=0;

	undef(@phist);

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

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

		my($b) = $ants_[$s][$tempF]*10;								# 10th of a degree histogram resolution

		$phist[$b]++;

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

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

	}

	for ($max=$modeBin; $phist[$max]||$phist[$max+1]; $max++) {}	# outliers after 2 gaps

	for ($min=$modeBin; $phist[$min]||$phist[$min-1]; $min--) {}

	$max /= 10; $min /= 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);

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

	# edit pressure spikes based on gradients

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

	for (my($s)=1; $s<$nscans; $s++) {								# calculated 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);

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

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

}

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_frequency',1/$opt_s);

&antsAddParams('sampling_rate',$opt_s);

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

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

$sp = $st = $sc = $np = $nt = $nc = 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;

	}

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

}

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

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

# Calculating derived quantities

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

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

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

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

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

	croak(sprintf("$CNVfile: unrealistic depth %d m at elapsed=%.1f s (r=$r)\n",

		$depth[$r],$elapsed[$r]))

			if numberp($depth[$r]) && ($depth[$r]<0 || $depth[$r]>6100);

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

	croak(sprintf("$CNVfile: unrealistic soundspeed %dm/s at elapsed=%.1fs & depth=%.1fm ($cond[$r],$temp[$r],$press[$r])\n",

		$sspd[$r],$elapsed[$r],$depth[$r]))

			if numberp($sspd[$r]) && ($sspd[$r]<1400 || $sspd[$r]>1600);

}

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

}

push(@w,nan);

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

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

# Low-pass filter velocity data

#	- interpolate missing vertical velocities first

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

if ($opt_l > 0) {

	print(STDERR "Low-pass filtering vertical package velocity...") if ($opt_v);

	&antsAddParams('w_lowpass_cutoff',$opt_l);

	my($trimmed) = 0;

	shift(@w),shift(@depth),shift(@elapsed),shift(@sspd),$trimmed++

		until numberp($w[0]);

	my($interpolated) = 0;

	for ($r=1; $r<@w; $r++) {

		next if numberp($w[$r]);

		my($lv) = $r-1;

		for ($nv=$r+1; $nv<@depth && !numberp($w[$nv]); $nv++) {}

		if ($nv < @depth) {

			while ($r < $nv) {

				$w[$r] = $w[$lv] + ($r-$lv)/($nv-$lv) * ($w[$nv]-$w[$lv]);

				$interpolated++;

				$r++;

			}

		} else {

			$trimmed += @w-$r;

			splice(@w,$r); splice(@depth,$r);

			splice(@elapsed,$r); splice(@sspd,$r);

		}

	}

	&antsAddParams('w_interpolated',$interpolated);

	printf(STDERR "\n\t%d/%d vertical velocities trimmed/interpolated",$trimmed,$interpolated) if ($opt_v > 1);

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

# Zero Pad Data

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

	for ($pot=1; $pot<@w; $pot<<=1) {}									# determine power of two

	for ($r=0; $r<@w; $r++) {											# copy data

		$fftbuf[2*$r] = $w[$r];

		$fftbuf[2*$r+1] = 0;

	}

	printf(STDERR "\n\t%d zero records added",$pot-$r) if ($opt_v > 1);

	while ($r < $pot) { 												# pad with zeroes

		$fftbuf[2*$r] = $fftbuf[2*$r+1] = 0;

		$r++;

	}

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

# Low-Pass Filter

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

	@fco = &FOUR1(-1,@fftbuf);											# forward FFT

	$n = @fco/2;

	for (my($ip)=2; $ip<=$n; $ip+=2) {									# +ve freq fco

		my($in) = 2*$n-$ip; 											# -ve freq fco

		my($f)	= $ip/2/$n*$opt_s; 										# frequency

		$fco[$ip] = $fco[$ip+1] = $fco[$in] = $fco[$in+1] = 0

			if ($f > 1/$opt_l); 										# low-pass filter

	}

	@w_lp = &FOUR1(1,@fco); 											# inverse FFT

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

} else {

	@w_lp = @w;

}

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

print(STDERR "Writing output...\n") if ($opt_v);

@antsNewLayout = ('elapsed','depth','sspd','w.raw','w');

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

	&antsOut($elapsed[$r],$depth[$r],$sspd[$r],$w[$r],$w_lp[2*$r]/@w_lp);

}

exit(0);															# don't flush @ants_