;+ ; Project : SOHO - CDS ; ; Name : CDS_GAUSS() ; ; Purpose : Fits Gaussian with constant, linear or quadratic background. ; ; Explanation : Uses the standard IDL curvefit to fit a gaussian but unlike ; IDL's gaussfit, it allows the user to choose the background ; type, ie zero, contant, linear or quadratic. ; ; ie Fit y=F(x) where: ; F(x) = a0*exp(-z^2/2) + a3 + a4*x + a5*x^2 ; and z=(x-a1)/a2 ; ; a0 = height of exp, ; a1 = center of exp, ; a2 = 1/e width, ; a3 = constant term, ; a4 = linear term, ; a5 = quadratic term. ; ; The parameters a0,a1,a2,a3 are estimated and then curvefit is ; called. ; Parameters a3, a4, and a5 are optional, depending on the value ; of k ; ; If the (max-avg) of Y is larger than (avg-min) then it is ; assumed the line is an emission line, otherwise it is assumed ; there is an absorption 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 ; extreme until a point is found < the 1/e value. ; ; Use : yfit = cds_gauss(x,y,a,k) ; ; Inputs : x = independent variable, must be a vector. ; y = dependent variable, must have the same number of points ; as x. ; ; Opt. Inputs : k = determines the order of polynomial to be fitted in ; addition to the gaussian part: ; ; k = 3 : gaussian + a3 + a4*X + a5*X^2 ; k = 2 : gaussian + a3 + a4*X ; k = 1 : gausian + a3 ; k = 0 or any other k: only gaussian ; ; If k is not defined, it is assumed k=0 ; ; ; Outputs : Function returns value yfit = fitted function. ; ; Opt. Outputs: a = coefficients. a (3+k) element vector as described above. ; sigmaa = vector of coefficient uncertainties ; ; Keywords : INTER - if given a value the output array is of that dimension ; ie the fit is interpolated (eg for a smooth curve to ; overplot). Note that the return value is then an ; array of dimension (inter,2) containing both the ; interpolated x and y values. ; ; INITIAL_A : Array, same content as A. Set to supply initial ; values of the fit parameters. Before CURVEFIT is ; called, the parameter array a is updated with: ; ; A(0) = INITIAL_A ; ; Thus this array may contain fewer values than A ; (so you can drop estimating the polynomial terms ; if you like) ; ; Calls : None ; ; Common : None ; ; Restrictions: None ; ; Side effects: None ; ; Category : Util, numerical ; ; Prev. Hist. : Original - gaussfit ; Background flexibility by A. Fludra (MSSL) ; ; Written : CDS version by C D Pike, RAL, 7-Jan-94 ; ; Modified : Added INTER keyword. CDP, 27-Mar-96 ; ; Version 3, S.V.H.Haugan, UiO/GSFC, 26 November 1996 ; Added INITIAL_A keyword. ; Version 4, Fix bug in INTERPOLATE mode. ; version 5, Added sigmaa optional output, J A Klimchuk, 97-10-07 ; ; Version : Version 5, 9-Oct-97 ;- ; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ; ;+ ; Project : SOHO - CDS ; ; Name : CDS_GFUNCT ; ; Purpose : Evaluate gaussian and its partial derivatives. ; ; Explanation : Used in call to curvefit from cds_gauss ; ; Use : Only as parameter to curvefit ; ; Inputs : X = VALUES OF INDEPENDENT VARIABLE. ; A = PARAMETERS OF EQUATION DESCRIBED BELOW. ; ; Opt. Inputs : None ; ; Outputs : Value of function at each x value ; ; Opt. Outputs: None ; ; Keywords : None ; ; Calls : None ; ; Common : None ; ; Restrictions: For use with cds_gauss and curvefit ; ; Side effects: None ; ; Category : Util, numerical ; ; Prev. Hist. : As below for cds_gauss ; ; Written : CDS version C D Pike, RAL, 7-Jan-94 ; ; Modified : Trap no data. CDP, 18-Jan-95 ; ; Version : Version 2, 18-Jan-95 ;- pro cds_gfunct,X,A,F,PDER k = N_elements(a) Z = (X-A(1))/A(2) ;GET Z EZ = EXP(-Z^2/2.)*(ABS(Z) LE 7.) ;GAUSSIAN PART IGNORE SMALL TERMS f = a(0)*ez if k eq 4 then f = f + a(3) if k eq 5 then f =f + a(3) + a(4)*X if k eq 6 then f = f + a(3) + a(4)*X + a(5)*X^2 IF N_PARAMS(0) LE 3 THEN RETURN ;NEED PARTIAL? PDER = FLTARR(N_ELEMENTS(X),k) ;YES, MAKE ARRAY. PDER(0,0) = EZ ;COMPUTE PARTIALS PDER(0,1) = A(0) * EZ * Z/A(2) PDER(0,2) = PDER(*,1) * Z if k eq 4 then PDER(*,3) = 1. if k eq 5 then begin PDER(*,3) = 1. PDER(0,4) = X endif if k eq 6 then begin PDER(*,3) = 1. PDER(0,4) = X PDER(0,5) = X^2 endif RETURN END function CDS_Gauss, x, y, a, k, inter=inter, initial_a=initial_a, $ sigmaa=sigmaa ; ; check have some data ; if n_elements(x) eq 0 or n_elements(y) eq 0 then return,0 ; ; if parameter k not defined, assume = 0 ; if N_params(0) le 3 then k = 0 ; ; number of data ; n = n_elements(y) ; ; if background fit wanted, fit line ; yd = y if k gt 0 then begin c = poly_fit(x,y,1,yf) yd = y-yf endif ; ; x,y subscripts of extreme ; ymax=max(yd) & xmax=x(!c) & imax=!c ymin=min(yd) & xmin=x(!c) & imin=!c ; ; coefficient vector ; a=fltarr(3+k) ; ; emission or absorption? ; if abs(ymax) gt abs(ymin) then i0=imax else i0=imin i0 = i0 > 1 < (n-2) ;never take edges dy=yd(i0) ;diff between extreme and mean del = dy/exp(1.) ;1/e value i=0 ; ; estimate width by creeping along ; while ((i0+i+1) lt n) and $ ((i0-i) gt 0) and $ (abs(yd(i0+i)) gt abs(del)) and $ (abs(yd(i0-i)) gt abs(del)) do i=i+1 ; ; initial guesstimates of gaussian parameters ; a = [yd(i0), x(i0), abs(x(i0)-x(i0+i))] ; ; add parameters as requested ; if k eq 1 then a = [a,c(0)] if k eq 2 then a = [a,c(0), c(1)] if k eq 3 then a = [a,c(0), c(1), 0.] ; ; reset cursor for plotting ; !c=0 ; ; and calculate fitted values ; IF keyword_set(initial_a) THEN a(0) = initial_a fit = curvefit(x,y,replicate(1.,n),a,sigmaa,funct='cds_gfunct') ; ; was interpolation requested? ; if keyword_set(inter) then begin step = (last_item(x)-x(0))/(float(inter)-1) xx = findgen(inter)*step+x(0) yy = gauss_put(xx,0,a(0),a(1),a(2)) aa = dblarr(3) if n_elements(a) gt 3 then begin for i=3,n_elements(a)-1 do aa(i-3) = a(i) yy = yy + aa(0) + aa(1)*xx + aa(2)*xx*xx endif return,[[xx],[yy]] endif else begin return,fit endelse end