% jednoduche pouziti: o=jhjamie(V,I);
% kde o.IsTeVf=[Is,Te,Vf]
% o=jhjamie(V,I[,plots,iprint,comp_err,lw,col])
% error interpretation: within the interval o.IsTeVf+-o.err
function o=jhjamie(V,I,plots,iprint,comp_err,lw,col)
if nargin<3 plots=0; end;
if nargin<4 iprint=0; end;
if nargin<5 comp_err=0; end;
if nargin<6 lw=1; end;
if nargin<7 col='r'; end;
Te=200.0; % the highest possible Te
pl=[];
ifail=0;
% the Jamie Gunn code for
%Non-linear least-squares fit to an exponential with vertical offset.
%I = b + c*exp(d*V) => b=I_sat, c=I_e, d=1/Te
%The saturation must occur at negative bias for this to work, although
%the sign of the current is not important.
%This method is published in my Rev. Sci. Instrum. paper (1995) on
%flush-mounted probes.
%plots=1;
a=[NaN,NaN,NaN]; err=a; Is2=NaN; Vf2=NaN; IJG=NaN; R2=NaN;
[o.IJG,o.R2,o.R2adj]=deal(NaN);
[o.IsTeVf,o.err,o.std]=deal([1,1,1]*NaN);
o.N=length(V);
if o.N<4 ifail=6; return; end;
ind=find(~isnan(V) & ~isnan(I)); V=V(ind); I=I(ind); % firstly removing all NaN pairs:
o.N=length(V);
if o.N<4 ifail=6; return; end;
A=sum(I);
IV=I.*V;
%d=1/Te;
dTe=Te;
ex=exp(V/Te);
B=sum(ex);
C=sum(I.*ex);
ex2=ex.^2;
D=sum(ex2);
E=sum(IV.*ex);
F=sum(V.*ex);
G=sum(V.*ex2);
f=C*G*o.N-A*B*G+A*D*F-D*E*o.N+B^2*E-B*C*F;
%The solution of this transcedental equation f(d)=0
%corresponds to the minimization of chi-squared.
%I use a simple bisection routine.
direction=-1;
Tpred=0; i=0; maxi=30; Tf=[];
while abs(Tpred-Te)>.005*Te
Tpred=Te;
i=i+1;
f0=f;
dTe=dTe/2;
Te=Te+dTe.*direction;
% d=1/Te;
ex=exp(V/Te);
B=sum(ex);
C=sum(I.*ex);
ex2=ex.^2;
D=sum(ex2);
E=sum(IV.*ex);
F=sum(V.*ex);
G=sum(V.*ex2);
f=C*G*o.N-A*B*G+A*D*F-D*E*o.N+B^2*E-B*C*F;
if sign(f.*f0)==-1 direction=-direction; end
% Tf(i)=Te+1i*f;
% this bisection routine can be highly improved by a linear estimation
if Te<0.2 | i>maxi ifail=7; Te=NaN; break; end;
end
%figure(138); plot(Tf);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Te=1/d;
warning off
c=(C*o.N-A*B)/(o.N*D-B^2);
b=(A-c*B)/o.N;
o.IJG=b + c*exp(V/Te);
ESS=sum((I-o.IJG).^2); % error sum of squares
RSS=sum((o.IJG-mean(I)).^2); % residual sum of squares
%TSS=sum((I-mean(I)).^2); % residual sum of squares ; it's the same
TSS=RSS+ESS; % Total sum of squares
o.R2=RSS/TSS; % relative fit goodness = The index of determination in statistics
o.R2adj=1-(1-o.R2)*(o.N-1)/(o.N-3);
if o.R2adj<0.1 ifail=4; end;
Is=b;
if -c/b<0 Vf=NaN; else Vf=-log(-c/b)*Te; end;
warning on
o.IsTeVf=[Is,Te,Vf];
if 1 % Horacek's method of error computation; written on 20/01/2011, published in http://server.ipp.cas.cz/~horacek/pub/Adamek_EPS2011.pdf
% based on page 181 in Statisticke metody , Jiri Andel, Matfyzpress 1993 (knihu ma Vladimir Weinzettl)
% F=[1+0*V exp(d*V) c*exp(d*V).*V];
warning off
F=[o.IJG/Is (o.IJG-Is).*(V-Vf)/Te^2 (o.IJG-Is)/Te]; % these are derivatives of IJG=Is*(1-exp(-(V-Vf)/Te)) by the parameters [Is, Te, Vf]
o.std=sqrt(diag(ESS*inv(F'*F))); % the diagonal terms yields the parameter error, ie. sqrt(FWF)
o.err=jhkvantil(0.05,o.N-3)*o.std/sqrt(o.N-3); % the diagonal terms yields the parameter error, ie. sqrt(FWF)
warning on
end
if plots
hh=ishold; hold on;
pl=plot(V,I,V,plots*o.IJG,'color',col,'LineWidth',lw);
if comp_err errorbar(mean(V(indIs)),plots*Is,plots*o.err(1),'r'); end;
if ~isnan(Vf2)
plot(V(ind),V(ind)*B(1)+B(2));
plot(Vf2+errVf2*[-1,1],[0,0],'b');
end
if hh hold on; else hold off; end;
end;
function u=jhkvantil(alpha,f)
% function u=jhkvantil(alpha,f)
% critical value of t-destribution at level 1-alpha/2 for f number of freedom
% at 1-a/2 credibility for normal distribution N(0,1) of the expected error
F= [NaN 0.05 0.01 0.001];
df=[ ...
1 12.706 63.656 636.578
2 4.303 9.925 31.600
3 3.182 5.841 12.924
4 2.776 4.604 8.610
5 2.571 4.032 6.869
6 2.447 3.707 5.959
7 2.365 3.499 5.408
8 2.306 3.355 5.041
9 2.262 3.250 4.781
10 2.228 3.169 4.587
11 2.201 3.106 4.437
12 2.179 3.055 4.318
13 2.160 3.012 4.221
14 2.145 2.977 4.140
15 2.131 2.947 4.073
16 2.120 2.921 4.015
17 2.110 2.898 3.965
18 2.101 2.878 3.922
19 2.093 2.861 3.883
20 2.086 2.845 3.850
21 2.080 2.831 3.819
22 2.074 2.819 3.792
23 2.069 2.807 3.768
24 2.064 2.797 3.745
25 2.060 2.787 3.725
30 2.042 2.750 3.646
35 2.030 2.724 3.591
40 2.021 2.704 3.551
45 2.014 2.690 3.520
50 2.009 2.678 3.496
60 2.000 2.660 3.460
70 1.994 2.648 3.435
80 1.990 2.639 3.416
90 1.987 2.632 3.402
100 1.984 2.626 3.390
300 1.960 2.576 3.290
];
u=df(closest(f,df(:,1)),alpha==F);
