/*----------------------------------------------------------------------------*/
/*--  Uses the Steinhart-Hart equation to compute thermistor coefficients,  --*/
/*--  then generates MegaSquirt transfer function tables as appropriate.    --*/
/*--  For full documentation on the derivation see thermistor.html.         --*/
/*--                                                                        --*/
/*--  Copyright (c) 2001 by Bruce Bowling.  All rights reserved.            --*/
/*--                                                                        --*/
/*--  This program is free software; you can redistribute it and/or         --*/
/*--  modify it under the terms of the GNU General Public License           --*/
/*--  as published by the Free Software Foundation; either version          --*/
/*--  2 of the License, or (at your option) any later version.              --*/
/*--  See http://www.gnu.org/licenses/gpl.txt                               --*/
/*----------------------------------------------------------------------------*/
/*--  History                                                               --*/
/*--  1.0 - Bruce Bowling's original code.                                  --*/
/*--  2.0 - initial rewrite by Eric Fahlgren <eric.fahlgren@mscsoftware.com>,--*/
/*--        incorporating Steinhart-Hart.                                   --*/
/*--  2.1 - efahl changes to MAT file numbers, get rid of temperature       --*/
/*--        limit; add some useful comments and a version number.           --*/
/*--  2.2 - Add -b and -s switches, use GM defaults, fix help.              --*/
/*----------------------------------------------------------------------------*/

#define VERSION "2.2"

#include <stdlib.h>
#include <math.h>
#include <ctype.h>

typedef struct {
   double t;
   double r;
} dataPoint;

int    userPoints = 0;
double A = 0.00149031826;
double B = 0.00022745952;
double C = 1.1639e-7;

void computeCoefficients(dataPoint p[3])
{
   if (userPoints) {
      double c11, c12, c13;
      double c21, c22, c23;
      double c31, c32, c33;

      c11 = log(p[0].r); c12 = pow(c11, 3.0); c13 = 1.0 / p[0].t;
      c21 = log(p[1].r); c22 = pow(c21, 3.0); c23 = 1.0 / p[1].t;
      c31 = log(p[2].r); c32 = pow(c31, 3.0); c33 = 1.0 / p[2].t;

      C = ((c23-c13) - (c33-c13)*(c21-c11)/(c31-c11)) / ((c22-c12) - (c32-c12)*(c21-c11)/(c31-c11));
      B = (c33-c13 - C*(c32-c12)) / (c31-c11);
      A = c13 - B*c11 - C*c12;
   }
}

double k2f(double t) { return (t * 9.0/5.0) - 459.67; }
double k2c(double t) { return t - 273.15; }
#define f2k(t) ((t+459.67) * 5.0/9.0)

double Tk(double R) /* Kelvins */
{
   return 1.0 / (A + B*log(R) + C*pow(log(R), 3));
}

double Tc(double R) /* Celsius    */ { return k2c(Tk(R)); }
double Tf(double R) /* Fahrenheit */ { return k2f(Tk(R)); }

/*============================================================================*/

#include <stdio.h>

void Usage(char *argv[], int iarg)
{
   fprintf(stderr,
      "usage: %s -l<label> -i -c -b<biasValue>\n"
      "          -r<resistances> -t<temperatures> [-F|-C]\n"
      "          -s<coefficents>\n", argv[0]);
   fprintf(stderr, "   as in\n"
                   "      %s -t0,70,100 -r6000,450,190 -C -i -c\n", argv[0]);
   fprintf(stderr,
      "\n"
      "   -i = Generate the manifold AIRDENFACTOR table in airdenfactor.inc.\n"
      "   -c = Generate the coolant THERMFACTOR table in thermfactor.inc.\n"
      "   -l = Supply an optional label, which is written to comments at the top\n"
      "        of the output files.\n"
      "   -b = Specifies the bias resistor value (R4 for MAT, R7 for CTS).\n"
      "\n"
      "   If you already have the three Steinhart-Hart coefficent values:\n"
      "   -s = Comma separated list of three coefficient values, A, B, C.\n"
      "\n"
      "   If you don't know them, then you must supply three measurements:\n"
      "   -r = Comma separated list of three resistance values.\n"
      "   -t = Comma separated list of three temperature values (default Kelvins).\n"
      "   -F = Above temperatures are in Fahrenheit.\n"
      "   -C = Above temperatures are in Celsius.\n"
   );
   exit(0);
}

void Error(char *argv[], int iarg, char *msg)
{
   fprintf(stderr, "ERROR: %s \"%s\"\n", msg, argv[iarg]);
   fprintf(stderr, "For help: %s -h\n", argv[0]);
   exit(1);
}

/*----------------------------------------------------------------------------*/

typedef enum { ctsOut = 1, matOut = 2, dbgOut = 4, allOut = 7 } outType;

char   *label        = "GM 121 463 12 sensor data.";
int     userLabel    = 0;
FILE   *cts          = NULL;
FILE   *mat          = NULL;
int     dbg          = 1;
double  biasResistor = 2490.0;

dataPoint p[3] = {
   /* GM defaults from original gentherm.c */
   {f2k( -4), 28681},
   {f2k( 86),  2238},
   {f2k(210),   177}
};

#define T(i) (p[i].t)
#define R(i) (p[i].r)

/*----------------------------------------------------------------------------*/

void parseValues(char *string, double Values[], int *count)
{
   *count = 0;

   while (isdigit(*string) || *string == '-') {
      Values[*count] = strtod(string, &string);
      if (*string == ',') string++;
      (*count)++;
   }
}

/*----------------------------------------------------------------------------*/

#define adjustI(n) {argv[iarg] += (n); if (argv[iarg][0] == '\0') ++iarg;}

void parseArgs(int argc, char *argv[])
{
   int    iarg;
   int    i, n;
   double v[3];

   enum { Kelvin, Celsius, Fahrenheit} tempsIn = Kelvin;

   for ( iarg = 1; iarg < argc; iarg++ ) {
      if ( argv[iarg][0] != '-' ) Error(argv, iarg, "Invalid argument.");

      switch ( argv[iarg][1] ) {
         case 'l':
            adjustI(2);
            label = argv[iarg];
            userLabel = 1;
            break;

         case 't':
            adjustI(2);
            parseValues(argv[iarg], v, &n);
            if (n != 3) Error(argv, iarg, "Wrong number of temperature values (should be 3)");
            for (i = 0; i < 3; i++) T(i) = v[i];
            userPoints = 1;
            if (!userLabel) label = "";
            break;

         case 'r':
            adjustI(2);
            parseValues(argv[iarg], v, &n);
            if (n != 3) Error(argv, iarg, "Wrong number of resistance values (should be 3)");
            for (i = 0; i < 3; i++) R(i) = v[i];
            userPoints = 1;
            if (!userLabel) label = "";
            break;

         case 'C':
            if (tempsIn != Kelvin) Error(argv, iarg, "Temperature already set to Fahrenheit.");
            tempsIn = Celsius;
            break;

         case 'F':
            if (tempsIn != Kelvin) Error(argv, iarg, "Temperature already set to Celsius");
            tempsIn = Fahrenheit;
            break;

         case 'i':
            if (mat) Error(argv, iarg, "MAT file already opened.");
            mat = fopen("airdenfactor.inc", "w");
            dbg = mat == NULL;
            break;

         case 'c':
            if (cts) Error(argv, iarg, "CTS file already opened.");
            cts = fopen("thermfactor.inc", "w");
            dbg = cts == NULL;
            break;

         case 'b':
            adjustI(2);
            biasResistor = strtod(argv[iarg], NULL);
            break;

         case 's':
            adjustI(2);
            parseValues(argv[iarg], v, &n);
            if (n != 3) Error(argv, iarg, "Wrong number of Steinhart-Hart coefficients (should be 3)");
            A = v[0];
            B = v[1];
            C = v[2];
            if (!userLabel) label = "";
            break;

         case 'h':
         default : Usage(argv, iarg);
      }
   }

   switch (tempsIn) { /* Convert input temps to Kelvins. */
      case Kelvin:
         break;
      case Celsius:
         for (i = 0; i < 3; i++) T(i) += 273.15;
         break;
      case Fahrenheit:
         for (i = 0; i < 3; i++) T(i) = (T(i) + 459.67) * 5.0/9.0;
         break;
   }
}

/*============================================================================*/

#include <stdarg.h>

void write(int files, char *format, ...)
{
   va_list  args;
   char     outStr[1000];

   va_start(args, format);
      vsprintf(outStr, format, args);
   va_end(args);

   /* assert (strlen (outStr) < 1000); */

   if (files & matOut && mat) fprintf(mat, "%s", outStr);
   if (files & ctsOut && cts) fprintf(cts, "%s", outStr);
   if (files & dbgOut && dbg) printf("%s", outStr);
}

/*============================================================================*/

int main(int argc, char *argv[])
{
   int    adcCount;
   int    i;

   parseArgs(argc, argv);
   computeCoefficients(p);

   write(allOut, "; %s v.%s\n;\n", argv[0], VERSION);
   if (label) write(allOut, "; %s\n", label);
   write(matOut, "; For use with %.0f ohm bias resistor installed at R4.\n", biasResistor);
   write(ctsOut, "; For use with %.0f ohm bias resistor installed at R7.\n", biasResistor);
   write(allOut, ";\n; %s Steinhart-Hart coefficients: A=%g  B=%g  C=%g\n", userPoints ? "Computed" : "Known", A, B, C);
   if (userPoints) {
      write(allOut, ";\n; Input data confirmation (2nd and 3rd columns had better be the same):\n;\n");
      write(allOut, ";    Resistance    Input Temp     Computed\n");
      for (i = 0; i < 3; i++) write(allOut, ";    % 10.2f    % 10.2f   % 10.2f\n", R(i), k2f(T(i)), Tf(R(i)));
   }
   write(allOut, "\n");

   write(matOut, "AIRDENFACTOR:\n\t\t\t; ADC - Temp (dF)\n");
   write(ctsOut, "THERMFACTOR:\n\t\t\t; ADC - Temp (dF)\n");
   write(dbgOut, "ADC -        Res       Temp    Density\n");

   /* Loop over all of the A/D range (0 to 255 for 8 bits). */
   for (adcCount = 0; adcCount < 256; adcCount++) {
      if (adcCount < 1 || adcCount > 254) {
         write(ctsOut, "\tDB\t210T\t; %3d - sensor failure.\n", adcCount);
         write(matOut, "\tDB\t100T\t; %3d - sensor failure.\n", adcCount);
         write(dbgOut, "%3d - sensor failure.\n", adcCount);
      }
      else {
         /* Determine resistance for corresponding adcCount. */
         double res  = biasResistor / (255.0/(double)adcCount - 1.0);
         double temp = Tf(res);
         double ad, adp, tCTS;

         /* Limit range and write temp as deg F plus 40 degree offset. */
         tCTS = temp;
         if (temp > 215.0) tCTS = 215.0;
         if (temp < -40.0) tCTS = -40.0;
         write(ctsOut, "\tDB\t%3dT\t; %3d - % 5.1f\n", (int) (tCTS + 40), adcCount, temp);

         /* Determine air density correction factor and write. */
         ad  = 37.943 / ((temp + 459.7) * 1728.0);
         adp = 100.0 * (ad / 4.157E-5);
         write(matOut, "\tDB\t%3dT\t; %3d - % 5.1f\n", (int) (adp), adcCount, temp);

         /* If not writing files, then spew all sorts of stuff. */
         write(dbgOut, "%3d - %10.2f   % 8.2f   % 8.2f\n", adcCount, res, temp, adp);
      }
   }

   return 0;
}