// http://www.sourcecodesworld.com/source/show.asp?ScriptID=1086

// aug 2024 - Jens Dalsgaard Nielsen
// code for calculating Steinhart-Hart 3rd order approx to NTC
// by three sets of temp(C) and resistance(Ohm) measurements
//
// compile by gcc  gcc findabc.c  -lm
// -lm links math library bq we do use nat log function
//
#include<stdio.h>
#include<stdlib.h>
#include <math.h>

#define matsize 3

float A[matsize][matsize];
float I[matsize][matsize];
float T[matsize];

float R1 = 25415.0,
      R2 = 10021.0,
      R3 = 6545.0,
      T1 = 5.0,
      T2 = 25.0,
      T3 = 35.0;


float ABC[3];

float temp;

void
fillMat ()
{
    A[0][0] = 1.0;
    A[1][0] = 1.0;
    A[2][0] = 1.0;
    A[0][1] = log (R1);  // natural log (ln) 2.71281828
    A[1][1] = log (R2);
    A[2][1] = log (R3);
    A[0][2] = pow (A[0][1], 3.0);
    A[1][2] = pow (A[1][1], 3.0);
    A[2][2] = pow (A[2][1], 3.0);
}

void
invT ()
{
    T[0] = 1.0 / (T1 + 273.15);
    T[1] = 1.0 / (T2 + 273.15);
    T[2] = 1.0 / (T3 + 273.15);
}

void
invA ()
{
    int i, j, k;
    // fill unit matr
    for (i = 0; i < matsize; i++) {
        for (j = 0; j < matsize; j++) {
            if (i == j) {
                I[i][j] = 1;
            } else {
                I[i][j] = 0;
            }
        }
    }
    for (k = 0; k < matsize; k++) {
        temp = A[k][k];
        for (j = 0; j < matsize; j++) {
            A[k][j] /= temp;
            I[k][j] /= temp;

        }
        for (i = 0; i < matsize; i++) {
            temp = A[i][k];
            for (j = 0; j < matsize; j++) {
                if (i == k) {
                    break;
                }
                A[i][j] -= A[k][j] * temp;
                I[i][j] -= I[k][j] * temp;
            }
        }
    }
}

void
calcABC ()
{
    for (int i = 0; i < 3; i++) {
        ABC[i] = 0.0;
        for (int j = 0; j < 3; j++) {
            ABC[i] += I[i][j] * T[j];
        }
    }
}

int
main ()
{
    int i, j, k;

    fillMat ();

    printf ("\n A matrix\n");
    for (i = 0; i < 3; i++) {
        for (j = 0; j < 3; j++) {
            printf ("%.10e	", A[i][j]);
        }
        printf ("\n");
    }

    invT ();
    invA ();
    calcABC ();

    printf ("\n A matrix as unit\n");
    for (i = 0; i < 3; i++) {
        for (j = 0; j < 3; j++) {
            printf ("%.10e	", A[i][j]);
        }
        printf ("\n");
    }

    printf ("\ninv A\n");
    for (i = 0; i < matsize; i++) {
        for (j = 0; j < matsize; j++) {
            printf ("%f	", I[i][j]);
        }
        printf ("\n");
    }

    printf ("\n A B C are\n");
    printf ("\n %.10e %.10e %.10e", ABC[0], ABC[1], ABC[2]);
    printf ("\n---\n");
    return 0;
}