// Fit the Steinhart-hart thermistor coefficients to measured temp-resistance data
// apply https://www.thinksrs.com/downloads/PDFs/ApplicationNotes/LDC%20Note%204%20NTC%20Calculatorold.pdf
// or https://www.thinksrs.com/downloads/programs/therm%20calc/ntccalibrator/ntccalculator.html
// to fitting an NTC thermistor curve with an Arduino
// and write a thermistor program
// for https://forum.arduino.cc/t/where-do-you-buy-your-ky-013-thermistor-module/1103253
// Code on Wokwi at https://wokwi.com/projects/359481663807896577

#include <BasicLinearAlgebra.h>

using namespace BLA;

// Measure and set these values:
float resistances[] = {135248.21, 10019.55, 1271.33}; // Ohms
float temperatures[] = {-24.0, 25.0, 80.0};  // °C
float c1, c2, c3;
constexpr int DEBUG = 0;

constexpr int NTC_pin = A1; // pin with NTC pull-down and 10Kohm pullup
constexpr int program_pin = 3; // for printing out program

void setup()
{
  Serial.begin(115200);
  pinMode(program_pin, INPUT_PULLUP);

  const auto nObs = sizeof(resistances)/sizeof(resistances[0]);

  // Let's go ahead and declare a system of equations
  Matrix<nObs, 3>  A ;

  Matrix<nObs> b ;
  A.Fill(0.0);
  b.Fill(0.0);

  // populate with the thermistor data points
  for (auto i = 0; i < nObs; ++i) {
    A(i, 0) = 1.0;
    A(i, 1) = log(resistances[i]);
    A(i, 2) = pow(log(resistances[i]), 3);
    b(i) = 1 / (temperatures[i] + 273.15);
  }

  // To solve for x you might be tempted to take the inverse of A then
  // calculate x like so: x = A^-1 * b

  // As it turns out though, actually taking the inverse of A to solve these kinds of equations is quite inefficient.
  // Textbooks often talk about x = A^-1 * b, but in practice we don't actually compute x this way.

  // Instead we use something called an LU decomposition (or if A has some special properties you can use some other
  // type of decomposition but we won't get into that here). LU decomposition factors an A matrix into a permutation
  // matrix, a lower and an upper triangular matrix:
  auto A_decomp = A;  // LUDecompose will destroy A here so we'll pass in a copy so we can refer back to A later
  auto decomp = LUDecompose(A_decomp);

  // You can take a look at these matrices if you like:
  decomp.P();  // P essentially rearranges the rows of the matrix that results when we multiply L by U
  decomp.L();  // Will have ones along its diagonal and zeros above it
  decomp.U();  // Will have zeros below its diagonal

  // And if we multiply them all together we'll recover the original A matrix:
  if (DEBUG)Serial << "reconstructed A: " << decomp.P() * decomp.L() * decomp.U() << "\n";

  // Once we've done the decomposition we can solve for x very efficiently:
  Matrix<nObs> x_lusolve = LUSolve(decomp, b);

  if (DEBUG)Serial << "x (via LU decomposition): " << x_lusolve << "\n";
  c1 = x_lusolve(0);
  c2 = x_lusolve(1);
  c3 = x_lusolve(2);


  //
  char buff[128];

  // 1/T = A+ B ln Rt + ⋅ ( ) C ln Rt ( ) ( ) 3

  // Print the coefficients:
  // ##########################################################################
  //###########################################################################
  Serial.println("// Three point calibration Steinhart-Hart coefficient fitting");
  Serial.println("//Measured Data:");
  for(auto &t:temperatures){
    int index = &t -temperatures;
    Serial.print("// data ");
    Serial.print(index);
    Serial.print(" T:");
    Serial.print(t);
    Serial.print(" R:");
    Serial.print(resistances[index]);
    Serial.println();

  }
  Serial.println("// Fitted Coefficients");
  Serial.print(F("// float c1 = "));
  dtostre(c1,buff,7,0);
  Serial.print(buff);
  Serial.println(F("; //coefficient A"));

  Serial.print(F("// float c2 = "));
  dtostre(c2,buff,7,0);
  Serial.print(buff);
  Serial.println(F("; //coefficient B"));

  Serial.print(F("// float c3 = "));
  dtostre(c3,buff,7,0);
  Serial.print(buff);
  Serial.println(F("; //coefficient C"));

}

void loop() {

  if (digitalRead(program_pin) == LOW) {
    printProgram();
    while (digitalRead(program_pin) == LOW) {
      ;
    }
  } else {
    getTemperature(NTC_pin);
  }

}

float getTemperature(int ThermistorPin) {
  constexpr float R1 = 10000.0;
  float Vo = analogRead(ThermistorPin) + 0.5;
  float R2 = R1 / (1024.0 / (float)Vo - 1.0); // VCC-R1-R_therm-GND
  //R2 = R1 * ( 1024.0/Vo -1 );         // VCC-R_therm-R-GND (opposite wiring)

  float logR2 = log(R2);
  float T = (1.0 / (c1 + c2 * logR2 + c3 * logR2 * logR2 * logR2));
  float Tc = T - 273.15;
  float Tf = (Tc * 9.0) / 5.0 + 32.0;

  Serial.print("Resistance:  ");
  Serial.print(R2);
  Serial.print("  Temperature:   ");
  Serial.print(Tc);
  Serial.println(" C");

  delay(500);

}


void printProgram() {

  char buff[20];
  Serial.print(F(R"(
//  Steinhart-Hart coefficients fitted to this data per 
//  https://wokwi.com/projects/359481663807896577
// 

)"));
  // Print the coefficients:
  // ##########################################################################
  //###########################################################################
  Serial.println("// Three point calibration Steinhart-Hart coefficient fitting");
  Serial.println("//Measured Data:");
  for(auto &t:temperatures){
    int index = &t -temperatures;
    Serial.print("// data ");
    Serial.print(index);
    Serial.print(" T:");
    Serial.print(t);
    Serial.print(" R:");
    Serial.print(resistances[index]);
    Serial.println();
  }
  Serial.println("// Fitted Coefficients");
  Serial.print(F("// float c1 = "));
  dtostre(c1,buff,7,0);
  Serial.print(buff);
  Serial.println(F("; //coefficient A"));

  Serial.print(F("// float c2 = "));
  dtostre(c2,buff,7,0);
  Serial.print(buff);
  Serial.println(F("; //coefficient B"));

  Serial.print(F("// float c3 = "));
  dtostre(c3,buff,7,0);
  Serial.print(buff);
  Serial.println(F("; //coefficient C"));


  Serial.print(
    F(R"(
// From https://forum.arduino.cc/t/problema-usando-termistores-pasate-aca/513763
// and https://www.circuitbasics.com/arduino-thermistor-temperature-sensor-tutorial/
// With modified values
//This is an example code on how to read a thermistor, the "Thermimistor.h" Lib out there only acepts Beta 
//coeficient and in my case yield to incorrects results, this a way more accuerrate way to read the 
//thermistor, in case you have odd or wrong meassurements please follow this steps: 
//
//For get the acurrate results for this code you will need; 
//a multymeter, a NTC thermistor, another accurrate themperature 
//probe meter. 

//Step 1.- Set multimiter on resistance meassurement mode 

//Step 2.- Read and anotate the actual resistance of the thermistor
//and the actual temperature (allow 1min to get stable meassurement).
//Some Hot water and a cup. 

//Step 3.- place both sensors (Thermistor and temperature probe in a 
//recipient containing water at ambient temperature).
//In another cup heat up some water. 
//Add hot water until you heat more than 10°C the temp probe, wait for
//stable meassurement and anotate the temperature and the resistance. 
//Add more water to heat up the element 20° from the first meassurement. 
//Take note of the temperature and resistance. 

//Step 4.- 
//Go to the website: 
//http://www.thinksrs.com/downloads/programs/Therm%20Calc/NTCCalibrator/NTCcalculator.htm

//and set your data on it. 
//The calculator will deliver three values we need on the code: A, B and C. 

//Step 5.- 
//Replace the values you get in the calculator on this code. 

//Step 6.- Upload and test it. 
//Place both sensors on ambien water, warm water and hot water, use the temperature
//probe to chek for accurracy. 

//Original code from: https://www.youtube.com/watch?v=-_XkGju35MI

//Procedure: Alex Santiago - 12/03/2018
//Tested on a 10K NTC B3450 -12/03/2018
//Arduino Mega at 5Vcc.

int ThermistorPin = A1;
float Vo;
float R1 = 10000;
float logR2, R2, T, Tc, Tf;
)"));

  Serial.print(F("float c1 = "));
  dtostre(c1,buff,7,0);
  Serial.print(buff);
  Serial.println(F("; //coefficient A"));

  Serial.print(F("float c2 = "));
  dtostre(c2,buff,7,0);
  Serial.print(buff);
  Serial.println(F("; //coefficient B"));

  Serial.print(F("float c3 = "));
  dtostre(c3,buff,7,0);
  Serial.print(buff);
  Serial.println(F("; //coefficient C"));

  Serial.print(F(R"(
void setup() {
Serial.begin(9600);
}

void loop() {

  Vo = analogRead(ThermistorPin)+0.5;
  // There are two ways to wire a thermistor as R2:
  R2 = R1 / (1024.0 / (float)Vo - 1.0); // VCC-R1-R2-GND (pullup+NTC)
  //R2 = R1 * ( 1024.0/Vo -1 );         // VCC-R2-R1-GND  (NTC+pulldown)

  logR2 = log(R2);
  T = (1.0 / (c1 + c2*logR2 + c3*logR2*logR2*logR2));
  Tc = T - 273.15;
  Tf = (Tc * 9.0)/ 5.0 + 32.0; 

  Serial.print("Resistance:  ");
  Serial.print(R2);
  Serial.print("  Temperature:   "); 
  Serial.print(Tc);
  Serial.println(" C"); 

  delay(500);
}
//
)"));


}