diff --git a/src/properties/thermo_fit.rs b/src/properties/thermo_fit.rs
index a5eba68..7ed1abb 100644
--- a/src/properties/thermo_fit.rs
+++ b/src/properties/thermo_fit.rs
@@ -22,7 +22,7 @@ pub struct SpeciesElement {
 }
 
 impl SpeciesThermoData {
-    pub fn cp_over_r(&self, temp: f64) -> f64 {
+    pub fn polynomial_at(&self, temp: f64) -> &ThermoPolynomial {
         //TODO: Not the most efficient. Can refactor to pre-compute tables
         //and do 1-d linear interpolation if needed
         //
@@ -33,35 +33,7 @@ impl SpeciesThermoData {
             .iter()
             .rposition(|polynomial| temp > polynomial.temp_range.0)
             .unwrap_or(0);
-        self.polynomials[i_polynomial].cp_over_r(temp)
-    }
-
-    pub fn h_over_rt(&self, temp: f64) -> f64 {
-        //TODO: Not the most efficient. Can refactor to pre-compute tables
-        //and do 1-d linear interpolation if needed
-        //
-        //TODO: I Think condensed species need to be treated differently. Verify how that works in
-        //the paper.
-        let i_polynomial = self
-            .polynomials
-            .iter()
-            .rposition(|polynomial| temp > polynomial.temp_range.0)
-            .unwrap_or(0);
-        self.polynomials[i_polynomial].h_over_rt(temp)
-    }
-
-    pub fn s_over_r(&self, temp: f64) -> f64 {
-        //TODO: Not the most efficient. Can refactor to pre-compute tables
-        //and do 1-d linear interpolation if needed
-        //
-        //TODO: I Think condensed species need to be treated differently. Verify how that works in
-        //the paper.
-        let i_polynomial = self
-            .polynomials
-            .iter()
-            .rposition(|polynomial| temp > polynomial.temp_range.0)
-            .unwrap_or(0);
-        self.polynomials[i_polynomial].s_over_r(temp)
+        &self.polynomials[i_polynomial]
     }
 }
 
@@ -108,24 +80,28 @@ impl ThermoPolynomial {
 #[cfg(test)]
 mod test {
     use super::ThermoPolynomial;
-    use crate::assert_delta;
-
-    fn poly(a: Vec<f64>) -> ThermoPolynomial {
-        ThermoPolynomial {
-            a,
-            temp_range: (0.0, 0.0),
-        }
-    }
+    use crate::{
+        assert_delta, assert_vec_delta,
+        properties::thermo_fit::{Phase, SpeciesThermoData},
+    };
 
     #[test]
     fn test_cp_over_r() {
         // At T=1: a[0]/T^2 + a[1]/T + a[2] + a[3]*T + a[4]*T^2 + a[5]*T^3 + a[6]*T^4
         // = 1 + 2 + 3 + 4 + 5 + 6 + 7 = 28
-        let result = poly(vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]).cp_over_r(1.0);
+        let result = ThermoPolynomial {
+            a: vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
+            temp_range: (0.0, 0.0),
+        }
+        .cp_over_r(1.0);
         assert_delta!(result, 28.0, 1e-10);
 
         // At T=2: 4/4 + 2/2 + 1 + 2*2 + 1*4 + 1*8 + 1*16 = 35
-        let result = poly(vec![4.0, 2.0, 1.0, 2.0, 1.0, 1.0, 1.0]).cp_over_r(2.0);
+        let result = ThermoPolynomial {
+            a: vec![4.0, 2.0, 1.0, 2.0, 1.0, 1.0, 1.0],
+            temp_range: (0.0, 0.0),
+        }
+        .cp_over_r(2.0);
         assert_delta!(result, 35.0, 1e-10);
     }
 
@@ -134,12 +110,20 @@ mod test {
         // At T=1: ln(1/T)=0, so the a[1] term vanishes
         // -a[0] + 0 + a[2] + a[3]/2 + a[4]/3 + a[5]/4 + a[6]/5 + a[7]
         // = -1 + 0 + 3 + 12/2 + 9/3 + 8/4 + 5/5 + 8 = -1+3+6+3+2+1+8 = 22
-        let result = poly(vec![1.0, 2.0, 3.0, 12.0, 9.0, 8.0, 5.0, 8.0]).h_over_rt(1.0);
+        let result = ThermoPolynomial {
+            a: vec![1.0, 2.0, 3.0, 12.0, 9.0, 8.0, 5.0, 8.0],
+            temp_range: (0.0, 0.0),
+        }
+        .h_over_rt(1.0);
         assert_delta!(result, 22.0, 1e-10);
 
         // At T=2: a[1]=0 so log term vanishes; a[4..6]=0 to avoid fractions
         // -4/4 + 0 + 3 + 4*2/2 + 0 + 0 + 0 + 2/2 = -1+3+4+1 = 7
-        let result = poly(vec![4.0, 0.0, 3.0, 4.0, 0.0, 0.0, 0.0, 2.0]).h_over_rt(2.0);
+        let result = ThermoPolynomial {
+            a: vec![4.0, 0.0, 3.0, 4.0, 0.0, 0.0, 0.0, 2.0],
+            temp_range: (0.0, 0.0),
+        }
+        .h_over_rt(2.0);
         assert_delta!(result, 7.0, 1e-10);
     }
 
@@ -148,12 +132,62 @@ mod test {
         // At T=1: ln(T)=0, so a[2] term vanishes
         // -a[0] - a[1] + 0 + a[3] + a[4]/2 + a[5]/3 + a[6]/4 + a[8]
         // = -1-2+0+4+6/2+12/3+8/4+5 = -1-2+4+3+4+2+5 = 15
-        let result = poly(vec![1.0, 2.0, 3.0, 4.0, 6.0, 12.0, 8.0, 0.0, 5.0]).s_over_r(1.0);
+        let result = ThermoPolynomial {
+            a: vec![1.0, 2.0, 3.0, 4.0, 6.0, 12.0, 8.0, 0.0, 5.0],
+            temp_range: (0.0, 0.0),
+        }
+        .s_over_r(1.0);
         assert_delta!(result, 15.0, 1e-10);
 
         // At T=2: a[2]=0 so log term vanishes; a[5..6]=0 to avoid fractions
         // -4/4 - 2/2 + 0 + 2*2 + 2*4/2 + 0 + 0 + 3 = -1-1+4+4+3 = 9
-        let result = poly(vec![4.0, 2.0, 0.0, 2.0, 2.0, 0.0, 0.0, 0.0, 3.0]).s_over_r(2.0);
+        let result = ThermoPolynomial {
+            a: vec![4.0, 2.0, 0.0, 2.0, 2.0, 0.0, 0.0, 0.0, 3.0],
+            temp_range: (0.0, 0.0),
+        }
+        .s_over_r(2.0);
         assert_delta!(result, 9.0, 1e-10);
     }
+
+    #[test]
+    fn test_polynomial_at() {
+        let polynomials = vec![
+            ThermoPolynomial {
+                a: vec![1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0],
+                temp_range: (1.0, 100.0),
+            },
+            ThermoPolynomial {
+                a: vec![4.0, 2.0, 1.0, 2.0, 1.0, 1.0, 1.0],
+                temp_range: (100.0, 300.0),
+            },
+        ];
+
+        let data = SpeciesThermoData {
+            name: "".to_string(),
+            elements: vec![],
+            phase: Phase::Gas,
+            polynomials,
+            molecular_weight: 1.0,
+            h_formation: 1.0,
+        };
+
+        assert!(std::ptr::eq(data.polynomial_at(0.5), &data.polynomials[0]));
+        assert!(std::ptr::eq(data.polynomial_at(50.0), &data.polynomials[0]));
+        assert!(std::ptr::eq(
+            data.polynomial_at(100.0),
+            &data.polynomials[0]
+        ));
+        assert!(std::ptr::eq(
+            data.polynomial_at(100.0 + 1e-12),
+            &data.polynomials[1]
+        ));
+        assert!(std::ptr::eq(
+            data.polynomial_at(150.0 + 1e-12),
+            &data.polynomials[1]
+        ));
+        assert!(std::ptr::eq(
+            data.polynomial_at(500.0 + 1e-12),
+            &data.polynomials[1]
+        ));
+    }
 }