From b7e294cc74c04e0d005cdab2312d374b78099e50 Mon Sep 17 00:00:00 2001
From: Gernot Bauer <bauer@itt.uni-stuttgart.de>
Date: Thu, 20 Jan 2022 16:37:34 +0100
Subject: [PATCH] Add partial derivatives of pressure wrt density and
 isothermal compressibility

---
 src/python/state.rs     | 44 +++++++++++++++++++++++++++++++++++++++++
 src/state/properties.rs | 17 ++++++++++++++++
 2 files changed, 61 insertions(+)

diff --git a/src/python/state.rs b/src/python/state.rs
index 669b293..0106b17 100644
--- a/src/python/state.rs
+++ b/src/python/state.rs
@@ -368,6 +368,23 @@ macro_rules! impl_state {
                 PySINumber::from(self.0.dp_dv(contributions.0))
             }
 
+            /// Return partial derivative of pressure w.r.t. density.
+            ///
+            /// Parameters
+            /// ----------
+            /// contributions: Contributions, optional
+            ///     the contributions of the helmholtz energy.
+            ///     Defaults to Contributions.Total.
+            ///
+            /// Returns
+            /// -------
+            /// SINumber
+            #[args(contributions = "PyContributions::Total()")]
+            #[pyo3(text_signature = "($self, contributions)")]
+            fn dp_drho(&self, contributions: PyContributions) -> PySINumber {
+                PySINumber::from(self.0.dp_drho(contributions.0))
+            }
+
             /// Return partial derivative of pressure w.r.t. temperature.
             ///
             /// Parameters
@@ -419,6 +436,23 @@ macro_rules! impl_state {
                 PySINumber::from(self.0.d2p_dv2(contributions.0))
             }
 
+            /// Return second partial derivative of pressure w.r.t. density.
+            ///
+            /// Parameters
+            /// ----------
+            /// contributions: Contributions, optional
+            ///     the contributions of the helmholtz energy.
+            ///     Defaults to Contributions.Total.
+            ///
+            /// Returns
+            /// -------
+            /// SINumber
+            #[args(contributions = "PyContributions::Total()")]
+            #[pyo3(text_signature = "($self, contributions)")]
+            fn d2p_drho2(&self, contributions: PyContributions) -> PySINumber {
+                PySINumber::from(self.0.d2p_drho2(contributions.0))
+            }
+
             /// Return molar volume of each component.
             ///
             /// Parameters
@@ -866,6 +900,16 @@ macro_rules! impl_state {
                 PySINumber::from(self.0.isentropic_compressibility())
             }
 
+            /// Return isothermal compressibility coefficient.
+            ///
+            /// Returns
+            /// -------
+            /// SINumber
+            #[pyo3(text_signature = "($self)")]
+            fn isothermal_compressibility(&self) -> PySINumber {
+                PySINumber::from(self.0.isothermal_compressibility())
+            }
+
             /// Return structure factor.
             ///
             /// Returns
diff --git a/src/state/properties.rs b/src/state/properties.rs
index 70ac8c1..7bd9bef 100644
--- a/src/state/properties.rs
+++ b/src/state/properties.rs
@@ -243,6 +243,11 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
         self.evaluate_property(Self::dp_dv_, contributions, true)
     }
 
+    /// Partial derivative of pressure w.r.t. density: $\left(\frac{\partial p}{\partial \rho}\right)_{T,N_i}$
+    pub fn dp_drho(&self, contributions: Contributions) -> QuantityScalar<U> {
+        -self.volume / self.density * self.dp_dv(contributions)
+    }
+
     /// Partial derivative of pressure w.r.t. temperature: $\left(\frac{\partial p}{\partial T}\right)_{V,N_i}$
     pub fn dp_dt(&self, contributions: Contributions) -> QuantityScalar<U> {
         self.evaluate_property(Self::dp_dt_, contributions, true)
@@ -258,6 +263,12 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
         self.evaluate_property(Self::d2p_dv2_, contributions, true)
     }
 
+    /// Second partial derivative of pressure w.r.t. density: $\left(\frac{\partial^2 p}{\partial \rho^2}\right)_{T,N_j}$
+    pub fn d2p_drho2(&self, contributions: Contributions) -> QuantityScalar<U> {
+        self.volume / (self.density * self.density)
+            * (self.volume * self.d2p_dv2(contributions) + 2.0 * self.dp_dv(contributions))
+    }
+
     /// Partial molar volume: $v_i=\left(\frac{\partial V}{\partial N_i}\right)_{T,p,N_j}$
     pub fn molar_volume(&self, contributions: Contributions) -> QuantityArray1<U> {
         let func = |s: &Self, evaluate: Evaluate| -s.dp_dni_(evaluate) / s.dp_dv_(evaluate);
@@ -444,6 +455,12 @@ impl<U: EosUnit, E: EquationOfState> State<U, E> {
         -self.c_v(c) / (self.c_p(c) * self.dp_dv(c) * self.volume)
     }
 
+    /// Isothermal compressibility: $\kappa_T=-\frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{T,N_i}$
+    pub fn isothermal_compressibility(&self) -> QuantityScalar<U> {
+        let c = Contributions::Total;
+        -1.0 / (self.dp_dv(c) * self.volume)
+    }
+
     /// Structure factor: $S(0)=k_BT\left(\frac{\partial\rho}{\partial p}\right)_{T,N_i}$
     pub fn structure_factor(&self) -> f64 {
         -(U::gas_constant() * self.temperature * self.density)