#!/usr/bin/env perl
#--------------------------------------------------------------------
#    The MB-system:	mbm_bpr.perl	5/13/2002
#
#    Copyright (c) 2012-2023 by
#    David W. Caress (caress@mbari.org)
#      Monterey Bay Aquarium Research Institute
#      Moss Landing, California, USA
#    Dale N. Chayes 
#      Center for Coastal and Ocean Mapping
#      University of New Hampshire
#      Durham, New Hampshire, USA
#    Christian dos Santos Ferreira
#      MARUM
#      University of Bremen
#      Bremen Germany
#      
#    MB-System was created by Caress and Chayes in 1992 at the
#      Lamont-Doherty Earth Observatory
#      Columbia University
#      Palisades, NY 10964
#
#    See README.md file for copying and redistribution conditions.
#--------------------------------------------------------------------
#
# Command:
#   mbm_bpr
#
# Purpose:
#   MB-System macro to process data from a Seabird SBE53 pressure
#   sensor into a tidal model for use by mbprocess
#
# Basic Usage:
#   mbm_bpr -Ifile -Ooutput [-Doffset -Fformat Rlongitude/latitude -Swindow -Tmode -H -V]
#
# Author:
#   David W. Caress
#   Monterey Bay Aquarium Research Institute
#   Moss Landing, CA 95039
#   May 13, 2002
#
#
#
use Getopt::Std;

my $Program_name = "mbm_bpr";

# some default values
$format = 0;
$tmode = 0;
$surfacepressure = 14.0;
$tidemin = 0.0;
$tidemax = 0.0;
$PI = 3.1415926;

# Deal with command line arguments using the "standard" perl library
if ( ! (getopts('D:d:F:f:HhI:i:O:o:R:r:S:s:TtVv'))) {
  print "$0: error on input\n";
  usage();
  exit 0;
}

my $help     = ($opt_H || $opt_h);
my $format   = ($opt_F || $opt_f);
my $ifile    = ($opt_I || $opt_i);
my $ofile    = ($opt_O || $opt_o);
my $location = ($opt_R || $opt_r);
my $offset   = ($opt_D || $opt_d);
my $tmode    = ($opt_T || $opt_t);
my $smooth   = ($opt_S || $opt_s);
my $verbose  = ($opt_V || $opt_v);

# print out help message if required
if ($help)
	{
	print "\n$ProgramName:\n";
	print "\nMB-System macro to process data from a Seabird SBE53 pressure\n";
	print "sensor into a tidal model for use by mbprocess\n";
	print "\nUsage: \n";
	print "\t$Program_name -Ibprfile -Otidefile [-Doffset -Fformat -Rbprlon/bprlat -T -H -V]\n";
	exit 0;
	}

# Get position if that has been specified
if ($location && $location =~ /^(\S+)\/(\S+)/) {
    ($longitude, $latitude) = $location =~ /^(\S+)\/(\S+)/;
}
else {
    $latitude = 0.0;
}

# Open the input and output files. Failure on either is fatal.
open(FI,$ifile) || die "Cannot open input file: $ifile\n$ProgramName aborted.\n";
open(FO,">".$ofile) || die "Cannot open output file: $ofile\n$ProgramName aborted.\n";

# Loop over the input records for SeaBird format
if ($format == 0) {
    $count = 0;
    while ($line=<FI>) {		# reading from the input
        ($tcount, $month, $day, $year, $hour, $minute, $second, $pressure, $temperature)
            = $line =~ /(\S+)\s+(\S+)\/(\S+)\/(\S+)\s+(\S+):(\S+):(\S+)\s+(\S+)\s+(\S+)/;

        # Use the first value as the surface pressure provided it is valid
        if ($tcount == 1 && $pressure < 99999.0) {
            $surfacepressure = $pressure;
        }

        # only use valid pressure values
        if ($pressure < 99999.0) {
            $time_d = `mbtime -T$year/$month/$day/$hour/$minute/$second`;
            chop($time_d);
            push (@time_ds, $time_d);

            push (@years, $year);
            push (@months, $month);
            push (@days, $day);
            push (@hours, $hour);
            push (@minutes, $minute);
            push (@seconds, $second);

            # Calculate the depth from pressure
            #
            # The following derives from Sea-Bird Electronics
            # Application Note number 69: Conversion of Pressure to Depth
            #	http://www.seabird.com/application_notes/AN69.htm
            #
            # Sea-Bird uses the formula in UNESCO Technical Papers in Marine Science No. 44.
            # This is an empirical formula that takes compressibility (that is, density) into account.
            # An ocean water column at 0 degree C (t = 0) and 35 PSU (s = 35) is assumed.
            #
            # The gravity variation with latitude and pressure is computed as:
            #
            # g (m/sec2) = 9.780318 * [ 1.0 + ( 5.2788x10 -3  + 2.36x10 -5  * x) * x ] + 1.092x10 -6  * p
            #
            # where
            # x = [sin (latitude / 57.29578) ] 2
            # p = pressure (decibars)
            #
            # Then, depth is calculated from pressure:
            #
            # depth (meters) = [(((-1.82x10 -15  * p + 2.279x10 -10 ) * p - 2.2512x10 -5 ) * p + 9.72659) * p] / g
            #
            # where
            # p = pressure (decibars)
            # g = gravity (m/sec2)
            #
            # N. P. Fofonoff and R. C. Millard, Jr., Algorithms for computation of fundamental properties of seawater, Unesco Tech. Papers in Mar. Sci., No. 44 1983.
            #
            $pressuredbar = 0.6894757293 * ($pressure - $surfacepressure);
            $x = (sin($latitude * $PI / 180.0))**2;
            $g = 9.780318 * ( 1.0 + ( 5.2788e-3  + 2.36e-5  * $x) * $x ) + 1.092e-6  * $pressuredbar;
            $depth = ((((-1.82e-15  * $pressuredbar + 2.279e-10 ) * $pressuredbar - 2.2512e-5 ) * $pressuredbar + 9.72659) * $pressuredbar) / $g;

            push (@depths, $depth);
            $count++;
        }
    }
}

# Loop over the input records for Sonardyne format
# Here the raw pressure values are in kPa == 10 dBar
elsif ($format == 1 || $format == 2) {
    $count = 0;
    $depthcount = 0;
    $depthsum = 0.0;
    $time_d_save = 0.0;
    if ($format == 1) {
      $tag = "PR2";
    } else {
      $tag = "PRS";
    }
    while ($line=<FI>) {		# reading from the input
        @fields = split(",", $line);
        if ($fields[0] eq $tag && $fields[1] ne "Record Time") {
            $numcolons = $fields[1] =~ tr/://;
            if ($numcolons == 2) {
                ($year,$month,$day,$hour,$minute,$second) = $fields[1] =~ /(\S+)\/(\S+)\/(\S+)\s+(\S+):(\S+):(\S+)/;
              }
            else {
              ($month,$day,$year,$hour,$minute) = $fields[1] =~ /(\S+)\/(\S+)\/(\S+)\s+(\S+):(\S+)/;
              $second = 0.0;
            }
            if ($year < 2000) {
              if ($year < 50) {
                $year += 2000;
              } else {
                $year += 1900;
              }
            }
            $time_d = `mbtime -T$year/$month/$day/$hour/$minute/$second`;
            chop($time_d);
            ($pressure) = $fields[4] =~ /(\S+)/;

            # Use the first value as the surface pressure
            if ($count == 0 && $depthcount == 0) {
                $surfacepressure = $pressure;
            }

            # Calculate the depth from pressure
            #
            # The following derives from Sea-Bird Electronics
            # Application Note number 69: Conversion of Pressure to Depth
            #	http://www.seabird.com/application_notes/AN69.htm
            #
            # Sea-Bird uses the formula in UNESCO Technical Papers in Marine Science No. 44.
            # This is an empirical formula that takes compressibility (that is, density) into account.
            # An ocean water column at 0 degree C (t = 0) and 35 PSU (s = 35) is assumed.
            #
            # The gravity variation with latitude and pressure is computed as:
            #
            # g (m/sec2) = 9.780318 * [ 1.0 + ( 5.2788x10 -3  + 2.36x10 -5  * x) * x ] + 1.092x10 -6  * p
            #
            # where
            # x = [sin (latitude / 57.29578) ] 2
            # p = pressure (decibars)
            #
            # Then, depth is calculated from pressure:
            #
            # depth (meters) = [(((-1.82x10 -15  * p + 2.279x10 -10 ) * p - 2.2512x10 -5 ) * p + 9.72659) * p] / g
            #
            # where
            # p = pressure (decibars)
            # g = gravity (m/sec2)
            #
            # N. P. Fofonoff and R. C. Millard, Jr., Algorithms for computation of fundamental properties of seawater, Unesco Tech. Papers in Mar. Sci., No. 44 1983.
            #
            $pressuredbar = 0.1 * ($pressure - $surfacepressure);
            $x = (sin($latitude * $PI / 180.0))**2;
            $g = 9.780318 * ( 1.0 + ( 5.2788e-3  + 2.36e-5  * $x) * $x ) + 1.092e-6  * $pressuredbar;
            $depth = ((((-1.82e-15  * $pressuredbar + 2.279e-10 ) * $pressuredbar - 2.2512e-5 ) * $pressuredbar + 9.72659) * $pressuredbar) / $g;
            if ($time_d != $time_d_save && $depthcount > 0) {
                $depth_save = $depthsum / $depthcount;
                push (@time_ds, $time_d_save);
                push (@years, $year_save);
                push (@months, $month_save);
                push (@days, $day_save);
                push (@hours, $hour_save);
                push (@minutes, $minute_save);
                push (@seconds, $second_save);
                push (@depths, $depth_save);
                $count++;
                $depthcount = 0;
                $depthsum = 0.0;

# if ($verbose && $count % 100 == 0) {
#  print "Output $count pressure values $time_d_save $depth_save\n";
# }
            }

            $time_d_save = $time_d;
            $year_save = $year;
            $month_save = $month;
            $day_save = $day;
            $hour_save = $hour;
            $minute_save = $minute;
            $second_save = $second;
            $depthsum += $depth;
            $depthcount++;

        }
    }

    if ($depthcount > 0) {
        $depth_save = $depthsum / $depthcount;
        push (@time_ds, $time_d_save);
        push (@years, $year_save);
        push (@months, $month_save);
        push (@days, $day_save);
        push (@hours, $hour_save);
        push (@minutes, $minute_save);
        push (@seconds, $second_save);
        push (@depths, $depth_save);
    }
}


# Loop over the input records for generic format
# Here the raw pressure values are in dBar
elsif ($format == 3) {
    $count = 0;
    while ($line=<FI>) {		# reading from the input
        ($time_d, $pressure) = $line =~ /(\S+)\s+(\S+)/;
        $time_str = `mbtime -M$time_d`;
        ($year, $month, $day, $hour, $minute, $second) = $time_str =~ /(\S+)\/(\S+)\/(\S+)\/(\S+)\/(\S+)\/(\S+)/;

        if ($year > 1970 && $year < 2100) {

            # Calculate the depth from pressure
            #
            # The following derives from Sea-Bird Electronics
            # Application Note number 69: Conversion of Pressure to Depth
            #	http://www.seabird.com/application_notes/AN69.htm
            #
            # Sea-Bird uses the formula in UNESCO Technical Papers in Marine Science No. 44.
            # This is an empirical formula that takes compressibility (that is, density) into account.
            # An ocean water column at 0 degree C (t = 0) and 35 PSU (s = 35) is assumed.
            #
            # The gravity variation with latitude and pressure is computed as:
            #
            # g (m/sec2) = 9.780318 * [ 1.0 + ( 5.2788x10 -3  + 2.36x10 -5  * x) * x ] + 1.092x10 -6  * p
            #
            # where
            # x = [sin (latitude / 57.29578) ] 2
            # p = pressure (decibars)
            #
            # Then, depth is calculated from pressure:
            #
            # depth (meters) = [(((-1.82x10 -15  * p + 2.279x10 -10 ) * p - 2.2512x10 -5 ) * p + 9.72659) * p] / g
            #
            # where
            # p = pressure (decibars)
            # g = gravity (m/sec2)
            #
            # N. P. Fofonoff and R. C. Millard, Jr., Algorithms for computation of fundamental properties of seawater, Unesco Tech. Papers in Mar. Sci., No. 44 1983.
            #
            $pressuredbar = 0.689475729 * $pressure;
            $x = (sin($latitude * $PI / 180.0))**2;
            $g = 9.780318 * ( 1.0 + ( 5.2788e-3  + 2.36e-5  * $x) * $x ) + 1.092e-6  * $pressuredbar;
            $depth = ((((-1.82e-15  * $pressuredbar + 2.279e-10 ) * $pressuredbar - 2.2512e-5 ) * $pressuredbar + 9.72659) * $pressuredbar) / $g;

            push (@time_ds, $time_d);
            push (@years, $year);
            push (@months, $month);
            push (@days, $day);
            push (@hours, $hour);
            push (@minutes, $minute);
            push (@seconds, $second);
            push (@depths, $depth);
            $count++;
        }
    }
}

close(FI);

# if requested smooth the tide data
if ($smooth && $count > 1) {
    $interval = int(($time_ds[$count-1] - $time_ds[0]) / ($count - 1));
    $dj = int(0.5 * $smooth / $interval);
    for ($i = 0; $i < $count; $i++) {
        $js = $i - $dj;
        if ($js < 0) {
            $js = 0;
        }
        $jf = $i + $dj;
        if ($jf > $count - 1) {
            $jf = $count - 1;
        }
        $depthsum = 0.0;
        $depthcount = 0;
        for ($j = $js; $j <= $jf; $j++) {
            $depthsum += $depths[$j];
            $depthcount++;
        }
        if ($depthcount > 0) {
            $depth = $depthsum / $depthcount;
        }
        push (@sdepths, $depth);
    }
} else {
    for ($i = 0; $i < $count; $i++) {
        push (@sdepths, $depths[$i]);
    }

}

# tell program status
if ($verbose) {
    print "\nProgram $ProgramName Status:\n";
    print "  $count pressure values read from $ifile\n";
    if ($offset) {
      print "  Depth offset: $offset m\n";
    }
    if ($smooth) {
        print "  Smoothing applied to tide data with a $smooth second window\n";
    }
    if ($location) {
      print "  Vertical reference to tidal model for position $longitude $latitude\n";
    } else {
      print "  Vertical reference will be the mean of the middle half of the data\n";
    }
    if ($tmode == 0) {
      print "  Tide will be output as <time_d tide> values\n";
      print "  A plot will be generated\n";
    }
   else {
      print "  Tide will be output as <yr mo da hr mi sc tide> values\n";
      print "  No plot will be generated\n";
    }
}

# If location specified then use mbotps to get a tidal model
# and then correct the bpr data to a tidal signal by
# subtracting the mean difference between the tidal model and
# the middle half of the data
if ($location && $count > 10) {
    # Execute mbotps to get tidal model for the specified location
    $i = $count / 4;
    $j = 3 * $count / 4;
    $interval = int(($time_ds[$j] - $time_ds[$i]) / ($j - $i));
    $start = "$years[0]/$months[0]/$days[0]/$hours[0]/$minutes[0]/$seconds[0]";
    $end = "$years[$count-1]/$months[$count-1]/$days[$count-1]/$hours[$count-1]/$minutes[$count-1]/$seconds[$count-1]";
    $modelfile = $ifile . "_tidemodel.txt";
    $cmd = "mbotps -A1 -D$interval -R$location -B$start -E$end -O$modelfile";
    if ($verbose) {
      print "  Executing: $cmd\n";
    }
    `$cmd`;

    $nmean = 0;
    $ntide = 0;
    $mean = 0.0;
    open(FI,$modelfile) || die "Cannot open temporary file: $modelfile\n$ProgramName aborted.\n";
    while ($line=<FI>) {		# reading from the input
      if ($line =~ /^\#.+/) {
      }
      else {
	($time_d, $tide) = $line =~ /(\S+)\s+(\S+)/;
	if ($ntide > $count / 4 && $ntide < 3 * $count / 4) {
	  $mean += $sdepths[$i + $nmean] - $tide;
	  $nmean++;
	}
	if ($tide < $tidemin) {
	  $tidemin = $tide;
	}
	if ($tide > $tidemax) {
	  $tidemax = $tide;
	}
	$ntide++;
      }
    }
    close(FI);
    if ($nmean > 0) {
        $mean /= $nmean;
    }
}

# Else if no location specified then correct to a tidal signal by
# subtracting the mean depth of the middle half of the data
# Calculate the mean depth for the middle half of the data
else {
    $nmean = 0;
    $mean = 0.0;
    for ($i = $count / 4; $i < 3 * $count / 4; $i++) {
	$mean += $sdepths[$i];
	$nmean++;
    }
    $mean /= $nmean;
}

# Output the tidal data corrected for the mean (or referenced to a model)
$ocount = 0;
for ($i = 0; $i < $count; $i++) {
  $tide = $sdepths[$i] - $mean + $offset;
  if (abs($tide) < 2.5 && ($i == 0 || $time_ds[$i] != $time_ds[$i-1])) {
    if ($tmode == 0) {
      printf FO "$time_ds[$i] $tide\n";
    }
    else {
      printf FO "%4.4d %2.2d %2.2d %2.2d %2.2d %2.2d %f\n",
    		$years[$i], $months[$i], $days[$i], $hours[$i],
		$minutes[$i], $seconds[$i], $tide;
    }
  $ocount++;
  if ($tide < $tidemin) {
    $tidemin = $tide;
  }
  if ($tide > $tidemax) {
    $tidemax = $tide;
  }
  }
}

close(FO);

# tell program status
if ($verbose) {
    print "  $ocount pressure values output to $ofile\n";
    print "  Vertical reference: $mean m\n";
    }

# generate plot of tide observations (and model if location specified)
if ($tmode == 0) {
  $plotfile = $ofile . "_tideplot";
  if (abs($tidemin) > $tidemax) {
    $tidemax = abs($tidemin);
  }
  $tidemax *= 1.1;
  if ($ntide > 0) {
    $cmd = "mbm_xyplot -R$time_ds[0]/$time_ds[$count-1]/-$tidemax/$tidemax -IWred:$modelfile -IWblack:$ofile -O$plotfile -L\"Tide Data from BPR <$ofile> (black) & Tide Model (red):Seconds:Tide (meters)\" -V";
  }
  else {
    $cmd = "mbm_xyplot -R$time_ds[0]/$time_ds[$count-1]/-$tidemax/$tidemax -IWblack:$ofile -O$plotfile -L\"Tide Data from BPR <$ofile>:Seconds:Tide (meters)\" -V";
  }
  print "  Executing $cmd\n\n";
  $line = `$cmd`;
  print $line;
}

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