From ad45ea71754582f03c8f0d0946f690f8e8e1362f Mon Sep 17 00:00:00 2001
From: Zachary Sunforge <zachary.sunforge@bfpower.io>
Date: Fri, 7 Jul 2023 14:41:50 -0700
Subject: [PATCH] Thermocouple type K conversion

---
 .../conversion/thermocouples/type_k.rs        | 91 ++++++++++++++-----
 1 file changed, 66 insertions(+), 25 deletions(-)

diff --git a/src/transducer/conversion/thermocouples/type_k.rs b/src/transducer/conversion/thermocouples/type_k.rs
index e59b39e..8bc742e 100644
--- a/src/transducer/conversion/thermocouples/type_k.rs
+++ b/src/transducer/conversion/thermocouples/type_k.rs
@@ -1,17 +1,21 @@
 use crate::transducer::InvalidValue;
-use uom::si::electric_potential::millivolt;
+use uom::si::electric_potential::{millivolt, volt};
 use uom::si::f32;
 use uom::si::thermodynamic_temperature::degree_celsius;
 
-/// Convert from a voltage produced by a type k thermocouple to a temperature.
+/// Convert from a voltage produced by a type k thermocouple to a temperature using polynomial with
+/// NIST coefficients, this method has a maximum error of 0.06°C in addition to the underlying
+/// thermocouple inaccuracy.
 ///
 /// This function uses the [NIST type K thermocouple linearisation polynomial](https://srdata.nist.gov/its90/type_k/kcoefficients_inverse.html).
 pub fn convert(
     voltage: f32::ElectricPotential,
-    cold_junction: f32::ThermodynamicTemperature,
+    r_junction: f32::ThermodynamicTemperature,
+    r_junction_offset: fn(
+        r_junction: f32::ThermodynamicTemperature,
+    ) -> Result<f32::ElectricPotential, InvalidValue>,
 ) -> Result<f32::ThermodynamicTemperature, InvalidValue> {
-    //TODO: Add cold junction correction
-    let mv = voltage.get::<millivolt>();
+    let mv = voltage.get::<millivolt>() + r_junction_offset(r_junction)?.get::<millivolt>();
     let mv_pow2 = mv * mv;
     let mv_pow3 = mv_pow2 * mv;
     let mv_pow4 = mv_pow3 * mv;
@@ -23,14 +27,14 @@ pub fn convert(
         let mv_pow8 = mv_pow7 * mv;
 
         let celsius = 1.0
-            + 2.5173462e1 * mv
+            + 2.5173462E+1 * mv
             + -1.1662878 * mv_pow2
             + -1.0833638 * mv_pow3
-            + -8.9773540e-1 * mv_pow4
-            + -3.7342377e-1 * mv_pow5
-            + -8.6632643e-2 * mv_pow6
-            + -1.0450598e-2 * mv_pow7
-            + -5.1920577e-4 * mv_pow8;
+            + -8.9773540E-1 * mv_pow4
+            + -3.7342377E-1 * mv_pow5
+            + -8.6632643E-2 * mv_pow6
+            + -1.0450598E-2 * mv_pow7
+            + -5.1920577E-4 * mv_pow8;
 
         Ok(f32::ThermodynamicTemperature::new::<degree_celsius>(celsius))
     } else if mv > 0.0 && mv < 20.644 {
@@ -39,28 +43,65 @@ pub fn convert(
         let mv_pow9 = mv_pow8 * mv;
 
         let celsius = 1.0
-            + 2.508355e1 * mv
-            + 7.860106e-2 * mv_pow2
-            + -2.503131e-1 * mv_pow3
-            + 8.315270e-2 * mv_pow4
-            + -1.228034e-2 * mv_pow5
-            + 9.804036e-4 * mv_pow6
-            + -4.413030e-5 * mv_pow7
-            + 1.057734e-6 * mv_pow8
-            + -1.052755e-8 * mv_pow9;
+            + 2.508355E+1 * mv
+            + 7.860106E-2 * mv_pow2
+            + -2.503131E-1 * mv_pow3
+            + 8.315270E-2 * mv_pow4
+            + -1.228034E-2 * mv_pow5
+            + 9.804036E-4 * mv_pow6
+            + -4.413030E-5 * mv_pow7
+            + 1.057734E-6 * mv_pow8
+            + -1.052755E-8 * mv_pow9;
 
         Ok(f32::ThermodynamicTemperature::new::<degree_celsius>(celsius))
     } else if mv >= 20.644 && mv <= 54.886 {
         let celsius = 1.318058e2
-            + 4.830222e1 * mv
+            + 4.830222E+1 * mv
             + -1.646031 * mv_pow2
-            + 5.464731e-2 * mv_pow3
-            + -9.650715e-4 * mv_pow4
-            + 8.802193e-6 * mv_pow5
-            + -3.110810e-8 * mv_pow6;
+            + 5.464731E-2 * mv_pow3
+            + -9.650715E-4 * mv_pow4
+            + 8.802193E-6 * mv_pow5
+            + -3.110810E-8 * mv_pow6;
 
         Ok(f32::ThermodynamicTemperature::new::<degree_celsius>(celsius))
     } else {
         Err(InvalidValue)
     }
 }
+
+pub fn r_junction_offset_poly(
+    temperature: f32::ThermodynamicTemperature,
+) -> Result<f32::ElectricPotential, InvalidValue> {
+    const T0: f32 = 2.5000000E+01;
+    const V0: f32 = 1.0003453E+00;
+    const P1: f32 = 4.0514854E-02;
+    const P2: f32 = -3.8789638E-05;
+    const P3: f32 = -2.8608478E-06;
+    const P4: f32 = -9.5367041E-10;
+    const Q1: f32 = -1.3948675E-03;
+    const Q2: f32 = -6.7976627E-05;
+
+    let celsius = temperature.get::<degree_celsius>();
+    if celsius >= -20.0 && celsius <= 70.0 {
+        let offset_tcj = celsius - T0;
+        let numerator = offset_tcj + (P1 + offset_tcj * (P2 + offset_tcj * (P3 + P4 * offset_tcj)));
+        let denominator = 1.0 + offset_tcj * (Q1 + Q2 * offset_tcj);
+        let mv = V0 + numerator / denominator;
+
+        Ok(f32::ElectricPotential::new::<millivolt>(mv))
+    } else {
+        Err(InvalidValue)
+    }
+}
+
+pub fn r_junction_offset_lin(
+    temperature: f32::ThermodynamicTemperature,
+) -> Result<f32::ElectricPotential, InvalidValue> {
+    let celsius = temperature.get::<degree_celsius>();
+    if celsius >= -2.0 && celsius <= 1000.0 {
+        let mv = 0.041 * celsius;
+        Ok(f32::ElectricPotential::new::<millivolt>(mv))
+    } else {
+        Err(InvalidValue)
+    }
+}