//1D heat transfer simulation through aluminum rod, heated on one end, 
// https://wokwi.com/projects/359193529297155073
// for https://forum.arduino.cc/t/pid-with-simulated-heater/1093539

#include <BasicLinearAlgebra.h> // https://github.com/tomstewart89/BasicLinearAlgebra
/*
  This models a 1-D heat transfer problem of flow through with n_ele elements

        convection--v                                          v- convection at ends to Tamb
     heater q_in -> [###round aluminum rod(rod_dia, rod_len)###]

  Code at https://wokwi.com/projects/359193529297155073

*/

// All the functions in BasicLinearAlgebra are wrapped up inside the namespace BLA, so specify that we're using it like
// so:
using namespace BLA;

void setup()
{
  Serial.begin(115200);
}

// https://wokwi.com/projects/359235477558390785
// 0-D heat transfer through 1x1x2sm aluminum block
//  T = T + Q * interval / 1000 / mass / Cps - (T - Tamb) * area * h;
float h = 5; // W/m2K thermal convection coefficient for Al passive
float Cps = 0.89; // J/g°C
float C_tcond = 237; // W/mK thermal conduction coefficient for Al
float Tmelt = 660; // °C
float rho = 2.710e6; // kg/m^3 density // Al: 2,710kg/m3
float rod_dia = 0.02; // m
float rod_len = 0.05; //m
float area = pow(rod_dia, 2) * PI / 4 + 1 * rod_len / 3 * PI * rod_dia; // m2 area for convection
float mass = rod_len * pow(rod_dia, 2) * PI / 4 * rho ; // g
float Tamb = 25; // °C


const int n_ele = 4;
const int n_force = 2;
// setup for three-element 1D heat transfer model
BLA::Matrix<n_ele> x; //initial
BLA::Matrix<n_force> u; // Forcings at ends Watts
BLA::Matrix<n_ele, 2> G; // map from forcings to elements, dT=f(W)

//
BLA::Matrix<n_ele, n_ele> Fa ; // transition matrix

float dt = 0.05;
int dtPrint = 500;

void loop() {
  // Or as a more real life example, here's how you might calculate an updated state estimate for a third order state
  // space model:
  static float t = 0;
  static uint32_t last = 0;
  Tamb = analogRead(A1) * 250.0 / 1024;
  float q = analogRead(A0) * 10.0 / 1024;  // 10W heater
  static bool oneShot = true;
  if (oneShot == true) {
    oneShot = false;
    for (int i = 0; i < n_ele; ++i) {
      x(i) = Tamb;
      //
      if (i > 0) { // conduction to the left
        Fa(i, i - 1) = C_tcond * pow(rod_dia, 2) * PI / 4 / (rod_len / n_ele);
        Fa(i, i) -= Fa(i, i - 1);
      }
      if (i < n_ele - 1) { // conduction to the right
        Fa(i, i + 1) = C_tcond * pow(rod_dia, 2) * PI / 4 / (rod_len / n_ele);
        Fa(i, i) -= Fa(i, i + 1);
      }
    }
    Fa = Fa * 1.0f ; // *20 extra conductive
    G(0, 0) = Cps;
    G(n_ele - 1, 1) = Cps;
    u(0) = 0;
    u(n_force - 1) = 0;
    Serial <<"F:"<< Fa << "\n";
    Serial <<"G:" << G << "\n";
    Serial.println("Oneshot initialization finished");
  }
  if (millis() - last >= dtPrint) {
    last += dtPrint;
    Serial << "t:" << t << " Tamb=" << Tamb << " q=" << q << " u=" << u << " x: " << x << "\n";
  }
  // convection + heater
  u(0) = ((Tamb - x(0)) * area * h) + q;
  // convection
  u(1) = (Tamb - x(3)) * area * h;
  x += (Fa * x + G * u) * dt ;
  t += dt;

  delay(1000 * dt); // realtime
}