# HG changeset patch
# User A.M. Thurnherr <athurnherr@yahoo.com>
# Date 1567029094 14400
# Node ID 4ccccbf69dfdbc21883058aaf1380aab5f7680a1
# Parent  0f6d9e64cc4fd2c6566ab4446dd70fc895a3a8ef
working?

diff -r 0f6d9e64cc4f -r 4ccccbf69dfd LADCP_w_CTD
--- a/LADCP_w_CTD	Tue Aug 27 19:11:54 2019 -0400
+++ b/LADCP_w_CTD	Wed Aug 28 17:51:34 2019 -0400
@@ -2,9 +2,9 @@
 #======================================================================
 #                    L A D C P _ W _ C T D 
 #                    doc: Mon Nov  3 17:34:19 2014
-#                    dlm: Mon Apr 29 18:00:44 2019
+#                    dlm: Wed Aug 28 17:47:08 2019
 #                    (c) 2014 A.M. Thurnherr
-#                    uE-Info: 111 15 NIL 0 0 72 2 2 4 NIL ofnI
+#                    uE-Info: 296 0 NIL 0 0 72 2 2 4 NIL ofnI
 #======================================================================
 
 $antsSummary = 'pre-process SBE 9plus CTD data for LADCP_w';
@@ -87,6 +87,8 @@
 #				  - BUG: ITS was not set. How is this possible?????
 #	Apr 21, 2019: - modified code to allow production of 24Hz files (previous code required
 #					min 2 samples per bin, allowing for max 12Hz sampling rate)
+#	Aug 27, 2019: - began adding correction for dropped CTD scans
+#	Aug 28, 2019: - made it work
 
 # NOTES:
 #	w_CTD is positive during the downcast to make the sign of the apparent
@@ -150,8 +152,6 @@
 	($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 s\n")
-#		unless (abs($sampint-1/24) < 1e-5);
 	_croak("$CNVfile: insufficient time resolution ($sampint s) for ${opt_s}Hz time series\n")
 		if (round(1/$sampint/$opt_s) < 1);
 
@@ -193,6 +193,15 @@
 			$condHistRes = 2;
 		}
 	}
+
+	$systimeF = &fnrNoErr('timeY');									# apply correction for dropped scans
+	$xmerrF  = &fnrNoErr('modError');
+	if (defined($systimeF) && defined($xmerrF)) {
+		$fill_gaps = 1;
+	} else {
+		print(STDERR "\n\n") if ($opt_v > 1);
+		print(STDERR "WARNING: timeY and/or modError missing from $CNVfile -- cannot correct for CTD modulo errors\n");
+    }	
 	
 	$latF = &fnrNoErr('lat');										# GPS data if available (to make files useful for u/v processing)
 	$lonF = &fnrNoErr('lon');
@@ -230,11 +239,123 @@
 	}
 }
 
-printf(STDERR "\n\t%d scans (%d seconds)",$nscans,int($nscans*$sampint))
+printf(STDERR "\n\t%d scans (%d minutes)",$nscans,round($nscans*$sampint/60))
 	if ($opt_v > 1);
 printf(STDERR "\n") if ($opt_v);
 
 #----------------------------------------------------------------------
+# Fill gaps in CTD time series due to dropped scans (modulo errors)
+#----------------------------------------------------------------------
+
+if ($fill_gaps) {
+	print(STDERR "Filling CTD time-series gaps...") if ($opt_v);
+
+	my($scans_filled) = 0;
+	my($scans_replaced) = 0;
+	my($scans_deleted) = 0;
+
+	for (my($scan)=30; $scan<@ants_; $scan++) {										# start a bit more than 1 second into the cast
+		next if ($ants_[$scan][$systimeF] == $ants_[$scan-1][$systimeF]);			# skip forward to next systime second
+
+#		while ($ants_[$scan][$systimeF] > $ants_[$scan-1][$systimeF]+1) {			# gap spans full second
+#			my(@splicescan);														# scan to splice in 
+#			$splicescan[$systimeF] = $ants_[$scan-1][$systimeF] + 1;
+#			$splicescan[$xmerrF] = $ants_[$scan-1][$xmerrF];
+#			for (my($i)=0; $i<24; $i++) {
+#				splice(@ants_,$scan,0,\@splicescan);
+#				$scans_filled++;
+#			}
+#		}
+			
+		my($second_start) = $scan - 1;												# backtrack to beginning of second
+		my($scans_this_sec) = 1;
+		while ($ants_[$second_start-1][$systimeF] == $ants_[$second_start][$systimeF]) {
+			$second_start--; $scans_this_sec++;
+		}
+		die("\nscan = $scan, second_start = $second_start, st = $ants_[$second_start][$systimeF]")
+			unless ($second_start > 0);
+
+		while ($scans_this_sec > 24) {												# CTD clock running fast => remove scans
+			splice(@ants_,$scan-1,1);
+			$scans_deleted++; $scans_this_sec--; $scan--;
+		}
+
+		my($gap_len) = 24 - $scans_this_sec;										# gap length in this second only
+		die("gap_len = $gap_len at scan#$scan") unless ($gap_len >= 0);
+		
+		my($ngaps) = 0;																# gaps between end of previous and beginning of next section
+		for (my($s)=$second_start; $s<=$scan; $s++) {
+			$ngaps++ if ($ants_[$s][$xmerrF] > $ants_[$s-1][$xmerrF]);
+		}
+
+		my(@splicescan);															# scan to splice in 
+		$splicescan[$systimeF] = $ants_[$second_start][$systimeF];
+		$splicescan[$xmerrF] = $ants_[$second_start][$xmerrF];
+
+		if ($ngaps==0 && $gap_len>0) {												# clock drift => add scan
+#			print(STDERR "adding $gap_len scans to make up for clock drift starting at scan $scan-1 (ss = $second_start) (st = $splicescan[$systimeF])\n");
+			for (my($i)=0; $i<$gap_len; $i++) {										# fill gap
+				splice(@ants_,$scan-1,0,\@splicescan);
+				$scans_filled++; $scan++;
+			}
+        } elsif ($ngaps == 1) {														# single gap => fill it
+			my($gap_location) = $second_start;										# find it
+			while ($ants_[$gap_location][$xmerrF] == $ants_[$second_start-1][$xmerrF]) {
+				$gap_location++;
+			}
+			die('INTERNAL ERROR') unless ($gap_location <= $scan);
+#			print(STDERR "filling $gap_len gap starting at scan $gap_location (ss = $second_start) (st = $splicescan[$systimeF])\n");
+			for (my($i)=0; $i<$gap_len; $i++) {										# fill gap
+				splice(@ants_,$gap_location,0,\@splicescan);
+				$scans_filled++; $scan++;
+			}
+		} else {																	# multiple gaps => delete data
+			my($last_gap);
+			for (my($s)=$second_start; $s<$scan; $s++) {							# replace ambiguous data scans
+				next if ($ants_[$s][$xmerrF] == $ants_[$second_start-1][$xmerrF]);
+				last if ($ants_[$s][$xmerrF] == $ants_[$scan][$xmerrF]);
+				splice(@ants_,$s,1,\@splicescan);									# replace
+				$scans_replaced++;
+				$last_gap = $s;
+            }		
+			for (my($i)=0; $i<$gap_len; $i++) {										# fill gaps
+				splice(@ants_,$last_gap+1,0,\@splicescan);
+				$scans_filled++; $scan++;
+			}
+		}
+
+		my($scans_this_sec) = 0;
+		die("$ants_[$second_start][$systimeF]==$ants_[$second_start-1][$systimeF]+1")
+			unless ($ants_[$second_start][$systimeF]==$ants_[$second_start-1][$systimeF]+1);
+		for (my($i)=0; $ants_[$second_start+$i][$systimeF]==$ants_[$second_start-1][$systimeF]+1; $i++) {
+			$scans_this_sec++;
+		}
+		die("\n$scans_this_sec scans in second $ants_[$second_start][$systimeF] at scans $second_start .. $scan-1 (gap_len $gap_len, $ngaps gaps)\n")
+			unless ($scans_this_sec == 24);
+	} # for ($scan
+
+#	my($sts) = 1;
+#	for (my($scan)=1; $scan<@ants_; $scan++) {
+#		$sts++,next if ($ants_[$scan][$systimeF] == $ants_[$scan-1][$systimeF]);
+#		printf(STDERR "sts = $sts beginning at scan %d (st = $ants_[$scan-1][$systimeF])\n",$scan-$sts)
+#			unless ($sts == 24);
+#		$sts = 1;
+#	}
+		
+
+	&antsAddParams('gaps_scans_filled',$scans_filled,
+				   'scans_removed',$scans_deleted,
+				   'gap_scans_cleared',$scans_replaced);
+
+	if ($opt_v>1 && $scans_filled) {
+		printf(STDERR "\n\t%d scans removed (clock drift)",$scans_deleted);
+		printf(STDERR "\n\t%d scans added (gaps & clock drift)",$scans_filled);
+		printf(STDERR "\n\t%d scans cleared (gap clusters)",$scans_replaced);
+	}
+	printf(STDERR "\n") if ($opt_v);
+} # if ($fill_gaps
+
+#----------------------------------------------------------------------
 # Redirect STDOUT to %.6Hz & create %_w_CTD.ps,%_sspd.ps if STDOUT is a tty
 #----------------------------------------------------------------------
 
@@ -309,6 +430,7 @@
 	my($press_rez) = 2;
 	my($outliers) = my($modeSamp) = 0; my($modeBin,$min,$max); local(@hist);
 	for (my($s)=0; $s<$nscans; $s++) {
+		$ants_[$s][$pressF] = nan unless defined($ants_[$s][$pressF]);
 		next unless ($ants_[$s][$pressF]>=-100 && $ants_[$s][$pressF]<=6500);
 		my($b) = ($ants_[$s][$pressF]+100) / $press_rez;
 		$hist[$b]++;
@@ -320,6 +442,7 @@
 	($min,$max) = validRange($modeBin);
 	$min = $press_rez*$min-100; $max = $press_rez*$max-100;
 	for (my($s)=0; $s<$nscans; $s++) {
+		next unless numberp($ants_[$s][$pressF]);
 		if ($ants_[$s][$pressF] > $max) { $outliers++; $ants_[$s][$pressF] = nan; }
 		if ($ants_[$s][$pressF] < $min) { $outliers++; $ants_[$s][$pressF] = nan; }
 	}
@@ -334,6 +457,7 @@
 	$outliers = $modeSamp = 0;
 	undef(@hist);
 	for (my($s)=0; $s<$nscans; $s++) {
+		$ants_[$s][$condF] = nan unless defined($ants_[$s][$condF]);
 		next unless ($ants_[$s][$condF] > 0);
 		my($b) = $ants_[$s][$condF]*$condHistRes;					# 1/10 S/m histogram resolution (1 mS/cm)
 		$hist[$b]++;
@@ -343,6 +467,7 @@
 	($min,$max) = validRange($modeBin);
 	$min /= $condHistRes; $max /= $condHistRes;
 	for (my($s)=0; $s<$nscans; $s++) {
+		next unless numberp($ants_[$s][$condF]);
 		if ($ants_[$s][$condF] > $max) { $outliers++; $ants_[$s][$condF] = nan; }
 		if ($ants_[$s][$condF] < $min) { $outliers++; $ants_[$s][$condF] = nan; }
 	}
@@ -362,6 +487,7 @@
 	$outliers = $modeSamp = 0;
 	undef(@hist);
 	for (my($s)=0; $s<$nscans; $s++) {
+		$ants_[$s][$tempF] = nan unless defined($ants_[$s][$tempF]);
 		next unless ($ants_[$s][$tempF] >= -10);
 		my($b) = ($ants_[$s][$tempF] + 10) * 0.5;
 		$hist[$b]++;
@@ -374,6 +500,7 @@
 	$min = ($min / 0.5) - 10;
 	$max = ($max / 0.5) - 10;
 	for (my($s)=0; $s<$nscans; $s++) {
+		next unless numberp($ants_[$s][$tempF]);
 		if ($ants_[$s][$tempF] > $max) { $outliers++; $ants_[$s][$tempF] = nan; }
 		if ($ants_[$s][$tempF] < $min) { $outliers++; $ants_[$s][$tempF] = nan; }
 	}