|
Hi Gerard,
The curvefit function is pretty confusing, so I went ahead and made an example of how to use it to determine what you are looking for.Basically, gfunct has as inputs:
x-value (time)
an array of the values for the fit (i.e. A = [m,b] for y = m*x + b)
The function value, which gets modified in place
The derivative of the function with respect to x. This is pder, for partial derivative, but in this case there is only X for the values to change with respect to
I went ahead and modified the gfunct equation to be exp(-t/tau) and here is an example of how it works:
PRO gfunct, X, A, F, pder
bx = EXP(-X/A[1])
F = A[0] * bx + A[2]
;If the procedure is called with four parameters, calculate the
;partial derivatives.
IF N_PARAMS() GE 4 THEN $
pder = [[bx], [A[0] * X * bx], [replicate(1.0, N_ELEMENTS(X))]]
END
pro exponential_fit_example
compile_opt idl2
ireset, /no_prompt
X = findgen(50)
;exponential decay constant, i.e. half life,
tau = 5
;starting amplitude
amp = 5
;add some random noise
Y = amp*EXP(-X/tau) + randomu(seed,n_elements(X))
;initial guess
A = [1., 4, 3]
;leave weights undefined, because we do not want them
yfit = CURVEFIT(X, Y, weights, A, SIGMA, FUNCTION_NAME='gfunct')
print, ''
print, 'Fit equation: y = ' +strtrim(A[0],1) + '*exp(-x/' +$
strtrim(A[1],1) + ') + ' + strtrim(A[2],1)
print, 'Decay Constant = ' + strtrim(A[1],1)
print, ''
;make our plots
p1 = scatterplot(X, Y, title = 'Input Data')
p2 = plot(X, yfit, color= 'red', /current, /overplot)
end
|