spcfun.m

ASUFIT－Matlab－全局拟合程序

function f = spcfun(params)

global data
global irf
global weight
global fit
global amps
global numTao
global numTaoFixed
global taoFixed
global t
global chisq
global textBox
global spcLeftEnd
global spcRightEnd
global absLeftEnd2
global absRightEnd2
%This is the function referred to in the call to "leastsq".
%It assumes the form:
%y(x)=a11*exp(-k1*t)+a12*exp(-k2*t)+. . .+a1n*exp(-kn*t)
%     a21*exp(-k1*t)+a22*exp(-k2*t)+. . .+a2n*exp(-kn*t)	
%     .
%     .
%     .
%     am1*exp(-k1*t)+am2*exp(-k2*t)+. . .+amn*exp(-kn*t)
%and makes available the vector of k values and the matrix
%of a values to the script file "fit.m".

tpch = t(2)-t(1);
dataTemp = data(spcLeftEnd:spcRightEnd,:);
irfTemp = irf(spcLeftEnd:spcRightEnd,:);
numPoints = size(dataTemp,1);
numCurves = size(dataTemp,2);
curve = zeros(numPoints,numTao + numTaoFixed);
fit = zeros(size(dataTemp));
a = zeros(numTao + numTaoFixed,numCurves);
taoFit = zeros(1,numTao + numTaoFixed);
tao = params(1:numTao);
shift = params(numTao + 1:numTao+numCurves);
shift = shift./tpch;
taoFit(1,1:numTaoFixed) = taoFixed;
taoFit(1,numTaoFixed+1:numTao + numTaoFixed) = tao;
taoFit = taoFit./tpch;
%a=params(numTao+2:length(params))';
bottom = floor(shift);
top = ceil(shift);
R = exp(-1*taoFit.^(-1));
for i= 1:numCurves,
   if (top(i) < 1)
      top(i) = 1;
   end
   curve(1,:) = irfshift(irfTemp(top(i),i),irfTemp(top(i)+1,i),shift(i));
   for count = 1:(numPoints-1),
      top(i) = bottom(i) + count;
      if (top(i) < 1) top(i) = 1; end
      if (top(i) > numPoints-2) top(i) = numPoints-2; end
      curve(count+1,:) = R.*(0.5*irfshift(irfTemp(top(i),i),irfTemp(top(i)+1,i),shift(i))+...
      curve(count,:)) + 0.5 * irfshift(irfTemp(top(i)+1,i),irfTemp(top(i)+2,i),shift(i));
   end
      a(:,i) = curve\dataTemp(:,i);
      z(:,i) = curve*a(:,i);
end
f = (z - dataTemp)./weight(spcLeftEnd:spcRightEnd,:);
amps = a;
fit=z;
chisq = trace(f'*f)/(numPoints*numCurves-(numTao + numTaoFixed)-(numTao + numTaoFixed)*numCurves-numCurves);
set(textBox,'string',['chisq = ' num2str(chisq)])
drawnow;

function G = irfshift(f1,f2,x0)
G = f1 + (x0 - floor(x0)) * (f2 - f1);