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

#                    L A D C P P R O C . U H C O D E 

#                    doc: Fri Sep 17 20:27:53 2010

#                    (c) 2010 A.M. Thurnherr & E. Firing

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

# PERLified functions from Eric's [merge.c]; with mods

# NOTES:

#	- velocity integration removed

#	- percent-good flag removed (no point for single-ping ensembles)

#	- need the following per-ensemble fields from previous steps:

#		CTD_DEPTH

#		CTD_SVEL

#	- negative depths are not allowed (should not happen given CTD_DEPTH)

# WEIRDNESSES:

#	- w reference layer in set_misc_flags is calculated without taking

#	  the ref layer bins into account

#	- u,v ref lr vels in set_wake_flag use w ref layer bins

#	- distance to bin 1 center is not sound-speed corrected [CORRECTED]

#	- $tilt calculation is wrong

# HISTORY:

#	Sep 17, 2010: - created

#	Oct 13, 2010: - first working version

#	Oct 14, 2010: - renamed from LADCPshear.UHcode

#	Oct 15, 2010: - added support for -pPI edit suppresion

#	Oct 19, 2010: - reversed semantics of -p

#	Oct 25, 2010: - added W_CTD_BIT, renamed W_BIT to W_OUTLIER_BIT

#	Oct 26, 2010: - added TILT_BIT

#	Dec 10, 2010: - modified assertion to allow processing of UL data

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

# VELOCITY EDITING

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

my($BADVEL_BIT)  	= 0x01;

my($ERRVEL_BIT)  	= 0x02;

my($CORREL_BIT)		= 0x04;

my($W_OUTLIER_BIT) 	= 0x08;

my($SHEAR_BIT) 		= 0x10;

my($SIDELOBE_BIT)	= 0x20;

my($WAKE_BIT)		= 0x40;

my($PPI_BIT)		= 0x80;

my($W_CTD_BIT)		= 0x100;

my($TILT_BIT)		= 0x200;

my($DELTA_TILT_BIT)	= 0x400;

my(%flag_count);

# apparently, in Eric's code, DISTANCE_TO_BIN1_CENTER is not sound-speed corrected

sub dzToBin($$)

{

	my($ens,$bin) = @_;

	my($sscorr) = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL} / $LADCP{ENSEMBLE}[$ens]->{SPEED_OF_SOUND};

	return $sscorr * ($LADCP{DISTANCE_TO_BIN1_CENTER} + $bin*$LADCP{BIN_LENGTH});

}

sub depthOfBin($$)

{

	my($ens,$bin) = @_;

	return $LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_UP} ?

		   $LADCP{ENSEMBLE}[$ens]->{DEPTH} - &dzToBin($ens,$bin) :

		   $LADCP{ENSEMBLE}[$ens]->{DEPTH} + &dzToBin($ens,$bin);

}

sub set_wake_flag($$)

{

	my($ens,$De) = @_;

	my($n) = 0;

	my($uref) = my($vref) = my($wref) = 0;

	for (my($bin)=$wbin_start-1; $bin<$wbin_end; $bin++) {		# calc crude ref lr vel from all data

		next if ($edit_flags[$ens][$bin]);

		$uref += $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$U];

		$vref += $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$V];

		$wref += $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W];

		$n++;

	}

	return if ($n==0);

	$uref /= $n;

	$vref /= $n;

	$wref /= $n;

	## if upward (=negative) velocity greater than minimum, calculate wake

	## 		heading and inclination

	if ($wref < -$min_wake_w) {

		my($wake_hd) = 180 / 3.14159265358979 * atan2($uref,$vref);

		my($speed) 	 = sqrt($uref*$uref + $vref*$vref);

		my($wake_ang)= abs(180 / 3.14159265358979 * atan($speed/$wref));

		$wake_hd += 360

			if ($wake_hd < 0);

		$LADCP{ENSEMBLE}[$ens]->{HEADING} += 360

			if ($LADCP{ENSEMBLE}[$ens]->{HEADING} < 0);

		my($wake_mod) = $wake_hd % 90;		# % returns integer part of remainder, but that's sufficient

		my($adcp_mod) = $LADCP{ENSEMBLE}[$ens]->{HEADING} % 90;

		if (((abs($wake_mod - $adcp_mod) < $wake_hd_dif) ||

             (abs($wake_mod - $adcp_mod) > (90 - $wake_hd_dif))) &&

			 ($wake_ang > $wake_ang_min)) {

			for (my($bin)=0; $bin<$n_wake_bins; $bin++) {

				$flag_count{$WAKE_BIT}[$bin]++;

				$edit_flags[$ens][$bin] |= $WAKE_BIT;

			}

		}

	} ## if ($wref < -min_wake_w)

	## This does not make a lot of sense, because it trims points

	## on only one side of the wake error, and that side depends

	## on whether the integration is forward or backward.

	if (($edit_flags[$ens+$De][0]&$WAKE_BIT) &&

		($edit_flags[$ens][0]&$WAKE_BIT == 0)) {

		for (my($bin)=0; $bin<$n_wake_bins; $bin++) {

			$flag_count{$WAKE_BIT}[$bin]++;

			$edit_flags[$ens][$bin] |= $WAKE_BIT;

		}

	}

}

sub set_misc_flags($$)

{

	my($ens,$De) = @_;

	my($ww) = my($n) = 0;

	my($SLIfac) = 1 - cos(rad($LADCP{BEAM_ANGLE}));

	for (my($bin)=0; $bin<$w_ref_bin; $bin++) {				# ref-lr w

   		if (!$edit_flags[$ens][$bin]) {

   			$ww += $LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W];

			$n++;

		}

	}

	$ww /= $n if ($n > 0);

	for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {

		next if ($edit_flags[$ens][$bin]);

		## We use the standard criterion for bottom interference; e.g. for

		## 30-degree beams, the last 15% of the profile is contaminated

		## by the sidelobe bounce.	1.5 bin length is added to allow for

		## the length of the bin and pulse, that is, contamination of part of a

		## bin.  Profiler tilt does not require a more stringent criterion.

		if ($LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_DOWN}) {

			if (numberp($water_depth) &&

				$water_depth - &depthOfBin($ens,$bin) <=

					$SLIfac * ($water_depth - $LADCP{ENSEMBLE}[$ens]->{DEPTH})

						+ 1.5 * $LADCP{BIN_LENGTH}) {

				$edit_flags[$ens][$bin] |= $SIDELOBE_BIT;

				$flag_count{$SIDELOBE_BIT}++;

			}

		} else { ## upward-looking

			if (&depthOfBin($ens,$bin) <=

				$SLIfac * $LADCP{ENSEMBLE}[$ens]->{DEPTH}

					+ 1.5 * $LADCP{BIN_LENGTH}) {

				$edit_flags[$ens][$bin] |= $SIDELOBE_BIT;

				$flag_count{$SIDELOBE_BIT}++;

			}

		}

		if ($ww != 0 &&

			abs($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W] - $ww) > $w_dif) {

			$flag_count{$W_OUTLIER_BIT}++;

			$edit_flags[$ens][$bin] |= $W_OUTLIER_BIT;

		}

		if (abs($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$E]) > $e_max) {

			$flag_count{$ERRVEL_BIT}++;

			$edit_flags[$ens][$bin] |= $ERRVEL_BIT;

		}

		for (my($beam)=0; $beam<=3; $beam++) {

			if ($LADCP{ENSEMBLE}[$ens]->{CORRELATION}[$bin][$beam] < $min_cor) {

				$flag_count{$CORREL_BIT}++;

				$edit_flags[$ens][$bin] |= $CORREL_BIT;

			}

		}

		if ($bin < $shbin_start-1 || $bin >= $shbin_end) {

			$edit_flags[$ens][$bin] |= $BADVEL_BIT;

			$flag_count{$BADVEL_BIT}++;

		}

	} # for ($bin=0...

}

## The following is for editing out the second bottom bounce.

sub set_PPI_flags($$)

{

	my($ens,$De) = @_;

	my($dt_ping) = $LADCP{ENSEMBLE}[$ens]->{UNIX_TIME} - $LADCP{ENSEMBLE}[$ens-1]->{UNIX_TIME};

	if ($LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_DOWN}) {

		if (numberp($water_depth)) {

			$clip_z1 = $water_depth

						- $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL}/2 * $dt_ping

							* cos(rad($LADCP{BEAM_ANGLE} + $LADCP{ENSEMBLE}[$ens]->{TILT}))

						+ $clip_margin;

			$clip_z0 = $water_depth

						- $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL}/2 * $dt_ping

							* cos(rad($LADCP{BEAM_ANGLE} - $LADCP{ENSEMBLE}[$ens]->{TILT}))

	                    - $clip_margin;

	    }

	} else { # upward-looking

		$clip_z1 = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL}/2 * $dt_ping

						* cos(rad($LADCP{BEAM_ANGLE} - $LADCP{ENSEMBLE}[$ens]->{TILT}))

					+ $clip_margin;

		$clip_z0 = $LADCP{ENSEMBLE}[$ens]->{CTD_SVEL}/2 * $dt_ping

						* cos(rad($LADCP{BEAM_ANGLE} + $LADCP{ENSEMBLE}[$ens]->{TILT}))

					- $clip_margin;

	}

	for (my($bin)=$first_clip_bin-1; $bin<$LADCP{N_BINS}; $bin++) {

		next if ($edit_flags[$ens][$bin]);

		my($dob) = depthOfBin($ens,$bin);

		if ($dob >= $clip_z0 && $dob <= $clip_z1) {

			$edit_flags[$ens][$bin] |= $PPI_BIT;

			$flag_count{$PPI_BIT}++;

		}

	}

}

sub edit_velocity($$)

{

	my($start,$end) = @_;												# ensemble indices

	my($De) = $start<$end ? 1 : -1;										# downcast: De = 1; upcast: De = -1

	$flag_count{$WAKE_BIT} = $flag_count{$W_OUTLIER_BIT} = $flag_count{$ERRVEL_BIT} =

	$flag_count{$CORREL_BIT} = $flag_count{$SHEAR_BIT} = $flag_count{$BADVEL_BIT} =

	$flag_count{$SIDELOBE_BIT} = $flag_count{$PPI_BIT} = $flag_count{$W_CTD_BIT} =

	$flag_count{$TILT_BIT} = $flag_count{$DELTA_TILT_BIT} = 0;

	for (my($ens)=$start; $ens!=$end+$De; $ens+=$De) {					# loop over all ens from start to end

		if (abs($LADCP{ENSEMBLE}[$ens]->{CTD_W}-$LADCP{ENSEMBLE}[$ens]->{W}) > $w_max_err) {	# get rid of aliased vels (ambiguity)

			for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {

				$edit_flags[$ens][$bin] |= $W_CTD_BIT;

				$flag_count{$W_CTD_BIT}++;

			}

			next;

		}

		if ($LADCP{ENSEMBLE}[$ens]->{TILT} > $max_tilt) {				# get rid ensembles with large tilt

			for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {

				$edit_flags[$ens][$bin] |= $TILT_BIT;

				$flag_count{$TILT_BIT}++;

			}

			next;

		}																# get rid ensembles after large rotation

		if (abs($LADCP{ENSEMBLE}[$ens]->{TILT}-$LADCP{ENSEMBLE}[$ens-$De]->{TILT}) > $max_delta_tilt) {

			for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) {

				$edit_flags[$ens][$bin] |= $DELTA_TILT_BIT;

				$flag_count{$DELTA_TILT_BIT}++;

			}

			next;

		}

		for (my($bin)=$shbin_start-1; $bin<$shbin_end; $bin++) {		# flag bad velocities

			$edit_flags[$ens][$bin] |= $BADVEL_BIT,$flag_count{$BADVEL_BIT}++

   				unless defined($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$bin][$W]);

		}

		set_wake_flag($ens,$De);

		set_misc_flags($ens,$De);

		set_PPI_flags($ens,$De)

    }

}

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

# CALCULATE VELOCITY SHEAR

#	- output in @ush_mu,@vsh_mu,@wsh_mu,@ush_sig,@vsh_sig,@wsh_sig

# 	- @sh_i0, @sh_i1, @dsh, @ush, @vsh, @wsh are defined "local" in calc_shear

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

sub uv_to_shear($)

{

	my($ens) = @_;

	my($nvel) = 0;

	@sh_i0 = @sh_i1 = ();

	for (my($bin)=0; $bin<$LADCP{N_BINS}; $bin++) { 

		next if ($edit_flags[$ens][$bin]);

		$nvel++;

		push(@sh_i1,$bin) if (@sh_i0);

		push(@sh_i0,$bin) if ($bin < $LADCP{N_BINS}-1);

    }

	@dsh = ();

	for (my($i)=0; $i<@sh_i1; $i++) {								# calc and bin shears

		my($d0) = &depthOfBin($ens,$sh_i0[$i]);

		my($d1) = &depthOfBin($ens,$sh_i1[$i]);

		die("$0: assertion failed (ens=$ens i=$i depth=$LADCP{ENSEMBLE}[$ens]->{DEPTH} sh_i0[$i]=$sh_i0[$i] sh_i1[$i]=$sh_i1[$i] d0=$d0 d1=$d1)")

			unless ($d0>=0 && $d1>=0);

		die("$0: assertion failed (ens=$ens i=$i sh_i0[$i]=$sh_i0[$i] sh_i1[$i]=$sh_i1[$i] d0=$d0 d1=$d1)")

			unless (($LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_DOWN} && $d1-$d0>0) ||

					($LADCP{ENSEMBLE}[$ens]->{XDUCER_FACING_UP}   && $d1-$d0<0));

		$dsh[$i] = ($d1 + $d0) / 2;

		$ush[$i] = ($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i1[$i]][$U] -

					$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i0[$i]][$U]) / ($d1-$d0);

		$vsh[$i] = ($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i1[$i]][$V] -

					$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i0[$i]][$V]) / ($d1-$d0);

		$wsh[$i] = ($LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i1[$i]][$W] -

					$LADCP{ENSEMBLE}[$ens]->{VELOCITY}[$sh_i0[$i]][$W]) / ($d1-$d0);

	}

	return $nvel;

}

# here ush_mu, sh_n, etc are still set from the pre-gridding pass

sub set_shear_flag($)

{

	my($ens) = @_;

	my(@ibad,@ush_dev,@vsh_dev);

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

		die("$0: assertion failed") unless numberp($dsh[$i]);

		my($bsgi) = int($dsh[$i] / $SHEAR_PREGRID_DZ);

		$ush_dev[$i] = $ush[$i] - $ush_mu[$bsgi];

		$vsh_dev[$i] = $vsh[$i] - $vsh_mu[$bsgi];

		push(@ibad,$i) if ($ush_sig[$i] > 0 &&

							(abs($ush_dev[$i]/$ush_sig[$bsgi]) > $max_shdev ||

               				 abs($vsh_dev[$i]/$vsh_sig[$bsgi]) > $max_shdev));

	} ## end of loop through shears

	## Look for internal glitches: a positive shear followed

	## immediately by a compensating negative shear, for

	## example.  When one is found, flag the common velocity

	## sample, and untag the two shears by setting each ibad

	## to -1.

	for (my($bi)=0; $bi<@ibad-1; $bi++) {

		next unless ($ibad[$bi]+1 == $ibad[$bi+1]);

		my($i) = $ibad[$bi];

		my($bsgi) = int($dsh[$i] / $SHEAR_PREGRID_DZ);

		if ($ush_sig[$bsgi] > 0 && $vsh_sig[$bsgi] > 0 && 

			sqrt(($ush_dev[$i]+$ush_dev[$i+1])**2/$ush_sig[$bsgi]**2) < $max_shdev_sum &&

			sqrt(($vsh_dev[$i]+$vsh_dev[$i+1])**2/$vsh_sig[$bsgi]**2) < $max_shdev_sum) {

				$flag_count{$SHEAR_BIT}++;

				$edit_flags[$ens][$sh_i1[$i]] |= $SHEAR_BIT;

				$ibad[$bi] = $ibad[$bi+1] = -1;

		}

	} ## end of first loop through bad shears

	## Now flag all remaining velocities involved in the shears

	## listed by ibad.

	for (my($bi)=0; $bi<@ibad; $bi++) {

		next if ($ibad[$bi] < 0);

		$flag_count{$SHEAR_BIT} += 2;

		$edit_flags[$ens][$sh_i0[$ibad[$bi]]] |= $SHEAR_BIT;

		$edit_flags[$ens][$sh_i1[$ibad[$bi]]] |= $SHEAR_BIT;

	}

}

sub calc_shear($$$$)

{

	my($start,$end,$grid_dz,$edit_shear) = @_;

	my($De) = $start<$end ? 1 : -1;										# downcast: De = 1; upcast: De = -1

	my($nvel,$nsh) = (0,0);

	for (my($ens)=$start; $ens!=$end+$De; $ens+=$De) {					# loop over all ens from start to end

		local(@sh_i0,@sh_i1);

		uv_to_shear($ens);

		if ($edit_shear) {

			set_shear_flag($ens);

			$nvel += uv_to_shear($ens);

			$nsh += @dsh;

		}

		for (my($i)=0; $i<@dsh; $i++) {									# save shears for binning calculations

			my($gi) = int($dsh[$i] / $grid_dz);

			push(@{$ush_vals[$gi]},$ush[$i]);

			push(@{$vsh_vals[$gi]},$vsh[$i]);

			push(@{$wsh_vals[$gi]},$wsh[$i]);

        }			

	} # $ens loop

	@ush_mu  = @vsh_mu  = @wsh_mu  = ();

	@ush_sig = @vsh_sig = @wsh_sig = ();

	for (my($gi)=0; $gi<@ush_vals; $gi++) {								# calc grid means & stddev

		my($sum_ush,$sum_vsh,$sum_wsh);

		$sh_n[$gi] = @{$ush_vals[$gi]};

		for (my($vi)=0; $vi<$sh_n[$gi]; $vi++) {

			$sum_ush += $ush_vals[$gi][$vi];

			$sum_vsh += $vsh_vals[$gi][$vi];

			$sum_wsh += $wsh_vals[$gi][$vi];

		}

		$ush_mu[$gi] = $sh_n[$gi] ? $sum_ush/$sh_n[$gi] : nan;

		$vsh_mu[$gi] = $sh_n[$gi] ? $sum_vsh/$sh_n[$gi] : nan;

		$wsh_mu[$gi] = $sh_n[$gi] ? $sum_wsh/$sh_n[$gi] : nan;

	}

	for (my($gi)=0; $gi<@ush_vals; $gi++) {								# calc & grid stddevs

		my($sumsq_ush,$sumsq_vsh,$sumsq_wsh);

		for (my($vi)=0; $vi<$sh_n[$gi]; $vi++) {

			$sumsq_ush += ($ush_vals[$gi][$vi] - $ush_mu[$gi])**2;

			$sumsq_vsh += ($vsh_vals[$gi][$vi] - $vsh_mu[$gi])**2;

			$sumsq_wsh += ($wsh_vals[$gi][$vi] - $wsh_mu[$gi])**2;

		}

		$ush_sig[$gi] = $sh_n[$gi]>1 ? sqrt($sumsq_ush/($sh_n[$gi]-1)) : nan;

		$vsh_sig[$gi] = $sh_n[$gi]>1 ? sqrt($sumsq_vsh/($sh_n[$gi]-1)) : nan;

		$wsh_sig[$gi] = $sh_n[$gi]>1 ? sqrt($sumsq_wsh/($sh_n[$gi]-1)) : nan;

	}

	if ($edit_shear && $opt_d) {

		print(STDERR "\n\t\t$nvel valid velocities");

		print(STDERR "\n\t\t$nsh valid shears");

		print(STDERR "\n\t\tflag counts:");

		print(STDERR "\n\t\t\tBADVEL_BIT     = $flag_count{$BADVEL_BIT}");

		print(STDERR "\n\t\t\tERRVEL_BIT     = $flag_count{$ERRVEL_BIT}");

		print(STDERR "\n\t\t\tCORREL_BIT     = $flag_count{$CORREL_BIT}");

		print(STDERR "\n\t\t\tW_OUTLIER_BIT  = $flag_count{$W_OUTLIER_BIT}");

		print(STDERR "\n\t\t\tSHEAR_BIT      = $flag_count{$SHEAR_BIT}");

	    print(STDERR "\n\t\t\tSIDELOBE_BIT   = $flag_count{$SIDELOBE_BIT}");

		print(STDERR "\n\t\t\tWAKE_BIT       = $flag_count{$WAKE_BIT}");

	    print(STDERR "\n\t\t\tPPI_BIT        = $flag_count{$PPI_BIT}");

	    printf(STDERR "\n\t\t\tW_CTD_BIT      = $flag_count{$W_CTD_BIT} (%d ensembles)",

														$flag_count{$W_CTD_BIT}/$LADCP{N_BINS});

	    printf(STDERR "\n\t\t\tTILT_BIT       = $flag_count{$TILT_BIT} (%d ensembles)",

														$flag_count{$TILT_BIT}/$LADCP{N_BINS});

	    printf(STDERR "\n\t\t\tDELTA_TILT_BIT = $flag_count{$DELTA_TILT_BIT} (%d ensembles)",

														$flag_count{$DELTA_TILT_BIT}/$LADCP{N_BINS});

	}

}

1;