RISC-V/V203F6P6/termistor/compute/termistor.py

#!/usr/bin/python3
# -*- coding: utf-8 -*-
import os
import cffi
import numpy as np
import matplotlib.pyplot as plt
from math import exp,log

c_header = '''double Compute (double x);
void plot ();
'''
def print_table (tab):
  return
  K  = 4095.0
  R0 = 10.0
  s  = 'static constexpr Pair measured [] = {\n'
  l  = len(tab)
  r  = range (0, l)
  for n in r:
    m  = tab [l - n - 1]
    R  = m[1]
    u  = K * R / (R + R0)
    s += '  {{ {:+6.1f},{:8.3f} }},\n'.format(u, m[0])
  s += '};\n'
  print (s)

def adc (R):
  K0 = 4095.0
  R0 = 10.0
  return K0 * R / (R + R0)

def R_NTC_T (t):
  T  = 273.15 + t
  T0 = 298.15
  R0 = 10000.0;
  B  = 3977.0
  R  = R0 * exp (B * (1.0 / T - 1.0 / T0))
  return 0.001 * R

def beta_diff (table):
  print ('t[°C] : R comp [kΩ] | R table [kΩ] | odchylka [%]')
  for m in table:
    t = m[0]
    R = R_NTC_T (t)
    D = 100.0 * (R - m[1]) / R
    print('{:+5g} : {:8.3f}    |  {:8.3f}    |   {:+.2g}%'.format(t, R, m[1], D))

def ComuputeStHaCoeff (table):
  e = table[ 5]; a = log (1000.0 * e[1]); f = 1.0 / (e[0] + 273.15) # 0 °C
  e = table[ 8]; b = log (1000.0 * e[1]); g = 1.0 / (e[0] + 273.15) # 25°C
  e = table[11]; c = log (1000.0 * e[1]); h = 1.0 / (e[0] + 273.15) # 50°C
  e = table[13]; d = log (1000.0 * e[1]); i = 1.0 / (e[0] + 273.15) # 70°C
  M = np.matrix ([[1.0, a, a**2, a**3],
                  [1.0, b, b**2, b**3],
                  [1.0, c, c**2, c**3],
                  [1.0, d, d**2, d**3]])
  I = M.getI(); # inverze matice
  V = np.matrix ([f, g, h, i])
  X = I * V.transpose()
  letters = ['A','B','C','D']; n=0;
  X = X.tolist()
  R = []
  print('Steinhart-Hart koeficienty :')
  for x in X:
    print('  {:s} = {:+20.14e}'.format(letters[n], x[0]))
    n += 1
    R.append (x[0])
  return R

def SteinhartHart (r, coe):
  R  = 1000.0 * r
  A  = coe[0]; B = coe[1]; C = coe[2]; D = coe[3];
  L  = log(R); K  = L;  # Ku podivu koeficienty vypočítané ze 4 bodů tabulky
  W  = A                # Steinhart-Hart metodou dávají lepší výsledky než
  W += B * K;  K *= L;  # výrobcem udávaný koeficient beta.
  W += C * K;  K *= L;
  W += D * K
  return 1.0 / W - 273.15   # [°C]
  
def stha_diff (table, c):
  print ('rozdíl teplot vypočtený podle Steinhart-Hart formule :')
  for m in table:
    t = SteinhartHart(m[1], c)
    print ('{:+5g} => {:g}'.format(m[0], t))
def SHCompute (u, c):
  K = 4095.0; R0 = 10.0
  R = u * R0 / (K - u)
  return SteinhartHart(R, c)
def graph_err (tab, C,c):
  #return
  xs = np.arange(134.0, 4010.0, 1.0); ys = []
  for x in xs: ys.append (SHCompute(x,c) - C.Compute(x))
  xc = []; yc = []
  for m in tab:
    x = adc (m[1])
    xc.append (x)
    yc.append (SHCompute(x,c) - m[0])
  fig,ap = plt.subplots (1, figsize=(6.0,4.0))
  fig.suptitle ('Chyba NTC', fontsize=15)
  ap.plot (xs, ys)
  ap.plot (xc, yc, '.')
  ap.legend (['Odchylka od vzorce', 'Body tabulky'])
  plt.grid()
  plt.xlabel('ADC [0-4095]')
  plt.ylabel('teplota [°C]')
  plt.savefig('err.png')
  #plt.show()
def graph(tab, c):
  #return
  xs = [];  ys = [];
  for m in tab:
    xs.append (adc (m[1]))
    ys.append (m[0])
  xc = np.arange(100.0, 4050.0, 10.0);  yc = [];
  for x in xc: yc.append(SHCompute(x,c))
  fig,ap = plt.subplots (1, figsize=(6.0,4.0))
  fig.suptitle ('NTC 10k', fontsize=15)
  ap.plot (xs, ys, linewidth=3)
  ap.plot (xc, yc)
  ap.legend (['Tabulka výrobce','Steinhart-Hart vzorce'])
  plt.grid()
  plt.xlabel('ADC [0-4095]')
  plt.ylabel('teplota [°C]')
  plt.savefig('ntc.png')
  #plt.show()

def main_func():
  table = [(-40,334.202),(-38,293.759),(-30,177.797),(-25,131.380),(-20,97.923),(-15,73.612),(-10,55.805),
           (-5,42.655),(0,32.869),(5,25.527),(10,19.977),(15,15.750),(20,12.507),(25,10.0),(30,8.049),
           (35,6.521),(40,5.315),(45,4.358),(50,3.594),(55,2.980),(60,2.483),(65,2.080),(70,1.750),(80,1.256),
           (90,0.916),(100,0.678),(110,0.509),(120,0.386),(125,0.338)]
  print_table (table)
  c = ComuputeStHaCoeff(table)
  #stha_diff   (table, c)
  ffi = cffi.FFI()
  ffi.cdef(c_header)
  C = ffi.dlopen("./spline.so")
  C.plot()
  graph_err   (table, C,c)
  ffi.dlclose (C)
  graph       (table, c)
  return
  beta_diff   (table)

if __name__ == '__main__':
  main_func()
  print ("END OK")