; NAME:
; rgfit
; PURPOSE:
; Evaluate the sum of a gaussian and a 2nd order polynomial.
; DESCRIPTION:
; Evaluate the sum of a gaussian and a 2nd order polynomial
; and optionally return the value of it's partial derivatives.
; normally, this function is used by curvefit to fit the
; sum of a line and a varying background to actual data.
;
; CATEGORY:
; E2 - curve and surface fitting.
; CALLING SEQUENCE:
; funct,r,a,f,pder
; INPUTS:
; r = values of independent variable.
; a = parameters of equation described below.
; OUTPUTS:
; f = value of function at each r(i).
;
; OPTIONAL OUTPUT PARAMETERS:
; pder = (n_elements(r),5) array containing the
; partial derivatives. p(i,j) = derivative
; at ith point w/respect to jth parameter.
; COMMON BLOCKS:
; SIDE EFFECTS:
; RESTRICTIONS:
; PROCEDURE:
; F = A(0)*EXP(-Z^2/2) + A(2) + A(3)*R + A(4)*R^2
; Z = R/A(1)
; MODIFICATION HISTORY:
; WRITTEN, DMS, RSI, SEPT, 1982.
; Modified, DMS, Oct 1990. Avoids divide by 0 if A(1) is 0.
; Added to Gauss_fit, when the variable function name to
; Curve_fit was implemented. DMS, Nov, 1990.
; Converted to a radial fitting program, MWB, Dec, 1992.
;
PRO RGAUSS_FUNCT,R,A,F,PDER
ON_ERROR,2 ;Return to caller if an error occurs
if a[1] ne 0.0 then Z = R/A[1] $ ;GET Z
else z= 10.
EZ = EXP(-Z^2/2.)*(ABS(Z) LE 7.) ;GAUSSIAN PART IGNORE SMALL TERMS
F = A[0]*EZ + A[2] + A[3]*R + A[4]*R^2 ;FUNCTIONS.
IF N_PARAMS(0) LE 3 THEN RETURN ;NEED PARTIAL?
;
PDER = FLTARR(N_ELEMENTS(R),5) ;YES, MAKE ARRAY.
PDER[0,0] = EZ ;COMPUTE PARTIALS
if a[1] ne 0. then PDER[0,1] = A[0] * EZ * Z^2/A[1]
PDER[*,2] = 1.
PDER[0,3] = R
PDER[0,4] = R^2
RETURN
END
function rgfit, r, y, a
;+
; NAME:
; rgfit
;
; PURPOSE:
; Fit a radial gaussian function and 2nd order polynominal to the input data.
;
; DESCRIPTION:
; Fit the equation y=f(r) where:
;
; F(r) = A0*EXP(-z^2/2) + A2 + A3*r + A4*r^2
; and
; z=r/A2
;
; A0 = height of exp, A1 = sigma (the width).
; A2 = constant term, A3 = linear term, A4 = quadratic term.
; The parameters A0, A1, A2 are estimated and then CURVEFIT is
; called. The gaussian is assumed to be centered at r=0.
;
; CATEGORY:
; Function fitting
;
; CALLING SEQUENCE:
; Result = GAUSSFIT(R, Y [, A])
;
; INPUTS:
; R: The independent variable. R must be a vector.
; Y: The dependent variable. Y must have the same number of points
; as R.
;
; OUTPUTS:
; The fitted function is returned.
;
; OPTIONAL OUTPUT PARAMETERS:
; A: The coefficients of the fit. A is a five-element vector as
; described under PURPOSE.
;
; COMMON BLOCKS:
; None.
;
; SIDE EFFECTS:
; None.
;
; RESTRICTIONS:
;
; PROCEDURE:
; The initial guess of the intensity is the value of the point
; with the smallest r value. The first point that drops below
; half of the first is the guess at the width.
;
; If the (MAX-AVG) of Y is larger than (AVG-MIN) then it is assumed
; that the line is an emission line, otherwise it is assumed there
; is an absorbtion line. The estimated center is the MAX or MIN
; element. The height is (MAX-AVG) or (AVG-MIN) respectively.
; The width is found by searching out from the extrema until
; a point is found less than the 1/e value.
;
; MODIFICATION HISTORY:
; DMS, RSI, Dec, 1983.
; Rewritten to change initial guess, MWB, Lowell Obs., Dec. 1992.
; 93/05/11, put in failsafe on gaussian width guess, MWB, Lowell Obs.
;-
;
on_error,2 ;Return to caller if an error occurs
n = n_elements(y) ;# of points.
idx=sort(r) ;Put all the values in increasing r order.
ymax=y[idx[0]] ;The max occurs closest to zero.
ymin=min(y) ;Take whatever min there is.
a=fltarr(5) ;coefficient vector
dy=ymax-ymin ;diff between extreme and mean
del = dy/exp(1.)+ymin ;1/e value
z=where( y[idx] lt del,count) ;Index for r closest to 1/e
if count eq 0 then z=n_elements(idx)/2 ;Failsafe in case of pathological data.
ewid = r[idx[min(z)]] ;The r closest to 1/e of peak.
a = [dy, ewid, ymin, 0., 0.] ;estimates
return,curvefit(r,y,replicate(1.,n),a,sigmaa, $
function_name = "RGAUSS_FUNCT") ;call curvefit
end