function thermistorResistanceCalibration(n,datadump)
% thermistorResistanceCalibration  Curve fit thermistor calibration data T = f(R)
%                    where R is the resistance of the thermistor
%
% Synopsis:  thermistorResistanceCalibration
%            thermistorResistanceCalibration(n)
%            thermistorResistanceCalibration(n,datadump)
%
% Input:  n = optional degree of polynomial to attempt for the curve fit
%         datadump = flag (true or false) determining whether the T = f(R)
%                    data is written to a text file.
%                    Default:  datadump = false

if nargin<1, n=3;             end
if nargin<2, datadump=false;  end

% --- Store data for thermistor 1 and perform the curve fit
%     Data taken at home with Fluke 77 multimeter.  R is resistance in Ohm
%     Reference temperature measured with CDN Pro Quick-read digital thermometer
R = [ 4.99 3.60 21.64 2.767 17.40 4.69 11.90 6.43  9.17  7.87  7.50 10.80 5.89]*1000;
T = [43.4 52.9   5.9 60.7   10.9 45.2  20.2  36.3 26.95 30.95 32.2  22.8 38.7];

f = figure;   set(f,'Name','Thermistor 1');
c1 = thermistorFit(R(:),T(:),n);   %  R(:) and T(:) guarantees that data sent
                                   %  to thermistorFit are column vectors
% --- Evaluate the curve fit
Rfit = linspace( min(R), max(R) );
Tfit1 = (1./ (c1(1) + c1(2)*log(Rfit) + c1(3)*(log(Rfit)).^3 ) ) - 273.15;

% --- Plot the two curve fits on the same axes
f = figure;  set(f,'Name','Thermistor curve fits');
plot(R/1000,T,'o',Rfit/1000,Tfit1,'r:');
xlabel('Resistance (k\Omega)');  ylabel('T (^{\circ}C)');
legend('Thermistor 1 data','SHH fit for thermistor 1');

R10 = 10e3;
TR10 = 1./ (c1(1) + c1(2)*log(R10) + c1(3)*(log(R10)).^3 ) - 273.15;
hold('on');
plot(R10*[1 1]/1000,[0 TR10],'k:',[0 R10]/1000,TR10*[1 1],'k:');
hold('off');



% --- Inverse calibration:  R = f(T)
ni = 2;
ci = polyfit(1./(T+273.15),log(R),ni);
Tifit = linspace(min(T),max(T));
Rifit = exp( polyval(ci,1./(Tifit+273.15)) );
fprintf('\nInverse fit coefficients:\n');
fprintf('%18.7e',ci);  fprintf('\n');

figure('Name','Inverse fit');
plot(T,R,'o',Tifit,Rifit,'r--');
xlabel('T (C)');  ylabel('R  (k\Omega)');
legend('Thermistor 1 data','Thermistor 1 fit');

% --- Save the raw data in a text file
if ~datadump, return;  end

fout = fopen('thermistor_RT_data.txt','wt');
for i=1:length(R)
  fprintf(fout,'%f\t%f\n',R(i),T(i));
end
fclose(fout);


% ===============================================
function c = thermistorFit(R,T,n)
% thermistorFit  Curve fit the Steinhart-Hart equation and a degree n polynomial
%                to the thermistor data
%
% Synopsis:  c = thermistorFit(R,T)
%            c = thermistorFit(R,T,n)
%
% Input:  R = column vector of delay times returned by the RCTIME command
%             R is proportional to the resistance of the thermistor
%         T = column vector of temperatures in C  (degrees C because the
%             curve fit coeficients will be obtained with T in kelvin
%         n = optional degree of the polynomial used in the alternative fit
%
% Output: c = vector of coefficients in the Steinhart-Hart equation

if nargin<3, n=3;  end   %  Provide defalt value for n

% --- Assemble the m-by-3 design matrix.  m rows of data.  Each column
%     is a basis function for the fit.  The first column is ones, which
%     corresponds the constant in the curve fit.
%     Solve the overdetermined system with the backslash operator, which
%     uses a QR factorization to obtain the least squares fit.  The right
%     hand side vector for the overdetermined system is a column of 1/Tk
%     values, where Tk is the temperature in kelvin.
A = [ones(size(R(:)))  log(R(:))   (log(R(:))).^3];
TKinv = 1 ./ (T(:) + 273.15);
c = A\TKinv;
fprintf('Steinhart-Hart coefficients\n');
fprintf('%17.9e\n',c)

% --- Use a polynomial fit T = f(R).  Converting the temperature scale
%     to kelvin does not improve the fit.
cp = polyfit(R,T,n);

% --- Evaluate and plot the curve fits over the range of R data
Rfit = linspace(min(R),max(R))';
Tfit = polyval(cp,Rfit);
TfitSH = 1./ (c(1) + c(2)*log(Rfit) + c(3)*(log(Rfit)).^3 );

plot(R/1000,T,'o',Rfit/1000,Tfit,'b:',Rfit/1000,TfitSH-273.15,'k--');
xlabel('Resistance (k\Omega)');  ylabel('T (^{\circ}C)');
legend('Data',sprintf('Degree %d polynomial',n),'Steinhart-Hart');

