#include "common.h"

#include <math.h>

// https://www.sciencedirect.com/science/article/abs/pii/S0263224111002594
// https://en.wikipedia.org/wiki/Steinhart%E2%80%93Hart_equation
// https://onlinelibrary.wiley.com/doi/abs/10.1002/andp.18280890511
// https://onlinelibrary.wiley.com/doi/abs/10.1002/andp.18441370202
// https://bmcnoldy.earth.miami.edu/Humidity.html

// Thermistor general constants
#define VIN   5.0f    // Source voltage (VCC)
#define GAIN  6.144f  // PGA gain setting 0x00
#define AMAX  0x7FFF  // 16-bit ADC resolution
#define R1    10000   // 10k Ohm resistor

// Thermistor coeficients
#define SC1   1.009249522e-3f
#define SC2   2.378405444e-4f
#define SC3   2.019202697e-7f

// Conversion constants
#define TH1   100.0f
#define TH2   243.04f
#define TH3   17.625f

float Resistance(int adc_raw)
{
	float av, r2;

	// [GND] - [10K R] - [ADC IN] - [THERMISTOR] - [VCC]

	// The ADC measurement is taken between a 10k Ohm
	// resistor and the thermistor. This is also called
	// a voltage divider circuit which has the following
	// equation:

	// Vout = Vin * (R2 / R1 + R2)

	// We can solve for resistance R2 by rearranging and
	// simplifying the equation:

	// R2 = R1 / (Vin / Vout - 1)

	// If the ADC reference voltage and voltage divider
	// source voltage (Vin) are the same, then the
	// following is true:

	// Vin / Vout = AMAX / AV

	// This means the ratio of voltage divider input
	// voltage to output voltage is the same as the ratio
	// of the ADC full range value (AMAX) to the value
	// returned by the ADC (AV). For a 16 bit ADC like
	// the ADS1115 the AMAX is 0x7FFF.

	av = adc_raw * (GAIN / VIN);
	r2 = R1 / (AMAX / av - 1.0f);

	return r2;
}

float SteinhartHart(float res)
{
	float logr, t;

	// Convert resistance to temperature

	// When the temperature increases, NTC thermistor
	// resistance will decrease:

	//  25C = 10000 Ohms
	// 100C =  6744 Ohms

	logr = log(res);
	t = 1.0f / (SC1 + SC2 * logr + SC3 * pow(logr, 3));

	return t - 273.15f; // Kelvin to celsius
}

float Dewpoint(float t, float rh)
{
	float a, b;

	// 243.04*(LN(RH/100)+((17.625*T)/(243.04+T)))/
	// (17.625-LN(RH/100)-((17.625*T)/(243.04+T)))

	a = log(rh / TH1);
	b = TH3 * t / (TH2 + t);

	return TH2 * (a + b) / (TH3 - a - b);
}

float Temperature(float td, float rh)
{
	float a, b;

	// 243.04*(((17.625*TD)/(243.04+TD))-LN(RH/100))/
	// (17.625+LN(RH/100)-((17.625*TD)/(243.04+TD)))

	a = log(rh / TH1);
	b = TH3 * td / (TH2 + td);

	return TH2 * (b - a) / (TH3 + a - b);
}

float RelHumidity(float t, float td)
{
	float a, b;

	// 100*(EXP((17.625*TD)/(243.04+TD))/
	//      EXP((17.625*T)/(243.04+T)))

	a = TH3 * td / (TH2 + td);
	b = TH3 * t / (TH2 + t);

	return TH1 * (exp(a) / exp(b));
}