/*----------------------------------------------------------------------------*/
/*--  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 Eric Fahlgren.  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 - Stole gentherm2 and genmap, mashed them into one program.       --*/
/*----------------------------------------------------------------------------*/

#define VERSION "1.0"

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

/*============================================================================*/
/*==  Begin Thermistor Code  =================================================*/

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

bool    userLabel  = false;
char   *tLabel     = "GM 121 463 12 sensor data.";
bool    userPoints = false;
double  A          = 0.00149031826;
double  B          = 0.00022745952;
double  C          = 1.1639e-7;

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

bool   doTherm       = false;
bool   doMap         = false;

FILE   *cts          = NULL;
FILE   *mat          = NULL;
bool    dbg          = false;
double  biasResistor = 2490.0;

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)

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

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

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 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)); }

/*==  End Thermistor Code  ===================================================*/
/*============================================================================*/
/*==  Begin MAP Code  ========================================================*/


typedef struct {
   /*--  Max number of points is 100, more than enough.  I'd like to use    --*/
   /*--  incomplete types, but can't in C.                                  --*/
   char   *pLabel;
   int     nPts;      // Length of the volt and pres arrays.
   double  volt[100]; // Voltages (0-5) corresponding to below pressures.
   double  pres[100]; // Pressure array must be same length as volt.
   int     pLo;       // Pressure limits for this sensor.
   int     pHi;
} MapData;

static MapData MPX4115 = { // Motorola MPX4115
   "; Motorola MPX4115A",
   2,
   {  0.204,   4.794 },
   { 15.000, 115.000 },
    11,
   118
};

static MapData MPX4250 = { // Motorola MPX4250
   "; Motorola MPX4250A",
   2,
   {  0.204,   4.896 },
   { 20.000, 250.000 },
    11,
   253
};

static MapData Bosch300 = {
   "ERROR, need to scale output to non-kPa because of values > 255!!!"
   "; Bosch 300 kPa MAP sensor, data from Robin Johansson",
   2,
   {  0.25,   4.8 },
   { 20.0 , 300.0 },
    15,
   310
};

static MapData GMold = {
   "; GM 100 kPa MAP sensor, data from Tom Dolsky",
   2,
   {  0.425,    4.025 },
   { 20.000,  100.000 },
    15, // Made up.
   105
};

static MapData GMnew = {
   "; GM 150? kPa MAP sensor, data from Tom Dolsky",
   4,
   { 4.46,   3.20,   2.12,    1.06  },
   { 6.90,  55.158, 96.527, 137.896 },
    15, // Made up.
   150
};

MapData *md = &MPX4250;

/* end user defined area here */
/******************************/

bool outOfRange(int pressure)
{
   return pressure < md->pLo || pressure > md->pHi;
}

void genMap()
{
   double   voltage;
   int      pressure, barocor;
   int      i, adcCount;
   FILE    *kpa;
   FILE    *bcf;
   bool     out;

   kpa = fopen("kpafactor.inc", "w");
   fprintf(kpa, "%s\n", md->pLabel);
   fprintf(kpa, "KPAFACTOR: ;\tKPA\t  ADC\tVolts (* = out of sensor range)\n");

   bcf = fopen("barofactor.inc", "w");
   fprintf(bcf, "%s\n", md->pLabel);
   fprintf(bcf, "BAROFAC: ;\tcorrFac\t  ADC\tVolts (* = out of sensor range)\n");

   /* Loop over all of the A/D range (0 to 255 for 8 bits) */
   i = 1; /* Since they are sorted by ascending values, we can do this once outside the loop. */
   for (adcCount = 0; adcCount < 256; adcCount++) {
      /* determine voltage for corresponding adcCount */
      voltage = (double)adcCount * 5.0 / 255.0;

      /* Find where voltage value lies in volt[] array */
      while (i < (md->nPts - 1) && voltage > md->volt[i]) i++;

      /* Linear interpolate */
      pressure = md->pres[i] + ((md->pres[i] - md->pres[i - 1]) * (voltage - md->volt[i])) / (md->volt[i] - md->volt[i - 1]);

      /* Limit range to specified sensor values and write */
      out = outOfRange(pressure);
      if (out) pressure = 100.0;
      fprintf(kpa, "\tDB\t%3dT\t; %3d - %5.3f%s\n", (int) (pressure), adcCount, voltage, out?"*":"");

      /* Generate barometer correction factor based on pressure value */
      barocor = 100.0 * (1.0 - 0.0047 * ((double)pressure - 100.0));
      fprintf(bcf, "\tDB\t%3dT\t; %3d - %5.3f%s\n", (int) (barocor), adcCount, voltage, out?"*":"");
   }

   fclose(kpa);
   fclose(bcf);
}

/*==  End MAP Code  ==========================================================*/
/*============================================================================*/
/*==  Argument Parser  =======================================================*/

char *progName = NULL;

void Usage()
{
   fprintf(stderr,
      "usage: %s -l<label> -mat -cts -map -b<biasValue>\n"
      "          -r<resistances> -t<temperatures> [-F|-C]\n"
      "          -s<coefficents>\n", progName);
   fprintf(stderr, "   as in\n"
                   "      %s -cts -t0,70,100 -r6000,450,190 -C\n", progName);
   fprintf(stderr,
      "\n"
      "For temperature sensor table generation, the following arguments apply:\n"
      "   -mat = Generate the manifold AIRDENFACTOR table in airdenfactor.inc.\n"
      "   -cts = Generate the coolant THERMFACTOR table in thermfactor.inc.\n"
      "   -dbg = generate some debug output in addition to any files.\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 the above, then you must supply three points:\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"
      "\n"
      "For pressure sensor table generation use the following parameters:\n"
      "   -map = generates both the kpafactor.inc and barofactor.inc files.\n"
      "   -Mn  = use canned table n, n=1 MPX4250, n=2 MPX4115, 3-5 try 'em.\n"
      "   -v   = Comma separated list of voltages, with fewer than 100 values.\n"
      "   -p   = Comma separated list of pressures (default units kPa).\n"
      "   -P   = above pressure values are in psig.\n"
      "   -B   = above pressure values are in BAR.\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", progName);
   exit(1);
}

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

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;

   progName = argv[0];

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

      switch (argv[iarg][1]) {
         case 'd':
            if (strncmp(argv[iarg]+2, "bg", 2) != 0)
               Error(argv, iarg, "Invalid argument.");
            else {
               dbg     = true;
               doTherm = true;
            }
            break;

         case 'l':
            adjustI(2);
            tLabel     = argv[iarg];
            md->pLabel = argv[iarg];
            userLabel  = true;
            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 = true;
            if (!userLabel) tLabel = "";
            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 = true;
            if (!userLabel) tLabel = "";
            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 'm':
            if (strncmp(argv[iarg]+2, "at", 2) == 0) {
               if (mat) Error(argv, iarg, "MAT file already opened.");
               mat = fopen("airdenfactor.inc", "w");
               doTherm = true;
            }
            else if (strncmp(argv[iarg]+2, "ap", 2) == 0) {
               doMap = true;
            }
            else
               Error(argv, iarg, "Invalid argument.");
            break;

         case 'c':
            if (strncmp(argv[iarg]+2, "ts", 2) != 0)
               Error(argv, iarg, "Invalid argument.");
            else {
               if (cts) Error(argv, iarg, "CTS file already opened.");
               cts = fopen("thermfactor.inc", "w");
               doTherm = true;
            }
            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) tLabel = "";
            break;

         case 'M':
            adjustI(2);
            i = strtol(argv[iarg], NULL, 10);
            switch(i) {
               case 1: md = &MPX4250;  break;
               case 2: md = &MPX4115;  break;
               case 3: md = &Bosch300; break;
               case 4: md = &GMold;    break;
               case 5: md = &GMnew;    break;
            }
            break;


         case 'h':
         default : Usage();
      }
   }

   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);
}

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

void genTherm()
{
   int    adcCount;
   int    i;

   computeCoefficients(p);

   write(allOut, "; %s v.%s\n;\n", progName, VERSION);
   if (tLabel) write(allOut, "; %s\n", tLabel);
   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);
      }
   }
}

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

int main(int argc, char *argv[])
{
   parseArgs(argc, argv);
   if (doTherm) genTherm();
   if (doMap  ) genMap  ();
}