diff --git a/src/properties/thermo_db.rs b/src/properties/thermo_db.rs
index 96cc957..7db0c05 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 {
-        name,
-        polynomials,
+    Ok(SpeciesThermoData::new(
+        &name,
         elements,
         phase,
+        polynomials,
         molecular_weight,
         h_formation,
-    })
+    ))
 }
 
 fn parse_polynomials_block<'a>(
@@ -295,9 +295,17 @@ mod test {
             -1.742171366e+01,
         ];
 
-        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);
+        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
+        );
 
         let real_coeff_2 = [
             -3.523782900e+05,
@@ -311,9 +319,17 @@ mod test {
             -2.695610360e+00,
         ];
 
-        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);
+        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
+        );
     }
 
     #[test]
@@ -363,10 +379,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.polynomials.len(), 2);
+        assert_eq!(alcl3.num_polynomials(), 2);
 
         assert_vec_delta!(
-            alcl3.polynomials[0].a,
+            alcl3.polynomial_at(650.0).unwrap().a,
             [
                 7.750600970e+04,
                 -1.440779717e+03,
@@ -380,11 +396,19 @@ END REACTANTS
             ],
             1e-9
         );
-        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_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_vec_delta!(
-            alcl3.polynomials[1].a,
+            alcl3.polynomial_at(3500.0).unwrap().a,
             [
                 -1.378630916e+05,
                 -5.579207290e+01,
@@ -398,8 +422,16 @@ END REACTANTS
             ],
             1e-9
         );
-        assert_delta!(alcl3.polynomials[1].temp_range.0, 1000.0, 1e-3);
-        assert_delta!(alcl3.polynomials[1].temp_range.1, 6000.0, 1e-3);
+        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
+        );
 
         // --- Air (reactant 0) ---
         let air = &thermo_db.reactants[0];
@@ -416,10 +448,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.polynomials.len(), 2);
+        assert_eq!(air.num_polynomials(), 2);
 
         assert_vec_delta!(
-            air.polynomials[0].a,
+            air.polynomial_at(650.0).unwrap().a,
             [
                 1.009950160e+04,
                 -1.968275610e+02,
@@ -433,11 +465,11 @@ END REACTANTS
             ],
             1e-9
         );
-        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_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_vec_delta!(
-            air.polynomials[1].a,
+            air.polynomial_at(3500.0).unwrap().a,
             [
                 2.415214430e+05,
                 -1.257874600e+03,
@@ -451,8 +483,16 @@ END REACTANTS
             ],
             1e-9
         );
-        assert_delta!(air.polynomials[1].temp_range.0, 1000.0, 1e-3);
-        assert_delta!(air.polynomials[1].temp_range.1, 6000.0, 1e-3);
+        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
+        );
 
         // --- n-Butanol (reactant 1) ---
         let butanol = &thermo_db.reactants[1];
@@ -467,6 +507,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.polynomials.len(), 0);
+        assert_eq!(butanol.num_polynomials(), 0);
     }
 }
diff --git a/src/properties/thermo_fit.rs b/src/properties/thermo_fit.rs
index 7ed1abb..893bbcf 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,
-    pub polynomials: Vec<ThermoPolynomial>,
+    polynomials: Vec<ThermoPolynomial>,
     pub molecular_weight: f64,
     pub h_formation: f64,
 }
@@ -22,18 +22,44 @@ pub struct SpeciesElement {
 }
 
 impl SpeciesThermoData {
-    pub fn polynomial_at(&self, temp: f64) -> &ThermoPolynomial {
+    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> {
         //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);
-        &self.polynomials[i_polynomial]
+        Some(&self.polynomials[i_polynomial])
     }
 }
 
@@ -55,7 +81,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 * inv_temp.ln()
+            + self.a[1] * inv_temp * temp.ln()
             + self.a[2]
             + self.a[3] * temp / 2.0
             + self.a[4] * temp * temp / 3.0
@@ -67,7 +93,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 - self.a[1] * inv_temp
+        -self.a[0] * inv_temp * inv_temp / 2.0 - self.a[1] * inv_temp
             + self.a[2] * temp.ln()
             + self.a[3] * temp
             + self.a[4] * temp * temp / 2.0
@@ -81,22 +107,19 @@ impl ThermoPolynomial {
 mod test {
     use super::ThermoPolynomial;
     use crate::{
-        assert_delta, assert_vec_delta,
+        assert_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(1.0);
-        assert_delta!(result, 28.0, 1e-10);
+        .cp_over_r(100.0);
+        assert_delta!(result, 706050403.0201, 1e-4);
 
-        // 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),
@@ -107,46 +130,36 @@ 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(1.0);
-        assert_delta!(result, 22.0, 1e-10);
+        .h_over_rt(2.0);
+        assert_delta!(result, 63.44314718055995, 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(2.0);
-        assert_delta!(result, 7.0, 1e-10);
+        .h_over_rt(100.0);
+        assert_delta!(result, 203.0196, 1e-4);
     }
 
     #[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![1.0, 2.0, 3.0, 4.0, 6.0, 12.0, 8.0, 0.0, 5.0],
+            a: vec![2.0, 3.0, 4.0, 1.0, 5.0, 2.0, 8.0, 1.0, 12.0],
             temp_range: (0.0, 0.0),
         }
-        .s_over_r(1.0);
-        assert_delta!(result, 15.0, 1e-10);
+        .s_over_r(100.0);
+        assert_delta!(result, 200691797.0572474, 1e-7);
 
-        // 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.0, 1e-10);
+        assert_delta!(result, 9.5, 1e-10);
     }
 
     #[test]
@@ -171,22 +184,28 @@ 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(100.0),
+            data.polynomial_at(0.5).unwrap(),
             &data.polynomials[0]
         ));
         assert!(std::ptr::eq(
-            data.polynomial_at(100.0 + 1e-12),
+            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.polynomials[1]
         ));
         assert!(std::ptr::eq(
-            data.polynomial_at(150.0 + 1e-12),
+            data.polynomial_at(150.0 + 1e-12).unwrap(),
             &data.polynomials[1]
         ));
         assert!(std::ptr::eq(
-            data.polynomial_at(500.0 + 1e-12),
+            data.polynomial_at(500.0 + 1e-12).unwrap(),
             &data.polynomials[1]
         ));
     }
diff --git a/src/properties/transport_db.rs b/src/properties/transport_db.rs
index 2ebb5f5..32bcc00 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 {
-        name,
+    Ok(SpeciesTransportData::new(
+        &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 8489d3c..6145a2d 100644
--- a/src/properties/transport_fit.rs
+++ b/src/properties/transport_fit.rs
@@ -1,7 +1,45 @@
 pub struct SpeciesTransportData {
     pub name: String,
-    pub viscosities: Vec<TransportFit>,
-    pub conductivities: Vec<TransportFit>,
+    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 struct TransportFit {
@@ -11,3 +49,148 @@ 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
+        );
+    }
+}