diff --git a/src/properties/thermo_db.rs b/src/properties/thermo_db.rs
index 7db0c05..96cc957 100644
--- a/src/properties/thermo_db.rs
+++ b/src/properties/thermo_db.rs
@@ -98,14 +98,14 @@ fn parse_species<'a>(
         lines.next().ok_or(PropertiesError::InvalidFile)?;
     }
 
-    Ok(SpeciesThermoData::new(
-        &name,
+    Ok(SpeciesThermoData {
+        name,
+        polynomials,
         elements,
         phase,
-        polynomials,
         molecular_weight,
         h_formation,
-    ))
+    })
 }
 
 fn parse_polynomials_block<'a>(
@@ -295,17 +295,9 @@ mod test {
             -1.742171366e+01,
         ];
 
-        assert_vec_delta!(species.polynomial_at(650.0).unwrap().a, real_coeff_1, 1e-9);
-        assert_delta!(
-            species.polynomial_at(650.0).unwrap().temp_range.0,
-            300.000,
-            1e-3
-        );
-        assert_delta!(
-            species.polynomial_at(650.0).unwrap().temp_range.1,
-            1000.000,
-            1e-3
-        );
+        assert_vec_delta!(species.polynomials[0].a, real_coeff_1, 1e-9);
+        assert_delta!(species.polynomials[0].temp_range.0, 300.000, 1e-3);
+        assert_delta!(species.polynomials[0].temp_range.1, 1000.000, 1e-3);
 
         let real_coeff_2 = [
             -3.523782900e+05,
@@ -319,17 +311,9 @@ mod test {
             -2.695610360e+00,
         ];
 
-        assert_vec_delta!(species.polynomial_at(3500.0).unwrap().a, real_coeff_2, 1e-9);
-        assert_delta!(
-            species.polynomial_at(3500.0).unwrap().temp_range.0,
-            1000.000,
-            1e-3
-        );
-        assert_delta!(
-            species.polynomial_at(3500.0).unwrap().temp_range.1,
-            6000.000,
-            1e-3
-        );
+        assert_vec_delta!(species.polynomials[1].a, real_coeff_2, 1e-9);
+        assert_delta!(species.polynomials[1].temp_range.0, 1000.000, 1e-3);
+        assert_delta!(species.polynomials[1].temp_range.1, 6000.000, 1e-3);
     }
 
     #[test]
@@ -379,10 +363,10 @@ END REACTANTS
         assert!(matches!(alcl3.phase, Phase::Gas));
         assert_delta!(alcl3.molecular_weight, 133.3405380, 1e-7);
         assert_delta!(alcl3.h_formation, -584678.863, 1e-3);
-        assert_eq!(alcl3.num_polynomials(), 2);
+        assert_eq!(alcl3.polynomials.len(), 2);
 
         assert_vec_delta!(
-            alcl3.polynomial_at(650.0).unwrap().a,
+            alcl3.polynomials[0].a,
             [
                 7.750600970e+04,
                 -1.440779717e+03,
@@ -396,19 +380,11 @@ END REACTANTS
             ],
             1e-9
         );
-        assert_delta!(
-            alcl3.polynomial_at(650.0).unwrap().temp_range.0,
-            300.0,
-            1e-3
-        );
-        assert_delta!(
-            alcl3.polynomial_at(650.0).unwrap().temp_range.1,
-            1000.0,
-            1e-3
-        );
+        assert_delta!(alcl3.polynomials[0].temp_range.0, 300.0, 1e-3);
+        assert_delta!(alcl3.polynomials[0].temp_range.1, 1000.0, 1e-3);
 
         assert_vec_delta!(
-            alcl3.polynomial_at(3500.0).unwrap().a,
+            alcl3.polynomials[1].a,
             [
                 -1.378630916e+05,
                 -5.579207290e+01,
@@ -422,16 +398,8 @@ END REACTANTS
             ],
             1e-9
         );
-        assert_delta!(
-            alcl3.polynomial_at(3500.0).unwrap().temp_range.0,
-            1000.0,
-            1e-3
-        );
-        assert_delta!(
-            alcl3.polynomial_at(3500.0).unwrap().temp_range.1,
-            6000.0,
-            1e-3
-        );
+        assert_delta!(alcl3.polynomials[1].temp_range.0, 1000.0, 1e-3);
+        assert_delta!(alcl3.polynomials[1].temp_range.1, 6000.0, 1e-3);
 
         // --- Air (reactant 0) ---
         let air = &thermo_db.reactants[0];
@@ -448,10 +416,10 @@ END REACTANTS
         assert!(matches!(air.phase, Phase::Gas));
         assert_delta!(air.molecular_weight, 28.9651159, 1e-7);
         assert_delta!(air.h_formation, -125.530, 1e-3);
-        assert_eq!(air.num_polynomials(), 2);
+        assert_eq!(air.polynomials.len(), 2);
 
         assert_vec_delta!(
-            air.polynomial_at(650.0).unwrap().a,
+            air.polynomials[0].a,
             [
                 1.009950160e+04,
                 -1.968275610e+02,
@@ -465,11 +433,11 @@ END REACTANTS
             ],
             1e-9
         );
-        assert_delta!(air.polynomial_at(650.0).unwrap().temp_range.0, 300.0, 1e-3);
-        assert_delta!(air.polynomial_at(650.0).unwrap().temp_range.1, 1000.0, 1e-3);
+        assert_delta!(air.polynomials[0].temp_range.0, 300.0, 1e-3);
+        assert_delta!(air.polynomials[0].temp_range.1, 1000.0, 1e-3);
 
         assert_vec_delta!(
-            air.polynomial_at(3500.0).unwrap().a,
+            air.polynomials[1].a,
             [
                 2.415214430e+05,
                 -1.257874600e+03,
@@ -483,16 +451,8 @@ END REACTANTS
             ],
             1e-9
         );
-        assert_delta!(
-            air.polynomial_at(3500.0).unwrap().temp_range.0,
-            1000.0,
-            1e-3
-        );
-        assert_delta!(
-            air.polynomial_at(3500.0).unwrap().temp_range.1,
-            6000.0,
-            1e-3
-        );
+        assert_delta!(air.polynomials[1].temp_range.0, 1000.0, 1e-3);
+        assert_delta!(air.polynomials[1].temp_range.1, 6000.0, 1e-3);
 
         // --- n-Butanol (reactant 1) ---
         let butanol = &thermo_db.reactants[1];
@@ -507,6 +467,6 @@ END REACTANTS
         assert!(matches!(butanol.phase, Phase::Condensed));
         assert_delta!(butanol.molecular_weight, 74.1216000, 1e-7);
         assert_delta!(butanol.h_formation, -278510.000, 1e-3);
-        assert_eq!(butanol.num_polynomials(), 0);
+        assert_eq!(butanol.polynomials.len(), 0);
     }
 }
diff --git a/src/properties/thermo_fit.rs b/src/properties/thermo_fit.rs
index 893bbcf..7ed1abb 100644
--- a/src/properties/thermo_fit.rs
+++ b/src/properties/thermo_fit.rs
@@ -7,7 +7,7 @@ pub struct SpeciesThermoData {
     pub name: String,
     pub elements: Vec<SpeciesElement>,
     pub phase: Phase,
-    polynomials: Vec<ThermoPolynomial>,
+    pub polynomials: Vec<ThermoPolynomial>,
     pub molecular_weight: f64,
     pub h_formation: f64,
 }
@@ -22,44 +22,18 @@ pub struct SpeciesElement {
 }
 
 impl SpeciesThermoData {
-    pub fn new(
-        name: &str,
-        elements: Vec<SpeciesElement>,
-        phase: Phase,
-        polynomials: Vec<ThermoPolynomial>,
-        molecular_weight: f64,
-        h_formation: f64,
-    ) -> Self {
-        Self {
-            name: name.to_string(),
-            elements,
-            phase,
-            polynomials,
-            molecular_weight,
-            h_formation,
-        }
-    }
-
-    pub fn num_polynomials(&self) -> usize {
-        self.polynomials.len()
-    }
-
-    pub fn polynomial_at(&self, temp: f64) -> Option<&ThermoPolynomial> {
+    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
         //
         //TODO: I Think condensed species need to be treated differently. Verify how that works in
         //the paper.
-        if self.polynomials.is_empty() {
-            return None;
-        }
-
         let i_polynomial = self
             .polynomials
             .iter()
             .rposition(|polynomial| temp > polynomial.temp_range.0)
             .unwrap_or(0);
-        Some(&self.polynomials[i_polynomial])
+        &self.polynomials[i_polynomial]
     }
 }
 
@@ -81,7 +55,7 @@ impl ThermoPolynomial {
     pub fn h_over_rt(&self, temp: f64) -> f64 {
         let inv_temp = 1.0 / temp;
         -self.a[0] * inv_temp * inv_temp
-            + self.a[1] * inv_temp * temp.ln()
+            + self.a[1] * inv_temp * inv_temp.ln()
             + self.a[2]
             + self.a[3] * temp / 2.0
             + self.a[4] * temp * temp / 3.0
@@ -93,7 +67,7 @@ impl ThermoPolynomial {
     /// NOTE: This is normalized and unitless
     pub fn s_over_r(&self, temp: f64) -> f64 {
         let inv_temp = 1.0 / temp;
-        -self.a[0] * inv_temp * inv_temp / 2.0 - self.a[1] * inv_temp
+        -self.a[0] * inv_temp * inv_temp - self.a[1] * inv_temp
             + self.a[2] * temp.ln()
             + self.a[3] * temp
             + self.a[4] * temp * temp / 2.0
@@ -107,19 +81,22 @@ impl ThermoPolynomial {
 mod test {
     use super::ThermoPolynomial;
     use crate::{
-        assert_delta,
+        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 = 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(100.0);
-        assert_delta!(result, 706050403.0201, 1e-4);
+        .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 = ThermoPolynomial {
             a: vec![4.0, 2.0, 1.0, 2.0, 1.0, 1.0, 1.0],
             temp_range: (0.0, 0.0),
@@ -130,36 +107,46 @@ mod test {
 
     #[test]
     fn test_h_over_rt() {
+        // 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 = 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(2.0);
-        assert_delta!(result, 63.44314718055995, 1e-10);
+        .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 = 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(100.0);
-        assert_delta!(result, 203.0196, 1e-4);
+        .h_over_rt(2.0);
+        assert_delta!(result, 7.0, 1e-10);
     }
 
     #[test]
     fn test_s_over_r() {
+        // 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 = ThermoPolynomial {
-            a: vec![2.0, 3.0, 4.0, 1.0, 5.0, 2.0, 8.0, 1.0, 12.0],
+            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(100.0);
-        assert_delta!(result, 200691797.0572474, 1e-7);
+        .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 = 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.5, 1e-10);
+        assert_delta!(result, 9.0, 1e-10);
     }
 
     #[test]
@@ -184,28 +171,22 @@ mod test {
             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(0.5).unwrap(),
+            data.polynomial_at(100.0),
             &data.polynomials[0]
         ));
         assert!(std::ptr::eq(
-            data.polynomial_at(50.0).unwrap(),
-            &data.polynomials[0]
-        ));
-        assert!(std::ptr::eq(
-            data.polynomial_at(100.0).unwrap(),
-            &data.polynomials[0]
-        ));
-        assert!(std::ptr::eq(
-            data.polynomial_at(100.0 + 1e-12).unwrap(),
+            data.polynomial_at(100.0 + 1e-12),
             &data.polynomials[1]
         ));
         assert!(std::ptr::eq(
-            data.polynomial_at(150.0 + 1e-12).unwrap(),
+            data.polynomial_at(150.0 + 1e-12),
             &data.polynomials[1]
         ));
         assert!(std::ptr::eq(
-            data.polynomial_at(500.0 + 1e-12).unwrap(),
+            data.polynomial_at(500.0 + 1e-12),
             &data.polynomials[1]
         ));
     }
diff --git a/src/properties/transport_db.rs b/src/properties/transport_db.rs
index 32bcc00..2ebb5f5 100644
--- a/src/properties/transport_db.rs
+++ b/src/properties/transport_db.rs
@@ -50,11 +50,11 @@ fn parse_species_transport_block<'a>(
             (ViscosityOrConductivity::Conductivity, fit) => conductivities.push(fit),
         }
     }
-    Ok(SpeciesTransportData::new(
-        &name,
+    Ok(SpeciesTransportData {
+        name,
         viscosities,
         conductivities,
-    ))
+    })
 }
 
 fn parse_species_header_line(line: &str) -> Result<(String, usize, usize), PropertiesError> {
diff --git a/src/properties/transport_fit.rs b/src/properties/transport_fit.rs
index 6145a2d..8489d3c 100644
--- a/src/properties/transport_fit.rs
+++ b/src/properties/transport_fit.rs
@@ -1,45 +1,7 @@
 pub struct SpeciesTransportData {
     pub name: String,
-    pub(crate) viscosities: Vec<TransportFit>,
-    pub(crate) conductivities: Vec<TransportFit>,
-}
-
-impl SpeciesTransportData {
-    pub fn new(
-        name: &str,
-        viscosities: Vec<TransportFit>,
-        conductivities: Vec<TransportFit>,
-    ) -> Self {
-        SpeciesTransportData {
-            name: name.to_string(),
-            viscosities,
-            conductivities,
-        }
-    }
-
-    pub fn viscosity_at(&self, temp: f64) -> Option<f64> {
-        if self.viscosities.is_empty() {
-            return None;
-        }
-        let i_viscosity = self
-            .viscosities
-            .iter()
-            .rposition(|viscosity| temp > viscosity.temp_range.0)
-            .unwrap_or(0);
-        Some(self.viscosities[i_viscosity].compute(temp))
-    }
-
-    pub fn conductivity_at(&self, temp: f64) -> Option<f64> {
-        if self.conductivities.is_empty() {
-            return None;
-        }
-        let i_conductivity = self
-            .conductivities
-            .iter()
-            .rposition(|conductivity| temp > conductivity.temp_range.0)
-            .unwrap_or(0);
-        Some(self.conductivities[i_conductivity].compute(temp))
-    }
+    pub viscosities: Vec<TransportFit>,
+    pub conductivities: Vec<TransportFit>,
 }
 
 pub struct TransportFit {
@@ -49,148 +11,3 @@ pub struct TransportFit {
     pub c: f64,
     pub d: f64,
 }
-
-impl TransportFit {
-    pub fn compute(&self, temp: f64) -> f64 {
-        let inv_temp = 1.0 / temp;
-        self.a * temp.ln() + self.b * inv_temp + self.c * inv_temp * inv_temp + self.d
-    }
-}
-
-#[cfg(test)]
-mod test {
-    use crate::{
-        assert_delta,
-        properties::transport_fit::{SpeciesTransportData, TransportFit},
-    };
-
-    #[test]
-    fn test_transport_fit_compute() {
-        let fit = TransportFit {
-            temp_range: (1.0, 2.0),
-            a: 10.0,
-            b: 20.0,
-            c: 30.0,
-            d: 40.0,
-        };
-
-        assert_delta!(fit.compute(4.0), 60.73794361119891, 1e-12);
-    }
-
-    #[test]
-    fn test_calc_transport_properties() {
-        let viscosities = vec![
-            TransportFit {
-                temp_range: (1.0, 100.0),
-                a: 10.0,
-                b: 20.0,
-                c: 30.0,
-                d: 40.0,
-            },
-            TransportFit {
-                temp_range: (100.0, 200.0),
-                a: 1.0,
-                b: 2.0,
-                c: 3.0,
-                d: 4.0,
-            },
-        ];
-
-        let conductivities = vec![
-            TransportFit {
-                temp_range: (1.0, 100.0),
-                a: 10.0,
-                b: 20.0,
-                c: 30.0,
-                d: 40.0,
-            },
-            TransportFit {
-                temp_range: (100.0, 200.0),
-                a: 1.0,
-                b: 2.0,
-                c: 3.0,
-                d: 4.0,
-            },
-        ];
-
-        let data = SpeciesTransportData {
-            viscosities,
-            conductivities,
-            name: "".to_string(),
-        };
-
-        assert_delta!(
-            data.conductivity_at(0.5),
-            data.conductivities[0].compute(0.5),
-            1e-12
-        );
-        assert_delta!(
-            data.conductivity_at(1.0),
-            data.conductivities[0].compute(1.0),
-            1e-12
-        );
-        assert_delta!(
-            data.conductivity_at(50.0),
-            data.conductivities[0].compute(50.0),
-            1e-12
-        );
-        assert_delta!(
-            data.conductivity_at(100.0),
-            data.conductivities[0].compute(100.0),
-            1e-12
-        );
-        assert_delta!(
-            data.conductivity_at(100.0 + 1e-12),
-            data.conductivities[1].compute(100.0 + 1e-12),
-            1e-12
-        );
-
-        assert_delta!(
-            data.conductivity_at(200.0),
-            data.conductivities[1].compute(200.0),
-            1e-12
-        );
-        assert_delta!(
-            data.conductivity_at(500.0),
-            data.conductivities[1].compute(500.0),
-            1e-12
-        );
-
-        assert_delta!(
-            data.viscosity_at(0.5),
-            data.viscosities[0].compute(0.5),
-            1e-12
-        );
-        assert_delta!(
-            data.viscosity_at(1.0),
-            data.viscosities[0].compute(1.0),
-            1e-12
-        );
-        assert_delta!(
-            data.viscosity_at(50.0),
-            data.viscosities[0].compute(50.0),
-            1e-12
-        );
-        assert_delta!(
-            data.viscosity_at(100.0),
-            data.viscosities[0].compute(100.0),
-            1e-12
-        );
-        assert_delta!(
-            data.viscosity_at(100.0 + 1e-12),
-            data.viscosities[1].compute(100.0 + 1e-12),
-            1e-12
-        );
-
-        assert_delta!(
-            data.viscosity_at(200.0),
-            data.viscosities[1].compute(200.0),
-            1e-12
-        );
-        assert_delta!(
-            data.viscosity_at(500.0),
-            data.viscosities[1].compute(500.0),
-            1e-12
-        );
-    }
-}