;+
; 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