gsw/gsw_internal_funcs/

mod.rs

//! Internal Functions
//!
//! Functions not intended to be used outside this library

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]
/// Non-dimensional pressure
///
/// The polynomial approximation solutions proposed by Roquet (2015) are based
/// on non-dimensional salinity, temperature, and pressure. Here we scale
/// pressure by p_u (1e4 [dbar]) to obtain the non-dimensional quantity \pi
/// (TEOS-10 Manual appendix K, or \zeta on Roquet (2015)).
///
/// # Argument
///
/// * `p`: sea pressure \[dbar\] (i.e. absolute pressure - 10.1325 dbar)
///
/// # Returns
///
/// * `\pi`: Non-dimensional pressure
///
/// # Notes
///
/// * The original formulation is a scaling of p by p_u. The MatLab and C
///   implementations of GSW operate as the product with 1e-4, which does make
///   sense since it is a lighter operation than a division is for computers.
///   The issue here is on the inability of f64 to fully represent certain
///   fractions. For instance, while 3812.0 can be perfectly represented,
///   0.3812 is rounded to 0.381200000000000038813. Similarly 1e4 is fine but
///   1e-4 on f64 is rounded to 0.000100000000000000004792, thus 3812/1e4 is
///   diffrent than 3812*1e-4.
///
/// # Example
/// ```
/// let p = 3812;
/// assert!((3_812.0 * 1e-4) > (3_812.0 / 1e4))
/// ```
pub(crate) fn non_dimensional_p(p: f64) -> f64 {
    if cfg!(feature = "compat") {
        p * 1e-4
    } else {
        p / GSW_PU
    }
}

/// Specific Volume of Standard Ocean Salinity and CT=0
///
/// This function calculates specific volume at the Standard Ocean Salinity,
/// SSO, and at a Conservative Temperature of zero degrees C, as a function
/// of pressure, p, in dbar, using a streamlined version of the 75-term CT
/// version of specific volume, that is, a streamlined version of the code
/// "specvol(SA,CT,p)".
///
/// version: 3.06.12
///
/// If using compat (truncated constants) there is a difference of O[1e-19],
/// which is negligible but enough to fail the validation tests.
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;

    /// specvol() at SSO & CT=0 should be identical to specvol_sso_0()
    #[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::conversions::t68_from_t90(t);
    let ft68: f64 = (t68 - 15.0) / (1.0 + K * (t68 - 15.0));

    // Find the initial estimates of Rtx (Rtx0) and of the derivative dSP_dRtx
    // at SP = 2.
    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);

    //  Begin a single modified Newton-Raphson iteration (McDougall and
    //  Wotherspoon, 2013) to find Rt at SP = 2.

    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;

    // This is the end of one full iteration of the modified Newton-Raphson
    // iterative equation solver.  The error in Rtx at this point is
    // equivalent to an error in SP of 9e-16 psu.

    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]
    // These test values were obtained by running GSW-Matlab. It was common
    // a difference of O(16), so we assumed that the values obtained with
    // Rust were correct.
    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)]
/// Gibbs energy and its derivatives
///
/// # Notes
/// No effort yet on optimizing, but just reproducing the results.
pub(crate) fn gibbs(ns: u8, nt: u8, np: u8, sa: f64, t: f64, p: f64) -> Result<f64> {
    // Temporary solution to avoid failure with powerpc64
    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;
    // Note.The input pressure (p) is sea pressure in units of dbar.
    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)
    //% Note. This pressure derivative of the gibbs function is in units of (J/kg) (Pa^-1) = m^3/kg
    } 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)
    // Note. This derivative of the Gibbs function is in units of (m^3/kg)/(g/kg) = m^3/g,
    // that is, it is the derivative of specific volume with respect to Absolute
    // Salinity measured in g/kg.
    } 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)
    // Note. This derivative of the Gibbs function is in units of (m^3/(K kg)),
    // that is, the pressure of the derivative in Pa.
    } 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)
    // Note. This is the second derivative of the Gibbs function with respect to
    //pressure, measured in Pa.  This derivative has units of (J/kg) (Pa^-2).
    } 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)
        // Note. This derivative of the Gibbs function is in units of (m^3/(K kg)), that is, the pressure of the derivative in Pa.
    } else {
        Err(Error::Undefined)
    }
}

#[cfg(test)]
mod test_gibbs {
    use super::gibbs;

    #[test]
    // NaN input results in NaN output.
    // Other libraries using GSW-rs might rely on this behavior to propagate
    // and handle invalid elements.
    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]
    /// A sanity check with a single set of (S,T,P).
    /// Check values match Matlab's output.
    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);
        }
    }
}

/// Ice specific Gibbs energy and derivatives up to order 2.
///
/// # Arguments
/// - `nt`: `u8` order of t derivative \[integers 0, 1, or 2\]
/// - `np`: `u8` order of p derivative \[integers 0, 1, or 2\]
/// - `t`: `f64` in-situ temoerature (ITS-90) \[deg C\]
/// - `p`: `f64` sea pressure \[dbar\]
///
/// # Returns
/// `gibbs_ice` : `Result<Ok(f64), Err(Error::OutOfBounds)>`
///     - specific gibbs energy of ice or its derivatives
///     - gibbs energy in \[J/kg\]
///     - temperature derivatives in: \[(J/kg) (K)^(-nt)\]
///     - pressure derivatives in: \[(J/kg) (Pa)^(-np)\]
///     - mixed derivatives in \[(J/kg) (K)^(-nt) (Pa)^(-np)\]
#[allow(dead_code)]
pub(crate) fn gibbs_ice(nt: u8, np: u8, t: f64, p: f64) -> Result<f64> {
    // use a complex number crate. eventually replaced by calculations by hand
    use num_complex::Complex;

    const GSW_T0: f64 = 273.15;
    const REC_TT: f64 = 3.660_858_105_139_845_e-3;
    const REC_PT: f64 = 1.634_903_221_903_779_e-3;
    const T1: Complex<f64> = Complex::new(3.680_171_128_550_51_e-2, 5.108_781_149_595_72_e-2);
    const T2: Complex<f64> = Complex::new(3.373_157_410_654_16_e-1, 3.354_494_159_193_09_e-1);
    const G00: f64 = -6.320_202_333_358_86_e5;
    const G01: f64 = 6.550_222_136_589_55_e-1;
    const G02: f64 = -1.893_699_293_261_31_e-8;
    const G03: f64 = 3.397_461_232_710_53_e-15; // in the GSW-C library this is 3.39_746_123_271_053_04e-15
    const G04: f64 = -5.564_648_690_589_91_e-22; // in the GSW-C library this is -5.56_464_869_058_990_9e-22
    const R21: Complex<f64> = Complex::new(-5.571_076_980_301_23_e-5, 4.645_786_345_808_06_e-5);
    const R22: Complex<f64> = Complex::new(2.348_014_092_159_13_e-11, -2.856_511_429_049_72_e-11);
    const TT: f64 = 273.16;
    // R20 already named
    const R20_GIBBS_ICE: Complex<f64> =
        Complex::new(-7.259_745_743_292_20_e1, -7.810_084_271_128_70_e1);
    // R1 already named
    const R1_GIBBS_ICE: Complex<f64> =
        Complex::new(4.470_507_162_853_88_e1, 6.568_768_474_634_81_e1);

    // declare these initial variables:
    let s0: f64 = -3.327_337_564_921_68_e3;
    let tau = Complex::new((t + GSW_T0) * REC_TT, 0.0);
    let dzi: f64 = DB2PA * p * REC_PT;

    // pulling these constants out of the if statements to not be redundant
    let tau_t1 = tau / T1;
    let sqtau_t1 = tau_t1 * tau_t1;
    let tau_t2 = tau / T2;
    let sqtau_t2 = tau_t2 * tau_t2;
    let r2 = R20_GIBBS_ICE + dzi * (R21 + R22 * dzi);
    let r2p = REC_PT * (R21 + 2.0 * R22 * dzi);

    if nt == 0 && np == 0 {
        let g0: f64 = G00 + dzi * (G01 + dzi * (G02 + dzi * (G03 + dzi * G04)));

        let g = R1_GIBBS_ICE
            * (tau * ((1.0 + tau_t1) / (1.0 - tau_t1)).ln()
                + T1 * ((1.0 - sqtau_t1).ln() - sqtau_t1))
            + r2 * ((tau * ((1.0 + tau_t2) / (1.0 - tau_t2)).ln())
                + T2 * ((1.0 - sqtau_t2).ln() - sqtau_t2));

        Ok((g0 - TT * (s0 * tau - g.re)).re)
    } else if nt == 1 && np == 0 {
        let g = R1_GIBBS_ICE * (((1.0 + tau_t1) / (1.0 - tau_t1)).ln() - 2.0 * tau_t1)
            + r2 * (((1.0 + tau_t2) / (1.0 - tau_t2)).ln() - 2.0 * tau_t2);

        Ok(-s0 + g.re)
    } else if nt == 0 && np == 1 {
        let g0p = REC_PT * (G01 + dzi * (2.0 * G02 + dzi * (3.0 * G03 + 4.0 * G04 * dzi)));

        let g = r2p
            * (tau * ((1.0 + tau_t2) / (1.0 - tau_t2)).ln()
                + T2 * ((1.0 - sqtau_t2).ln() - sqtau_t2));

        Ok(g0p + TT * g.re)
    } else if nt == 1 && np == 1 {
        let g = r2p * (((1.0 + tau_t2) / (1.0 - tau_t2)).ln() - 2.0 * tau_t2);

        Ok(g.re)
    } else if nt == 2 && np == 0 {
        // start to rewrite without num-complex

        let g = R1_GIBBS_ICE * (1.0 / (T1 - tau) + 1.0 / (T1 + tau) - 2.0 / T1)
            + r2 * (1.0 / (T2 - tau));

        Ok(REC_TT * g.re)
    } else if nt == 0 && np == 2 {
        let sqrec_pt = REC_PT * REC_PT;

        let g0pp = sqrec_pt * (2.0 * G02 + dzi * (6.0 * G03 + 12.0 * G04 * dzi));

        let r2pp = 2.0 * R22 * sqrec_pt;

        let g = r2pp
            * (tau * ((1.0 + tau_t2) / (1.0 - tau_t2)).ln()
                + T2 * ((1.0 - sqtau_t2).ln() - sqtau_t2));

        Ok(g0pp + TT * g.re)
    } else {
        Err(Error::OutOfBounds)
    }
}

#[cfg(test)]
mod test_gibbs_ice {

    use crate::gsw_internal_funcs::gibbs_ice;
    use crate::Error;

    #[test]
    fn check_values() {
        // answers have not been checked for correctness, but are outputs of the function using num-complex
        // (nt, np, ans)
        let answers = [
            (0, 0, 98.267_598_402_919_25),
            (1, 0, 1_220.788_661_299_952_8),
            (0, 1, 0.001_090_843_442_926_435),
            (1, 1, 1.743_608_249_608_496_2E-7),
            (2, 0, -9.078_428_407_946_305),
            (0, 2, -1.284_848_246_397_617_7E-13),
        ];
        for (nt, np, ans) in answers.iter() {
            assert!(
                (gibbs_ice(*nt, *np, 0.0, 0.0).unwrap() - *ans).abs() < f64::EPSILON,
                "nt = {nt}, np = {np}, ans = {ans}"
            );
        }
        //assert!((gibbs_ice(0, 0, 4.0, 10.0).unwrap() - 5.029179813607595e3).abs() < f64::EPSILON); // this is super close
    }

    #[test]
    fn out_of_bounds() {
        match gibbs_ice(0, 3, 0.0, 0.0) {
            Err(Error::OutOfBounds) => (),
            _ => panic!("error should be OutOfBounds"),
        }
    }

    #[test]
    // NaN input results in NaN output.
    // Other libraries using GSW-rs might rely on this behavior to propagate
    // and handle invalid elements.
    fn nan() {
        let v = gibbs_ice(0, 0, 0.0, f64::NAN);
        assert!(v.unwrap().is_nan());

        let v = gibbs_ice(0, 0, f64::NAN, 0.0);
        assert!(v.unwrap().is_nan());

        let v = gibbs_ice(0, 0, f64::NAN, f64::NAN);
        assert!(v.unwrap().is_nan());
    }
}

/*
gsw_SAAR
gsw_Fdelta
gsw_deltaSA_atlas
gsw_SA_from_SP_Baltic
gsw_SP_from_SA_Baltic
gsw_infunnel
gsw_entropy_part
gsw_entropy_part_zerop
gsw_quadprog
gsw_wiggliness
gsw_data_interp
gsw_interp_ref_cast
gsw_linear_interp_SA_CT
gsw_pchip_interp_SA_CT
gsw_rr68_interp_SA_CT
gsw_spline_interp_SA_CT
gsw_gibbs_pt0_pt0
gsw_gibbs_ice_part_t
gsw_gibbs_ice_pt0
*/