swstate.h

Go to the documentation of this file.

// Copyright 2013-2021:
//   GobySoft, LLC (2013-)
//   Massachusetts Institute of Technology (2007-2014)
//   Community contributors (see AUTHORS file)
// File authors:
//   Toby Schneider <toby@gobysoft.org>
//
//
// This file is part of the Goby Underwater Autonomy Project Libraries
// ("The Goby Libraries").
//
// The Goby Libraries are free software: you can redistribute them and/or modify
// them under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 2.1 of the License, or
// (at your option) any later version.
//
// The Goby Libraries are distributed in the hope that they will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with Goby.  If not, see <http://www.gnu.org/licenses/>.
 
// modified for C++ by s. petillo spetillo@mit.edu
// ocean engineering graduate student - mit / whoi joint program
// massachusetts institute of technology (mit)
// laboratory for autonomous marine sensing systems (lamss)
 
#ifndef GOBY_UTIL_SEAWATER_SWSTATE_H
#define GOBY_UTIL_SEAWATER_SWSTATE_H
 
#include <cmath>
 
#include <boost/units/quantity.hpp>
#include <boost/units/systems/si.hpp>
#include <boost/units/systems/si/prefixes.hpp>
#include <boost/units/systems/temperature/celsius.hpp>
 
#include "units.h"
 
namespace goby
{
namespace util
{
namespace seawater
{
template <typename DimensionlessUnit = boost::units::si::dimensionless,
          typename TemperatureUnit = boost::units::celsius::temperature,
          typename PressureUnit = decltype(boost::units::si::deci* bar)>
boost::units::quantity<boost::units::si::mass_density>
density_anomaly(boost::units::quantity<DimensionlessUnit> salinity,
                boost::units::quantity<boost::units::absolute<TemperatureUnit> > temperature,
                boost::units::quantity<PressureUnit> pressure)
{
    using namespace boost::units;
 
    double S = salinity;
    double T = quantity<absolute<celsius::temperature> >(temperature).value();
    double P0 = quantity<decltype(si::deci * bar)>(pressure).value();
 
    /*
 
      SIGMA = density_anomaly(S,T,P) returns the density anomaly SIGMA (kg/m^3)
      given the salinity S (ppt), temperature T (deg C) and pressure
      P (dbars).
    
      Notes: RP (WHOI 2/dec/91)
 
      This stuff is directly copied from the UNESCO algorithms, with some
      minor changes to make it Matlab compatible (like adding ";" and changing
      "*" to "*" when necessary.
 
      RP  (WHOI 3/dec/91)
 
      ******************************************************
      SPECIFIC VOLUME ANOMALY (STERIC ANOMALY) BASED ON 1980 EQUATION
      OF STATE FOR SEAWATER AND 1978 PRACTICAL SALINITY SCALE.
      REFERENCES
      MILLERO, ET AL (1980) DEEP-SEA RES.,27A,255-264
      MILLERO AND POISSON 1981,DEEP-SEA RES.,28A PP 625-629.
      BOTH ABOVE REFERENCES ARE ALSO FOUND IN UNESCO REPORT 38 (1981)
      
      UNITS:      
      PRESSURE        P0       DECIBARS
      TEMPERATURE     T        DEG CELSIUS (IPTS-68)
      SALINITY        S        (IPSS-78)
      SPEC. VOL. ANA. SVAN     M**3/KG *1.0E-8
      DENSITY ANA.    SIGMA    KG/M**3
      ******************************************************************
      CHECK VALUE: SVAN=981.3021 E-8 M**3/KG.  FOR S = 40 (IPSS-78) ,
      T = 40 DEG C, P0= 10000 DECIBARS.
      CHECK VALUE: SIGMA = 59.82037  KG/M**3 FOR S = 40 (IPSS-78) ,
      T = 40 DEG C, P0= 10000 DECIBARS.
      CHECK VALUE: FOR S = 40 (IPSS-78) , T = 40 DEG C, P0= 10000 DECIBARS.
      DR/DP                  DR/DT                 DR/DS
      DRV(1,7)              DRV(2,3)             DRV(1,8)
 
      FINITE DIFFERENCE WITH 3RD ORDER CORRECTION DONE IN DOUBLE PRECSION
 
      3.46969238E-3       -.43311722           .705110777
 
      EXPLICIT DIFFERENTIATION SINGLE PRECISION FORMULATION EOS80 
 
      3.4696929E-3        -.4331173            .7051107
 
      (RP...I think this ---------^^^^^^ should be -.4431173!);
    */
 
    // *******************************************************
    using namespace std;
 
    // DATA
    double R3500 = 1028.1063;
    double R4 = 4.8314E-4;
    double DR350 = 28.106331;
 
    // CONVERT PRESSURE TO BARS AND TAKE SQUARE ROOT SALINITY.
    double P = P0 / 10.;
    //double SAL=S;
    double SR = sqrt(abs(S));
    // *********************************************************
    // PURE WATER DENSITY AT ATMOSPHERIC PRESSURE
    //   BIGG P.H.,(1967) BR. J. APPLIED PHYSICS 8 PP 521-537.
 
    double R1 = ((((6.536332E-9 * T - 1.120083E-6) * T + 1.001685E-4) * T - 9.095290E-3) * T +
                 6.793952E-2) *
                    T -
                28.263737;
    // SEAWATER DENSITY ATM PRESS.
    //  COEFFICIENTS INVOLVING SALINITY
    double R2 = (((5.3875E-9 * T - 8.2467E-7) * T + 7.6438E-5) * T - 4.0899E-3) * T + 8.24493E-1;
    double R3 = (-1.6546E-6 * T + 1.0227E-4) * T - 5.72466E-3;
    //  INTERNATIONAL ONE-ATMOSPHERE EQUATION OF STATE OF SEAWATER
    double SIG = (R4 * S + R3 * SR + R2) * S + R1;
    // SPECIFIC VOLUME AT ATMOSPHERIC PRESSURE
    double V350P = 1.0 / R3500;
    double SVA = -SIG * V350P / (R3500 + SIG);
    // double SIGMA = SIG + DR350;
    // double V0 = 1.0/(1000.0 + SIGMA);
    // SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS
    // double SVAN=SVA*1.0E+8;
 
    // ******************************************************************
    // ******  NEW HIGH PRESSURE EQUATION OF STATE FOR SEAWATER ********
    // ******************************************************************
    //        MILLERO, ET AL , 1980 DSR 27A, PP 255-264
    //               CONSTANT NOTATION FOLLOWS ARTICLE
    // ********************************************************
    // COMPUTE COMPRESSION TERMS
    double E = (9.1697E-10 * T + 2.0816E-8) * T - 9.9348E-7;
    double BW = (5.2787E-8 * T - 6.12293E-6) * T + 3.47718E-5;
    double B = BW + E * S; // Bulk Modulus (almost)
                           //  CORRECT B FOR ANAMOLY BIAS CHANGE
                           //    double Bout = B + 5.03217E-5;
 
    double D = 1.91075E-4;
    double C = (-1.6078E-6 * T - 1.0981E-5) * T + 2.2838E-3;
    double AW = ((-5.77905E-7 * T + 1.16092E-4) * T + 1.43713E-3) * T - 0.1194975;
    double A = (D * SR + C) * S + AW;
    //  CORRECT A FOR ANAMOLY BIAS CHANGE
    //    double Aout = A + 3.3594055;
 
    double B1 = (-5.3009E-4 * T + 1.6483E-2) * T + 7.944E-2;
    double A1 = ((-6.1670E-5 * T + 1.09987E-2) * T - 0.603459) * T + 54.6746;
    double KW = (((-5.155288E-5 * T + 1.360477E-2) * T - 2.327105) * T + 148.4206) * T - 1930.06;
    double K0 = (B1 * SR + A1) * S + KW;
 
    // EVALUATE PRESSURE POLYNOMIAL
    // ***********************************************
    //   K EQUALS THE SECANT BULK MODULUS OF SEAWATER
    //   DK=K(S,T,P)-K(35,0,P)
    //  K35=K(35,0,P)
    // ***********************************************
    double DK = (B * P + A) * P + K0;
    double K35 = (5.03217E-5 * P + 3.359406) * P + 21582.27;
    double GAM = P / K35;
    double PK = 1.0 - GAM;
    SVA = SVA * PK + (V350P + SVA) * P * DK / (K35 * (K35 + DK));
    //  SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS
    //SVAN=SVA*1.0E+8;      // Volume anomaly
    V350P = V350P * PK;
    // ****************************************************
    // COMPUTE DENSITY ANAMOLY WITH RESPECT TO 1000.0 KG/M**3
    //  1) DR350: DENSITY ANAMOLY AT 35 (IPSS-78), 0 DEG. C AND 0 DECIBARS
    //  2) DR35P: DENSITY ANAMOLY 35 (IPSS-78), 0 DEG. C ,  PRES. VARIATION
    //  3) DVAN : DENSITY ANAMOLY VARIATIONS INVOLVING SPECFIC VOL. ANAMOLY
    // ********************************************************************
    // CHECK VALUE: SIGMA = 59.82037  KG/M**3 FOR S = 40 (IPSS-78),
    // T = 40 DEG C, P0= 10000 DECIBARS.
    // *******************************************************
    double DR35P = GAM / V350P;
    double DVAN = SVA / (V350P * (V350P + SVA));
    double SIGMA = DR350 + DR35P - DVAN; // Density anomaly
 
    return SIGMA * si::kilograms_per_cubic_meter;
}
 
template <typename TemperatureUnit = boost::units::celsius::temperature,
          typename PressureUnit = decltype(boost::units::si::deci* bar)>
boost::units::quantity<boost::units::si::mass_density>
density_anomaly(double salinity,
                boost::units::quantity<boost::units::absolute<TemperatureUnit> > temperature,
                boost::units::quantity<PressureUnit> pressure)
{
    return density_anomaly(boost::units::quantity<boost::units::si::dimensionless>(salinity),
                           temperature, pressure);
}
 
} // namespace seawater
} // namespace util
} // namespace goby
 
#endif