use crate::gsw_internal_const::{DB2PA, GSW_PU, GSW_SFAC};
use crate::gsw_sp_coefficients::*;
use crate::gsw_specvol_coefficients::{V005, V006};
use crate::{Error, Result};
const G0: f64 = 2.641_463_563_366_498e-1;
const G1: f64 = 2.007_883_247_811_176e-4;
const G2: f64 = -4.107_694_432_853_053e-6;
const G3: f64 = 8.401_670_882_091_225e-8;
const G4: f64 = -1.711_392_021_989_21e-9;
const G5: f64 = 3.374_193_893_377_38e-11;
const G6: f64 = -5.923_731_174_730_784e-13;
const G7: f64 = 8.057_771_569_962_299e-15;
const G8: f64 = -7.054_313_817_447_962e-17;
const G9: f64 = 2.859_992_717_347_235e-19;
#[inline]
pub(crate) fn non_dimensional_p(p: f64) -> f64 {
if cfg!(feature = "compat") {
p * 1e-4
} else {
p / GSW_PU
}
}
pub fn specvol_sso_0(p: f64) -> f64 {
const VXX0: f64 = if cfg!(feature = "compat") {
9.726_613_854_843_87e-4
} else {
9.726_613_854_843_871e-4
};
const VXX1: f64 = if cfg!(feature = "compat") {
-4.505_913_211_160_929e-5
} else {
-4.505_913_211_160_931e-5
};
const VXX2: f64 = if cfg!(feature = "compat") {
7.130_728_965_927_127e-6
} else {
7.130_728_965_927_128e-6
};
const VXX3: f64 = if cfg!(feature = "compat") {
-6.657_179_479_768_312e-7
} else {
-6.657_179_479_768_313e-7
};
const VXX4: f64 = if cfg!(feature = "compat") {
-2.994_054_447_232_88e-8
} else {
-2.994_054_447_232_877_6e-8
};
let p = non_dimensional_p(p);
VXX0 + p * (VXX1 + p * (VXX2 + p * (VXX3 + p * (VXX4 + p * (V005 + V006 * p)))))
}
#[cfg(test)]
mod test_specvol_sso_0 {
use super::specvol_sso_0;
use crate::gsw_internal_const::GSW_SSO;
use crate::volume::specvol;
#[test]
fn specvol_vs_specvol_sso_0() {
let p_to_test: [f64; 5] = [0., 10., 100., 1_000., 5_000.];
for p in p_to_test.iter().cloned() {
let specvol = specvol(GSW_SSO, 0., p).unwrap();
let specvol_sso_0 = specvol_sso_0(p);
assert!((specvol - specvol_sso_0).abs() < f64::EPSILON);
}
}
}
pub(crate) fn enthalpy_sso_0(p: f64) -> f64 {
const H006: f64 = -2.107_876_881_00e-9;
const H007: f64 = 2.801_929_132_90e-10;
let z = p * 1.0e-4;
let dynamic_enthalpy_sso_0_p = z
* (9.726_613_854_843_87e-4
+ z * (-2.252_956_605_630_465e-5
+ z * (2.376_909_655_387_404e-6
+ z * (-1.664_294_869_986_011e-7
+ z * (-5.988_108_894_465_758e-9 + z * (H006 + H007 * z))))));
dynamic_enthalpy_sso_0_p * DB2PA * 1.0e4
}
pub(crate) fn hill_ratio_at_sp2(t: f64) -> f64 {
let sp2 = 2.0;
let t68: f64 = crate::practical_salinity::t68_from_t90(t);
let ft68: f64 = (t68 - 15.0) / (1.0 + K * (t68 - 15.0));
let rtx0: f64 = G0
+ t68
* (G1
+ t68
* (G2
+ t68
* (G3
+ t68
* (G4
+ t68
* (G5
+ t68
* (G6
+ t68 * (G7 + t68 * (G8 + t68 * G9))))))));
let dsp_drtx: f64 = A1
+ (2.0 * A2 + (3.0 * A3 + (4.0 * A4 + 5.0 * A5 * rtx0) * rtx0) * rtx0) * rtx0
+ ft68 * (B1 + (2.0 * B2 + (3.0 * B3 + (4.0 * B4 + 5.0 * B5 * rtx0) * rtx0) * rtx0) * rtx0);
let sp_est: f64 = A0
+ (A1 + (A2 + (A3 + (A4 + A5 * rtx0) * rtx0) * rtx0) * rtx0) * rtx0
+ ft68 * (B0 + (B1 + (B2 + (B3 + (B4 + B5 * rtx0) * rtx0) * rtx0) * rtx0) * rtx0);
let rtx: f64 = rtx0 - (sp_est - sp2) / dsp_drtx;
let rtxm: f64 = 0.5 * (rtx + rtx0);
let dsp_drtx: f64 = A1
+ (2.0 * A2 + (3.0 * A3 + (4.0 * A4 + 5.0 * A5 * rtxm) * rtxm) * rtxm) * rtxm
+ ft68 * (B1 + (2.0 * B2 + (3.0 * B3 + (4.0 * B4 + 5.0 * B5 * rtxm) * rtxm) * rtxm) * rtxm);
let rtx: f64 = rtx0 - (sp_est - sp2) / dsp_drtx;
let x: f64 = 400.0 * rtx * rtx;
let sqrty: f64 = 10.0 * rtx;
let part1: f64 = 1.0 + x * (1.5 + x);
let part2: f64 = 1.0 + sqrty * (1.0 + sqrty * (1.0 + sqrty));
let sp_hill_raw_at_sp2: f64 = sp2 - A0 / part1 - B0 * ft68 / part2;
2.0 / sp_hill_raw_at_sp2
}
#[cfg(test)]
mod test_hill_ratio_at_sp2 {
use super::hill_ratio_at_sp2;
#[test]
fn example_values() {
let ratio = hill_ratio_at_sp2(0.);
assert!((ratio - 0.999_834_364_412_733_2).abs() <= f64::EPSILON);
let ratio = hill_ratio_at_sp2(10.);
assert!((ratio - 0.999_958_676_175_969_7).abs() <= f64::EPSILON);
let ratio = hill_ratio_at_sp2(25.);
assert!((ratio - 1.000_076_483_096_847_2).abs() <= f64::EPSILON);
let ratio = hill_ratio_at_sp2(30.);
assert!((ratio - 1.000_104_961_868_661_5).abs() <= f64::EPSILON,);
let ratio = hill_ratio_at_sp2(40.);
assert!((ratio - 1.000_151_589_022_029_8).abs() <= f64::EPSILON);
let ratio = hill_ratio_at_sp2(50.);
assert!((ratio - 1.000_188_151_345_468_5).abs() <= f64::EPSILON);
}
}
#[allow(dead_code, clippy::excessive_precision, clippy::manual_range_contains)]
pub(crate) fn gibbs(ns: u8, nt: u8, np: u8, sa: f64, t: f64, p: f64) -> Result<f64> {
if sa.is_nan() | t.is_nan() | p.is_nan() {
return Ok(f64::NAN);
}
let sa: f64 = if sa < 0.0 {
if cfg!(feature = "compat") {
0.0
} else if cfg!(feature = "invalidasnan") {
return Ok(f64::NAN);
} else {
return Err(Error::NegativeSalinity);
}
} else {
sa
};
if (sa > 120.0) | (t < -12.0) | (t > 80.0) | (p > 12_000.0) {
if cfg!(feature = "invalidasnan") {
return Ok(f64::NAN);
} else {
return Err(Error::Undefined);
}
}
let x2 = GSW_SFAC * sa;
let x = libm::sqrt(x2);
let y = t * 0.025;
let z = non_dimensional_p(p);
if (ns == 0) & (nt == 0) & (np == 0) {
let g03 = 101.342_743_139_674
+ z * (100_015.695_367_145
+ z * (-2_544.576_542_036_3
+ z * (284.517_778_446_287
+ z * (-33.314_675_425_361_1
+ (4.202_631_088_030_84 - 0.546_428_511_471_039 * z) * z))))
+ y * (5.905_783_479_094_02
+ z * (-270.983_805_184_062
+ z * (776.153_611_613_101
+ z * (-196.512_550_881_22
+ (28.979_652_629_417_5 - 2.132_900_835_183_27 * z) * z)))
+ y * (-12_357.785_933_039
+ z * (1_455.036_454_046_8
+ z * (-756.558_385_769_359
+ z * (273.479_662_323_528
+ z * (-55.560_406_381_721_8 + 4.344_206_719_171_97 * z))))
+ y * (736.741_204_151_612
+ z * (-672.507_783_145_07
+ z * (499.360_390_819_152
+ z * (-239.545_330_654_412
+ (48.801_251_859_387_2 - 1.663_071_062_089_05 * z) * z)))
+ y * (-148.185_936_433_658
+ z * (397.968_445_406_972
+ z * (-301.815_380_621_876
+ (152.196_371_733_841 - 26.374_837_723_280_2 * z) * z))
+ y * (58.025_912_584_257_1
+ z * (-194.618_310_617_595
+ z * (120.520_654_902_025
+ z * (-55.272_305_234_015_2
+ 6.481_906_680_772_21 * z)))
+ y * (-18.984_384_651_417_2
+ y * (3.050_816_464_879_67 - 9.631_081_193_930_62 * z)
+ z * (63.511_393_664_178_5
+ z * (-22.289_731_714_045_9
+ 8.170_605_418_181_12 * z))))))));
let mut g08 = x2
* (1_416.276_484_841_97
+ z * (-3_310.491_540_448_39
+ z * (384.794_152_978_599
+ z * (-96.532_432_010_745_8
+ (15.840_817_276_682_4 - 2.624_801_565_909_92 * z) * z)))
+ x * (-2_432.146_623_817_94
+ x * (2_025.801_156_036_97
+ y * (543.835_333_000_098
+ y * (-68.557_250_920_449_1
+ y * (49.366_769_485_625_4
+ y * (-17.139_757_741_978_8 + 2.496_970_095_695_08 * y)))
- 22.668_355_851_282_9 * z)
+ x * (-1_091.668_410_429_67 - 196.028_306_689_776 * y
+ x * (374.601_237_877_84 - 48.589_106_902_540_9 * x
+ 36.757_162_299_580_5 * y)
+ 36.028_419_561_108_6 * z)
+ z * (-54.791_913_353_288_7
+ (-4.081_939_789_122_61 - 30.175_511_197_116_1 * z) * z))
+ z * (199.459_603_073_901
+ z * (-52.294_090_928_133_5
+ (68.044_494_272_645_9 - 3.412_519_324_412_82 * z) * z))
+ y * (-493.407_510_141_682
+ z * (-175.292_041_186_547
+ (83.192_392_780_181_9 - 29.483_064_349_429 * z) * z)
+ y * (-43.066_467_597_804_2
+ z * (383.058_066_002_476
+ z * (-54.191_726_251_711_2 + 25.639_848_738_991_4 * z))
+ y * (-10.022_737_086_187_5 - 460.319_931_801_257 * z
+ y * (0.875_600_661_808_945 + 234.565_187_611_355 * z)))))
+ y * (168.072_408_311_545
+ z * (729.116_529_735_046
+ z * (-343.956_902_961_561
+ z * (124.687_671_116_248
+ z * (-31.656_964_386_073 + 7.046_588_033_154_49 * z))))
+ y * (880.031_352_997_204
+ y * (-225.267_649_263_401
+ y * (91.426_044_775_125_9
+ y * (-21.660_324_087_531_1 + 2.130_169_708_471_83 * y)
+ z * (-297.728_741_987_187
+ (74.726_141_138_756 - 36.487_291_900_158_8 * z) * z))
+ z * (694.244_814_133_268
+ z * (-204.889_641_964_903
+ (113.561_697_840_594 - 11.128_273_432_641_3 * z) * z)))
+ z * (-860.764_303_783_977
+ z * (337.409_530_269_367
+ z * (-178.314_556_207_638
+ (44.204_035_830_8 - 7.920_015_472_116_82 * z) * z))))));
if x > 0.0 {
g08 += x2 * (5_812.814_566_267_32 + 851.226_734_946_706 * y) * libm::log(x);
}
Ok(g03 + g08)
} else if (ns == 1) & (nt == 0) & (np == 0) {
let mut g08 = 8_645.367_535_951_26
+ z * (-6_620.983_080_896_78
+ z * (769.588_305_957_198
+ z * (-193.064_864_021_491_6
+ (31.681_634_553_364_8 - 5.249_603_131_819_84 * z) * z)))
+ x * (-7_296.439_871_453_82
+ x * (8_103.204_624_147_88
+ y * (2_175.341_332_000_392
+ y * (-274.229_003_681_796_4
+ y * (197.467_077_942_501_6
+ y * (-68.559_030_967_915_2 + 9.987_880_382_780_32 * y)))
- 90.673_423_405_131_6 * z)
+ x * (-5_458.342_052_148_35 - 980.141_533_448_88 * y
+ x * (2_247.607_427_267_04 - 340.123_748_317_786_3 * x
+ 220.542_973_797_483 * y)
+ 180.142_097_805_543 * z)
+ z * (-219.167_653_413_154_8
+ (-16.327_759_156_490_44 - 120.702_044_788_464_4 * z) * z))
+ z * (598.378_809_221_703
+ z * (-156.882_272_784_400_5
+ (204.133_482_817_937_7 - 10.237_557_973_238_46 * z) * z))
+ y * (-1_480.222_530_425_046
+ z * (-525.876_123_559_641
+ (249.577_178_340_545_71 - 88.449_193_048_287 * z) * z)
+ y * (-129.199_402_793_412_6
+ z * (1_149.174_198_007_428
+ z * (-162.575_178_755_133_6 + 76.919_546_216_974_2 * z))
+ y * (-30.068_211_258_562_5 - 1_380.959_795_403_770_8 * z
+ y * (2.626_801_985_426_835 + 703.695_562_834_065 * z)))))
+ y * (1_187.371_551_569_795_9
+ z * (1_458.233_059_470_092
+ z * (-687.913_805_923_122
+ z * (249.375_342_232_496
+ z * (-63.313_928_772_146 + 14.093_176_066_308_98 * z))))
+ y * (1_760.062_705_994_408
+ y * (-450.535_298_526_802
+ y * (182.852_089_550_251_8
+ y * (-43.320_648_175_062_2 + 4.260_339_416_943_66 * y)
+ z * (-595.457_483_974_374
+ (149.452_282_277_512 - 72.974_583_800_317_6 * z) * z))
+ z * (1_388.489_628_266_536
+ z * (-409.779_283_929_806
+ (227.123_395_681_188 - 22.256_546_865_282_6 * z) * z)))
+ z * (-1_721.528_607_567_954
+ z * (674.819_060_538_734
+ z * (-356.629_112_415_276
+ (88.408_071_661_6 - 15.840_030_944_233_64 * z) * z)))));
if x > 0.0 {
g08 += (11_625.629_132_534_64 + 1_702.453_469_893_412 * y) * libm::log(x);
}
if x == 0.0 {
g08 = f64::NAN;
}
Ok(0.5 * GSW_SFAC * g08)
} else if (ns == 0) & (nt == 1) & (np == 0) {
let g03 = 5.905_783_479_094_02
+ z * (-270.983_805_184_062
+ z * (776.153_611_613_101
+ z * (-196.512_550_881_22
+ (28.979_652_629_417_5 - 2.132_900_835_183_27 * z) * z)))
+ y * (-24_715.571_866_078
+ z * (2_910.072_908_093_6
+ z * (-1_513.116_771_538_718
+ z * (546.959_324_647_056
+ z * (-111.120_812_763_443_6 + 8.688_413_438_343_94 * z))))
+ y * (2_210.223_612_454_836_3
+ z * (-2_017.523_349_435_21
+ z * (1_498.081_172_457_456
+ z * (-718.635_991_963_235_9
+ (146.403_755_578_161_6 - 4.989_213_186_267_150_5 * z) * z)))
+ y * (-592.743_745_734_632
+ z * (1_591.873_781_627_888
+ z * (-1_207.261_522_487_504
+ (608.785_486_935_364 - 105.499_350_893_120_8 * z) * z))
+ y * (290.129_562_921_285_47
+ z * (-973.091_553_087_975
+ z * (602.603_274_510_125
+ z * (-276.361_526_170_076 + 32.409_533_403_861_05 * z)))
+ y * (-113.906_307_908_503_21
+ y * (21.355_715_254_157_69 - 67.417_568_357_514_34 * z)
+ z * (381.068_361_985_070_96
+ z * (-133.738_390_284_275_4
+ 49.023_632_509_086_724 * z)))))));
let mut g08 = x2
* (168.072_408_311_545
+ z * (729.116_529_735_046
+ z * (-343.956_902_961_561
+ z * (124.687_671_116_248
+ z * (-31.656_964_386_073 + 7.046_588_033_154_49 * z))))
+ x * (-493.407_510_141_682
+ x * (543.835_333_000_098
+ x * (-196.028_306_689_776 + 36.757_162_299_580_5 * x)
+ y * (-137.114_501_840_898_2
+ y * (148.100_308_456_876_18
+ y * (-68.559_030_967_915_2 + 12.484_850_478_475_4 * y)))
- 22.668_355_851_282_9 * z)
+ z * (-175.292_041_186_547
+ (83.192_392_780_181_9 - 29.483_064_349_429 * z) * z)
+ y * (-86.132_935_195_608_4
+ z * (766.116_132_004_952
+ z * (-108.383_452_503_422_4 + 51.279_697_477_982_8 * z))
+ y * (-30.068_211_258_562_5 - 1_380.959_795_403_770_8 * z
+ y * (3.502_402_647_235_78 + 938.260_750_445_42 * z))))
+ y * (1_760.062_705_994_408
+ y * (-675.802_947_790_203
+ y * (365.704_179_100_503_6
+ y * (-108.301_620_437_655_52 + 12.781_018_250_830_98 * y)
+ z * (-1_190.914_967_948_748
+ (298.904_564_555_024 - 145.949_167_600_635_2 * z) * z))
+ z * (2_082.734_442_399_804_3
+ z * (-614.668_925_894_709
+ (340.685_093_521_782 - 33.384_820_297_923_9 * z) * z)))
+ z * (-1_721.528_607_567_954
+ z * (674.819_060_538_734
+ z * (-356.629_112_415_276
+ (88.408_071_661_6 - 15.840_030_944_233_64 * z) * z)))));
if x > 0.0 {
g08 += 851.226_734_946_706 * x2 * libm::log(x);
}
Ok((g03 + g08) * 0.025)
} else if (ns == 0) & (nt == 0) & (np == 1) {
let g03 = 100_015.695_367_145
+ z * (-5_089.153_084_072_6
+ z * (853.553_335_338_861_1
+ z * (-133.258_701_701_444_4
+ (21.013_155_440_154_2 - 3.278_571_068_826_234 * z) * z)))
+ y * (-270.983_805_184_062
+ z * (1_552.307_223_226_202
+ z * (-589.537_652_643_66
+ (115.918_610_517_67 - 10.664_504_175_916_349 * z) * z))
+ y * (1_455.036_454_046_8
+ z * (-1_513.116_771_538_718
+ z * (820.438_986_970_584
+ z * (-222.241_625_526_887_2 + 21.721_033_595_859_85 * z)))
+ y * (-672.507_783_145_07
+ z * (998.720_781_638_304
+ z * (-718.635_991_963_235_9
+ (195.205_007_437_548_8 - 8.315_355_310_445_25 * z) * z))
+ y * (397.968_445_406_972
+ z * (-603.630_761_243_752
+ (456.589_115_201_523 - 105.499_350_893_120_8 * z) * z)
+ y * (-194.618_310_617_595
+ y * (63.511_393_664_178_5 - 9.631_081_193_930_62 * y
+ z * (-44.579_463_428_091_8
+ 24.511_816_254_543_362 * z))
+ z * (241.041_309_804_05
+ z * (-165.816_915_702_045_6
+ 25.927_626_723_088_84 * z)))))));
let g08 = x2
* (-3_310.491_540_448_39
+ z * (769.588_305_957_198
+ z * (-289.597_296_032_237_4
+ (63.363_269_106_729_6 - 13.124_007_829_549_6 * z) * z))
+ x * (199.459_603_073_901
+ x * (-54.791_913_353_288_7 + 36.028_419_561_108_6 * x
- 22.668_355_851_282_9 * y
+ (-8.163_879_578_245_22 - 90.526_533_591_348_31 * z) * z)
+ z * (-104.588_181_856_267
+ (204.133_482_817_937_7 - 13.650_077_297_651_28 * z) * z)
+ y * (-175.292_041_186_547
+ (166.384_785_560_363_8 - 88.449_193_048_287 * z) * z
+ y * (383.058_066_002_476
+ y * (-460.319_931_801_257 + 234.565_187_611_355 * y)
+ z * (-108.383_452_503_422_4 + 76.919_546_216_974_2 * z))))
+ y * (729.116_529_735_046
+ z * (-687.913_805_923_122
+ z * (374.063_013_348_744
+ z * (-126.627_857_544_292 + 35.232_940_165_772_45 * z)))
+ y * (-860.764_303_783_977
+ y * (694.244_814_133_268
+ y * (-297.728_741_987_187
+ (149.452_282_277_512 - 109.461_875_700_476_41 * z) * z)
+ z * (-409.779_283_929_806
+ (340.685_093_521_782 - 44.513_093_730_565_2 * z) * z))
+ z * (674.819_060_538_734
+ z * (-534.943_668_622_914
+ (176.816_143_323_2 - 39.600_077_360_584_095 * z) * z)))));
Ok((g03 + g08) * 1e-8)
} else if (ns == 1) & (nt == 1) & (np == 0) {
let mut g08 = 1_187.371_551_569_795_9
+ z * (1_458.233_059_470_092
+ z * (-687.913_805_923_122
+ z * (249.375_342_232_496
+ z * (-63.313_928_772_146 + 14.093_176_066_308_98 * z))))
+ x * (-1_480.222_530_425_046
+ x * (2_175.341_332_000_392
+ x * (-980.141_533_448_88 + 220.542_973_797_483 * x)
+ y * (-548.458_007_363_592_9
+ y * (592.401_233_827_504_7
+ y * (-274.236_123_871_660_8 + 49.939_401_913_901_6 * y)))
- 90.673_423_405_131_6 * z)
+ z * (-525.876_123_559_641
+ (249.577_178_340_545_71 - 88.449_193_048_287 * z) * z)
+ y * (-258.398_805_586_825_2
+ z * (2_298.348_396_014_856
+ z * (-325.150_357_510_267_2 + 153.839_092_433_948_4 * z))
+ y * (-90.204_633_775_687_5 - 4_142.879_386_211_312_5 * z
+ y * (10.507_207_941_707_34 + 2_814.782_251_336_26 * z))))
+ y * (3_520.125_411_988_816
+ y * (-1_351.605_895_580_406
+ y * (731.408_358_201_007_2
+ y * (-216.603_240_875_311_03 + 25.562_036_501_661_96 * y)
+ z * (-2_381.829_935_897_496
+ (597.809_129_110_048 - 291.898_335_201_270_4 * z) * z))
+ z * (4_165.468_884_799_608_5
+ z * (-1_229.337_851_789_418
+ (681.370_187_043_564 - 66.769_640_595_847_8 * z) * z)))
+ z * (-3_443.057_215_135_908
+ z * (1_349.638_121_077_468
+ z * (-713.258_224_830_552
+ (176.816_143_323_2 - 31.680_061_888_467_28 * z) * z))));
if sa == 0.0 {
g08 = f64::NAN;
} else if x > 0.0 {
g08 += 1_702.453_469_893_412 * libm::log(x);
}
Ok(0.5 * GSW_SFAC * 0.025 * g08)
} else if (ns == 1) & (nt == 0) & (np == 1) {
let g08 = -6_620.983_080_896_78
+ z * (1_539.176_611_914_396
+ z * (-579.194_592_064_474_8
+ (126.726_538_213_459_2 - 26.248_015_659_099_2 * z) * z))
+ x * (598.378_809_221_703
+ x * (-219.167_653_413_154_8 + 180.142_097_805_543 * x
- 90.673_423_405_131_6 * y
+ (-32.655_518_312_980_88 - 362.106_134_365_393_25 * z) * z)
+ z * (-313.764_545_568_801
+ (612.400_448_453_813_2 - 40.950_231_892_953_84 * z) * z)
+ y * (-525.876_123_559_641
+ (499.154_356_681_091_43 - 265.347_579_144_861 * z) * z
+ y * (1_149.174_198_007_428
+ y * (-1_380.959_795_403_770_8 + 703.695_562_834_065 * y)
+ z * (-325.150_357_510_267_2 + 230.758_638_650_922_6 * z))))
+ y * (1_458.233_059_470_092
+ z * (-1_375.827_611_846_244
+ z * (748.126_026_697_488
+ z * (-253.255_715_088_584 + 70.465_880_331_544_9 * z)))
+ y * (-1_721.528_607_567_954
+ y * (1_388.489_628_266_536
+ y * (-595.457_483_974_374
+ (298.904_564_555_024 - 218.923_751_400_952_82 * z) * z)
+ z * (-819.558_567_859_612
+ (681.370_187_043_564 - 89.026_187_461_130_4 * z) * z))
+ z * (1_349.638_121_077_468
+ z * (-1_069.887_337_245_828
+ (353.632_286_646_4 - 79.200_154_721_168_19 * z) * z))));
Ok(g08 * GSW_SFAC * 0.5e-8)
} else if (ns == 0) & (nt == 1) & (np == 1) {
let g03 = -270.983_805_184_062
+ z * (1_552.307_223_226_202
+ z * (-589.537_652_643_66
+ (115.918_610_517_67 - 10.664_504_175_916_349 * z) * z))
+ y * (2_910.072_908_093_6
+ z * (-3_026.233_543_077_436
+ z * (1_640.877_973_941_168
+ z * (-444.483_251_053_774_4 + 43.442_067_191_719_7 * z)))
+ y * (-2_017.523_349_435_21
+ z * (2_996.162_344_914_912
+ z * (-2_155.907_975_889_708
+ (585.615_022_312_646_4 - 24.946_065_931_335_752 * z) * z))
+ y * (1_591.873_781_627_888
+ z * (-2_414.523_044_975_008
+ (1_826.356_460_806_092 - 421.997_403_572_483_2 * z) * z)
+ y * (-973.091_553_087_975
+ z * (1_205.206_549_020_25
+ z * (-829.084_578_510_228 + 129.638_133_615_444_2 * z))
+ y * (381.068_361_985_070_96 - 67.417_568_357_514_34 * y
+ z * (-267.476_780_568_550_8 + 147.070_897_527_260_17 * z))))));
let g08 = x2
* (729.116_529_735_046
+ z * (-687.913_805_923_122
+ z * (374.063_013_348_744
+ z * (-126.627_857_544_292 + 35.232_940_165_772_45 * z)))
+ x * (-175.292_041_186_547 - 22.668_355_851_282_9 * x
+ (166.384_785_560_363_8 - 88.449_193_048_287 * z) * z
+ y * (766.116_132_004_952
+ y * (-1_380.959_795_403_770_8 + 938.260_750_445_42 * y)
+ z * (-216.766_905_006_844_8 + 153.839_092_433_948_4 * z)))
+ y * (-1_721.528_607_567_954
+ y * (2_082.734_442_399_804_3
+ y * (-1_190.914_967_948_748
+ (597.809_129_110_048 - 437.847_502_801_905_65 * z) * z)
+ z * (-1_229.337_851_789_418
+ (1_022.055_280_565_346 - 133.539_281_191_695_6 * z) * z))
+ z * (1_349.638_121_077_468
+ z * (-1_069.887_337_245_828
+ (353.632_286_646_4 - 79.200_154_721_168_19 * z) * z))));
Ok((g03 + g08) * 2.5e-10)
} else if (ns == 2) & (nt == 0) & (np == 0) {
let mut g08 = 5_812.814_566_267_320
+ 851.226_734_946_706_0 * y
+ x * (-3_648.219_935_726_910
+ x * (8_103.204_624_147_880
+ x * (-8_187.513_078_222_526
+ x * (4_495.214_854_534_080 - 850.309_370_794_465_7 * x)))
+ y * (-740.111_265_212_523_0
+ x * (2_175.341_332_000_392
+ x * (-1_470.212_300_173_320 + 441.085_947_594_966_0 * x))
+ y * (-64.599_701_396_706_29 - 274.229_003_681_796_4 * x
+ y * (-15.034_105_629_281_25
+ 197.467_077_942_501_6 * x
+ y * (1.313_400_992_713_418 - 68.559_030_967_915_21 * x
+ 9.987_880_382_780_312 * x * y))))
+ z * (299.189_404_610_851_5
+ x * (-219.167_653_413_154_8 + 270.213_146_708_314_5 * x)
+ y * (-262.938_061_779_820_5 - 90.673_423_405_131_60 * x
+ y * (574.587_099_003_714_0
+ y * (-690.479_897_701_885_4 + 351.847_781_417_032_5 * y)))
+ z * (-78.441_136_392_200_25 - 16.327_759_156_490_44 * x
+ y * (124.788_589_170_272_9 - 81.287_589_377_566_80 * y)
+ z * (102.066_741_408_968_9 - 120.702_044_788_464_4 * x
+ y * (-44.224_596_524_143_50 + 38.459_773_108_487_10 * y)
- 5.118_778_986_619_230 * z))));
if x > 0.0 {
g08 /= x2;
}
if x == 0.0 {
g08 = f64::NAN;
}
Ok(0.5 * GSW_SFAC * GSW_SFAC * g08)
} else if (ns == 0) & (nt == 2) & (np == 0) {
let g03 = -24_715.571_866_078
+ z * (2_910.072_908_093_6
+ z * (-1_513.116_771_538_718
+ z * (546.959_324_647_056
+ z * (-111.120_812_763_443_6 + 8.688_413_438_343_94 * z))))
+ y * (4_420.447_224_909_672_5
+ z * (-4_035.046_698_870_42
+ z * (2_996.162_344_914_912
+ z * (-1_437.271_983_926_471_9
+ (292.807_511_156_323_2 - 9.978_426_372_534_301 * z) * z)))
+ y * (-1_778.231_237_203_896
+ z * (4_775.621_344_883_664
+ z * (-3_621.784_567_462_512
+ (1_826.356_460_806_092 - 316.498_052_679_362_44 * z) * z))
+ y * (1_160.518_251_685_141_9
+ z * (-3_892.366_212_351_9
+ z * (2_410.413_098_040_5
+ z * (-1_105.446_104_680_304 + 129.638_133_615_444_2 * z)))
+ y * (-569.531_539_542_516
+ y * (128.134_291_524_946_15 - 404.505_410_145_086_05 * z)
+ z * (1_905.341_809_925_355
+ z * (-668.691_951_421_377 + 245.118_162_545_433_62 * z))))));
let g08 = x2
* (1_760.062_705_994_408
+ x * (-86.132_935_195_608_4
+ x * (-137.114_501_840_898_2
+ y * (296.200_616_913_752_36
+ y * (-205.677_092_903_745_63 + 49.939_401_913_901_6 * y)))
+ z * (766.116_132_004_952
+ z * (-108.383_452_503_422_4 + 51.279_697_477_982_8 * z))
+ y * (-60.136_422_517_125 - 2_761.919_590_807_541_7 * z
+ y * (10.507_207_941_707_34 + 2_814.782_251_336_26 * z)))
+ y * (-1_351.605_895_580_406
+ y * (1_097.112_537_301_510_9
+ y * (-433.206_481_750_622_06 + 63.905_091_254_154_904 * y)
+ z * (-3_572.744_903_846_243_7
+ (896.713_693_665_072 - 437.847_502_801_905_65 * z) * z))
+ z * (4_165.468_884_799_608_5
+ z * (-1_229.337_851_789_418
+ (681.370_187_043_564 - 66.769_640_595_847_8 * z) * z)))
+ z * (-1_721.528_607_567_954
+ z * (674.819_060_538_734
+ z * (-356.629_112_415_276
+ (88.408_071_661_6 - 15.840_030_944_233_64 * z) * z))));
Ok((g03 + g08) * 0.000_625)
} else if (ns == 0) & (nt == 0) & (np == 2) {
let g03 = -5_089.153_084_072_6
+ z * (1_707.106_670_677_722_1
+ z * (-399.776_105_104_333_2
+ (84.052_621_760_616_8 - 16.392_855_344_131_17 * z) * z))
+ y * (1_552.307_223_226_202
+ z * (-1_179.075_305_287_32
+ (347.755_831_553_01 - 42.658_016_703_665_396 * z) * z)
+ y * (-1_513.116_771_538_718
+ z * (1_640.877_973_941_168
+ z * (-666.724_876_580_661_5 + 86.884_134_383_439_4 * z))
+ y * (998.720_781_638_304
+ z * (-1_437.271_983_926_471_9
+ (585.615_022_312_646_4 - 33.261_421_241_781 * z) * z)
+ y * (-603.630_761_243_752
+ (913.178_230_403_046 - 316.498_052_679_362_44 * z) * z
+ y * (241.041_309_804_05
+ y * (-44.579_463_428_091_8 + 49.023_632_509_086_724 * z)
+ z * (-331.633_831_404_091_2 + 77.782_880_169_266_52 * z))))));
let g08 = x2
* (769.588_305_957_198
+ z * (-579.194_592_064_474_8
+ (190.089_807_320_188_78 - 52.496_031_318_198_4 * z) * z)
+ x * (-104.588_181_856_267
+ x * (-8.163_879_578_245_22 - 181.053_067_182_696_62 * z)
+ (408.266_965_635_875_4 - 40.950_231_892_953_84 * z) * z
+ y * (166.384_785_560_363_8 - 176.898_386_096_574 * z
+ y * (-108.383_452_503_422_4 + 153.839_092_433_948_4 * z)))
+ y * (-687.913_805_923_122
+ z * (748.126_026_697_488
+ z * (-379.883_572_632_876 + 140.931_760_663_089_8 * z))
+ y * (674.819_060_538_734
+ z * (-1_069.887_337_245_828
+ (530.448_429_969_6 - 158.400_309_442_336_38 * z) * z)
+ y * (-409.779_283_929_806
+ y * (149.452_282_277_512 - 218.923_751_400_952_82 * z)
+ (681.370_187_043_564 - 133.539_281_191_695_6 * z) * z))));
Ok((g03 + g08) * 1e-16)
} else if (ns == 0) & (nt == 3) & (np == 0) {
let g03 = 4_420.447_224_909_672_51
+ y * (-3_556.462_474_407_791_89
+ y * (3_481.554_755_055_425_63
+ y * (-2_278.126_158_170_064_17 + 640.671_457_624_730_749 * y)))
+ z * (-4_035.046_698_870_420_19
+ y * (9_551.242_689_767_328_56
+ y * (-11_677.098_637_055_699_3
+ y * (7_621.367_239_701_419_75 - 2_022.527_050_725_430_23 * y)))
+ z * (2_996.162_344_914_911_96
+ y * (-7_243.569_134_925_023_72
+ y * (7_231.239_294_121_499_37 - 2_674.767_805_685_507_94 * y))
+ z * (-1_437.271_983_926_471_88
+ y * (3_652.712_921_612_183_89
+ y * (-3_316.338_314_040_911_88 + 980.472_650_181_734_480 * y))
+ z * (292.807_511_156_323_187
+ y * (-632.996_105_358_724_890 + 388.914_400_846_332_569 * y)
- 9.978_426_372_534_300_98 * z))));
let g08 = x2
* (-1_351.605_895_580_406
+ x * (-60.136_422_517_125_0 + 296.200_616_913_752 * x)
+ y * (2_194.225_074_603_02
+ y * (-1_299.619_445_251_87 + 255.620_365_016_620 * y)
+ x * (21.014_415_883_414_7
+ x * (-411.354_185_807_491 + 149.818_205_741_705 * y)))
+ z * (4_165.468_884_799_61 - 2_761.919_590_807_54 * x
+ y * (-7_145.489_807_692_49 + 5_629.564_502_672_52 * x)
+ z * (-1_229.337_851_789_42
+ 1_793.427_387_330_14 * y
+ z * (681.370_187_043_564
- 875.695_005_603_811 * y
- 66.769_640_595_847_8 * z))));
Ok((g03 + g08) * 0.000_625 * 0.025)
} else if (ns == 0) & (nt == 2) & (np == 1) {
let g03 = 2_910.072_908_093_60
+ y * (-4_035.046_698_870_42
+ y * (4_775.621_344_883_66
+ y * (-3_892.366_212_351_90
+ y * (1_905.341_809_925_35 - 404.505_410_145_086 * y))))
+ z * (-3_026.233_543_077_44
+ y * (5_992.324_689_829_82
+ y * (-7_243.569_134_925_02
+ y * (4_820.826_196_081_00 - 1_337.383_902_842_75 * y)))
+ z * (1_640.877_973_941_17
+ y * (-4_311.815_951_779_42
+ y * (5_479.069_382_418_28
+ y * (-3_316.338_314_040_91 + 735.354_487_636_301 * y)))
+ z * (-444.483_251_053_774
+ y * (1_171.230_044_625_29
+ y * (-1_265.992_210_717_45 + 518.552_534_461_777 * y))
+ z * (43.442_067_191_719_7 - 49.892_131_862_671_5 * y))));
let g08 = x2
* (-1_721.528_607_567_95 + 766.116_132_004_952 * x - 2_761.919_590_807_54 * x * y
+ y * (4_165.468_884_799_61
+ y * (-3_572.744_903_846_24 + 2_814.782_251_336_26 * x))
+ z * (1_349.638_121_077_47 - 216.766_905_006_845 * x
+ y * (-2_458.675_703_578_84 + 1_793.427_387_330_14 * y)
+ z * (-1_069.887_337_245_83
+ 153.839_092_433_949 * x
+ y * (2_044.110_561_130_69 - 1_313.542_508_405_72 * y)
+ z * (353.632_286_646_400
- 267.078_562_383_391 * y
- 79.200_154_721_168_2 * z))));
Ok((g03 + g08) * 2.5e-10 * 0.025)
} else if (ns == 1) & (nt == 1) & (np == 1) {
let g08 = 1_458.233_059_470_09
+ x * (-525.876_123_559_641 - 90.673_423_405_131_6 * x
+ y * (2_298.348_396_014_86
+ y * (-4_142.879_386_211_31 + 2_814.782_251_336_26 * y)))
+ y * (-3_443.057_215_135_91 + y * (4_165.468_884_799_61 - 2_381.829_935_897_50 * y))
+ z * (-1_375.827_611_846_24
+ x * (499.154_356_681_091 - 650.300_715_020_534 * y)
+ y * (2_699.276_242_154_94
+ y * (-2_458.675_703_578_84 + 1_195.618_258_220_10 * y))
+ z * (748.126_026_697_488
+ x * (-265.347_579_144_861 + 461.517_277_301_845 * y)
+ y * (-2_139.774_674_491_66
+ y * (2_044.110_561_130_69 - 875.695_005_603_811 * y))
+ z * (-253.255_715_088_584
+ y * (707.264_573_292_800 - 267.078_562_383_391 * y)
+ z * (70.465_880_331_544_9 - 158.400_309_442_336 * y))));
Ok(g08 * GSW_SFAC * 0.5e-8 * 0.025)
} else if (ns == 2) & (nt == 0) & (np == 1) {
let mut g08 = 299.189_404_610_851_5
+ x * (-219.167_653_413_154_8 + 270.213_146_708_314_5 * x)
+ y * (-262.938_061_779_820_5 - 90.673_423_405_131_60 * x
+ y * (574.587_099_003_714_0
+ y * (-690.479_897_701_885_4 + 351.847_781_417_032_5 * y)))
+ z * (-156.882_272_784_400_5 - 32.655_518_312_980_88 * x
+ y * (249.577_178_340_545_7 - 162.575_178_755_133_6 * y)
+ z * (306.200_224_226_906_6 - 362.106_134_365_393_2 * x
+ y * (-132.673_789_572_430_5 + 115.379_319_325_461_3 * y)
- 20.475_115_946_476_92 * z));
if x > 0.0 {
g08 /= x;
}
if x == 0.0 {
g08 = f64::NAN;
}
Ok(0.5 * GSW_SFAC * GSW_SFAC * g08 * 1e-8)
} else if (ns == 1) & (nt == 2) & (np == 0) {
let mut g08 = 3_520.125_411_988_816
+ x * (-258.398_805_586_825_2 - 548.458_007_363_592_9 * x)
+ y * (-2_703.211_791_160_812
+ x * (-180.409_267_551_375 + 1_184.802_467_655_009_4 * x)
+ y * (2_194.225_074_603_021_7
+ x * (31.521_623_825_122_02
+ x * (-822.708_371_614_982_47 + 199.757_607_655_606_4 * y))
+ y * (-866.412_963_501_244_1 + 127.810_182_508_309_8 * y)))
+ z * (-3_443.057_215_135_908
+ 2_298.348_396_014_856 * x
+ y * (8_330.937_769_599_217 - 8_285.758_772_422_625 * x
+ y * (-7_145.489_807_692_487_8 + 8_444.346_754_008_78 * x))
+ z * (1_349.638_121_077_468 - 325.150_357_510_267_2 * x
+ y * (-2_458.675_703_578_836 + 1_793.427_387_330_144 * y)
+ z * (-713.258_224_830_552
+ 153.839_092_433_948_4 * x
+ y * (1_362.740_374_087_128 - 875.695_005_603_811_2 * y)
+ z * (176.816_143_323_2
- 133.539_281_191_695_6 * y
- 31.680_061_888_467_3 * z))));
if x == 0.0 {
g08 = f64::NAN;
}
Ok(0.5 * GSW_SFAC * 0.025 * g08 * 0.025)
} else if (ns == 2) & (nt == 1) & (np == 0) {
let mut g08 = 851.226_734_946_705_96
+ x * (-740.111_265_212_522_99
+ x * (2_175.341_332_000_392_0
+ x * (-1_470.212_300_173_319_9 + 441.085_947_594_966_00 * x))
+ y * (-129.199_402_793_412_59 - 548.458_007_363_592_86 * x
+ y * (-45.102_316_887_843_749
+ 592.401_233_827_504_77 * x
+ y * (5.253_603_970_853_670_4
+ x * (-274.236_123_871_660_82 + 49.939_401_913_901_600 * y))))
+ z * (-262.938_061_779_820_49 - 90.673_423_405_131_601 * x
+ y * (1_149.174_198_007_428_0
+ y * (-2_071.439_693_105_656_3 + 1_407.391_125_668_129_9 * y))
+ z * (124.788_589_170_272_86 - 162.575_178_755_133_61 * y
+ z * (-44.224_596_524_143_500 + 76.919_546_216_974_197 * y))));
if x > 0.0 {
g08 /= x2;
}
if x == 0.0 {
g08 = f64::NAN;
}
Ok(0.5 * GSW_SFAC * GSW_SFAC * g08 * 0.025)
} else if (ns == 1) & (nt == 0) & (np == 2) {
let g08 = 1_539.176_611_914_396_06
+ x * (-313.764_545_568_801 - 32.655_518_312_980_9 * x)
+ y * (-1_375.827_611_846_244
+ 499.154_356_681_091_43 * x
+ y * (1_349.638_121_077_468 - 325.150_357_510_267 * x
+ y * (-819.558_567_859_612 + 298.904_564_555_024 * y)))
+ z * (-1_158.389_184_128_949_6
+ x * (1_224.800_896_907_626_3 - 724.212_268_730_786_5 * x)
+ y * (1_496.252_053_394_976 - 530.695_158_289_722 * x
+ y * (-2_139.774_674_491_656
+ 461.517_277_301_845_2 * x
+ y * (1_362.740_374_087_128 - 437.847_502_801_906 * y)))
+ z * (380.179_614_640_378 - 122.850_695_678_862 * x
+ y * (-759.767_145_265_752
+ y * (1_060.896_859_939_2 - 267.078_562_383_391 * y))
+ z * (-104.992_062_636_397
+ y * (281.863_521_326_18 - 316.800_618_884_673 * y))));
Ok(g08 * GSW_SFAC * 0.5e-8 * 1e-8)
} else if (ns == 0) & (nt == 1) & (np == 2) {
let g03 = 1_552.307_223_226_202
+ y * (-3_026.233_543_077_436
+ y * (2_996.162_344_914_912
+ y * (-2_414.523_044_975_008
+ y * (1_205.206_549_020_25 - 267.476_780_568_550_8 * y))))
+ z * (-1_179.075_305_287_32
+ y * (3_281.755_947_882_336
+ y * (-4_311.815_951_779_416
+ y * (3_652.712_921_612_184
+ y * (-1_658.169_157_020_456 + 294.141_795_054_520 * y))))
+ z * (347.755_831_553_010
+ y * (-1_333.449_753_161_323
+ y * (1_756.845_066_937_94
+ y * (-1_265.992_210_717_45 + 388.914_400_846_332_6 * y)))
+ z * (-42.658_016_703_665_4
+ y * (173.768_268_766_878_8 - 99.784_263_725_343_0 * y))));
let g08 = x2
* (-687.913_805_923_122_0
+ 166.384_785_560_363_8 * x
+ y * (1_349.638_121_077_468 - 216.766_905_006_844_8 * x
+ y * (-1_229.337_851_789_418 + 597.809_129_110_048_0 * y))
+ z * (748.126_026_697_488_0 - 176.898_386_096_574_0 * x
+ y * (-2_139.774_674_491_656
+ 307.678_184_867_896_8 * x
+ y * (2_044.110_561_130_692 - 875.695_005_603_811_3 * y))
+ z * (-379.883_572_632_876
+ y * (1_060.896_859_939_20 - 400.617_843_575_086_8 * y)
+ z * (140.931_760_663_089_8 - 316.800_618_884_672_8 * y))));
Ok((g03 + g08) * 2.5e-10 * 1e-8)
} else {
Err(Error::Undefined)
}
}
#[cfg(test)]
mod test_gibbs {
use super::gibbs;
#[test]
fn nan() {
for ns in 0..=1 {
for nt in 0..=1 {
for np in 0..=1 {
let v = gibbs(ns, nt, np, f64::NAN, 1.0, 1.0);
assert!(v.unwrap().is_nan());
let v = gibbs(ns, nt, np, 1.0, f64::NAN, 1.0);
assert!(v.unwrap().is_nan());
let v = gibbs(ns, nt, np, 1.0, 1.0, f64::NAN);
assert!(v.unwrap().is_nan());
}
}
}
}
#[test]
fn checking_values() {
let sa = 32.0;
let t = 10.0;
let p = 100.0;
let check_values = [
(0, 0, 0, 50.249_433_022_084_645),
(0, 0, 1, 9.756_573_803_380_614e-4),
(0, 0, 2, -4.319_267_169_958_817e-13),
(0, 1, 0, -144.761_662_219_825_35),
(0, 1, 1, 1.579_419_118_574_033e-7),
(0, 1, 2, 1.731_843_128_190_281e-15),
(0, 2, 0, -14.141_773_253_302_12),
(0, 2, 1, 9.947_007_722_211_213e-9),
(0, 3, 0, 4.830_012_068_790_147e-2),
(1, 0, 0, 60.388_852_305_031_49),
(1, 0, 1, -7.383_641_712_245_458e-7),
(1, 0, 2, 1.306_882_540_566_367e-15),
(1, 1, 0, 0.471_143_208_002_190_2),
(1, 1, 1, 1.849_778_736_567_631e-9),
(1, 2, 0, 1.904_117_940_743_395e-2),
(2, 0, 0, 2.269_679_195_264_845),
(2, 0, 1, 8.949_706_277_002_558e-10),
(2, 1, 0, 1.013_083_889_176_643e-2),
];
for (ns, nt, np, ans) in check_values.iter() {
assert!((gibbs(*ns, *nt, *np, sa, t, p).unwrap() - *ans).abs() < f64::EPSILON);
}
}
}