#!/usr/bin/env python

#############################################################################
## 
#  \file ntc.py
#  \brief Interpolation table calculation for NTC sensors
#  \version 0.1
#  \author Dimitri Denk
#  \date 14.11.2014
#
# This file implements the NTC sensor lookup table (LUT) calculation for
# fast interpolation algorithm with output in fix7_8 format.
# As measurement circuit a resistive voltage divider with NTC sensor
# referenced to ground is used.
#
# Copyright (c) 2014 
# Dimitri Denk
#
# Permission to use, copy, modify, distribute and sell this software
# and its documentation for any purpose is hereby granted without fee,
# provided that the above copyright notice appear in all copies and
# that both that copyright notice and this permission notice appear
# in supporting documentation. Author makes no representations about 
# the suitability of this software for any purpose. 
# It is provided "as is" without express or implied warranty.
#
# The fast interpolation algorithm is based on linear interpolation equation:
# f(x) = f0 + (f1 - f0) / (x1 - x0)
#
# It uses of aligned lookup table, which allows easy determinate search area
# and replace division by arithmetic binary shift operation:
# f(x) = LUT[i] + (LUT[i + 1] - LUT[i]) * (x - (i << (N - n))) >> (N - n))
# where:
# N - Resolution of input value in bits
# n - Number of bits allocated in lut entry
# i = x >> (N - n) - Index of actual element in lookup table
# The necessary size of lookup table can be calculated as (1 << n) + 1
#
# This implementation minimizes of calculation time and code size.
#
# The calculation script transforms known temperature values into
# n-bits aligned lookup table and generates C structure, which can be
# directly integrated into C code.
#
# Script parameters:
# R - Upper resistor value
# N - ADC resolution in bits
# n - Width of LUT in bits
# ADCref - ADC reference voltage
# NTCref - resistive divisor reference voltage
# Rntc - Temperature and resistance values from datasheet must be placed in tuple or list:
# ((temperature 1, resistance 1), (temperature 2, resistance 2),  ...)
#
#############################################################################

import numpy as np
import scipy.interpolate as inter
import sys

# ADC resolution
N = 12
# width of LUT in bits
n = 4
# ADC reference voltage
ADCref = 3.3
# NTC sensor reference voltage
NTCref = 3.3
# upper resistor value in kOhm
R = 10.
# Resistance values in kOhm an NTCLE100E3103 from datasheet
Rntc = (
	(-40, 332.094),
	(-35, 239.900),
	(-30, 175.200),
	(-25, 129.287),
	(-20, 96.358),
	(-15, 72.500),
	(-10, 55.046),
	(-5, 42.157),
	(0, 32.554),
	(5, 25.339),
	(10, 19.872),
	(15, 15.698),
	(20, 12.488),
	(25, 10.0),
	(30, 8.059),
	(35, 6.535),
	(40, 5.330),
	(45, 4.372),
	(50, 3.605),
	(55, 2.989),
	(60, 2.490),
	(65, 2.084),
	(70, 1.753),
	(75, 1.481),
	(80, 1.256),
	(85, 1.070),
	(90, 0.9154),
	(95, 0.786),
	(100, 0.6773),
	(105, 0.5857),
	(110, 0.5083),
	(115, 0.4426),
	(120, 0.3866),
	(125, 0.3387),
)

# resistance to temperature conversion function
Ty, Rx = np.transpose(Rntc)
T = inter.InterpolatedUnivariateSpline(Rx[::-1], Ty[::-1], k = 3)

# fast interpolation function based on aligned lookup table
def fast_interpolate(v, lut, N, n):
	i = (v >> (N - n))
	return lut[i] + (lut[i + 1] - lut[i]) * (v - (i << (N - n))) / (1 << (N - n))

# calculating of reference points for LUT
X = np.linspace(0, (1 << N), 2**n + 1)
# array for temperature values
Y = []
# calcule and generate C code
print "static const signed short lut[" + str(2**n + 1) + "] = {"
for cadc in X:
	# voltage on ADC
	vadc = (cadc * ADCref) / (1 << N)
	# check if value is valid
	if vadc < NTCref:
		# calculate NTC resistance
		r = (R * vadc) / (NTCref - vadc)
	else:
		# open circuit
		r = R * 1000
	# calculate temperature form resistance value
	t = T(r)
	# calculate a fixed point value
	fixt = int(t * 256)
	# range check
	if fixt > 32767:
		fixt = 32767
		t = float(32767) / 256
		print "\t" + str(fixt) + ",\t\t// Upper limit (temperature value is not correct), short circuit\tADC value:" + str(int(cadc))
	elif fixt < -32768:
		fixt = -32768
		t = float(-32768) / 256
		print "\t" + str(fixt) + ",\t\t// Lower limit (temperature value is not correct), open circuit\tADC value:" + str(int(cadc))
	else:
		print "\t" + str(fixt) + ",\t\t// T:" + str(round(t, 2)) + "'C,\tRt:" + str(round(r, 3)) + "k\tADC value:" + str(int(cadc))
	Y.append(t)
print "};"
print ""
print "static inline signed short ntcvalue(unsigned short adcval) {"
print "	unsigned short i = (adcval >> " + str(N - n) + ");"
print "	return lut[i] + (signed short)((signed long)(lut[i + 1] - lut[i]) * (adcval - (i << " + str(N - n) + ")) >> " + str(N - n) + ");"
print "}"

# comment this line if graphical plot of results is required
sys.exit(0)

# visualization of results
import matplotlib.pyplot as plt
import matplotlib.cm
import matplotlib.mlab as mlab
from matplotlib.figure import SubplotParams
# ADC code
U = np.arange(1, 2**N, 1)
# Temperature values
T1 = []
T2 = []
for u in U:
	T1.append(T(u * R / (2**N - u)))
	T2.append(fast_interpolate(u, Y, N, n))
# plot
sp = SubplotParams(left=0.15, bottom=0.1, right=0.9, top=0.92, wspace=0., hspace=0.)
fig = plt.figure(frameon=False, subplotpars=sp, figsize=(7., 5.))
ax1 = fig.add_subplot(1, 1, 1)

ax1.plot(U, T1, 'r', label='Temperature response')
ax1.plot(X, Y, 'bo', label='LUT points')
ax1.plot(U, T2, 'g', label='Interpolated values')
#ax1.plot(U, (np.array(T2) - np.array(T1)) * 100, 'k', label='Error * 100')

ax1.set_title("NTCLE100E3103 interpolation, N = " + str(N) + ", n = " + str(n))
ax1.grid(True)
ax1.set_ylim(-128, 127)
ax1.set_xlim(0., 1 << N)
ax1.set_xlabel("ADC code")
ax1.set_ylabel("Temperature 'C")
ax1.legend()

fig.savefig("NTCLE100E3103_" + str(N) + "_" + str(n) + ".png", format='png', dpi=75)
plt.show()