Deriving Steinhart-Hart (http://en.wikipedia.org/wiki/Steinhart–Hart_equation) coefficients given a known thermistor response curve.

steinhart-test-results.txt

CL-USER> (test-steinhart-hart)          
Steinhart-hart coefficients for the ACI/10K-CP curve are:
A: 0.0011212672
B: 2.3534849E-4
C: 8.3802405E-8

Resistance: 336450.0 Expected: -40.0 Calculated: -39.999985
Resistance: 242660.0 Expected: -35.0 Calculated: -35.009888
Resistance: 176960.0 Expected: -30.0 Calculated: -30.018707
Resistance: 130410.0 Expected: -25.0 Calculated: -25.02591
Resistance: 97072.0 Expected: -20.0 Calculated: -20.031494
Resistance: 72951.0 Expected: -15.0 Calculated: -15.0355835
Resistance: 55326.0 Expected: -10.0 Calculated: -10.037994
Resistance: 42326.0 Expected: -5.0 Calculated: -5.0384216
Resistance: 32650.0 Expected: 0.0 Calculated: -0.036224365
Resistance: 25391.0 Expected: 5.0 Calculated: 4.9650574
Resistance: 19889.0 Expected: 10.0 Calculated: 9.976471
Resistance: 15711.0 Expected: 15.0 Calculated: 14.966095
Resistance: 12492.0 Expected: 20.0 Calculated: 19.966736
Resistance: 10000.0 Expected: 25.0 Calculated: 24.967407
Resistance: 8057.0 Expected: 30.0 Calculated: 29.96875
Resistance: 6531.0 Expected: 35.0 Calculated: 34.974
Resistance: 5326.0 Expected: 40.0 Calculated: 39.97827
Resistance: 4368.0 Expected: 45.0 Calculated: 44.984497
Resistance: 3602.0 Expected: 50.0 Calculated: 49.99167
Resistance: 2986.0 Expected: 55.0 Calculated: 55.00003
Resistance: 2488.0 Expected: 60.0 Calculated: 60.00833
Resistance: 2083.0 Expected: 65.0 Calculated: 65.02008
Resistance: 1752.0 Expected: 70.0 Calculated: 70.03491
Resistance: 1479.0 Expected: 75.0 Calculated: 75.07663
Resistance: 1255.0 Expected: 80.0 Calculated: 80.09567
Resistance: 1070.0 Expected: 85.0 Calculated: 85.09799
Resistance: 915.4 Expected: 90.0 Calculated: 90.120514
Resistance: 786.6 Expected: 95.0 Calculated: 95.127045
Resistance: 678.6 Expected: 100.0 Calculated: 100.12717
Resistance: 587.6 Expected: 105.0 Calculated: 105.12454
Resistance: 510.6 Expected: 110.0 Calculated: 110.12079
Resistance: 445.2 Expected: 115.0 Calculated: 115.11606
Resistance: 389.6 Expected: 120.0 Calculated: 120.0961
Resistance: 341.9 Expected: 125.0 Calculated: 125.08865
Resistance: 301.0 Expected: 130.0 Calculated: 130.07501
Resistance: 265.8 Expected: 135.0 Calculated: 135.05676
Resistance: 235.4 Expected: 140.0 Calculated: 140.0351
Resistance: 209.0 Expected: 145.0 Calculated: 145.02298
Resistance: 186.1 Expected: 150.0 Calculated: 150.0

steinhart.lisp

;; Figure out the Steinhart-Hart coefficients to use for a given thermistor curve.

; The ACI/10K-CP thermistor curve from here:
; http://workaci.com/sites/default/files/additional_information/thertempcurve.pdf
(defparameter +ACI/10K-CP+
  '((336450.0 -40.0)
    (242660.0 -35.0)
    (176960.0 -30.0)
    (130410.0 -25.0)
    (97072.0 -20.0)
    (72951.0 -15.0)
    (55326.0 -10.0)
    (42326.0 -5.0)
    (32650.0 0.0)
    (25391.0 5.0)
    (19889.0 10.0)
    (15711.0 15.0)
    (12492.0 20.0)
    (10000.0 25.0)
    (8057.0 30.0)
    (6531.0 35.0)
    (5326.0 40.0)
    (4368.0 45.0)
    (3602.0 50.0)
    (2986.0 55.0)
    (2488.0 60.0)
    (2083.0 65.0)
    (1752.0 70.0)
    (1479.0 75.0)
    (1255.0 80.0)
    (1070.0 85.0)
    (915.4 90.0)
    (786.6 95.0)
    (678.6 100.0)
    (587.6 105.0)
    (510.6 110.0)
    (445.2 115.0)
    (389.6 120.0)
    (341.9 125.0)
    (301.0 130.0)
    (265.8 135.0)
    (235.4 140.0)
    (209.0 145.0)
    (186.1 150.0)))

(defun celcius->kelvin (C)
  (+ C 273.15))

(defun kelvin->celcius (K)
  (- K 273.15))

(defun find-steinhart-hart-coefficients (R1 C1 R2 C2 R3 C3)
  "Given 3 pairs of resistance/temperature (in degrees Celcius), derive the Steinhart-Hart
coefficients that will generate the curve used by a thermistor. Derived from Wikipedia:
http://en.wikipedia.org/wiki/Steinhart–Hart_equation#Steinhart.E2.80.93Hart_coefficients"
  (let* ((K1 (celcius->kelvin C1))
         (K2 (celcius->kelvin C2))
         (K3 (celcius->kelvin C3))
         (L1 (log R1))
         (L2 (log R2))
         (L3 (log R3))
         (Y1 (/ 1 K1))
         (Y2 (/ 1 K2))
         (Y3 (/ 1 K3))
         (gamma2 (/ (- Y2 Y1) (- L2 L1)))
         (gamma3 (/ (- Y3 Y1) (- L3 L1)))
         (C (* (/ (- gamma3 gamma2)
                  (- L3 L2))
               (expt (+ L1 L2 L3) -1)))
         (B (- gamma2 (* C (+ (expt L1 2)
                              (* L1 L2)
                              (expt L2 2)))))
         (A (- Y1 (* L1 (+ B (* (expt L1 2) C))))))
    (list A B C)))

(defun steinhart-hart-temperature (A B C R)
  "Return the temperature (in degrees Celcius) given by the Steinhart-Hart curve given three
coefficients and a resistance value."
  (let* ((lnR (log R))
         (lnR3 (expt lnR 3))
         (degreesK (+ A (* B lnR) (* C lnR3))))
    (kelvin->celcius (expt degreesK -1))))

(defun inverse-steinhart-hart (A B C degreesC)
  "Given three Steinhart-Hart coefficients and a temperature in degrees Celcius, return
the resistance predicted by the Steinhart-Hart equation."
  (let* ((degreesK (celcius->kelvin degreesC))
         (y (/ (- A (/ 1 degreesK))
               (* 2 C)))
         (x (sqrt (+ (expt (/ B (* 2 C)) 3) (expt y 2)))))
    (exp (- (expt (- x y) 1/3)
            (expt (+ x y) 1/3)))))

(defun coefficients-from-table (table)
  "Given a table like the +ACI/10K-CP+ one defined at the top of this file, return
the Steinhart-Hart coefficients that will generate the curve that table describes."
  (let* ((len (length table))
         (rt1 (first table))
         (rt2 (elt table (floor (/ len 2))))
         (rt3 (car (last table)))
         (results (find-steinhart-hart-coefficients (first rt1) (second rt1)
                                                    (first rt2) (second rt2)
                                                    (first rt3) (second rt3))))
    (values (first results) (second results) (third results))))

(defun test-steinhart-hart ()
  (multiple-value-bind (A B C) (coefficients-from-table +ACI/10K-CP+)
    (format t "Steinhart-hart coefficients for the ACI/10K-CP curve are:~%")
    (format t "A: ~a~%B: ~a~%C: ~a~%~%" A B C)
    (loop for rt-pair in +ACI/10K-CP+ do
         (format t "Resistance: ~a Expected: ~a Calculated: ~a~%"
                 (first rt-pair)
                 (second rt-pair)
                 (steinhart-hart-temperature A B C (first rt-pair))))))