1 #!/usr/bin/perl |
1 #!/usr/bin/perl |
2 #====================================================================== |
2 #====================================================================== |
3 # L A D C P _ W |
3 # L A D C P _ W |
4 # doc: Fri Dec 17 18:11:13 2010 |
4 # doc: Fri Dec 17 18:11:13 2010 |
5 # dlm: Thu Apr 16 10:24:42 2015 |
5 # dlm: Thu Apr 16 14:27:45 2015 |
6 # (c) 2010 A.M. Thurnherr |
6 # (c) 2010 A.M. Thurnherr |
7 # uE-Info: 159 70 NIL 0 0 72 2 2 4 NIL ofnI |
7 # uE-Info: 162 60 NIL 0 0 72 2 2 4 NIL ofnI |
8 #====================================================================== |
8 #====================================================================== |
9 |
9 |
10 # TODO: |
10 # TODO: |
11 # detection so that editing stats make sense |
11 # detection so that editing stats make sense |
12 # - own seabed detection (P403) |
12 # - own seabed detection (P403) |
144 # - added editPPI() |
144 # - added editPPI() |
145 # Jul 6, 2014: - BUG: nan water depth had been interpreted as known |
145 # Jul 6, 2014: - BUG: nan water depth had been interpreted as known |
146 # - BUG: sVelProf[] was not allowed to have any gaps |
146 # - BUG: sVelProf[] was not allowed to have any gaps |
147 # Jul 9, 2014: - BUG: Jul 6 bug fixes had been applied to older |
147 # Jul 9, 2014: - BUG: Jul 6 bug fixes had been applied to older |
148 # version |
148 # version |
149 # - BUG: code meant to ensure gap-freee svel profiles did not work correctly |
149 # - BUG: code meant to ensure gap-free svel profiles did not work correctly |
150 # Jul 12, 2014: - finally made output files executable |
150 # Jul 12, 2014: - finally made output files executable |
151 # Apr 5, 2015: - added check for required software |
151 # Apr 5, 2015: - added check for required software |
152 # - BUG: removed dc/uc mean w fields from .prof again |
152 # - BUG: removed dc/uc mean w fields from .prof again |
153 # Apr 7, 2015: - made LADCP_w callable from installation directory |
153 # Apr 7, 2015: - made LADCP_w callable from installation directory |
154 # - BUG: -v default was wrong in usage message |
154 # - BUG: -v default was wrong in usage message |
155 # - replaced 'ens' in output files by 'ensemble' |
155 # - replaced 'ens' in output files by 'ensemble' |
156 # Apr 16, 2015: - turned output specifies into lists (re-design of |
156 # Apr 16, 2015: - turned output specifies into lists (re-design of |
157 # plotting sub-system) |
157 # plotting sub-system) |
158 # - removed 30s sleep from PostProcess.sh call |
158 # - removed 30s sleep from PostProcess.sh call |
159 # - disabled active output when ANTS are not available |
159 # - disabled active output when ANTS are not available |
|
160 # - removed /bin/ksh requirement |
|
161 # - BUG: error messages were not reported in the log file |
|
162 # - made seabed detection code more flexible |
160 |
163 |
161 # CTD REQUIREMENTS |
164 # CTD REQUIREMENTS |
162 # - elapsed elapsed seconds; see note below |
165 # - elapsed elapsed seconds; see note below |
163 # - depth |
166 # - depth |
164 # - ss sound speed |
167 # - ss sound speed |
213 ($ADCP_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$}); |
216 ($ADCP_TOOLS) = (`which mkProfile` =~ m{^(.*)/[^/]*$}); |
214 ($WCALC) = ($0 =~ m{^(.*)/[^/]*$}); |
217 ($WCALC) = ($0 =~ m{^(.*)/[^/]*$}); |
215 $WCALC = '.' if ($WCALC eq ''); |
218 $WCALC = '.' if ($WCALC eq ''); |
216 $ANTS_TOOLS_AVAILABLE = (`which list` ne ''); |
219 $ANTS_TOOLS_AVAILABLE = (`which list` ne ''); |
217 |
220 |
218 die("$0: Korn shell (/bin/ksh) required but not found\n") |
|
219 unless (-x '/bin/ksh'); |
|
220 die("$0: Generic Mapping Tools (GMT) required but not found (bad \$PATH?)\n") |
221 die("$0: Generic Mapping Tools (GMT) required but not found (bad \$PATH?)\n") |
221 unless (`which psxy` ne ''); |
222 unless (`which psxy` ne ''); |
222 die("$0: ANTSlib required but not found (bad \$PATH?)\n") |
223 die("$0: ANTSlib required but not found (bad \$PATH?)\n") |
223 unless ($ANTS ne ''); |
224 unless ($ANTS ne ''); |
224 die("$0: ADCP Tools required but not found (bad \$PATH?)\n") |
225 die("$0: ADCP Tools required but not found (bad \$PATH?)\n") |
298 @length_of_timelag_windows = split(',',$opt_w); |
299 @length_of_timelag_windows = split(',',$opt_w); |
299 croak("$0: cannot decode -w $opt_w\n") |
300 croak("$0: cannot decode -w $opt_w\n") |
300 unless numberp($length_of_timelag_windows[0]) && numberp($length_of_timelag_windows[1]); |
301 unless numberp($length_of_timelag_windows[0]) && numberp($length_of_timelag_windows[1]); |
301 |
302 |
302 |
303 |
303 |
|
304 |
|
305 croak("$0: \$out_basename undefined\n") # all plotting routines use this |
304 croak("$0: \$out_basename undefined\n") # all plotting routines use this |
306 unless defined($out_basename); |
305 unless defined($out_basename); |
307 &antsAddParams('out_basename',$out_basename); |
306 &antsAddParams('out_basename',$out_basename); |
308 &antsAddParams('run_label',$RUN); |
307 &antsAddParams('run_label',$RUN); |
309 |
308 |
357 |
362 |
358 progress("Reading LADCP data from <$LADCP_file>...\n"); |
363 progress("Reading LADCP data from <$LADCP_file>...\n"); |
359 readData($LADCP_file,\%LADCP); |
364 readData($LADCP_file,\%LADCP); |
360 progress("\t%d ensembles\n",scalar(@{$LADCP{ENSEMBLE}})); |
365 progress("\t%d ensembles\n",scalar(@{$LADCP{ENSEMBLE}})); |
361 |
366 |
362 croak("$LADCP_file: not enough LADCP bins ($LADCP{N_BINS}) for choice of -r\n") |
367 error("$LADCP_file: not enough LADCP bins ($LADCP{N_BINS}) for choice of -r\n") |
363 unless ($LADCP{N_BINS} >= $refLr_lastBin); |
368 unless ($LADCP{N_BINS} >= $refLr_lastBin); |
364 |
369 |
365 croak("$0: first reference-layer bin outside valid range\n") |
370 error("$0: first reference-layer bin outside valid range\n") |
366 unless ($refLr_firstBin>=1 && $refLr_firstBin<=$LADCP{N_BINS}); |
371 unless ($refLr_firstBin>=1 && $refLr_firstBin<=$LADCP{N_BINS}); |
367 croak("$0: last reference-layer bin outside valid range\n") |
372 error("$0: last reference-layer bin outside valid range\n") |
368 unless ($refLr_lastBin>=1 && $refLr_lastBin<=$LADCP{N_BINS}); |
373 unless ($refLr_lastBin>=1 && $refLr_lastBin<=$LADCP{N_BINS}); |
369 croak("$0: first reference-layer bin > last reference-layer bin\n") |
374 error("$0: first reference-layer bin > last reference-layer bin\n") |
370 unless ($refLr_firstBin <= $refLr_lastBin); |
375 unless ($refLr_firstBin <= $refLr_lastBin); |
371 |
376 |
372 $LADCP_lastBin = $LADCP{N_BINS}-1 |
377 $LADCP_lastBin = $LADCP{N_BINS}-1 |
373 if ($LADCP_lastBin eq '*'); |
378 if ($LADCP_lastBin eq '*'); |
374 croak("$0: first valid LADCP bin outside valid range\n") |
379 error("$0: first valid LADCP bin outside valid range\n") |
375 unless ($LADCP_firstBin>=1 && $LADCP_firstBin<=$LADCP{N_BINS}); |
380 unless ($LADCP_firstBin>=1 && $LADCP_firstBin<=$LADCP{N_BINS}); |
376 croak("$0: last valid LADCP bin outside valid range\n") |
381 error("$0: last valid LADCP bin outside valid range\n") |
377 unless ($LADCP_lastBin>=1 && $LADCP_lastBin<=$LADCP{N_BINS}); |
382 unless ($LADCP_lastBin>=1 && $LADCP_lastBin<=$LADCP{N_BINS}); |
378 croak("$0: first valid LADCP bin > last valid LADCP bin\n") |
383 error("$0: first valid LADCP bin > last valid LADCP bin\n") |
379 unless ($LADCP_firstBin <= $LADCP_lastBin); |
384 unless ($LADCP_firstBin <= $LADCP_lastBin); |
380 |
385 |
381 warning(0,"first reference-layer bin < first valid LADCP bin\n") |
386 warning(0,"first reference-layer bin < first valid LADCP bin\n") |
382 unless ($refLr_firstBin >= $LADCP_firstBin); |
387 unless ($refLr_firstBin >= $LADCP_firstBin); |
383 warning(0,"last reference-layer bin > last valid LADCP bin\n") |
388 warning(0,"last reference-layer bin > last valid LADCP bin\n") |
430 correctAttitude($ens,$pitch_bias,$roll_bias,$heading_bias); |
435 correctAttitude($ens,$pitch_bias,$roll_bias,$heading_bias); |
431 $nvv += countValidBeamVels($ens); |
436 $nvv += countValidBeamVels($ens); |
432 $cte += editCorr($ens,$opt_c); |
437 $cte += editCorr($ens,$opt_c); |
433 # $pte += editTilt($ens,$opt_t); |
438 # $pte += editTilt($ens,$opt_t); |
434 } |
439 } |
435 croak("$LADCP_file: no valid data\n") unless ($nvv > 0); |
440 error("$LADCP_file: no valid data\n") unless ($nvv > 0); |
436 progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$opt_c,$cte,round(100*$cte/$nvv)); |
441 progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$opt_c,$cte,round(100*$cte/$nvv)); |
437 # progress("\tattitude threshold (-t %d deg): %d velocites removed (%d%% of total)\n",$opt_t,$pte,round(100*$pte/$nvv)); |
442 # progress("\tattitude threshold (-t %d deg): %d velocites removed (%d%% of total)\n",$opt_t,$pte,round(100*$pte/$nvv)); |
438 } else { |
443 } else { |
439 progress("Editing velocity data...\n"); |
444 progress("Editing velocity data...\n"); |
440 croak("$LADCP_file: cannot apply beamvel-mask $opt_m to earth-coordinate data\n") |
445 error("$LADCP_file: cannot apply beamvel-mask $opt_m to earth-coordinate data\n") |
441 if defined($opt_m); |
446 if defined($opt_m); |
442 $nvv = $cte = 0; |
447 $nvv = $cte = 0; |
443 for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) { |
448 for ($ens=0; $ens<=$#{$LADCP{ENSEMBLE}}; $ens++) { |
444 $nvv += countValidBeamVels($ens); |
449 $nvv += countValidBeamVels($ens); |
445 $cte += editCorr_Earthcoords($ens,$opt_c); |
450 $cte += editCorr_Earthcoords($ens,$opt_c); |
446 # $pte += editTilt($ens,$opt_t); |
451 # $pte += editTilt($ens,$opt_t); |
447 } |
452 } |
448 croak("$LADCP_file: no valid data\n") unless ($nvv > 0); |
453 error("$LADCP_file: no valid data\n") unless ($nvv > 0); |
449 progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$opt_c,$cte,round(100*$cte/$nvv)); |
454 progress("\tcorrelation threshold (-c %d counts): %d velocites removed (%d%% of total)\n",$opt_c,$cte,round(100*$cte/$nvv)); |
450 # progress("\tattitude threshold (-t %d deg): %d velocites removed (%d%% of total)\n",$opt_t,$pte,round(100*$pte/$nvv)); |
455 # progress("\tattitude threshold (-t %d deg): %d velocites removed (%d%% of total)\n",$opt_t,$pte,round(100*$pte/$nvv)); |
451 } |
456 } |
452 |
457 |
453 #------------------------------------------------------------------- |
458 #------------------------------------------------------------------- |
521 } |
526 } |
522 } |
527 } |
523 progress("\t$nvw valid velocities in bins $LADCP_firstBin-$LADCP_lastBin\n"); |
528 progress("\t$nvw valid velocities in bins $LADCP_firstBin-$LADCP_lastBin\n"); |
524 } |
529 } |
525 |
530 |
526 croak("$LADCP_file: insufficient valid velocities\n") unless ($nvw > 1000); |
531 error("$LADCP_file: insufficient valid velocities\n") unless ($nvw > 1000); |
527 |
532 |
528 #---------------------------------------------- |
533 #---------------------------------------------- |
529 # STEP: Edit earth-coordinate -velocity data |
534 # STEP: Edit earth-coordinate -velocity data |
530 # 1) error-velocity threshold |
535 # 1) error-velocity threshold |
531 # 2) vertical-velocity outliers |
536 # 2) vertical-velocity outliers |
554 |
559 |
555 progress("Calculating LADCP time-series...\n"); |
560 progress("Calculating LADCP time-series...\n"); |
556 |
561 |
557 ($firstGoodEns,$lastGoodEns,$LADCP_atbottom,$LADCP_w_gap_time) = |
562 ($firstGoodEns,$lastGoodEns,$LADCP_atbottom,$LADCP_w_gap_time) = |
558 calcLADCPts(\%LADCP,$opt_s,$refLr_firstBin,$refLr_lastBin,$opt_g); |
563 calcLADCPts(\%LADCP,$opt_s,$refLr_firstBin,$refLr_lastBin,$opt_g); |
559 croak("$LADCP_file: no good ensembles\n") |
564 error("$LADCP_file: no good ensembles\n") |
560 unless defined($firstGoodEns) && ($lastGoodEns-$firstGoodEns > 0); |
565 unless defined($firstGoodEns) && ($lastGoodEns-$firstGoodEns > 0); |
561 |
566 |
562 my($cast_duration) = $LADCP{ENSEMBLE}[$lastGoodEns]->{ELAPSED} - |
567 my($cast_duration) = $LADCP{ENSEMBLE}[$lastGoodEns]->{ELAPSED} - |
563 $LADCP{ENSEMBLE}[$firstGoodEns]->{ELAPSED}; |
568 $LADCP{ENSEMBLE}[$firstGoodEns]->{ELAPSED}; |
564 croak("$0: implausibly short cast ($cast_duration seconds)\n") |
569 error("$0: implausibly short cast ($cast_duration seconds)\n") |
565 unless ($cast_duration > 600); |
570 unless ($cast_duration > 600); |
566 |
571 |
567 $LADCP{MEAN_DT} = $cast_duration / ($lastGoodEns-$firstGoodEns-1); |
572 $LADCP{MEAN_DT} = $cast_duration / ($lastGoodEns-$firstGoodEns-1); |
568 |
573 |
569 progress("\tStart of cast : %s (#%5d)\n", |
574 progress("\tStart of cast : %s (#%5d)\n", |
588 @antsNewLayout = ('ensemble','elapsed','reflr_w','reflr_w.stddev','reflr_w.nsamp','depth'); |
593 @antsNewLayout = ('ensemble','elapsed','reflr_w','reflr_w.stddev','reflr_w.nsamp','depth'); |
589 |
594 |
590 foreach my $of (@out_LADCP) { |
595 foreach my $of (@out_LADCP) { |
591 progress("<$of> "); |
596 progress("<$of> "); |
592 $of = ">$of" unless ($of =~ /^$|^\s*\|/); |
597 $of = ">$of" unless ($of =~ /^$|^\s*\|/); |
593 open(STDOUT,$of) || croak("$of: $!\n"); |
598 open(STDOUT,$of) || error("$of: $!\n"); |
594 undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE); |
599 undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE); |
595 for (my($ens)=$firstGoodEns; $ens<=$lastGoodEns; $ens++) { |
600 for (my($ens)=$firstGoodEns; $ens<=$lastGoodEns; $ens++) { |
596 &antsOut($LADCP{ENSEMBLE}[$ens]->{NUMBER}, |
601 &antsOut($LADCP{ENSEMBLE}[$ens]->{NUMBER}, |
597 $LADCP{ENSEMBLE}[$ens]->{ELAPSED}, |
602 $LADCP{ENSEMBLE}[$ens]->{ELAPSED}, |
598 $LADCP{ENSEMBLE}[$ens]->{REFLR_W}, |
603 $LADCP{ENSEMBLE}[$ens]->{REFLR_W}, |
603 &antsOut('EOF'); open(STDOUT,'>&2'); |
608 &antsOut('EOF'); open(STDOUT,'>&2'); |
604 } |
609 } |
605 progress("\n"); |
610 progress("\n"); |
606 } |
611 } |
607 |
612 |
|
613 error("deepest depth is at end of cast (no upcast data)\n") |
|
614 if ($lastGoodEns-$LADCP_atbottom < 100); |
|
615 |
608 #---------------------------------------------------------------------- |
616 #---------------------------------------------------------------------- |
609 # More editing |
617 # More editing |
610 # - this requires ${first,last}GoodEns to be known |
618 # - this requires ${first,last}GoodEns to be known |
611 # - TILT field is set as a side-effect |
619 # - TILT field is set as a side-effect |
612 #---------------------------------------------------------------------- |
620 #---------------------------------------------------------------------- |
633 #-------------- |
641 #-------------- |
634 # Read CTD data |
642 # Read CTD data |
635 #-------------- |
643 #-------------- |
636 |
644 |
637 progress("Reading CTD data from <$CTD_file>...\n"); |
645 progress("Reading CTD data from <$CTD_file>...\n"); |
638 open(STDIN,$CTD_file) || croak("$CTD_file: $!\n"); |
646 open(STDIN,$CTD_file) || error("$CTD_file: $!\n"); |
639 croak("$CTD_file: no data\n") unless (&antsIn()); |
647 error("$CTD_file: no data\n") unless (&antsIn()); |
640 undef($antsOldHeaders); |
648 undef($antsOldHeaders); |
641 ($CTD_elapsed,$CTD_depth,$CTD_svel,$CTD_w) = &fnr('elapsed','depth','ss','w'); |
649 ($CTD_elapsed,$CTD_depth,$CTD_svel,$CTD_w) = &fnr('elapsed','depth','ss','w'); |
642 $CTD_temp = &fnrNoErr('temp'); |
650 $CTD_temp = &fnrNoErr('temp'); |
643 &antsAddParams('lat',$P{lat}) if defined($P{lat}); |
651 &antsAddParams('lat',$P{lat}) if defined($P{lat}); |
644 &antsAddParams('lon',$P{lon}) if defined($P{lon}); |
652 &antsAddParams('lon',$P{lon}) if defined($P{lon}); |
645 |
653 |
646 $CTD_maxdepth = -1; |
654 $CTD_maxdepth = -1; |
647 |
655 |
648 do { # read data |
656 do { # read data |
649 croak("$0: cannot deal with non-numeric CTD elapsed time\n") |
657 error("$0: cannot deal with non-numeric CTD elapsed time\n") |
650 unless &antsNumbers($CTD_elapsed); |
658 unless &antsNumbers($CTD_elapsed); |
651 push(@{$CTD{ELAPSED}},$ants_[0][$CTD_elapsed]); |
659 push(@{$CTD{ELAPSED}},$ants_[0][$CTD_elapsed]); |
652 push(@{$CTD{DEPTH}}, $ants_[0][$CTD_depth]); |
660 push(@{$CTD{DEPTH}}, $ants_[0][$CTD_depth]); |
653 push(@{$CTD{SVEL}}, $ants_[0][$CTD_svel]); |
661 push(@{$CTD{SVEL}}, $ants_[0][$CTD_svel]); |
654 push(@{$CTD{W}}, $ants_[0][$CTD_w]); |
662 push(@{$CTD{W}}, $ants_[0][$CTD_w]); |
667 for (my($s)=1; $s<@{$CTD{ELAPSED}}-1; $s++) { # centered differences for both |
675 for (my($s)=1; $s<@{$CTD{ELAPSED}}-1; $s++) { # centered differences for both |
668 $CTD{W_t} [$s] = ($CTD{W}[$s+1] - $CTD{W}[$s-1]) / (2*$CTD{DT}); |
676 $CTD{W_t} [$s] = ($CTD{W}[$s+1] - $CTD{W}[$s-1]) / (2*$CTD{DT}); |
669 $CTD{W_tt}[$s] = ($CTD{W}[$s+1] + $CTD{W}[$s-1] - 2*$CTD{W}[$s]) / $CTD{DT}**2; |
677 $CTD{W_tt}[$s] = ($CTD{W}[$s+1] + $CTD{W}[$s-1] - 2*$CTD{W}[$s]) / $CTD{DT}**2; |
670 } |
678 } |
671 |
679 |
672 croak("$0: CTD start depth must be numeric\n") |
680 error("$0: CTD start depth must be numeric\n") |
673 unless numberp($CTD{DEPTH}[0]); |
681 unless numberp($CTD{DEPTH}[0]); |
674 if (($CTD{DEPTH}[0] < -$opt_d) && !defined($opt_d)) { |
682 if (($CTD{DEPTH}[0] < -$opt_d) && !defined($opt_d)) { |
675 $opt_d = -1 * round($CTD{DEPTH}[0]); |
683 $opt_d = -1 * round($CTD{DEPTH}[0]); |
676 } |
684 } |
677 if ($opt_d > 0) { |
685 if ($opt_d > 0) { |
745 #------------------------ |
753 #------------------------ |
746 |
754 |
747 $CTD{TIME_LAG} = |
755 $CTD{TIME_LAG} = |
748 calc_lag($number_of_timelag_windows[0],$length_of_timelag_windows[0], |
756 calc_lag($number_of_timelag_windows[0],$length_of_timelag_windows[0], |
749 int(1/$CTD{DT}+0.5),$firstGoodEns,$lastGoodEns); |
757 int(1/$CTD{DT}+0.5),$firstGoodEns,$lastGoodEns); |
750 croak("$0: Cannot proceed without valid lag!\n") unless defined($CTD{TIME_LAG}); |
758 error("$0: Cannot proceed without valid lag!\n") unless defined($CTD{TIME_LAG}); |
751 progress("\telapsed(CTD) ~ elapsed(LADCP) + %.2fs\n",$CTD{TIME_LAG}); |
759 progress("\telapsed(CTD) ~ elapsed(LADCP) + %.2fs\n",$CTD{TIME_LAG}); |
752 |
760 |
753 #--------------------------------- |
761 #--------------------------------- |
754 # stage 2: piece-wise time lagging |
762 # stage 2: piece-wise time lagging |
755 #--------------------------------- |
763 #--------------------------------- |
787 push(@CTD_time_lag,nan); |
795 push(@CTD_time_lag,nan); |
788 } |
796 } |
789 shift(@splits); |
797 shift(@splits); |
790 } |
798 } |
791 |
799 |
792 croak("$0: Cannot proceed without at least one lag!\n") # fill failed lag with surrounding data |
800 error("$0: Cannot proceed without at least one lag!\n") # fill failed lag with surrounding data |
793 unless defined($valid_lag); |
801 unless defined($valid_lag); |
794 while ($valid_lag < $#CTD_time_lag) { # forward |
802 while ($valid_lag < $#CTD_time_lag) { # forward |
795 $CTD_time_lag[$valid_lag+1] = $CTD_time_lag[$valid_lag]; |
803 $CTD_time_lag[$valid_lag+1] = $CTD_time_lag[$valid_lag]; |
796 $valid_lag++; |
804 $valid_lag++; |
797 } |
805 } |
849 |
857 |
850 if (defined($LADCP{ENSEMBLE}[$ens]->{REFLR_W}) # not a gap |
858 if (defined($LADCP{ENSEMBLE}[$ens]->{REFLR_W}) # not a gap |
851 && numberp($CTD{DEPTH}[$scan])) { |
859 && numberp($CTD{DEPTH}[$scan])) { |
852 $LADCP{ENSEMBLE}[$ens]->{REFLR_W_NOSSCORR} = $LADCP{ENSEMBLE}[$ens]->{REFLR_W}; |
860 $LADCP{ENSEMBLE}[$ens]->{REFLR_W_NOSSCORR} = $LADCP{ENSEMBLE}[$ens]->{REFLR_W}; |
853 $LADCP{ENSEMBLE}[$ens]->{REFLR_W} *= $CTD{SVEL}[$scan]/1500; # correct for sound-speed variations at source |
861 $LADCP{ENSEMBLE}[$ens]->{REFLR_W} *= $CTD{SVEL}[$scan]/1500; # correct for sound-speed variations at source |
854 croak(sprintf("\n$0: negative depth (%.1fm) in CTD file at elapsed(CTD) = %.1fs (use -d?)\n", |
862 error(sprintf("\n$0: negative depth (%.1fm) in CTD file at elapsed(CTD) = %.1fs (use -d?)\n", |
855 $CTD{DEPTH}[$scan],$CTD{ELAPSED}[$scan])) |
863 $CTD{DEPTH}[$scan],$CTD{ELAPSED}[$scan])) |
856 unless ($CTD{DEPTH}[$scan] >= 0); |
864 unless ($CTD{DEPTH}[$scan] >= 0); |
857 $LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} = $CTD{DEPTH}[$scan]; |
865 $LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH} = $CTD{DEPTH}[$scan]; |
858 $LADCP{ENSEMBLE}[$ens]->{CTD_SCAN} = $scan; |
866 $LADCP{ENSEMBLE}[$ens]->{CTD_SCAN} = $scan; |
859 my($reflr_ocean_w) = $LADCP{ENSEMBLE}[$ens]->{REFLR_W} - $CTD{W}[$scan]; |
867 my($reflr_ocean_w) = $LADCP{ENSEMBLE}[$ens]->{REFLR_W} - $CTD{W}[$scan]; |
878 &antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',sqrt($sumWsqI/$nWsqI)); |
886 &antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',sqrt($sumWsqI/$nWsqI)); |
879 progress("\t%.2f cm/s rms reference-layer w_ocean, %.2f cm/s away from boundaries\n", |
887 progress("\t%.2f cm/s rms reference-layer w_ocean, %.2f cm/s away from boundaries\n", |
880 100*sqrt($sumWsq/$nWsq),100*sqrt($sumWsqI/$nWsqI)); |
888 100*sqrt($sumWsq/$nWsq),100*sqrt($sumWsqI/$nWsqI)); |
881 warning(0,"%.2f cm/s reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI)) |
889 warning(0,"%.2f cm/s reference-layer w_ocean away from boundaries\n",100*sqrt($sumWsqI/$nWsqI)) |
882 if (sqrt($sumWsqI/$nWsqI) > 0.05); |
890 if (sqrt($sumWsqI/$nWsqI) > 0.05); |
883 # croak("$0: rms reference-layer w_ocean is too large\n") |
891 # error("$0: rms reference-layer w_ocean is too large\n") |
884 # unless (sqrt($sumWsqI/$nWsqI) < 0.07); |
892 # unless (sqrt($sumWsqI/$nWsqI) < 0.07); |
885 } elsif ($nWsq > 0) { |
893 } elsif ($nWsq > 0) { |
886 &antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',nan); |
894 &antsAddParams('rms_w_reflr_err',sqrt($sumWsq/$nWsq),'rms_w_reflr_err_interior',nan); |
887 progress("\t%.2f cm/s rms reference-layer w_ocean\n",100*sqrt($sumWsq/$nWsq)); |
895 progress("\t%.2f cm/s rms reference-layer w_ocean\n",100*sqrt($sumWsq/$nWsq)); |
888 } else { |
896 } else { |
889 croak("$0: no valid vertical velocities\n"); |
897 error("$0: no valid vertical velocities\n"); |
890 } |
898 } |
891 |
899 |
892 #---------------------------------------------------------------------- |
900 #---------------------------------------------------------------------- |
893 # Calculate Volume Scattering Coefficients |
901 # Calculate Volume Scattering Coefficients |
894 #---------------------------------------------------------------------- |
902 #---------------------------------------------------------------------- |
913 progress("Finding seabed...\n"); |
921 progress("Finding seabed...\n"); |
914 ($water_depth,$sig_water_depth) = |
922 ($water_depth,$sig_water_depth) = |
915 find_backscatter_seabed($LADCP{ENSEMBLE}[$LADCP_atbottom]->{CTD_DEPTH}); |
923 find_backscatter_seabed($LADCP{ENSEMBLE}[$LADCP_atbottom]->{CTD_DEPTH}); |
916 ($water_depth_BT,$sig_water_depth_BT) = |
924 ($water_depth_BT,$sig_water_depth_BT) = |
917 find_seabed(\%LADCP,$LADCP_atbottom,$LADCP{BEAM_COORDINATES}); |
925 find_seabed(\%LADCP,$LADCP_atbottom,$LADCP{BEAM_COORDINATES}); |
918 if (defined($water_depth_BT)) { |
926 if (defined($water_depth) && defined($water_depth_BT)) { |
919 my($dd) = abs($water_depth_BT - $water_depth); |
927 my($dd) = abs($water_depth_BT - $water_depth); |
920 warning(2,sprintf("Large RDI vs. own water-depth difference (%.1fm)\n",$dd)) |
928 warning(2,sprintf("Large instrument vs. backscatter-derived water-depth difference (%.1fm)\n",$dd)) |
921 if ($dd > 5); |
929 if ($dd > 5); |
922 } |
930 } |
|
931 if (!$SS_use_BT && !defined($water_depth) && defined($water_depth_BT)) { |
|
932 warning(1,"using water_depth from BT data\n"); |
|
933 $SS_use_BT = 1; |
|
934 } |
|
935 if ($SS_use_BT) { |
|
936 $water_depth = $water_depth_BT; |
|
937 $sig_water_depth = $sig_water_depth_BT; |
|
938 } |
923 } |
939 } |
924 |
940 |
925 &antsAddParams('water_depth',$water_depth,'water_depth.sig',$sig_water_depth); |
941 &antsAddParams('water_depth',$water_depth,'water_depth.sig',$sig_water_depth) |
|
942 if defined($water_depth); |
926 |
943 |
927 if (numberp($water_depth)) { |
944 if (numberp($water_depth)) { |
928 if (numberp($water_depth_BT)) { |
945 if (numberp($water_depth_BT)) { |
929 progress("\t%.1f(%.1f) m water depth (%.1f(%.1f)m from BT)\n", |
946 progress("\t%.1f(%.1f) m water depth (%.1f(%.1f)m from BT)\n", |
930 $water_depth,$sig_water_depth,$water_depth_BT,$sig_water_depth_BT); |
947 $water_depth,$sig_water_depth,$water_depth_BT,$sig_water_depth_BT); |
1264 &antsAddParams('BR_max_bin',max(scalar(@dc_bres),scalar(@uc_bres))); |
1281 &antsAddParams('BR_max_bin',max(scalar(@dc_bres),scalar(@uc_bres))); |
1265 |
1282 |
1266 foreach my $of (@out_BR) { |
1283 foreach my $of (@out_BR) { |
1267 progress("<$of> "); |
1284 progress("<$of> "); |
1268 $of = ">$of" unless ($of =~ /^$|^\s*\|/); |
1285 $of = ">$of" unless ($of =~ /^$|^\s*\|/); |
1269 open(STDOUT,$of) || croak("$of: $!\n"); |
1286 open(STDOUT,$of) || error("$of: $!\n"); |
1270 undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE); |
1287 undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE); |
1271 for (my($bin)=0; $bin<max(scalar(@dc_bres),scalar(@uc_bres)); $bin++) { |
1288 for (my($bin)=0; $bin<max(scalar(@dc_bres),scalar(@uc_bres)); $bin++) { |
1272 my($dc_avg) = avg(@{$dc_bres[$bin]}); |
1289 my($dc_avg) = avg(@{$dc_bres[$bin]}); |
1273 my($uc_avg) = avg(@{$uc_bres[$bin]}); |
1290 my($uc_avg) = avg(@{$uc_bres[$bin]}); |
1274 &antsOut($bin+1,$dc_avg,stddev2($dc_avg,@{$dc_bres[$bin]}),scalar(@{$dc_bres[$bin]}), |
1291 &antsOut($bin+1,$dc_avg,stddev2($dc_avg,@{$dc_bres[$bin]}),scalar(@{$dc_bres[$bin]}), |
1316 'pitch','roll','tilt','heading','3_beam','svel'); |
1333 'pitch','roll','tilt','heading','3_beam','svel'); |
1317 |
1334 |
1318 foreach my $of (@out_w) { |
1335 foreach my $of (@out_w) { |
1319 progress("<$of> "); |
1336 progress("<$of> "); |
1320 $of = ">$of" unless ($of =~ /^$|^\s*\|/); |
1337 $of = ">$of" unless ($of =~ /^$|^\s*\|/); |
1321 open(STDOUT,$of) || croak("$of: $!\n"); |
1338 open(STDOUT,$of) || error("$of: $!\n"); |
1322 undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE); |
1339 undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE); |
1323 |
1340 |
1324 for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) { # downcast |
1341 for ($ens=$firstGoodEns; $ens<$LADCP_atbottom; $ens++) { # downcast |
1325 next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}); |
1342 next unless numberp($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}); |
1326 my(@bindepth) = calc_binDepths($ens); |
1343 my(@bindepth) = calc_binDepths($ens); |
1405 'BT_w','BT_w.mad','BT_w.nsamp'); |
1422 'BT_w','BT_w.mad','BT_w.nsamp'); |
1406 |
1423 |
1407 foreach my $of (@out_profile) { |
1424 foreach my $of (@out_profile) { |
1408 progress("<$of> "); |
1425 progress("<$of> "); |
1409 $of = ">$of" unless ($of =~ /^$|^\s*\|/); |
1426 $of = ">$of" unless ($of =~ /^$|^\s*\|/); |
1410 open(STDOUT,$of) || croak("$of: $!\n"); |
1427 open(STDOUT,$of) || error("$of: $!\n"); |
1411 undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE); |
1428 undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE); |
1412 |
1429 |
1413 for (my($bi)=0; $bi<=max($#{$DNCAST{ENSEMBLE}},$#{$UPCAST{ENSEMBLE}},$#{$BT{NSAMP}}); $bi++) { |
1430 for (my($bi)=0; $bi<=max($#{$DNCAST{ENSEMBLE}},$#{$UPCAST{ENSEMBLE}},$#{$BT{NSAMP}}); $bi++) { |
1414 &antsOut(($bi+0.5)*$opt_o, # nominal depth |
1431 &antsOut(($bi+0.5)*$opt_o, # nominal depth |
1415 $DNCAST{MEAN_DEPTH}[$bi],$DNCAST{MEAN_ELAPSED}[$bi], |
1432 $DNCAST{MEAN_DEPTH}[$bi],$DNCAST{MEAN_ELAPSED}[$bi], |
1441 'reflr_ocean_w'); |
1458 'reflr_ocean_w'); |
1442 |
1459 |
1443 foreach my $of (@out_timeseries) { |
1460 foreach my $of (@out_timeseries) { |
1444 progress("<$of> "); |
1461 progress("<$of> "); |
1445 $of = ">$of" unless ($of =~ /^$|^\s*\|/); |
1462 $of = ">$of" unless ($of =~ /^$|^\s*\|/); |
1446 open(STDOUT,$of) || croak("$of: $!\n"); |
1463 open(STDOUT,$of) || error("$of: $!\n"); |
1447 undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE); |
1464 undef($antsActiveHeader) unless ($ANTS_TOOLS_AVAILABLE); |
1448 |
1465 |
1449 for ($ens=$firstGoodEns; $ens<=$realLastGoodEns; $ens++) { |
1466 for ($ens=$firstGoodEns; $ens<=$realLastGoodEns; $ens++) { |
1450 next unless defined($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}); |
1467 next unless defined($LADCP{ENSEMBLE}[$ens]->{CTD_DEPTH}); |
1451 my($reflr_oc_w) = defined($LADCP{ENSEMBLE}[$ens]->{REFLR_W}) |
1468 my($reflr_oc_w) = defined($LADCP{ENSEMBLE}[$ens]->{REFLR_W}) |