;+
; NAME:
;        FITHEXP
;
; PURPOSE:
;        This routine fits the waiting times for a hyperexponential distribution
;        distorted by a non-extended deadtime,
;
;        f(t) = P1*r1*exp(-r1*(t-tau)) + (1-P1)*r2*exp(-r2*(t-tau))
;
; CATEGORY:
;        Math
;
; CALLING SEQUENCE:
;
;        Result = FITHEXP( Wtbin, Noccur, Norm, Init )
;
; INPUTS:
;
;        Wtbin:    The binned waiting times of the distribution.
;
;        Noccur:   The number of waiting times occurring within a given waiting
;                  time bin.
;
;        Norm:     The total number of waiting times of the entire distribution.
;                  (NOTE: TOTAL(Noccur) < Norm, because Wtbin < infinity).
;
;        Init:     A 2-element array giving the initial values for r1 and r2.
;
; OPTIONAL INPUT KEYWORD PARAMETERS:
;
;        VERBOSE:  Set this keyword to display the current fit parameters before
;                  each iteration, (0=Default).
;
; OUTPUTS:
;
;        This function returns a structure holding the results of the fit.  Tags are:
;
;                  param : Structure holding the fitted parameters with tags:
;                    rate1  : Fitted count rate of source 1 without dead time.
;                    rate2  : Fitted count rate of source 2 without dead time.
;                    offset : Fitted non-extended dead time.
;
;                  sigma : Structure holding the 1 sigma error of the fitted
;                          parameters with tags:
;                    rate1  : 1 sigma error on the fitted count rate1 parameter.
;                    rate2  : 1 sigma error on the fitted count rate2 parameter.
;                    offset : 1 sigma error on the offset parameter.
;
;                  wtbin:    Waiting time bins of the fitted distribution.
;                  noccur:   Number of occurrences at each waiting time bin.
;                  prob:     Probability distribution of the fit.
;                  nwt:      Total number of waiting times (normalization).
;                  chi2:     Total chi-squared of the fit.
;                  dof:      Number of degrees of freedom of the fit.
;                  CL:       Associated confidence level of the chi-squared.
;
;
; COMMON BLOCKS:
;        FITHEXP:   Pass extra parameters to the USER function of CURVEFIT.
;
; MODIFICATION HISTORY:
;        Written by:    Han Wen, November 1996. (Adapted from FITWTIME)
;        06-DEC-1996:   Changed name from FITWTBUG -> FITHEXP. Changed OFFSET1
;                       tag to OFFSET, rate1->rate2, rate2->rate1, frac->Prob1,
;                       Param (output from WTHIST) -> Init. Convert from pro
;                       -> function.
;-
pro FITHEXP_, n, A, Prob

         common FITHEXP, par2_

         Mean1     = ABS(A(0))
         Mean2     = ABS(A(1))
         Offset    = ABS(A(2))
         Frac      = ABS(A(3))

         Norm      = par2_.Norm
         Bin       = par2_.Bin
         Verbose   = par2_.Verbose
         if KEYWORD_SET(Verbose) then print,'params:',STRCOMPRESS(A) & wait,1
         Prob = Norm*PROBHEXP(n, Mean1, Mean2, Offset, Frac, BIN=Bin)
end

function FITHEXP, Wtbin, Noccur, Norm, Init_Param, VERBOSE=Verbose

         common FITHEXP, par2_

         NP   = N_PARAMS()
         if (NP ne 4) then message, $
              'Must be called with 5 parameters: Wtbin, Noccur, Norm, Init'

         Binsize   = Wtbin(2) - Wtbin(1)
         if (N_ELEMENTS(Verbose) eq 0) then Verbose = 0

         nbin = N_ELEMENTS(wtbin)
         i0   = -1
         repeat i0=i0+1 until (noccur(i0) ne 0)
         tau1 = wtbin(i0)

         mean1= 1/Init_Param(0)
         mean2= 1/Init_Param(1)
         Init = [Mean1, Mean2, tau1, 0.5]

         A    =( A0 = Init)
         x    = wtbin(i0:nbin-1)
         y    = noccur(i0:nbin-1)
         w    = 1./y
         par2_= { BIN:Binsize, NORM:norm, VERBOSE:Verbose }
         yfit = CURVEFIT(x,y,w,A,sigA,FUNCTION_NAME='FITHEXP_', $
                        /NODERIVATIVE, CHI2=Chi2, ITER=Iter, ITMAX=50 )
         Prob = PROBHEXP(wtbin(0:nbin-1), A(0), A(1), A(2), A(3), $
                             BIN=Binsize)
         if (A(2) lt (wtbin(i0)-Binsize)) then message, /INFORMATIONAL, $
              'WARNING: Dead time from Source 1 < Smallest Waiting Time '+ $
              'with Counts.'

         param  = { rate1    :1/A(0),            $
                    rate2    :1/A(1),            $
                    offset   :A(2),              $
                    prob1    :A(3)               }
         sigma  = { rate1    :sigA(0)/A(0)^2,    $
                    rate2    :sigA(1)/A(1)^2,    $
                    offset   :sigA(2),           $
                    prob1    :sigA(3)            }

         dof  = N_ELEMENTS(y) - N_ELEMENTS(A)
         CL   = 1. - CHI_SQR1( Chi2*dof, dof )
         minima = { param    :param,             $
                    sigma    :sigma,             $
                    prob     :prob,              $
                    chi2     :chi2*dof,          $
                    dof      :dof,               $
                    CL       :CL                 }
         return, minima
end