#!/usr/bin/perl

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

#                    L A D C P _ W S P E C 

#                    doc: Thu Jun 11 12:02:49 2015

#                    dlm: Thu Mar 31 11:08:44 2016

#                    (c) 2012 A.M. Thurnherr

#                    uE-Info: 154 41 NIL 0 0 72 10 2 4 NIL ofnI

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

$antsSummary = 'calculate VKE window spectra from LADCP profiles';

# HISTORY:

#   Jun 11, 2015: - adapted from [binpgrams]

#   Jun 12, 2015: - renamed %PARAM prefixes

#   Jun 15, 2015: - BUG: de-meaning did not respect _gap variables

#                 - added %output_depth_resolution

#                 - reversed semantics of -d/-u

#   Jun 16, 2015: - reversed semantics of -t

#                 - re-added pwrdens.0 to make output consistent with [binpgrams]

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

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

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

#   Mar  1, 2016: - made trailing message much less frequent

#                 - BUG: program croaked gracelessly or entered infinite loop

#                        when presented with insufficient (no valid) input

#                        data

#   Mar 27, 2016: - added -z)ap

#   Mar 28, 2016: - removed -z

#                 - renamed nsamp to nspec

#                 - added w.nsamp.avg

#   Mar 31, 2016: - changed version %PARAM

#				  - BUG: nspec was nan insted of 0

#				  - replaced wspec:: %PARAM-prefix with LADCP_wspec::

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

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

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

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

require "$WCALC/version.pl";

require "$ANTSLIB/ants.pl";

require "$ANTSLIB/antsfilters.pl";

require "$ANTSLIB/libstats.pl";

require "$ANTSLIB/fft.pl";

require "$ANTSLIB/lfit.pl";

require "$ANTSLIB/nrutil.pl";

require "$ANTSBIN/.lsfit.poly";

require "$ANTSBIN/.nminterp.linear";

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

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

# Usage

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

&antsUsage('bc:dg:o:s:tuw:',0,

            '[poly-o)rder <n[0]> to de-mean data; -1 to disable>] [suppress cosine-t)aper]',

            '[-d)own/-u)pcast-only] [exclude -b)ottom window]',

            '[shortwave -c)utoff <kz or lambda>]',

            '[-s)urface <layer depth to exclude[150m]>',

            '[-g)ap <max depth layer to fill with interpolation[40m]>]',

            '[-w)indow <power-of-two input-records[32]>]',

            '[LADCP-profile(s)]');

&antsIntOpt(\$opt_o,0);                                     # polynomial order to remove

if ($opt_o >= 0) {                                              # init model

    &modelUsage();

    matrix(\@covar,1,$modelNFit,1,$modelNFit);

    vector(\@afunc,1,$modelNFit);

    &antsAddParams('LADCP_wspec::demean_poly_order',$opt_o);

}

croak("$0: cannot ignore both down- and upcast\n")

    unless ($opt_d+$opt_u < 2);

if ($opt_d) {

    &antsAddParams('LADCP_wspec::input_data','dc');

} elsif ($opt_u) {

    &antsAddParams('LADCP_wspec::input_data','uc');

} else {

    &antsAddParams('LADCP_wspec::input_data','dc/uc');

}

&antsAddParams('LADCP_wspec::cos_taper_applied',$opt_t ? 'no' : 'yes');

&antsAddParams('LADCP_wspec::btm_window_included',$opt_b ? 'no' : 'yes');

if (defined($opt_c)) {                                          # shortwave cutoff

    $kzlim = ($opt_c < 1) ? $opt_c : 2*$PI/$opt_c;

    &antsAddParams('LADCP_wspec::shortwave_cutoff',$kzlim);

} else {

    $kzlim = 9e99;

}

&antsCardOpt($opt_w);                                           # window size

&antsCardOpt(\$opt_g,40);                                       # gap length [m]

&antsAddParams('LADCP_wspec::min_gap_thickness',$opt_g);

&antsCardOpt(\$opt_s,150);                                      # surface layer

&antsAddParams('LADCP_wspec::surface_layer',$opt_s);

&ISUsage;                                                       # interpolation model

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

# Read Data & Define Layout

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

croak("LADCP_VKE: spectral-input mode must be selected manually (-f)\n")

    if defined(fnrNoErr('pwrdens.0')) && !defined(fnrNoErr('dc_w'));

$zfnr = fnr('depth');                                           # required fields

$dcwfnr = fnr('dc_w'); $dcsfnr = fnr('dc_w.nsamp');

$ucwfnr = fnr('uc_w'); $ucsfnr = fnr('uc_w.nsamp');

$habfnr = fnrNoErr('hab');                                      # optional fields

&antsInstallBufFull(0);                                         # read entire file

&antsIn();

while (@ants_ && $ants_[0][$zfnr] < $opt_s) {					# remove surface layer

	shift(@ants_);

}

croak("$0: insufficient data (no valid records found)\n")

	unless (@ants_);

for ($trimmed=0; @ants_ && !numberp($ants_[$#ants_][$dcwfnr]) && !numberp($ants_[$#ants_][$ucwfnr]); $trimmed++) {

	pop(@ants_);

}

&antsInfo("$trimmed trailing non-numeric records trimmed")

	if ($trimmed > 1);											# 1 is very common

croak("$0: insufficient data (no valid records found)\n")

	unless (@ants_);

$dz = &antsXCheck($zfnr,0,$#ants_,1.01);						# calc dT; 1% jitter

&antsAddParams("LADCP_wspec::input_depth_resolution",$dz);

$opt_g = int(($opt_g - 1) / $dz);								# [m] -> [records]

unless (defined($opt_w)) {										# window size

	for ($opt_w=32; $opt_w*$dz>600; $opt_w/=2) {}

	&antsInfo("%d-m windows ($opt_w records)",$opt_w*$dz);

}

&antsAddParams('LADCP_wspec::window_size',$opt_w,'LADCP_wspec::output_depth_resolution',$dz*$opt_w);

croak(sprintf("$0: insufficient data (%d records found, %d required)\n",scalar(@ants_),$opt_w))

	unless (@ants_ >= $opt_w);

$zrange = $opt_w * $dz;											# NB: not equal to max-min!!!

$resolution_bandwidth = 1 / $zrange;

$resolution_bandwidth *= 2*$PI;

&antsAddParams('LADCP_wspec::resolution_bandwidth',$resolution_bandwidth);

push(@antsNewLayout,'widx','depth','depth.min','depth.max','hab','nspec','w.nsamp.avg');

for (my($i)=0; $i<$opt_w/2+1; $i++) {

	my($kz) = 2*$PI*$i/$zrange;

	last if ($kz > $kzlim);

	&antsAddParams(sprintf('LADCP_wspec::k.%d',$i),$kz);

	&antsAddParams(sprintf('LADCP_wspec::lambda.%d',$i),$i ? $zrange/$i : inf);

	push(@antsNewLayout,sprintf('pwrdens.%d',$i));

}

push(@antsNewLayout,'pwr.tot');

&antsActivateOut();

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

# interpolate short gaps

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

&ISInit($dcwfnr,$zfnr);

&ISInit($ucwfnr,$zfnr);

my($dcLastValid,$ucLastValid,$interp);

for (my($r)=0; $r<=$#ants_; $r++) {

	if (numberp($ants_[$r][$dcwfnr])) {								# number => no interp.

		$dcLastValid = $r;

	} elsif (defined($dcLastValid)) {

		$ants_[$r][$dcwfnr] = &interpolate($zfnr,$opt_g,$dcwfnr,$dcLastValid,$r);

		$interp++ if numberp($ants_[$r][$dcwfnr]);

	}

	if (numberp($ants_[$r][$ucwfnr])) {								# number => no interp.

		$ucLastValid = $r;

	} elsif (defined($ucLastValid)) {

		$ants_[$r][$ucwfnr] = &interpolate($zfnr,$opt_g,$ucwfnr,$ucLastValid,$r);

		$interp++ if numberp($ants_[$r][$ucwfnr]);

	}

}

&antsInfo("$interp non-numeric values interpolated")

	if ($interp);

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

# loop over windows

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

sub avgF($$)														# average field over window

{

	my($f,$r) = @_;

	my(@vals);

	push(@vals,$ants_[$r++][$f])

		while (@vals < $opt_w);

	return avg(@vals);

}

sub medianF($$)														# average field over window

{

	my($f,$r) = @_;

	my(@vals);

	push(@vals,$ants_[$r++][$f])

		while (@vals < $opt_w);

	return median(@vals);

}

unless ($opt_t) {													# compile taper function only if needed

	sub cosTaperWeight($)

	{

		my($z) = $_[0] - $ants_[$fromR][$zfnr];						# elapsed time

		return ($z<0.1*$zrange || $z>0.9*$zrange) ?

					0.5*(1+cos(10*$PI*$z/$zrange-$PI)) : 1;

	}

}

WINDOW: for (my($widx)=1; 1; $widx++) {

	undef(@out);

	$out[0] = $widx;

	local($fromR) = round(($widx-1)*($opt_w/2));					# local, cuz it's used in cosTaperWeight

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

	# partial bottom window

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

	if ($fromR+$opt_w-1 > $#ants_) {

		last if ($opt_b);

		$out[0] = -1;

		$fromR = @ants_ - $opt_w;

		$opt_b = 1;

	}

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

	# output nan on non-numeric y values

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

	my($dc_gap) = $opt_u;											# exclude dc with -d, uc with -u

	my($uc_gap) = $opt_d;

	for (my($r)=$fromR; $r<$fromR+$opt_w; $r++) {

		$dc_gap |= !numberp($ants_[$r][$dcwfnr]);

		$uc_gap |= !numberp($ants_[$r][$ucwfnr]);

	}

	if ($dc_gap && $uc_gap) {

		push(@out,$ants_[$fromR+$opt_w/2][$zfnr]);		

		if ($ants_[0][$zfnr] > $ants_[1][$zfnr]) {		

			push(@out,$ants_[$fromR+$opt_w-1][$zfnr]);	

			push(@out,$ants_[$fromR][$zfnr]); 			

		} else {										

			push(@out,$ants_[$fromR][$zfnr]);			

			push(@out,$ants_[$fromR+$opt_w-1][$zfnr]);	

	    }

	    push(@out,defined($habfnr) ?								# hab

						avgF($habfnr,$fromR) : nan);

		push(@out,0);												# nspec

		push(@out,nan);												# w.nsamp.avg

		for ($i=0; $i<=$opt_w/2; $i++) {							# power

			push(@out,nan);

		}

		&antsOut(@out);												# output nan record and go to next window

		next WINDOW;

	}

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

	# save current values

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

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

		$dcwbuf[$i] = $ants_[$fromR+$i][$dcwfnr];

		$ucwbuf[$i] = $ants_[$fromR+$i][$ucwfnr];

	}

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

	# polynomial de-"mean"ing

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

	if ($opt_o >= 0) {

		my($calc);

		unless ($dc_gap) {

			&modelInit();												# dc

			for (my($a)=1; $a<=$modelNFit; $a++) { $iA[$a] = 1; }

			&lfit($zfnr,$dcwfnr,-1,\@A,\@iA,\@covar,\&modelEvaluate,$fromR,$fromR+$opt_w-1);

			&modelCleanup();

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

				modelEvaluate($i,$zfnr,\@afunc);

				for ($calc=0,my($p)=1; $p<=$modelNFit; $p++) {

					$calc += $A[$p] * $afunc[$p];

				}

				$ants_[$fromR+$i][$dcwfnr] -= $calc;

	        }

	    }

	    unless ($uc_gap) {

			&modelInit();												# uc

			for (my($a)=1; $a<=$modelNFit; $a++) { $iA[$a] = 1; }

			&lfit($zfnr,$ucwfnr,-1,\@A,\@iA,\@covar,\&modelEvaluate,$fromR,$fromR+$opt_w-1);

			&modelCleanup();

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

				modelEvaluate($i,$zfnr,\@afunc);

				for ($calc=0,my($p)=1; $p<=$modelNFit; $p++) {

					$calc += $A[$p] * $afunc[$p];

				}

				$ants_[$fromR+$i][$ucwfnr] -= $calc;

	        }

	    }

	}

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

	# taper data

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

	unless ($opt_t) {

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

			$ants_[$fromR+$i][$dcwfnr] *= &cosTaperWeight($ants_[$fromR+$i][$zfnr]);

			$ants_[$fromR+$i][$ucwfnr] *= &cosTaperWeight($ants_[$fromR+$i][$zfnr]);

		}

		$taper_correction = 1/0.875;

	} else {

		$taper_correction = 1;

	}

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

	# PSD Estimate

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

#	for (my($r)=0; $r<$opt_w; $r++) {

#		print(STDERR "$ants_[$fromR+$r][$dcwfnr], $ants_[$fromR+$r][$ucwfnr]\n");

#	}

	@dc_coeff = &cFFT($dcwfnr,nan,$opt_w,$fromR)				# FFT

		unless ($dc_gap);

	@uc_coeff = &cFFT($ucwfnr,nan,$opt_w,$fromR)

		unless ($uc_gap);

	croak("$0: -n $opt_w not a power-of-two\n")

		unless (@dc_coeff/2==$opt_w || @uc_coeff/2==$opt_w);

	@dc_pwr = &pgram_onesided($opt_w,@dc_coeff)					# total power

		unless ($dc_gap);

	@uc_pwr = &pgram_onesided($opt_w,@uc_coeff)

		unless ($uc_gap);

	push(@out,$ants_[$fromR+$opt_w/2][$zfnr]);					# middle z

	if ($ants_[0][$zfnr] > $ants_[1][$zfnr]) {					# input descending

		push(@out,$ants_[$fromR+$opt_w-1][$zfnr]);				# z.min

		push(@out,$ants_[$fromR][$zfnr]); 						# z.max

	} else {													# input ascending

		push(@out,$ants_[$fromR][$zfnr]);						# z.min

		push(@out,$ants_[$fromR+$opt_w-1][$zfnr]);				# z.max

    }

    push(@out,defined($habfnr) ?								# hab

					avgF($habfnr,$fromR) : nan);

	my($nspec) = !$dc_gap + !$uc_gap;							# nspec

	push(@out,$nspec);

	my($nsamp_sum) = my($nsn) = 0;								# w.nsamp.avg

	$nsamp_sum+=medianF($dcsfnr,$fromR),$nsn++ unless ($dc_gap);	# median to avoid biasing by short bottle stops

	$nsamp_sum+=medianF($ucsfnr,$fromR),$nsn++ unless ($uc_gap);

	push(@out,$nsamp_sum/$nsn);

	my($totP) = 0;												# power

	my($i);

	for ($i=0; $i<$opt_w/2+1; $i++) {				

		my($sumP) = 0;

		$sumP += $dc_pwr[$i] * $taper_correction unless ($dc_gap);

		$sumP += $uc_pwr[$i] * $taper_correction unless ($uc_gap);

		push(@out,$sumP/$nspec/$resolution_bandwidth)

			unless (antsParam("k.$i") > $kzlim);

		$totP += $sumP;

	}

	push(@out,$totP);											# total power

	&antsOut(@out);

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

	# undo tapering and/or de-meaning

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

	for ($i=0; $i<$opt_w; $i++) {

		$ants_[$fromR+$i][$dcwfnr] = $dcwbuf[$i];

		$ants_[$fromR+$i][$ucwfnr] = $ucwbuf[$i];

    }	

}

&antsInfo("$skipped windows with non-numeric values skipped")

	if ($skipped);

antsExit();