diff --git a/src/properties/thermo_db.rs b/src/properties/thermo_db.rs
index 96cc957..4468874 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 {
+    Ok(SpeciesThermoData::new(
         name,
-        polynomials,
         elements,
         phase,
+        polynomials,
         molecular_weight,
         h_formation,
-    })
+    ))
 }
 
 fn parse_polynomials_block<'a>(
@@ -295,9 +295,9 @@ 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).a, real_coeff_1, 1e-9);
+        assert_delta!(species.polynomial_at(650.0).temp_range.0, 300.000, 1e-3);
+        assert_delta!(species.polynomial_at(650.0).temp_range.1, 1000.000, 1e-3);
 
         let real_coeff_2 = [
             -3.523782900e+05,
@@ -311,9 +311,9 @@ 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).a, real_coeff_2, 1e-9);
+        assert_delta!(species.polynomial_at(3500.0).temp_range.0, 1000.000, 1e-3);
+        assert_delta!(species.polynomial_at(3500.0).temp_range.1, 6000.000, 1e-3);
     }
 
     #[test]
@@ -363,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.polynomials.len(), 2);
+        assert_eq!(alcl3.num_polynomials(), 2);
 
         assert_vec_delta!(
-            alcl3.polynomials[0].a,
+            alcl3.polynomial_at(650.0).a,
             [
                 7.750600970e+04,
                 -1.440779717e+03,
@@ -380,11 +380,11 @@ 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).temp_range.0, 300.0, 1e-3);
+        assert_delta!(alcl3.polynomial_at(650.0).temp_range.1, 1000.0, 1e-3);
 
         assert_vec_delta!(
-            alcl3.polynomials[1].a,
+            alcl3.polynomial_at(3500.0).a,
             [
                 -1.378630916e+05,
                 -5.579207290e+01,
@@ -398,8 +398,8 @@ 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).temp_range.0, 1000.0, 1e-3);
+        assert_delta!(alcl3.polynomial_at(3500.0).temp_range.1, 6000.0, 1e-3);
 
         // --- Air (reactant 0) ---
         let air = &thermo_db.reactants[0];
@@ -416,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.polynomials.len(), 2);
+        assert_eq!(air.num_polynomials(), 2);
 
         assert_vec_delta!(
-            air.polynomials[0].a,
+            air.polynomial_at(650.0).a,
             [
                 1.009950160e+04,
                 -1.968275610e+02,
@@ -433,11 +433,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).temp_range.0, 300.0, 1e-3);
+        assert_delta!(air.polynomial_at(650.0).temp_range.1, 1000.0, 1e-3);
 
         assert_vec_delta!(
-            air.polynomials[1].a,
+            air.polynomial_at(3500.0).a,
             [
                 2.415214430e+05,
                 -1.257874600e+03,
@@ -451,8 +451,8 @@ 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).temp_range.0, 1000.0, 1e-3);
+        assert_delta!(air.polynomial_at(3500.0).temp_range.1, 6000.0, 1e-3);
 
         // --- n-Butanol (reactant 1) ---
         let butanol = &thermo_db.reactants[1];
@@ -467,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.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..8a4c3fd 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,6 +22,28 @@ pub struct SpeciesElement {
 }
 
 impl SpeciesThermoData {
+    pub fn new(
+        name: String,
+        elements: Vec<SpeciesElement>,
+        phase: Phase,
+        polynomials: Vec<ThermoPolynomial>,
+        molecular_weight: f64,
+        h_formation: f64,
+    ) -> Self {
+        Self {
+            name,
+            elements,
+            phase,
+            polynomials,
+            molecular_weight,
+            h_formation,
+        }
+    }
+
+    pub fn num_polynomials(&self) -> usize {
+        self.polynomials.len()
+    }
+
     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
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 0f5ea20..60a7213 100644
--- a/src/properties/transport_fit.rs
+++ b/src/properties/transport_fit.rs
@@ -1,10 +1,22 @@
 pub struct SpeciesTransportData {
     pub name: String,
-    pub viscosities: Vec<TransportFit>,
-    pub conductivities: Vec<TransportFit>,
+    viscosities: Vec<TransportFit>,
+    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) -> f64 {
         let i_viscosity = self
             .viscosities