ct_from_pt

src/conversions.rs

@@ -3,21 +3,154 @@
 //! Other conversions between temperatures, salinities, entropy, pressure
 //! and height.
 
-use crate::gsw_internal_const::{DB2PA, DEG2RAD, GAMMA, GSW_P0};
+use crate::gsw_internal_const::{DB2PA, DEG2RAD, GAMMA, GSW_CP0, GSW_P0, GSW_SFAC, GSW_UPS};
 use crate::gsw_internal_funcs::enthalpy_sso_0;
+use crate::{Error, Result};
+
+/*
+/// Absolute Salinity Anomaly from Practical Salinity
+///
+pub fn deltasa_from_sp(sp: f64, p: f64, lon: f64, lat: f64) -> Result<f64> {
+    // Remove out of range values
+    if ((p < 100.0) & (sp > 120.0)) | ((p >= 100.0) & (sp > 42.0)) {
+        if cfg!(feature = "invalidasnan") {
+            return Ok(f64::NAN);
+        } else {
+            return Err(Error::Undefined);
+        }
+    }
+
+    if (p < -1.5) | (p > 12000.0) | (lon < 0.0) | (lon > 360.0) | (lat < -90.0) | (lat > 90.0) {
+        if cfg!(feature = "compat") {
+            return Err(Error::Undefined);
+        } else if cfg!(feature = "invalidasnan") {
+            return Ok(f64::NAN);
+        } else {
+            return Err(Error::Undefined);
+        }
+    }
+
+    let sp: f64 = if sp < 0.0 {
+        if cfg!(feature = "compat") {
+            0.0
+        } else if cfg!(feature = "invalidasnan") {
+            return Ok(f64::NAN);
+        } else {
+            return Err(Error::NegativeSalinity);
+        }
+    } else {
+        sp
+    };
+
+    let sa = sa_from_sp(sp, p, lon, lat);
+    let sr = sr_from_sp(sp);
+    let dsa = sa - sr;
+    Ok(dsa)
+}
+*/
 
 /*
-gsw_deltaSA_from_SP
 gsw_SA_Sstar_from_SP
-gsw_SR_from_SP
-gsw_SP_from_SR
+*/
+
+/// Reference Salinity from Practical Salinity
+///
+/// # Examples
+/// ```
+/// use gsw::conversions::sr_from_sp;
+///
+/// let sr = sr_from_sp(32.0);
+/// assert!((sr - 32.150893714285715).abs() <= f64::EPSILON);
+/// ```
+pub fn sr_from_sp(sp: f64) -> f64 {
+    if cfg!(feature = "compat") {
+        sp * 1.004715428571429
+    } else {
+        sp * GSW_UPS
+    }
+}
+
+/// Practical Salinity from Reference Salinity
+///
+/// # Examples
+/// ```
+/// use gsw::conversions::sp_from_sr;
+///
+/// let sp = sp_from_sr(32.0);
+/// assert!((sp - 31.849814474830684).abs() <= f64::EPSILON);
+/// ```
+pub fn sp_from_sr(sr: f64) -> f64 {
+    if cfg!(feature = "compat") {
+        sr * 0.995306702338459
+    } else {
+        sr / GSW_UPS
+    }
+}
+
+/*
 gsw_SP_from_SA
 gsw_Sstar_from_SA
 gsw_SA_from_Sstar
 gsw_SP_from_Sstar
 gsw_pt_from_CT
 gsw_t_from_CT
-gsw_CT_from_pt
+*/
+
+/// Conservative Temperature from potential temperature
+///
+/// # Examples
+/// ```
+/// use gsw::conversions::ct_from_pt;
+///
+/// let ct = ct_from_pt(32.0, 10.0).unwrap();
+/// assert!((ct - 10.047455620469973).abs() <= f64::EPSILON);
+/// ```
+pub fn ct_from_pt(sa: f64, pt: f64) -> Result<f64> {
+    // Doesn't apply the offset so can't use non_dimensional_sa
+    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
+    };
+
+    let x2 = GSW_SFAC * sa;
+    let x = libm::sqrt(x2);
+    let y = pt * 0.025;
+    let pot_enthalpy = 61.01362420681071e0
+        + y * (168776.46138048015e0
+            + y * (-2735.2785605119625e0
+                + y * (2574.2164453821433e0
+                    + y * (-1536.6644434977543e0
+                        + y * (545.7340497931629e0
+                            + (-50.91091728474331e0 - 18.30489878927802e0 * y) * y)))))
+        + x2 * (268.5520265845071e0
+            + y * (-12019.028203559312e0
+                + y * (3734.858026725145e0
+                    + y * (-2046.7671145057618e0
+                        + y * (465.28655623826234e0
+                            + (-0.6370820302376359e0 - 10.650848542359153e0 * y) * y))))
+            + x * (937.2099110620707e0
+                + y * (588.1802812170108e0
+                    + y * (248.39476522971285e0
+                        + (-3.871557904936333e0 - 2.6268019854268356e0 * y) * y))
+                + x * (-1687.914374187449e0
+                    + x * (246.9598888781377e0
+                        + x * (123.59576582457964e0 - 48.5891069025409e0 * x))
+                    + y * (936.3206544460336e0
+                        + y * (-942.7827304544439e0
+                            + y * (369.4389437509002e0
+                                + (-33.83664947895248e0 - 9.987880382780322e0 * y) * y))))));
+
+    Ok(pot_enthalpy / GSW_CP0)
+}
+
+/*
 gsw_pot_enthalpy_from_pt
 gsw_pt_from_t
 gsw_pt0_from_t
@@ -99,8 +232,57 @@ pub fn z_from_p(
     -2.0 * c / (b + libm::sqrt(b * b - 4.0 * a * c))
 }
 
+/// Pressure from height (75-term polynomial approximation)
+///
+/// # Notes
+///
+/// - The GSW-C does not allow dynamic height anomaly or geopotential at zero
+///   sea pressure.
+///
+/// # Examples
+/// ```
+/// use gsw::conversions::p_from_z;
+///
+/// let z = p_from_z(-1000.0, 15., None, None).unwrap();
+/// dbg!(&z);
+/// assert!((z - 1008.321764487538).abs() <= f64::EPSILON);
+/// ```
+pub fn p_from_z(
+    z: f64,
+    lat: f64,
+    geo_strf_dyn_height: Option<f64>,
+    sea_surface_geopotental: Option<f64>,
+) -> Result<f64> {
+    if z > 5.0 {
+        return Err(Error::Undefined);
+    }
+
+    let sinlat = libm::sin(lat * DEG2RAD);
+    let sin2 = sinlat * sinlat;
+    let gs = 9.780327 * (1.0 + (5.2792e-3 + (2.32e-5 * sin2)) * sin2);
+
+    // get the first estimate of p from Saunders (1981)
+    let c1 = 5.25e-3 * sin2 + 5.92e-3;
+    let p = -2.0 * z / ((1.0 - c1) + libm::sqrt((1.0 - c1) * (1.0 - c1) + 8.84e-6 * z));
+    // end of the first estimate from Saunders (1981)
+
+    // initial value of the derivative of f
+    let df_dp = DB2PA * crate::gsw_internal_funcs::specvol_sso_0(p);
+
+    let f = crate::gsw_internal_funcs::enthalpy_sso_0(p)
+        + gs * (z - 0.5 * GAMMA * (z * z))
+        + geo_strf_dyn_height.unwrap_or(0.0)
+        + sea_surface_geopotental.unwrap_or(0.0);
+    let p_old = p;
+    let p = p_old - f / df_dp;
+    let p_mid = 0.5 * (p + p_old);
+    let df_dp = DB2PA * crate::gsw_internal_funcs::specvol_sso_0(p_mid);
+    let p = p_old - f / df_dp;
+
+    Ok(p)
+}
+
 /*
-gsw_p_from_z
 gsw_z_from_depth
 gsw_depth_from_z
 */