001: /*
002: * GeoTools - OpenSource mapping toolkit
003: * http://geotools.org
004: * (C) 2003-2006, GeoTools Project Managment Committee (PMC)
005: * (C) 2001, Institut de Recherche pour le Développement
006: * (C) 1999, Pêches et Océans Canada
007: *
008: * This library is free software; you can redistribute it and/or
009: * modify it under the terms of the GNU Lesser General Public
010: * License as published by the Free Software Foundation;
011: * version 2.1 of the License.
012: *
013: * This library is distributed in the hope that it will be useful,
014: * but WITHOUT ANY WARRANTY; without even the implied warranty of
015: * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
016: * Lesser General Public License for more details.
017: *
018: * NOTE: permission has been given to the JScience project (http://www.jscience.org)
019: * to distribute this file under BSD-like license.
020: */
021: package org.geotools.nature;
022:
023: /**
024: * Sea water properties as a function of salinity, temperature and pressure.
025: * Density is computed using the 1980 definition of Equation of State (EOS80).
026: * Units are:
027: *
028: * <ul>
029: * <li>Salinity: Pratical Salinity Scale 1978 (PSS-78).</li>
030: * <li>Temperature: Celsius degrees according International Temperature Scale 1968 (ITS-68).</li>
031: * <li>Pressure: decibars (1 dbar = 10 kPa).
032: * </ul>
033: *
034: * @source $URL: http://svn.geotools.org/geotools/tags/2.4.1/modules/library/referencing/src/main/java/org/geotools/nature/SeaWater.java $
035: * @version $Id: SeaWater.java 24765 2007-03-15 03:50:56Z desruisseaux $
036: * @author Bernard Pelchat
037: * @author Martin Desruisseaux
038: *
039: * @since 2.1
040: */
041: public final class SeaWater {
042: /*
043: * Note: Les algorithmes originaux de l'UNESCO recevait en entrés
044: * des pressions en décibars. Les algorithmes écrites par
045: * Bernard Pelchat recevaient en entrés des pressions en
046: * MegaPascal. La première ligne de code des algorithmes
047: * de Bernard Pelchat multipliait donc les pressions par
048: * 100, afin de les convertir en decibars.
049: */
050:
051: /**
052: * Conductivity (in mS/cm) of a standard sea water sample.
053: * S is for <cite>Siemens</cite> (or Mho, its the same...).
054: */
055: public static final double STANDARD_CONDUCTIVITY = 42.914;
056:
057: /**
058: * Coéfficients de l'équation d'état EOS-80. La densité
059: * calculée par ces coéfficients est la densité Sigma-T.
060: */
061: private static final double EOS80_A[] = { -28.263737E+0,
062: 6.793952E-2, -9.095290E-3, 1.001685E-4, -1.120083E-6,
063: 6.536332E-9 }, EOS80_B[] = { 8.24493E-1, -4.0899E-3,
064: 7.6438E-5, -8.2467E-7, 5.3875E-9 }, EOS80_C[] = {
065: -5.72466E-3, 1.0227E-4, -1.6546E-6 }, EOS80_D = 4.8314E-4,
066: EOS80_E[] = { -1930.06E+0, 148.4206E+0, -2.327105E+0,
067: 1.360477E-2, -5.155288E-5 }, EOS80_F[] = {
068: 54.6746E+0, -6.03459E-1, 1.09987E-2, -6.1670E-5 },
069: EOS80_G[] = { 7.944E-2, 1.6483E-2, -5.3009E-4 },
070: EOS80_H[] = { -1.194975E-1, 1.43713E-3, 1.16092E-4,
071: -5.77905E-7 }, EOS80_I[] = { 2.2838E-3, -1.0981E-5,
072: -1.6078E-6 }, EOS80_J = 1.91075E-4, EOS80_K[] = {
073: 3.47718E-5, -6.12293E-6, 5.2787E-8 }, EOS80_M[] = {
074: -9.9348E-7, 2.0816E-8, 9.1697E-10 }, EOS80_N[] = {
075: 21582.27, 3.359406, 5.03217E-5 },
076: RHO_35_0_0 = 1028.1063, DR_35_0_0 = 28.106331;
077:
078: /**
079: * Coéfficients de l'équation d'état EOS-80. La densité
080: * calculée par ces coéfficients est la densité "vrai".
081: */
082: private static final double EOS80_At[] = { 999.842594, 6.793952E-2,
083: -9.095290E-3, 1.001685E-4, -1.120083E-6, 6.536332E-9 },
084: EOS80_Et[] = { 19652.21, 148.4206, -2.327105, 1.360477E-2,
085: -5.155288E-5 }, EOS80_Ht[] = { 3.239908,
086: 1.43713E-3, 1.16092E-4, -5.77905E-7 },
087: EOS80_Kt[] = { 8.50935E-5, -6.12293E-6, 5.2787E-8 };
088:
089: /**
090: * Coéfficients de l'équation de la salinité PSS-78.
091: */
092: private static final double PSS78_A[] = { 0.0080, -0.1692, 25.3851,
093: 14.0941, -7.0261, 2.7081 }, PSS78_B[] = { 0.0005, -0.0056,
094: -0.0066, -0.0375, 0.0636, -0.0144 },
095: PSS78_C[] = { 0.6766097, 2.00564E-2, 1.104259E-4,
096: -6.9698E-7, 1.0031E-9 }, PSS78_D[] = { 3.426E-2,
097: 4.464E-4, 4.215E-1, -3.107E-3 }, PSS78_E[] = {
098: 2.070E-5, -6.370E-10, 3.989E-15 }, PSS78_G[] = {
099: -0.1692, 50.7702, 42.2823, -28.1044, 13.5405 },
100: PSS78_H[] = { -0.0056 - 0.0132, -0.1125, 0.2544, -0.0720 },
101: PSS78_K = 0.0162;
102:
103: /**
104: * Coéfficients pour les salinités élevées,
105: */
106: private static final double PSS78_AR[] = { 7.737, -9.819, 8.663,
107: -2.625 }, PSS78_AT[] = { 3.473E-2, 3.188E-3, -4.655E-5 },
108: PSS78_CR[] = { -10.01E-2, 4.82E-2, -6.682E-4 };
109:
110: /**
111: * Constantes nécessaires au calcul de la chaleur spécifique.
112: *
113: * @see #specificHeat
114: */
115: private static final double HEAT_AA[] = { -7.643575, 0.1072763,
116: -1.38385E-3 }, HEAT_BB[] = { 0.1770383, -4.07718E-3,
117: 5.148E-5 }, HEAT_CC[] = { 4217.4, -3.720283, 0.1412855,
118: -2.654387E-3, 2.093236E-5 }, HEAT_A[] = { -4.9592E-1,
119: 1.45747E-2, -3.13885E-4, 2.0357E-6, 1.7168E-8 },
120: HEAT_B[] = { 2.4931E-4, -1.08645E-5, 2.87533E-7,
121: -4.0027E-9, 2.2956E-11 }, HEAT_C[] = { -5.422E-8,
122: 2.6380E-9, -6.5637E-11, 6.136E-13 }, HEAT_D[] = {
123: 4.9247E-3, -1.28315E-4, 9.802E-7, 2.5941E-8,
124: -2.9179E-10 }, HEAT_E[] = { -1.2331E-4, -1.517E-6,
125: 3.122E-8 }, HEAT_F[] = { -2.9558E-6, 1.17054E-7,
126: -2.3905E-9, 1.8448E-11 }, HEAT_G = 9.971E-8,
127: HEAT_H[] = { 5.540E-10, -1.7682E-11, 3.513E-13 },
128: HEAT_J = -1.4300E-12;
129:
130: /**
131: * Constantes nécessaires au calcul de la température adiabétique.
132: *
133: * @see #adiabeticTemperatureGradient
134: */
135: private static final double GRAD_A[] = { 3.5803E-05, 8.5258E-06,
136: -6.8360E-08, 6.6228E-10 }, GRAD_B[] = { 1.8932E-06,
137: -4.2393E-08 }, GRAD_C[] = { 1.8741E-08, -6.7795E-10,
138: 8.7330E-12, -5.4481E-14 }, GRAD_D[] = { -1.1351E-10,
139: 2.7759E-12 }, GRAD_E[] = { -4.6206E-13, 1.8676E-14,
140: -2.1687E-16 };
141:
142: /**
143: * Constantes nécessaires au calcul de la profondeur.
144: *
145: * @see #depth
146: */
147: private static final double DEPTH_C[] = { 9.72659, -2.2512E-5,
148: 2.279E-10, -1.82E-15 };
149:
150: /**
151: * Constantes nécessaires au calcul de la vitesse du son.
152: *
153: * @see #soundVelocity
154: */
155: private static final double SOUND_A0[] = { 1.389, -1.262E-2,
156: 7.164E-5, 2.006E-6, -3.21E-8 }, SOUND_A1[] = { 9.4742E-5,
157: -1.2580E-5, -6.4885E-8, 1.0507E-8, -2.0122E-10 },
158: SOUND_A2[] = { -3.9064E-7, 9.1041E-9, -1.6002E-10,
159: 7.988E-12 }, SOUND_A3[] = { 1.100E-10, 6.649E-12,
160: -3.389E-13 }, SOUND_B0[] = { -1.922E-2, -4.42E-5 },
161: SOUND_B1[] = { 7.3637E-5, 1.7945E-7 }, SOUND_C0[] = {
162: 1402.388, 5.03711, -5.80852E-2, 3.3420E-4,
163: -1.47800E-6, 3.1464E-9 }, SOUND_C1[] = { 0.153563,
164: 6.8982E-4, -8.1788E-6, 1.3621E-7, -6.1185E-10 },
165: SOUND_C2[] = { 3.1260E-5, -1.7107E-6, 2.5974E-8,
166: -2.5335E-10, 1.0405E-12 }, SOUND_C3[] = {
167: -9.7729E-9, 3.8504E-10, -2.3643E-12 },
168: SOUND_D0 = 1.727E-3, SOUND_D1 = -7.9836E-6;
169:
170: /**
171: * Constantes nécessaires au calcul de la saturation en oxygène dissous.
172: *
173: * @see #saturationO2
174: */
175: private static final double O2_AT[] = { -135.29996, 1.572288E+5,
176: -6.637149E+7, 1.243678E+10, -8.621061E+11 }, O2_AS[] = {
177: 0.020573, -12.142, 2363, 1 };
178:
179: /**
180: * Do not allow instantiation of this class.
181: */
182: private SeaWater() {
183: }
184:
185: /**
186: * Computes density as a function of salinity, temperature and pressure.
187: *
188: * @param S Salinity PSS-78 (0 to 42)
189: * @param T Temperature ITS-68 (-2 to 40°C)
190: * @param P Pressure in decibars (0 to 10<sup>5</sup> dbar), not including atmospheric pressure.
191: * @return Density (kg/m³).
192: */
193: public static double density(final double S, final double T,
194: double P) {
195: P /= 10.0;
196:
197: // Pure water density at atmospheric pressure
198: final double RHO_0_T_0 = polynome(T, EOS80_At);
199:
200: // Sea water density at atmospheric pressure
201: final double SR = Math.sqrt(S);
202: final double RHO_S_T_0 = (EOS80_D * S + polynome(T, EOS80_C)
203: * SR + polynome(T, EOS80_B))
204: * S + RHO_0_T_0;
205:
206: // Compression terms
207: final double K_S_T_0 = (polynome(T, EOS80_F) + polynome(T,
208: EOS80_G)
209: * SR)
210: * S + polynome(T, EOS80_Et);
211: final double K_S_T_P = K_S_T_0
212: + ((EOS80_J * SR + polynome(T, EOS80_I)) * S
213: + polynome(T, EOS80_Ht) + (polynome(T, EOS80_Kt) + polynome(
214: T, EOS80_M)
215: * S)
216: * P) * P;
217: return RHO_S_T_0 / (1.0 - P / K_S_T_P);
218: }
219:
220: /**
221: * Computes density sigma-T as a function of salinity, temperature and pressure.
222: * Density Sigma-T is equivalent to the true density minus 1000 kg/m³, and
223: * has typical values around 35. This computation avoid some rouding errors
224: * occuring in the true density computation.
225: *
226: * @param S Salinity PSS-78 (0 to 42)
227: * @param T Temperature ITS-68 (-2 to 40°C)
228: * @param P Pressure in decibars (0 to 10<sup>5</sup> dbar), not including atmospheric pressure.
229: * @return Density Sigma-T (kg/m³).
230: */
231: public static double densitySigmaT(final double S, final double T,
232: double P) {
233: P /= 10.0;
234: // Sea water density at atmospheric pressure
235: final double SR = Math.sqrt(S);
236: final double RHO = (EOS80_D * S + polynome(T, EOS80_C) * SR + polynome(
237: T, EOS80_B))
238: * S + polynome(T, EOS80_A);
239:
240: // Specific volume at atmospheric pressure
241: final double V_35_0_0 = 1.0 / RHO_35_0_0;
242: final double SVAN_S_T_0 = -RHO * V_35_0_0 / (RHO + RHO_35_0_0);
243: if (P <= 0) {
244: return RHO + DR_35_0_0;
245: }
246: // Compression terms, DK = K(S,T,P) - K(35,0,P)
247: final double K0 = (polynome(T, EOS80_F) + polynome(T, EOS80_G)
248: * SR)
249: * S + polynome(T, EOS80_E);
250: final double DK = K0
251: + (((EOS80_J * SR + polynome(T, EOS80_I)) * S + polynome(
252: T, EOS80_H)) + (polynome(T, EOS80_K) + polynome(
253: T, EOS80_M)
254: * S)
255: * P) * P;
256:
257: final double K_35_0_P = polynome(P, EOS80_N);
258: final double V_S_T_0 = SVAN_S_T_0 + V_35_0_0;
259: final double SVANS = SVAN_S_T_0 * (1.0 - P / K_35_0_P)
260: + V_S_T_0 * P * DK / (K_35_0_P * (K_35_0_P + DK));
261:
262: // Compute density anomaly
263: final double V_35_0_P = V_35_0_0 * (1.0 - P / K_35_0_P);
264: final double DR_35_0_P = P / (K_35_0_P * V_35_0_P);
265: final double DVAN = SVANS / (V_35_0_P * (V_35_0_P + SVANS));
266: return DR_35_0_0 + DR_35_0_P - DVAN;
267: }
268:
269: /**
270: * Computes volume as a function of salinity, temperature and pressure.
271: * This quantity if the inverse of density. This method is equivalent
272: * to <code>1/{@link #density density}(S,T,P)</code>.
273: *
274: * @param S Salinity PSS-78 (0 to 42)
275: * @param T Temperature ITS-68 (-2 to 40°C)
276: * @param P Pressure in decibars (0 to 10<sup>5</sup> dbar), not including atmospheric pressure.
277: * @return Volume (m³/kg).
278: */
279: public static double volume(final double S, final double T, double P) {
280: P /= 10.0;
281: // Sea water density at atmospheric pressure
282: final double SR = Math.sqrt(S);
283: final double RHO = (EOS80_D * S + polynome(T, EOS80_C) * SR + polynome(
284: T, EOS80_B))
285: * S + polynome(T, EOS80_A);
286:
287: // Specific volume at atmospheric pressure
288: final double V_35_0_0 = 1.0 / RHO_35_0_0;
289: final double SVAN_S_T_0 = -RHO * V_35_0_0 / (RHO + RHO_35_0_0);
290: if (P <= 0) {
291: return SVAN_S_T_0 + V_35_0_0;
292: }
293: // Compression terms, DK = K(S,T,P) - K(35,0,P)
294: final double K0 = (polynome(T, EOS80_F) + polynome(T, EOS80_G)
295: * SR)
296: * S + polynome(T, EOS80_E);
297: final double DK = K0
298: + (((EOS80_J * SR + polynome(T, EOS80_I)) * S + polynome(
299: T, EOS80_H)) + (polynome(T, EOS80_K) + polynome(
300: T, EOS80_M)
301: * S)
302: * P) * P;
303:
304: final double K_35_0_P = polynome(P, EOS80_N);
305: final double V_S_T_0 = SVAN_S_T_0 + V_35_0_0;
306: return (SVAN_S_T_0 + V_35_0_0) * (1.0 - P / K_35_0_P) + V_S_T_0
307: * P * DK / (K_35_0_P * (K_35_0_P + DK));
308: }
309:
310: /**
311: * Computes volumic anomaly as a function of salinity, temperature and pressure.
312: * Volumic anomaly is defined as the sea water sample's volume minus a standard
313: * sample's volume, where the standard sample is a sample of salinity 35, temperature
314: * 0°C and the same pressure. In pseudo-code, {@code volumeAnomaly} is equivalent
315: * to <code>{@link #volume volume}(S,T,P)-{@link #volume volume}(35,0,P)</code>.
316: *
317: * @param S Salinity PSS-78 (0 to 42)
318: * @param T Temperature ITS-68 (-2 to 40°C)
319: * @param P Pressure in decibars (0 to 10<sup>5</sup> dbar), not including atmospheric pressure.
320: * @return Volumic anomaly (m³/kg).
321: */
322: public static double volumeAnomaly(final double S, final double T,
323: double P) {
324: P /= 10.0;
325: // Sea water density at atmospheric pressure
326: final double SR = Math.sqrt(S);
327: final double RHO = (EOS80_D * S + polynome(T, EOS80_C) * SR + polynome(
328: T, EOS80_B))
329: * S + polynome(T, EOS80_A);
330:
331: // Specific volume at atmospheric pressure
332: final double V_35_0_0 = 1.0 / RHO_35_0_0;
333: final double SVAN_S_T_0 = -RHO * V_35_0_0 / (RHO + RHO_35_0_0);
334: if (P <= 0) {
335: return SVAN_S_T_0;
336: }
337: // Compression terms, DK = K(S,T,P) - K(35,0,P)
338: final double K0 = (polynome(T, EOS80_F) + polynome(T, EOS80_G)
339: * SR)
340: * S + polynome(T, EOS80_E);
341: final double DK = K0
342: + (((EOS80_J * SR + polynome(T, EOS80_I)) * S + polynome(
343: T, EOS80_H)) + (polynome(T, EOS80_K) + polynome(
344: T, EOS80_M)
345: * S)
346: * P) * P;
347:
348: final double K_35_0_P = polynome(P, EOS80_N);
349: final double V_S_T_0 = SVAN_S_T_0 + V_35_0_0;
350: return (SVAN_S_T_0 * (1.0 - P / K_35_0_P) + V_S_T_0 * P * DK
351: / (K_35_0_P * (K_35_0_P + DK)));
352: }
353:
354: /**
355: * Practical salinity scale 1978 definition
356: * with temperature correction, XR = SQRT( Rt )
357: */
358: private static double sal(double RT, double XT) {
359: return polynome(RT, PSS78_A) + (XT / (1.0 + PSS78_K * XT))
360: * polynome(RT, PSS78_B);
361: }
362:
363: /**
364: * {@code dsal(RT,XT)} function for derivative
365: * of {@code sal(RT,XT)} with <var>RT</var>.
366: */
367: private static double dsal(double RT, double XT) {
368: return polynome(RT, PSS78_G) + (XT / (1.0 + PSS78_K * XT))
369: * polynome(RT, PSS78_H);
370: }
371:
372: /**
373: * Computes salinity as a function of conductivity, temperature and pressure.
374: *
375: * @param C Conductivity in mS/cm (millisiemens by centimeters). Multiply
376: * par {@link #STANDARD_CONDUCTIVITY} if {@code C} is not a
377: * real conductivity, but instead the ratio between the sample's
378: * conductivity and the standard sample's conductivity.
379: * @param T Temperature ITS-68 (-2 to 40°C)
380: * @param P Pressure in decibars (0 to 10<sup>5</sup> dbar), not including atmospheric pressure.
381: * @return Salinity PSS-78.
382: *
383: * @todo What to do with pression!?! Check the equation of state.
384: */
385: public static double salinity(double C, final double T,
386: final double P) {
387: C /= STANDARD_CONDUCTIVITY;
388: if (!(C < 5E-4)) { // use '!' in order to accept NaN
389: final double XR = Math
390: .sqrt(C
391: / (polynome(T, PSS78_C) * (1.0 + polynome(
392: P, PSS78_E)
393: * P
394: / ((PSS78_D[1] * T + PSS78_D[0])
395: * T + 1.0 + (PSS78_D[3] * T + PSS78_D[2])
396: * C))));
397: final double S = sal(XR, T - 15.0); // Do not use an 'assert' statement invoking 'cond'.
398: if (!(S >= 42))
399: return S; // use '!' to accept NaN
400: /*
401: * Calcule la salinité pour une eau de conductivité,
402: * de température et de pression données. Cet algorithme
403: * doit être utilisé lorsque l'on s'attend à une salinité
404: * entre 42 et 50.
405: */
406: return 35
407: * C
408: + C
409: * (C - 1)
410: * (polynome(C, PSS78_AR) + T
411: * (polynome(T, PSS78_AT) + C
412: * (PSS78_CR[0] + PSS78_CR[1] * C + PSS78_CR[2]
413: * T)));
414: // TODO: VERIFIER CE QUE DEVIENT LA PRESSION ET IMPLEMENTER L'EQUATION D'ETAT.
415: } else {
416: return 0; // Zero conductivity trap
417: }
418: }
419:
420: /**
421: * Computes conductivity as a function of salinity, temperature and pressure.
422: *
423: * @param S Salinity PSS-78 (0 to 42)
424: * @param T Temperature ITS-68 (-2 to 40°C)
425: * @param P Pressure (0 to 10<sup>5</sup> dbar), not including atmospheric pressure.
426: * @return Conductivity in mS/cm.
427: */
428: public static double conductivity(final double S, final double T,
429: final double P) {
430: if (!(S < 0.02)) { // use '!' in order to accept NaN
431: double XT = T - 15.0;
432: double RT = Math.sqrt(S / 35.0); // First approximation
433: double SI = sal(RT, XT);
434: for (int n = 0; n < 10; n++) { // Iteration loop begin here with a maximum of 10 cycles
435: RT += (S - SI) / dsal(RT, XT);
436: SI = sal(RT, XT);
437: if (Math.abs(SI - S) < 1E-4)
438: break;
439: }
440: double RTT = polynome(T, PSS78_C) * (RT * RT);
441: double AT = PSS78_D[3] * T + PSS78_D[2];
442: double BT = (PSS78_D[1] * T + PSS78_D[0]) * T + 1.0;
443: double CP = RTT * (BT + polynome(P, PSS78_E) * P);
444: BT -= RTT * AT;
445: // Solve quadratic equation for C = RT35*RT*(1+C/AR+b)
446: double cnd = 0.5
447: * (Math.sqrt(Math.abs((BT * BT) + 4.0 * AT * CP)) - BT)
448: / AT;
449: return cnd * STANDARD_CONDUCTIVITY;
450: } else {
451: return 0; // Zero salinity trap
452: }
453: }
454:
455: /**
456: * Computes specific heat as a function of salinity, temperature and pressure.
457: *
458: * @param S Salinity PSS-78.
459: * @param T Temperature (°C).
460: * @param P Pressure (dbar), not including atmospheric pressure.
461: * @return Specific heat (J/(kg×°C)).
462: */
463: public static double specificHeat(final double S, final double T,
464: double P) {
465: P /= 10.0;
466: final double SR = Math.sqrt(S);
467: return (polynome(T, HEAT_CC)
468: + (polynome(T, HEAT_BB) * SR + polynome(T, HEAT_AA))
469: * S
470: + (((polynome(T, HEAT_C) * P + polynome(T, HEAT_B)) * P + polynome(
471: T, HEAT_A)) * P) + ((((HEAT_J * SR + polynome(
472: T, HEAT_H))
473: * S * P + (HEAT_G * SR + polynome(T, HEAT_F)) * S)
474: * P + (polynome(T, HEAT_E) * SR + polynome(T, HEAT_D))
475: * S) * P));
476: }
477:
478: /**
479: * Computes fusion temperature (melting point) as a function of salinity and pressure.
480: *
481: * @param S Salinity PSS-78.
482: * @param P Pressure (dbar), not including atmospheric pressure.
483: * @return Melting point (°C).
484: */
485: public static double fusionTemperature(final double S,
486: final double P) {
487: return (-0.0575 + 1.710523E-3 * Math.sqrt(S) + -2.154996E-4 * S)
488: * S + -7.53E-4 * P;
489: }
490:
491: /**
492: * Computes adiabetic temperature gradient as a function of salinity, temperature and pressure.
493: *
494: * @param S Salinity PSS-78.
495: * @param T Temperature (°C).
496: * @param P Pressure (dbar), not including atmospheric pressure.
497: * @return Adiabetic temperature gradient (°C/dbar).
498: */
499: public static double adiabeticTemperatureGradient(double S,
500: final double T, final double P) {
501: S -= 35.0;
502: return (polynome(T, GRAD_A) + polynome(T, GRAD_B) * S + (polynome(
503: T, GRAD_C)
504: + polynome(T, GRAD_D) * S + polynome(T, GRAD_E) * P)
505: * P);
506: }
507:
508: /**
509: * Computes depth as a function of pressure and latitude.
510: *
511: * @param P Pressure (dbar), not including atmospheric pressure.
512: * @param lat Latitude in degrees (-90 to 90°)
513: * @return Depth (m).
514: */
515: public static double depth(final double P, double lat) {
516: lat = Math.sin(lat);
517: lat *= lat;
518: lat = 9.780318 * (1.0 + 5.2788E-3 * lat + 2.36E-5 * (lat * lat));
519: return polynome(P, DEPTH_C) * P / (lat + (0.5 * 2.184E-6) * P);
520: }
521:
522: /**
523: * Computes sound velocity as a function of salinity, temperature and pressure.
524: *
525: * @param S Salinity PSS-78.
526: * @param T Temperature (°C).
527: * @param P Pressure (dbar), not including atmospheric pressure.
528: * @return Sound velocity (m/s).
529: */
530: public static double soundVelocity(final double S, final double T,
531: final double P) {
532: // S^0 terms
533: final double CW = ((polynome(T, SOUND_C3) * P + polynome(T,
534: SOUND_C2))
535: * P + polynome(T, SOUND_C1))
536: * P + polynome(T, SOUND_C0);
537: // S^1 terms
538: final double A = ((polynome(T, SOUND_A3) * P + polynome(T,
539: SOUND_A2))
540: * P + polynome(T, SOUND_A1))
541: * P + polynome(T, SOUND_A0);
542: // S^3/2 terms
543: final double B = polynome(T, SOUND_B0) + polynome(T, SOUND_B1)
544: * P;
545:
546: // S^2 terms
547: final double D = SOUND_D0 + SOUND_D1 * P;
548:
549: // sound speed return
550: return CW + (D * S + B * Math.sqrt(S) + A) * S;
551: }
552:
553: /**
554: * Computes saturation in disolved oxygen as a function of salinity and temperature.
555: *
556: * @param S Salinity PSS-78.
557: * @param T Temperature (°C).
558: * @return Saturation in disolved oxygen (µmol/kg).
559: */
560: public static double saturationO2(final double S, double T) {
561: T += 273.15;
562: return Math.exp(polynome_neg(T, O2_AT) + S
563: * polynome_neg(T, O2_AS));
564: }
565:
566: /**
567: * Calcule la valeur d'un polynôme.
568: * Cette fonction calcule la valeur de:
569: *
570: * <blockquote><pre>
571: * y = C[0] + C[1]*x + C[2]*x² + C[3]*x³
572: * </pre></blockquote>
573: *
574: * où C est un vecteur de coéfficients transmis en argument.
575: * Une exception sera levée si ce tableau ne contient pas
576: * au moins 1 élément.
577: *
578: * @param x Valeur x à laquelle calculer le polynôme.
579: * @param c Coéfficients C du polynôme.
580: * @return La valeur du polynôme au x spécifié.
581: *
582: * @see #poly_inv(double,double[])
583: */
584: private static double polynome(final double x, final double c[]) {
585: int n = c.length - 1;
586: double y = c[n];
587: while (n > 0) {
588: y = y * x + c[--n];
589: }
590: return y;
591: }
592:
593: /**
594: * Calcule la valeur de:
595: *
596: * <blockquote><pre>
597: * y = C[0] + C[1]/x + C[2]/x² + C[3]/x³
598: * </pre></blockquote>
599: *
600: * où C est un vecteur de coéfficients transmis en argument.
601: * Une exception sera levée si ce tableau ne contient pas
602: * au moins 1 élément.
603: *
604: * @param x Valeur x à laquelle calculer le polynôme.
605: * @param C Coéfficients C du polynôme.
606: * @return La valeur du polynôme au x spécifié.
607: *
608: * @see #polynome(double,double[])
609: */
610: private static double polynome_neg(final double x, final double c[]) {
611: int n = c.length - 1;
612: double y = c[n];
613: while (n > 0) {
614: y = y / x + c[--n];
615: }
616: return y;
617: }
618: }
|