Goby3  3.1.5a
2024.05.23
swstate.h
Go to the documentation of this file.
1 // Copyright 2013-2021:
2 // GobySoft, LLC (2013-)
3 // Massachusetts Institute of Technology (2007-2014)
4 // Community contributors (see AUTHORS file)
5 // File authors:
6 // Toby Schneider <toby@gobysoft.org>
7 //
8 //
9 // This file is part of the Goby Underwater Autonomy Project Libraries
10 // ("The Goby Libraries").
11 //
12 // The Goby Libraries are free software: you can redistribute them and/or modify
13 // them under the terms of the GNU Lesser General Public License as published by
14 // the Free Software Foundation, either version 2.1 of the License, or
15 // (at your option) any later version.
16 //
17 // The Goby Libraries are distributed in the hope that they will be useful,
18 // but WITHOUT ANY WARRANTY; without even the implied warranty of
19 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
20 // GNU Lesser General Public License for more details.
21 //
22 // You should have received a copy of the GNU Lesser General Public License
23 // along with Goby. If not, see <http://www.gnu.org/licenses/>.
24 
25 // modified for C++ by s. petillo spetillo@mit.edu
26 // ocean engineering graduate student - mit / whoi joint program
27 // massachusetts institute of technology (mit)
28 // laboratory for autonomous marine sensing systems (lamss)
29 
30 #ifndef GOBY_UTIL_SEAWATER_SWSTATE_H
31 #define GOBY_UTIL_SEAWATER_SWSTATE_H
32 
33 #include <cmath>
34 
35 #include <boost/units/quantity.hpp>
36 #include <boost/units/systems/si.hpp>
37 #include <boost/units/systems/si/prefixes.hpp>
38 #include <boost/units/systems/temperature/celsius.hpp>
39 
40 #include "units.h"
41 
42 namespace goby
43 {
44 namespace util
45 {
46 namespace seawater
47 {
55 template <typename DimensionlessUnit = boost::units::si::dimensionless,
56  typename TemperatureUnit = boost::units::celsius::temperature,
57  typename PressureUnit = decltype(boost::units::si::deci* bar)>
58 boost::units::quantity<boost::units::si::mass_density>
59 density_anomaly(boost::units::quantity<DimensionlessUnit> salinity,
60  boost::units::quantity<boost::units::absolute<TemperatureUnit> > temperature,
61  boost::units::quantity<PressureUnit> pressure)
62 {
63  using namespace boost::units;
64 
65  double S = salinity;
66  double T = quantity<absolute<celsius::temperature> >(temperature).value();
67  double P0 = quantity<decltype(si::deci * bar)>(pressure).value();
68 
69  /*
70 
71  SIGMA = density_anomaly(S,T,P) returns the density anomaly SIGMA (kg/m^3)
72  given the salinity S (ppt), temperature T (deg C) and pressure
73  P (dbars).
74 
75  Notes: RP (WHOI 2/dec/91)
76 
77  This stuff is directly copied from the UNESCO algorithms, with some
78  minor changes to make it Matlab compatible (like adding ";" and changing
79  "*" to "*" when necessary.
80 
81  RP (WHOI 3/dec/91)
82 
83  ******************************************************
84  SPECIFIC VOLUME ANOMALY (STERIC ANOMALY) BASED ON 1980 EQUATION
85  OF STATE FOR SEAWATER AND 1978 PRACTICAL SALINITY SCALE.
86  REFERENCES
87  MILLERO, ET AL (1980) DEEP-SEA RES.,27A,255-264
88  MILLERO AND POISSON 1981,DEEP-SEA RES.,28A PP 625-629.
89  BOTH ABOVE REFERENCES ARE ALSO FOUND IN UNESCO REPORT 38 (1981)
90 
91  UNITS:
92  PRESSURE P0 DECIBARS
93  TEMPERATURE T DEG CELSIUS (IPTS-68)
94  SALINITY S (IPSS-78)
95  SPEC. VOL. ANA. SVAN M**3/KG *1.0E-8
96  DENSITY ANA. SIGMA KG/M**3
97  ******************************************************************
98  CHECK VALUE: SVAN=981.3021 E-8 M**3/KG. FOR S = 40 (IPSS-78) ,
99  T = 40 DEG C, P0= 10000 DECIBARS.
100  CHECK VALUE: SIGMA = 59.82037 KG/M**3 FOR S = 40 (IPSS-78) ,
101  T = 40 DEG C, P0= 10000 DECIBARS.
102  CHECK VALUE: FOR S = 40 (IPSS-78) , T = 40 DEG C, P0= 10000 DECIBARS.
103  DR/DP DR/DT DR/DS
104  DRV(1,7) DRV(2,3) DRV(1,8)
105 
106  FINITE DIFFERENCE WITH 3RD ORDER CORRECTION DONE IN DOUBLE PRECSION
107 
108  3.46969238E-3 -.43311722 .705110777
109 
110  EXPLICIT DIFFERENTIATION SINGLE PRECISION FORMULATION EOS80
111 
112  3.4696929E-3 -.4331173 .7051107
113 
114  (RP...I think this ---------^^^^^^ should be -.4431173!);
115  */
116 
117  // *******************************************************
118  using namespace std;
119 
120  // DATA
121  double R3500 = 1028.1063;
122  double R4 = 4.8314E-4;
123  double DR350 = 28.106331;
124 
125  // CONVERT PRESSURE TO BARS AND TAKE SQUARE ROOT SALINITY.
126  double P = P0 / 10.;
127  //double SAL=S;
128  double SR = sqrt(abs(S));
129  // *********************************************************
130  // PURE WATER DENSITY AT ATMOSPHERIC PRESSURE
131  // BIGG P.H.,(1967) BR. J. APPLIED PHYSICS 8 PP 521-537.
132 
133  double R1 = ((((6.536332E-9 * T - 1.120083E-6) * T + 1.001685E-4) * T - 9.095290E-3) * T +
134  6.793952E-2) *
135  T -
136  28.263737;
137  // SEAWATER DENSITY ATM PRESS.
138  // COEFFICIENTS INVOLVING SALINITY
139  double R2 = (((5.3875E-9 * T - 8.2467E-7) * T + 7.6438E-5) * T - 4.0899E-3) * T + 8.24493E-1;
140  double R3 = (-1.6546E-6 * T + 1.0227E-4) * T - 5.72466E-3;
141  // INTERNATIONAL ONE-ATMOSPHERE EQUATION OF STATE OF SEAWATER
142  double SIG = (R4 * S + R3 * SR + R2) * S + R1;
143  // SPECIFIC VOLUME AT ATMOSPHERIC PRESSURE
144  double V350P = 1.0 / R3500;
145  double SVA = -SIG * V350P / (R3500 + SIG);
146  // double SIGMA = SIG + DR350;
147  // double V0 = 1.0/(1000.0 + SIGMA);
148  // SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS
149  // double SVAN=SVA*1.0E+8;
150 
151  // ******************************************************************
152  // ****** NEW HIGH PRESSURE EQUATION OF STATE FOR SEAWATER ********
153  // ******************************************************************
154  // MILLERO, ET AL , 1980 DSR 27A, PP 255-264
155  // CONSTANT NOTATION FOLLOWS ARTICLE
156  // ********************************************************
157  // COMPUTE COMPRESSION TERMS
158  double E = (9.1697E-10 * T + 2.0816E-8) * T - 9.9348E-7;
159  double BW = (5.2787E-8 * T - 6.12293E-6) * T + 3.47718E-5;
160  double B = BW + E * S; // Bulk Modulus (almost)
161  // CORRECT B FOR ANAMOLY BIAS CHANGE
162  // double Bout = B + 5.03217E-5;
163 
164  double D = 1.91075E-4;
165  double C = (-1.6078E-6 * T - 1.0981E-5) * T + 2.2838E-3;
166  double AW = ((-5.77905E-7 * T + 1.16092E-4) * T + 1.43713E-3) * T - 0.1194975;
167  double A = (D * SR + C) * S + AW;
168  // CORRECT A FOR ANAMOLY BIAS CHANGE
169  // double Aout = A + 3.3594055;
170 
171  double B1 = (-5.3009E-4 * T + 1.6483E-2) * T + 7.944E-2;
172  double A1 = ((-6.1670E-5 * T + 1.09987E-2) * T - 0.603459) * T + 54.6746;
173  double KW = (((-5.155288E-5 * T + 1.360477E-2) * T - 2.327105) * T + 148.4206) * T - 1930.06;
174  double K0 = (B1 * SR + A1) * S + KW;
175 
176  // EVALUATE PRESSURE POLYNOMIAL
177  // ***********************************************
178  // K EQUALS THE SECANT BULK MODULUS OF SEAWATER
179  // DK=K(S,T,P)-K(35,0,P)
180  // K35=K(35,0,P)
181  // ***********************************************
182  double DK = (B * P + A) * P + K0;
183  double K35 = (5.03217E-5 * P + 3.359406) * P + 21582.27;
184  double GAM = P / K35;
185  double PK = 1.0 - GAM;
186  SVA = SVA * PK + (V350P + SVA) * P * DK / (K35 * (K35 + DK));
187  // SCALE SPECIFIC VOL. ANAMOLY TO NORMALLY REPORTED UNITS
188  //SVAN=SVA*1.0E+8; // Volume anomaly
189  V350P = V350P * PK;
190  // ****************************************************
191  // COMPUTE DENSITY ANAMOLY WITH RESPECT TO 1000.0 KG/M**3
192  // 1) DR350: DENSITY ANAMOLY AT 35 (IPSS-78), 0 DEG. C AND 0 DECIBARS
193  // 2) DR35P: DENSITY ANAMOLY 35 (IPSS-78), 0 DEG. C , PRES. VARIATION
194  // 3) DVAN : DENSITY ANAMOLY VARIATIONS INVOLVING SPECFIC VOL. ANAMOLY
195  // ********************************************************************
196  // CHECK VALUE: SIGMA = 59.82037 KG/M**3 FOR S = 40 (IPSS-78),
197  // T = 40 DEG C, P0= 10000 DECIBARS.
198  // *******************************************************
199  double DR35P = GAM / V350P;
200  double DVAN = SVA / (V350P * (V350P + SVA));
201  double SIGMA = DR350 + DR35P - DVAN; // Density anomaly
202 
203  return SIGMA * si::kilograms_per_cubic_meter;
204 }
205 
213 template <typename TemperatureUnit = boost::units::celsius::temperature,
214  typename PressureUnit = decltype(boost::units::si::deci* bar)>
215 boost::units::quantity<boost::units::si::mass_density>
216 density_anomaly(double salinity,
217  boost::units::quantity<boost::units::absolute<TemperatureUnit> > temperature,
218  boost::units::quantity<PressureUnit> pressure)
219 {
220  return density_anomaly(boost::units::quantity<boost::units::si::dimensionless>(salinity),
221  temperature, pressure);
222 }
223 
224 } // namespace seawater
225 } // namespace util
226 } // namespace goby
227 
228 #endif
goby
The global namespace for the Goby project.
Definition: acomms_constants.h:33
prefixes.hpp
goby::util::seawater::density_anomaly
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)
Definition: swstate.h:59
goby::util::seawater::bar
static const boost::units::metric::bar_base_unit::unit_type bar
Definition: units.h:36
goby::util::units::rpm::dimensionless
boost::units::unit< boost::units::dimensionless_type, system > dimensionless
Definition: system.hpp:47
units.h
goby::util::seawater::salinity
boost::units::quantity< boost::units::si::dimensionless > salinity(boost::units::quantity< ConductivityUnit > conductivity, boost::units::quantity< boost::units::absolute< TemperatureUnit > > temperature, boost::units::quantity< PressureUnit > pressure)
Calculates salinity from conductivity, temperature, and pressure Adapted from "Algorithms for computa...
Definition: salinity.h:58
goby::util::seawater::pressure
boost::units::quantity< decltype(boost::units::si::deci *bar)> pressure(boost::units::quantity< DepthUnit > depth, boost::units::quantity< LatitudeUnit > latitude)
Calculates pressure from depth and latitude.
Definition: pressure.h:55
boost::units
Definition: time.hpp:18
B