;+ ; NAME: ; FIT2EXP ; ; PURPOSE: ; Fits data to a distribution of 2 exponential functions of the form: ; ; y = aL * exp( -bL*(x-cL) ), x < x0 ; = aR * exp( -bR*(x-cR) ), x > x0 ; ; CATEGORY: ; Math. ; ; CALLING SEQUENCE: ; ; Result = FIT2EXP( X, Y [, Param, Sigma_, X0 ] ) ; ; INPUTS: ; X,Y: Data to be fitted. ; ; OPTIONAL INPUT KEYWORD PARAMETERS: ; ; HELP: Set this keyword to display help messages on how to draw ; a line to determine the initial parameters. ; ; Any other keywords will be passed on as _EXTRA keywords to the plot ; procedure called in this routine. ; ; OUTPUTS: ; This function returns a vector of fitted Y values for the given ; parameters, Param. ; ; OPTIONAL OUTPUTS: ; ; Param: A vector containing the fitted parameters: [aL,bL,cL,aR,bR,cR] ; ; Sigma: A vector of standard deviations for the elements of Param. ; ; X0: The position along the abscissa at the intersection of ; the two exponentials. ; ; OPTIONAL OUTPUT KEYWORD PARAMETERS: ; ; CHI2: Set this keyword equal to a named variable that will contain ; the value of chi-squared. ; ; DOF: Number of degrees of freedom associated with the chi-squared. ; ; HERE: A longword vector that contains the one-dimensional subscripts ; of the fitted elements of Y. ; ; MODIFICATION HISTORY: ; Written by: Han Wen, September 1996. ;- function FIT2EXP, X, Y, Param, Sigma_, X0, HELP=Help, CHI2=Chi2, DOF=DOF, \$ HERE=Here, TITLE=Title, XRANGE=Xrange, YRANGE=Yrange, _EXTRA=Plot_keys if (N_ELEMENTS(Xrange) eq 0) then Xrange = [MIN(X,MAX=mx),mx] if (N_ELEMENTS(Help) eq 0) then Help = 0 ; Fit the "LEFT" exponential fitL = FITEXP( X, Y, paramL, sigmaL, HELP=Help, HERE=HereL, \$ TITLE= 'Fit LEFT exponential', XRANGE=xrange, \$ /RESTRICT, _EXTRA=Plot_keys ) if (N_ELEMENTS(paramL) eq 0) then return,-1 \$ else begin print,'Click on plot to continue...' cursor,dum,dum endelse ; Fit the "RIGHT" exponential fitR = FITEXP( X, Y, paramR, sigmaR, HELP=Help, HERE=HereR, \$ TITLE='Fit RIGHT exponential',XRANGE=xrange, \$ /RESTRICT, _EXTRA=Plot_keys ) if (N_ELEMENTS(paramR) eq 0) then return,-1 \$ else begin print,'Click on plot to continue...' cursor,dum,dum endelse Here = [HereL, HereR] Param = [paramL, paramR] Sigma_ = [sigmaL, sigmaR] nL = paramL(0) & nR = paramR(0) ; amplitude pL = paramL(1) & pR = paramR(1) ; counts/sec TL = paramL(2) & TR = paramR(2) ; offset ; Set X0 at the intersection of the two exponentials X0 = (ALOG(nL/nR) + pL*TL - pR*TR)/(pL - pR) ; Determine the chi-squared for the entire fit hlt = WHERE( X lt X0, nlt ) hge = WHERE( X ge X0, nge ) fit = [fitL(hlt),fitR(hge)] hhL = WHERE( X(hlt) gt X(HereL(0)), nptsL ) hhR = WHERE( X(hge) lt X(HereR(N_ELEMENTS(HereR)-1)), nptsR ) iL = hlt(hhL) iR = hge(hhR) i = [iL, iR] yexp = [fitL(iL),fitR(iR)] Chi2 = TOTAL((Y(i) - yexp)^2/yexp) npar = 2*3 ; constraints: nL, pL, TL, nR, pR, TR DOF = nptsL + nptsR - npar ; Make plot of data and fit loadct, 12 & red = 160 if (N_ELEMENTS(Title) eq 0) then Title = '' plot, X, Y>1, PSYM=1,SYMSIZE=0.5,/YTYPE,TITLE=Title,_EXTRA=Plot_keys oplot,X, fit, COLOR=Red return, fit end