#!/usr/bin/perl

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

#                    L A D C P _ V K E 

#                    doc: Tue Oct 14 11:05:16 2014 

#                    (c) 2012 A.M. Thurnherr

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

$antsSummary = 'calculate VKE from LADCP-derived vertical-velocity profiles';

# NOTES:

#	- LADCP_VKE takes either one .wspec file (produced by [LADCP_wspec]

#	  and, possibly, post-processed) or any number of .wprof files (the

#	  normal case) as input

#	- when .wprof files are supplied, [LADCP_wspec] is used to calculate

#	  the corresponding VKE spectra

#	- when multiple .wprof files are supplied, the resulting spectra

#	  are averaged

#	- the averaged spectra will only cover the windows that are present in

#	  the final(!) .wprof input files

# TODO:

#	! verify that p0fit.slope.sig is correct (-x scale factor)

# HISTORY:

#	Oct 14, 2014: - created from [LADCPfs]

#	Oct 15, 2014: - added parameterization output

#	Oct 17, 2014: - changed parameterization constant $c to 0.021

#	Nov  6, 2014: - restored from backup and adapted to ANTS V6.0

#	Nov  7, 2014: - changed parameterization constant $c to 0.0215

#	Apr 16, 2015: - disabled output activation unless ANTS tools are available

#				  - removed superfluous $ANTSBIN definition

#	May 18, 2015: - added -p)ulse <length>

#	Jun 11, 2015: - removed w_z code (yfname param requirement)

#	Jun 12, 2015: - adapted to &antsParam()

#				  - BUG: %k.0 and %lambda.0 had been required erroneously (.1 are first ones to be used)

#				  - made finescale limits optional

#				  - renamed -b=>-l, -c=>-o, -s=>-i

#	Jun 14, 2015: - renamed -p=>-a

#				  - added -p

#				  - removed weird evals

#	Jun 15, 2015: - added plot & other mods

#	Jun 16, 2015: - define default outputs when STDOUT is a tty

#				  - adapted to re-added pwrdens.0 in LADCP_wspec output

#				  - modified plot label

#	Jun 17, 2015: - added eps.w to plot

#	Jul 30, 2015: - moved main label and made it blue

#	Sep 28, 2015: - added a p0fit r test to ensure that correlation has the correct sign

#	Oct 12, 2015: - BUG: plot was in landscape mode

#				  - added input run label(s) to plot

#				  - require ANTSlibs V6.2 for release

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

#	Nov 30, 2015: - BUG: -a was not allowed in usage

#				  - added -f

#	Dec 27, 2015: - reduced minlim for eps to 1e-13 W/kg

#	Dec 29, 2015: - added 3rd consistency check (p0 limit)

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

#	Mar  7, 2016: - removed unused spectral normalization code

#				  - added latitude constraint for calculation of eps.w

#	Mar  8, 2016: - renamed fs_pwr to pwr.fs for consistency

#				  - renamed eps.w to eps.VKE for consistency

#				  - changed default output filenames for -d and -u

#				  - removed ./ from figure label

#				  - update ANTSlib to 6.3

#	Mar 16, 2016: - adapted to gmt5

#	Mar 17, 2016: - added version to plot

#				  - added -k)e dissipation

#				  - reduced eps scale

#	Mar 26, 2016: - added slope

#	Mar 27, 2016: - added slope stddev

#				  - added opt_x

#				  - re-designed QC checks

#				  - added -z)ap

#	Mar 28, 2016: - made QC tests consistent with I08S observations based on p0fit.rms == 0.4

#				  - added support for w.nsamp.avg

#				  - removed -z from LDCP_wspec

#				  - reduced -z default to 1, because 2 noise should affect only shotest scales

#	Mar 29, 2016: - changed default of -x to 1 & removed from usage

#				  - BUG: -l/-a did not work

#				  - removed p0fit.nsamp from output (is constant)

#	Mar 30, 2016: - added support for multiple wprof input files

#	Mar 31, 2016: - standardized %PARAMs (got rid of old binpgrams:: params)

#				  - added spectral averaging

#				  - disabled -l/-a code (override bin size & pulse length)

#				  - added removal of params on multiple files

#				  - disabled code to create different plot/output names on -d/-u

#				  - added -q)uatorial cutoff lat

#	Apr  1, 2016: - cosmetics

#	Apr  2, 2016: - added low-p0 power-law fit as alternative

#	Apr  3, 2016: - removed low-p0 power-law fit as it was not well thought out

#				  - changed -l default from 1e-7 to 5e-8, which agrees better with GHP

#					solution without latiudinal dependency in case of 2016 I08S

#				  - define own LADCP_wspec defaults here

#				  - changed default of -q from 5 to 3

#	Apr  6, 2016: - cosmetics

#	Apr 16, 2016: - assertion failed on I98S#153 => modify code to allow nans in LADCP_wspec

#					output

#	Apr 27, 2016: - cosmetics (usage message)

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

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

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

$ANTS_TOOLS_AVAILABLE = (`which list` ne '');

require "$WCALC/version.pl";

require "$ANTSLIB/ants.pl";

require "$ANTSLIB/libLADCP.pl";

require "$ANTSLIB/libGMT.pl";

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

use FileHandle;

use IPC::Open2;

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

# Empirical constants & defaults

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

#my($c) = 0.0215;						# Thurnherr et al. (GRL 2015)

my($c) = 0.026;							# increased by 20% for V1.2beta7

$opt_q = 3;								# Equatorial band: little more than a guess based on 2015 P16N

$opt_z = 1;								# number of w_ocean samples to require

$opt_o = 0;								# remove mean before calculating spectra

$opt_s = 150;							# surface layer to exclude from spectra

$opt_g = 40;							# max gap to interpolate over

$opt_c_default = 100;					# short-wavelength cutoff

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

# Usage

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

		    "[poly-o)rder <n[$opt_o]> to de-mean data; -1 to disable>] [apply cosine-t)aper]",

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

			"[-s)urface <layer depth to exclude[${opt_s}m]>",

            "[-g)ap <max depth layer to fill with interpolation[${opt_g}m]>]",

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

			"[shortwave -c)utoff <kz or lambda[${opt_c_default}m]>]",						# LADCP_VKE options

			"[-z) ignore velocities derived from fewer than <N[$opt_z]> samples]",

			'[o-m)it spectral correction] [spectral-tilt-correction -r)ange <max[0m]>]',

			"[-e)ps-parameterization <constant[${c}s^-0.5]>",

			'[file...]');

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

# Calculate VKE spectra with [LADCP_wspec] if input is a set of w_ocean files

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

my($widx_min) =  99;																# sentinels

my($widx_max) = -99;

if (defined(fnrNoErr('dc_w'))) {													# pre-process with LADCP_wspec when handed vertical-velocity input

	&antsInstallBufFull('eof');														# read first file

	&antsIn();

	my($opts);																		# set up options to pass

	$opts .= ' -d' if ($opt_d);														

	$opts .= ' -u' if ($opt_u);

	$opts .= ' -b' if ($opt_b);

	$opts .= ' -t' if ($opt_t);

	$opts .= " -s $opt_s" if defined($opt_s);

	$opts .= " -g $opt_g" if defined($opt_g);

	$opts .= " -w $opt_w" if defined($opt_w);

	$opts .= " -o $opt_o" if defined($opt_o);

	open2(\*FROMCLD,\*TOCLD,"LADCP_wspec $opts") ||									# spawn sub-process

		croak("LADCP_wspec $opts: $!\n");

	print(TOCLD $antsOldHeaders); 													# feed already gobbled header 

	for (my($r)=0; $r<@ants_; $r++) {												# feed all .wprof records to LADCP_wspec

		print(TOCLD "@{$ants_[$r]}\n");

    }

	for (my($bufi)=0; defined($ARGV[0]); $bufi++) {									# multiple input files: loop until 2nd last

		$input_list .= "$P{profile_id}($P{run_label}) ";

		if ($bufi == 0) {															# do once for mulitple files

			&antsAddParams('ADCP_bin_length','','ADCP_blanking_distance','',		# delete most %PARAMs, leaving

						   'ADCP_frequency','','ADCP_orientation','',				# 	only potentially useful ones (from last file):

						   'ADCP_pulse_length','','BT_max_bin_range_diff','',		#		%profile_id

						   'BT_max_range','','BT_max_w_error','',					#		%water_depth

						   'BT_rms_w_discrepancy','','CTD_time_lags','',			#		%lat

						   'LADCP_firstBin','','LADCP_lastBin','',					#		%lon

						   'LADCP_w_ocean::version','','SS_min_samp','',			#		%dnXX

						   'SS_max_allowed_depth_range','','SS_min_signal','',

						   'Sv_ref_bin','','TL_max_allowed_three_lag_spread','',

						   'dc_rms_accel_pkg','','dc_rms_tilt','',

						   'uc_rms_accel_pkg','','uc_rms_tilt','',

						   'max_depth','','max_elapsed','','max_ens','',

						   'min_depth','','min_elapsed','','min_ens','',

						   'out_basename','','outgrid_dz','','run_label','',

						   'outgrid_firstbin','','outgrid_lastbin','',

						   'outgrid_minsamp','','per_bin_valid_frac_lim','',

						   'processing_options','','refLr_firstBin','',

						   'refLr_lastBin','','rms_w_reflr_err','',

						   'rms_w_reflr_err_interior','',

						   'sidelobe_editing','','surface_layer_depth','',

						   'vessel_draft','','w_max_lim','',

						   'water_depth.sig','','water_depth_from','',

			);

		}

		close(TOCLD);																# close LADCP_VKE input

    	my(@specrec);

	    while (@specrec = &antsFileIn(FROMCLD)) {

	    	my($i) = ($specrec[$widxf]>0) ? $specrec[$widxf] : 0;

	    	@{$specbuf[$bufi][$i]} = @specrec;

#			print(STDERR "specrec[$i] = @specrec\n");

    	}

    	close(FROMCLD);

#		printf(STDERR "windows up to %d buffered\n",$#{$specbuf[$bufi]});

		$antsBufSkip = @ants_;

    	&antsIn();																	# read next .wprof file

		open2(\*FROMCLD,\*TOCLD,"LADCP_wspec $opts") || 							# process it

			croak("LADCP_wspec $opts: $!\n");

		print(TOCLD $antsOldHeaders);

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

			print(TOCLD "@{$ants_[$r]}\n");

	    }

	}

	close(TOCLD);																	# connect stdout from LADCP_VKE to stdin

	open(STDIN,"<&FROMCLD") || croak("dup(FROMCLD): $!\n");

	close(FROMCLD);

	<>;																				# clear EOF condition

	undef(%P);																		# shouldn't matter, because we'll get the same %PARAMs back

	undef(@antsLayout);																# shouldn't matter, because it will get overwritten

	undef($antsOldHeaders);															# forget those

	undef(@ants_);

} elsif (defined(fnrNoErr('pwrdens.0'))) {

	croak("$0: -d, -u, -b, -w, -s meaningless when $0 used with spectral input\n")

		if ($opt_d || $opt_u || $opt_b || defined($opt_w) || defined($opt_s) || defined($opt_g));

} else {

	if ($ARGV[0]) {

		croak("$ARGV[0]: no such file or directory\n");

	} else {

		croak("$0: empyt input\n");

	}

}

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

# Handle LADCP_VKE usage & read spectra from final file

#	- spectra from previous files are in @specbuf

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

$n_input_files = 1 + @specbuf;									# number of input files provided

&antsAddParams('LADCP_VKE::input_files.n',$n_input_files,

&antsFloatOpt(\$opt_e,$c);										# default parameterization

&antsFloatOpt(\$opt_x,1);										# spectral fit stddev scale factor

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

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

	&antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin);	# ensure eps.VKE is calculated below

} elsif (defined(antsParam('shortwave_cutoff'))) {				# cutoff already applied

	$lmin = 2*$PI/antsParam('shortwave_cutoff');

} else {														# use 100m default cutoff

	$lmin = $opt_c_default;

	&antsAddParams('LADCP_VKE::shortwave_cutoff',2*$PI/$lmin);	# ensure eps.VKE is calculated below

}

$lmax = 9e99;													# no longwave cutoff implemented yet

&antsInstallBufFull(0);											# load entire file

&antsIn();

$P{run_label} = "$input_list$P{profile_id}($P{run_label})"

	if ($n_input_files > 1);

my($Hbuf) = $antsOldHeaders;									# save for later (used on -i)

&antsRequireParam('profile_id');

&antsRequireParam('lambda.1');

&antsRequireParam('k.1');

&antsRequireParam('resolution_bandwidth');

&antsRequireParam('input_depth_resolution');

&antsRequireParam('output_depth_resolution');

&antsRequireParam('ADCP_bin_length');

&antsAddParams('ADCP_pulse_length',antsParam('ADCP_bin_length'))

	unless defined(antsParam('ADCP_pulse_length'));

$imin = 0;														# find frequency bin limits

for ($nfreq=1; defined(antsParam("lambda.$nfreq")); $nfreq++) {

	$imin = $nfreq if ($imin==0 && antsParam("lambda.$nfreq")<=$lmax);

	$imax = $nfreq if (antsParam("lambda.$nfreq") >= $lmin);

}

croak("$0: <lambda.min=$lmin> < min(lambda)")

	unless defined($imax);

$pg_fmin = fnr('pwrdens.0');								# first power field in spectra

$fs_fmin = $pg_fmin + $imin;								# first power field in finescale range

$fs_fmax = $pg_fmin + $imax;								# last power field in finescale range

$widxf = fnr('widx');

$df	   = fnr('depth');

$mindf = fnr('depth.min');

$maxdf = fnr('depth.max');

$wsf   = fnr('w.nsamp.avg');

$doff  = fnr('nspec');

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

# Redirect STDOUT & create plot if STDOUT is a tty

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

if (-t STDOUT) {

	$opt_f = '.' unless defined($opt_f);					# output directory

	croak("$opt_f: not a directory\n") unless (-d $opt_f);

	my($id) = &antsRequireParam('profile_id');

#	if ($opt_d) {

#		$opt_p = sprintf('%s/%03d_dc_VKE.ps',$opt_f,$id) unless defined($opt_p);

#		$outfile = sprintf('%s/%03d_dc.VKE',$opt_f,$id);

#	} elsif ($opt_u) {

#		$opt_p = sprintf('%s/%03d_uc_VKE.ps',$opt_f,$id) unless defined($opt_p);

#		$outfile = sprintf('%s/%03d_uc.VKE',$opt_f,$id);

#	} else {

		$opt_p = sprintf('%s/%03d_VKE.ps',$opt_f,$id) unless defined($opt_p);

		$outfile = sprintf('%s/%03d.VKE',$opt_f,$id);

#	}

	$outfile =~ s@^./@@;

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

} elsif (defined($opt_f)) {

	croak("-f can only be used without STDOUT redirection\n");

}

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

# Library

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

sub average_spectra($)

{

	my($r) = @_;

	for (my($f)=$fs_fmin-1; $f<=$fs_fmax; $f++) {								# average w.nsamp.avg & spectral densities

		my($sum) = my($ns) = 0;

		if (numberp($ants_[$r][$f])) {											# final file has valid spectra

			$sum = $ants_[$r][$f]; $ns = 1;										# NB: nspec is correct even if it doesn't

		}

		my($wi) = ($ants_[$r][$widxf]>0) ? $ants_[$r][$widxf] : 0;				# adjust for -1 

		for (my($bi)=0; $bi<@specbuf; $bi++) {									# loop over all buffered files

			next unless @{$specbuf[$bi][$wi]};									# skip input files w/o valid spectra

			if (abs($specbuf[$bi][$wi][$df] - $ants_[$r][$df]) > 0) {			# depth mismatch

				die("assertion failed") unless ($wi == 0);						# only allowed in bottom window

				if (abs($specbuf[$bi][$wi][$df] - $ants_[$r][$df]) >

					abs($ants_[$r][$maxdf] - $ants_[$r][$mindf])) {

						printf(STDERR "WARNING: ignoring bottom window from input file #%d because of depth mismatch\n",$bi+1)

							if ($f == $fs_fmin);

						next;

				}

			}

#			print(STDERR "specbuf[$bi][$wi][$f] = $specbuf[$bi][$wi][$f]\n");

			if (numberp($specbuf[$bi][$wi][$f])) {

				$sum += $specbuf[$bi][$wi][$f]; $ns++;

            }

			$ants_[$r][$doff] += $specbuf[$bi][$wi][$doff]						# update nspec once per input record

				if ($f == $fs_fmax);											#	... but for all files

		}

		$ants_[$r][$f] = ($ns > 0) ? $sum/$ns : nan;							# update averaged spectral density

	}

}	

sub integrate_fs_power($)														# integrate fs spectrum

{

	my($r) = @_;

	$ants_[$r][$fspwrf] = 0;

	for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {

		$ants_[$r][$fspwrf] += $ants_[$r][$f];

	}

	$ants_[$r][$fspwrf] *= antsParam('resolution_bandwidth');

}

sub fit_universal_w_spec($)														# vertical velocity => p0

{

	my($r) = @_;

	my($nsamp) = $fs_fmax - $fs_fmin + 1;

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

	# fit slope-2 line in log-log space (main estimator)

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

	if ($nsamp >= 2) {															# require min 2 wavenumber samples

		my($DOF) = 0;

		my($sumd,$sumx,$sumy) = (0,0,0);										# fit kz^-2 power law

		for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {

			my($i) = $f - $pg_fmin;

			$sumx += log10(antsParam("k.$i"));

			$sumy += log10($ants_[$r][$f]);

			$sumd += log10($ants_[$r][$f]) + 2*log10(antsRequireParam("k.$i"));

		}

		my($p0) = $sumd/$nsamp;

		$ants_[$r][$p0f] = 10**$p0;

		my($avgx) = $sumx/$nsamp;												# avg for r calc

		my($avgy) = $sumy/$nsamp;

		my($sumsqerr,$sxx,$syy,$sxy,$sumsqxt,$sumx,$sumy,$sumwt) =				# r, rms error, pwrlaw slope

			(0,0,0,0,0,0,0,0);

		for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {								

			my($i)  	= $f - $pg_fmin;

			my($x)  	= log10(&antsParam("k.$i"));

			my($y)  	= log10($ants_[$r][$f]);

			my($ysig) 	= $opt_x * $y / sqrt($ants_[$r][$doff]);

			my($xt) 	= $x - $avgx; $sxx += &SQR($xt);						# correlation coeff (r)

            my($yt) 	= $y - $avgy; $syy += &SQR($yt); $sxy += $xt * $yt;

			my($wt) 	= 1 / &SQR($ysig); $sumwt += $wt;						# slope (linear fit in log-log space)

			$sumx 	   += $x * $wt; 

			$sumy 	   += $y * $wt;

			$sumsqerr  += ($p0 - 2*$x - $y)**2;									# rms error

        }

        my($midx) = $sumx / $sumwt;

#		print(STDERR "$sumx:$sumy:$sumwt\n");

#		print(STDERR "$midx\n");

		my($sumsqdx,$sumslp) = (0,0);

		for (my($f)=$fs_fmin; $f<=$fs_fmax; $f++) {

			my($i)  	= $f - $pg_fmin;

			my($x)      = log10(&antsParam("k.$i"));

			my($y) 		= log10($ants_[$r][$f]);

			my($ysig) 	= $opt_x * $y / sqrt($ants_[$r][$doff]);

	        my($dx)   	= ($x - $midx) / $ysig;

			$sumsqdx   += &SQR($dx);

	        $sumslp    += $dx * $y / $ysig;

        }

        $ants_[$r][$rmsf] 	= sqrt($sumsqerr/$nsamp);

        $ants_[$r][$rf]   	= $sxy/(sqrt($sxx * $syy) + $SMALL_AMOUNT);

        $ants_[$r][$slpf] 	= $sumslp / $sumsqdx;

        $ants_[$r][$sslpf] 	= sqrt(1 / $sumsqdx);

	} else {

		&antsInfo("WARNING: no fit --- need min 2 samples");

	}

}

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

# Load Dissipation Data

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

my(@eps_ms,@depth_ms);										# output variables

if (defined($opt_k)) {

	my($file,$field) = split(':',$opt_k);

	open(ADDF,"$file") || croak("$file: $!\n");				# open file

	my(@afl) = &antsFileLayout(ADDF);						# read layout

	my($akf,$aef);

	for (my($f)=0; $f<=$#afl; $f++) {						# find depth & eps fields

		$akf = $f if ($afl[$f] eq 'depth');

		$aef = $f if ($afl[$f] eq $field);

	}

	croak("$file: fields 'depth' and '$field' required\n")

		unless defined($akf) && defined($aef);

	while (1) {												# load entire profile

		my(@ar);

		last unless @ar = &antsFileIn(ADDF);

		next unless numberp($ar[$akf]) && numberp($ar[$aef]);

		push(@depth_ms,$ar[$akf]); push(@eps_ms,$ar[$aef]);

	}

	close(ADDF);

}

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

# Process File

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

$fspwrf  = &antsNewField('pwr.fs');							# derived fields

$p0f	 = &antsNewField('p0');								# VKE density

$rmsf	 = &antsNewField('p0fit.rms');						# rms misfit

$rf      = &antsNewField('p0fit.r');						# correlation coefficient

$slpf	 = &antsNewField('p0fit.slope');					# power-law slope

$sslpf	 = &antsNewField('p0fit.slope.sig');				# power-law slope stddev

$wepsf   = &antsNewField('eps.VKE');						# epsilon from VKE

$msepsf  = &antsNewField('eps.ms')							# externally supplied microstructure eps

	if defined($opt_k);

my(@outLayout) = @antsNewLayout;							# save for later