SUBROUTINE clcvec(phacur, phasig, phader, nvpar, alpha, beta,
     &                  ier, ier2)

      INTEGER nvpar
      DOUBLE PRECISION alpha(nvpar, nvpar), beta(nvpar)
      REAL phacur(*), phasig(*), phader(*)

      INTEGER ier, ier2

c This is wrap-up for the the routines that calculate the alpha and beta
c arrays. Beta is the minus half times the vector of first derivatives of the 
c statistic wrt to the fit parameters and alpha half times the matrix of 
c second derivatives. Note that to avoid numerically calculating the second 
c derivatives analytic calculations are used to express the derivatives in 
c terms of the derivatives of the model (or (y_obs-y_mod)/sigma in the case 
c of chi-squared). The terms involving second derivatives of the model are 
c assumed to be small.

c	arguments:
c	phacur	r	 i: vector of T=(y_obs-y_mod)/sigma at current parameter 
c	phasig	r	 i: vector of sigma at current parameter 
c	phader	r	 i: vector of derivatives of T 
c	alpha	r*8	 r: array of derivatives
c	beta	r*8	 r: vector required by crfitx
c	ier	i	 r: error flag (see below)
c	ier2	i	 r: error info flag: (see below)
c	error conditions:
c	ier=	ier2=
c	0	-	normal return
c	-1	ivpar	the indicated variable parameter shows no variation 
c                       of the model when calculating the derivative with 
c                       respect to this parameter 
c       -2              failed to grab memory

      INCLUDE '../../inc/xspec.inc'

c Pointers for local dynamic arrays

      INTEGER idc1, idc2

c Other local variables

      DOUBLE PRECISION dsum

      INTEGER nbsave, ivpar, jvpar, ipha, iphaoff, jphaoff

      CHARACTER wrtstr*255

      DOUBLE PRECISION dlprior, d2lprior
      INTEGER dgdpnt, dgndst, dgnbne
      LOGICAL cgmerr 
      CHARACTER xgstat*11
      EXTERNAL xgstat, dgndst, dgnbne, dgdpnt, cgmerr
      EXTERNAL dlprior, d2lprior

      SAVE idc1, idc2, nbsave

      DATA idc1, idc2, nbsave/-1,-1,-1/

c Grab the memory for summation arrays

      IF ( dgnbne(dgndst()) .NE. nbsave ) THEN
         nbsave = dgnbne(dgndst())
         CALL udmget(nbsave, 7, idc1, ier2)
         CALL udmget(nbsave, 7, idc2, ier2)
         IF ( ier2 .NE. 0 ) THEN
            CALL xwrite('Memory allocation failed in CLCVEC', 2)
            ier = -2
            RETURN
         ENDIF
      ENDIF

c Do the case of chi-squared

      IF ( xgstat() .EQ. 'Chi-Squared') THEN

          CALL AB_CHI2(phacur, MEMD(idc1), MEMD(idc2))

c Do the case of C-statistic. 

      ELSEIF ( xgstat() .EQ. 'C-statistic' ) THEN

          IF(cgmerr()) call xwrite(
     &   'Warning: model error ignored used in C-statistics',5) 

          CALL AB_CSTAT(MEMR(dgdpnt('src')), MEMR(dgdpnt('bkg')),
     &                  MEMR(dgdpnt('mod')), MEMR(dgdpnt('reg')),
     &                  MEMR(dgdpnt('brg')), MEMR(dgdpnt('are')),
     &                  MEMR(dgdpnt('bre')), MEMD(idc1), MEMD(idc2))

c Do the case of L-statistic

      ELSEIF ( xgstat() .EQ. 'L-statistic' ) THEN

          IF(cgmerr()) call xwrite(
     &   'Warning: model error ignored used in L-statistics',5) 

          CALL AB_LSTAT(MEMR(dgdpnt('src')), MEMR(dgdpnt('bkg')),
     &                  MEMR(dgdpnt('mod')), MEMR(dgdpnt('reg')),
     &                  MEMR(dgdpnt('brg')), MEMR(dgdpnt('are')),
     &                  MEMR(dgdpnt('bre')), MEMD(idc1), MEMD(idc2))

      END IF

c If the statistic is either C or L then we reset phader so it is the
c derivative of M not of (O-M)/sigma.

      IF ( xgstat() .EQ. 'C-statistic' .OR. 
     &     xgstat() .EQ. 'L-statistic') THEN

          DO ivpar = 1, nvpar
             iphaoff = (ivpar-1)*nbsave
             DO ipha = 1, nbsave
                phader(ipha+iphaoff) = - phasig(ipha) 
     &                                 * phader(ipha+iphaoff)
             ENDDO
          ENDDO

       ENDIF

c Now loop over the variable parameters loading the alpha and beta vectors

      DO ivpar = 1, nvpar

         dsum = 0.D0
         iphaoff = (ivpar-1)*nbsave

         DO ipha = 1, nbsave
            dsum = dsum + MEMD(idc1+ipha-1)*phader(ipha+iphaoff)
         ENDDO

         beta(ivpar) = dsum

         DO jvpar = ivpar, nvpar

            dsum = 0.D0
            jphaoff = (jvpar-1)*nbsave

            DO ipha = 1, nbsave
               dsum = dsum + MEMD(idc2+ipha-1)*phader(ipha+iphaoff)
     &                                        *phader(ipha+jphaoff)
            ENDDO

            alpha(ivpar, jvpar) = dsum - d2lprior(ivpar, jvpar)
            alpha(jvpar, ivpar) = dsum - d2lprior(jvpar, ivpar)

         ENDDO


c  if there is no variation in the model for this parameter then return
c  with an error 4.

         IF ( alpha(ivpar,ivpar) .EQ. 0. ) THEN
            ier = -1
            ier2 = ivpar
            WRITE(wrtstr, '(a,i8)') 
     &  'CLCVEC: No variation in the model for variable parameter ',
     &  ivpar
            CALL xwrite(wrtstr,20)
            RETURN
         ENDIF

      ENDDO

      RETURN
      END 



c============================================================

      SUBROUTINE AB_CSTAT(phaobs, phabkg, phamod, region, brgion, 
     &                    areasc, breasc, dc1, dc2)

      IMPLICIT NONE

      DOUBLE PRECISION dc1(*), dc2(*)
      REAL phaobs(*), phamod(*), phabkg(*), region(*), brgion(*)
      REAL areasc(*), breasc(*)

c  calculate the summation arrays used in the alpha matrix and beta vector.

c       arguments:
c       phaobs  r        i: observed rate in the ith combined pha channel.
c       phabkg  r        i: background rate in the ith combined pha channel.
c       phamod  r        i: the current models' expected rate of the ith pha
c                           channel.
c       region  r        i: source backscal array
c       brgion  r        i: background backscal array
c       dc1     r*8      r: array used in first derivative summation
c       dc2     r*8      r: array used in second derivative summation

      DOUBLE PRECISION dconv, dbconv, dtconv, dbfact, dfact
      DOUBLE PRECISION dmod, dobs, dbobs
      DOUBLE PRECISION da, db, dc, dd, dq
      DOUBLE PRECISION df, dg, dh
      INTEGER ipha, idset, isign
      LOGICAL qbkgf

      INTEGER dgndst, dgnbne, dgnbnb
      REAL dgtime, dgbtme
      EXTERNAL dgndst, dgnbne, dgnbnb, dgbtme, dgtime

c Initialize the summation arrays
 
      DO ipha = 1, dgnbne(dgndst())
         dc1(ipha) = 0.d0
         dc2(ipha) = 0.d0
      ENDDO

c loop around the pha bins

      idset = 1
      dfact = DBLE(dgtime(idset))
      dbfact = DBLE(dgbtme(idset))
      qbkgf = ( dgbtme(idset) .NE. 0. )

      DO ipha = 1, dgnbne(dgndst())

c if this is the first bin or we are moving onto another dataset
c then set the dataset stuff

         IF ( ipha .GT. dgnbne(idset) ) THEN

            idset = idset + 1
            dfact = DBLE(dgtime(idset))
            dbfact = DBLE(dgbtme(idset)) 
            qbkgf = ( dgbtme(idset) .NE. 0. )

         ENDIF

         dconv = dfact * DBLE(areasc(ipha))
         dobs = dconv*DBLE(phaobs(ipha))

         dbconv = 0.0
         dbobs = 0.0
         IF ( qbkgf ) THEN
            IF ( region(ipha) .NE. 0. ) THEN
               dbconv = dbfact*DBLE(breasc(ipha))
     &                     *DBLE(brgion(ipha))/DBLE(region(ipha))
            ELSE
               dbconv = dbfact
            ENDIF
            dbobs = dbconv*DBLE(phabkg(ipha))
         ENDIF

         dtconv = dconv + dbconv

c Put a floor under the model of 1.e-5 counts/bin

         dmod = MAX(DBLE(phamod(ipha)), 1.D-5/dconv)

c First do case of no background spectrum

         IF ( .NOT. qbkgf ) THEN

            dc1(ipha) = dobs/dmod - dconv
            dc2(ipha) = dobs/dmod**2

c Now have background spectrum: special case of no data counts

         ELSEIF ( dobs .EQ. 0.0d0 ) THEN

            dc1(ipha) = -dconv
            dc2(ipha) = 0.0

c If there are no background counts...

         ELSEIF ( dbobs .EQ. 0.0d0 ) THEN

            IF ( dmod .LE. dobs/dtconv ) THEN
               dc1(ipha) = dbconv
               dc2(ipha) = 0.0
            ELSE
               dc1(ipha) = (dobs/dmod)-dconv
               dc2(ipha) = dobs/dmod**2
            ENDIF

c else if there are both data and background counts...

         ELSE

            da = dtconv
            db = dtconv*dmod - dobs - dbobs
            dc = -dbobs*dmod
            dd = sqrt(db**2 - 4*da*dc)
            isign = 1
            IF ( db .LT. 0.0d0 ) isign = -1
            dq = -0.5*(db + isign*dd)
            df = dq/da
            IF ( df .LT. 0.0d0 ) THEN
               df = dc/dq
               dg = (dtconv*dmod - dobs - dbobs - dd)/(2*dd)
               dh = 2*dtconv*(dobs/dd)*(dbobs/dd)/dd
            ELSE
               dg = (-dtconv*dmod + dobs - dbobs - dd)/(2*dd)
               dh = -2*dtconv*(dobs/dd)*(dbobs/dd)/dd
            ENDIF

cc            dobs1 = dobs/dtconv
cc            dbobs1 = dbobs/dtconv
cc            da1 =  SQRT( (dmod-dobs1-dbobs1)**2
cc     &                   +4*dbobs1*dmod )
cc            df = (dobs1 + dbobs1 - dmod + da1)/2.
cc            dg = (dmod - dobs1 + dbobs1 - da1)/2./da1
cc            dh = 2*(dobs1/da1)*(dbobs1/da1)/da1

            IF ( (dmod+df) .NE. 0.d0 .AND. df .NE. 0.d0 ) THEN

               dc1(ipha) = dobs*(1+dg)/(dmod+df) 
     &                        + dbobs*dg/df - dtconv*dg - dconv
               dc2(ipha) = dtconv*dh
     &                - dobs*(dh-(1+dg)**2/(dmod+df))/(dmod+df)
     &                - dbobs*(dh-dg**2/df)/df

            ENDIF

         ENDIF

      ENDDO

      RETURN
      END

c============================================================

      SUBROUTINE AB_LSTAT(phaobs, phabkg, phamod, region, brgion, 
     &                    areasc, breasc, dc1, dc2)

      DOUBLE PRECISION dc1(*), dc2(*)
      REAL phaobs(*), phamod(*), phabkg(*), region(*), brgion(*)
      REAL areasc(*), breasc(*)

c  calculate the summation arrays used in the alpha matrix and beta vector.

c       arguments:
c       phaobs  r        i: observed rate in the ith combined pha channel.
c       phabkg  r        i: background rate in the ith combined pha channel.
c       phamod  r        i: the current models' expected rate of the ith pha
c                           channel.
c       region  r        i: source backscal values
c       brgion  r        i: background backscal values
c       areasc  r        i: source areascal values
c       breasc  r        i: background areascal values
c       dc1     r*8      r: array used in first derivative summation
c       dc2     r*8      r: array used in second derivative summation

      DOUBLE PRECISION dconv, dbconv, dtconv, dfact, dbfact
      DOUBLE PRECISION dmod, dobs, dbobs
      DOUBLE PRECISION df, dg, dh, dt1, dt2
      INTEGER ipha, idset, nobs, nbobs
      LOGICAL qbkgf

      DOUBLE PRECISION lorsum
      INTEGER dgndst, dgnbne, dgnbnb
      REAL dgtime, dgbtme
      EXTERNAL dgndst, dgnbne, dgnbnb, dgbtme
      EXTERNAL dgtime, lorsum

c initialize the work arrays
 
      DO ipha = 1, dgnbne(dgndst())
         dc1(ipha) = 0.d0
         dc2(ipha) = 0.d0
      ENDDO

c loop around the pha bins

      idset = 1
      dfact = DBLE(dgtime(idset))
      dbfact = DBLE(dgbtme(idset))
      qbkgf = ( dgbtme(idset) .NE. 0. )

      DO ipha = 1, dgnbne(dgndst())

c if this is the first bin or we are moving onto another dataset
c then set the dataset stuff

         IF ( ipha .GT. dgnbne(idset) ) THEN

            idset = idset + 1
            dfact = DBLE(dgtime(idset))
            dbfact = DBLE(dgbtme(idset))
            qbkgf = ( dgbtme(idset) .NE. 0. )

         ENDIF

         dconv = DBLE(areasc(ipha)) * dfact
         dobs = DBLE(phaobs(ipha))
         dbobs = DBLE(phabkg(ipha))
         dmod = DBLE(phamod(ipha))

c If there is not a background file then this defaults to the C-statistic

         IF ( .NOT.qbkgf ) THEN

            IF ( dmod .NE. 0.d0 ) THEN
               dc1(ipha) = dconv*(dobs/dmod-1.d0)
               dc2(ipha) = dconv*dobs/dmod**2
            ENDIF

c If there is a background file...

         ELSE

            IF ( region(ipha) .NE. 0. ) THEN
               dbconv = dbfact*DBLE(breasc(ipha))
     &                        *DBLE(brgion(ipha))/DBLE(region(ipha))
            ELSE
               dbconv = dbfact
            ENDIF
            dtconv = dconv + dbconv

            dt1 = LOG(dtconv)
            dt2 = 2 * dt1

            nobs = NINT(dobs*dconv)
            nbobs = NINT(dbobs*dbconv)

c we use df for the log of the sum over dmod^j, dg for j dmod^(j-1) and 
c dh for j(j-1) dmod^(j-2).

            df = lorsum(dtconv*dmod, nobs, nbobs, 0)
            dg = lorsum(dtconv*dmod, nobs, nbobs, 1) + dt1
            dh = lorsum(dtconv*dmod, nobs, nbobs, 2) + dt2

c Now pull all the terms together and attempt to avoid numerical overflows

            dc1(ipha) = EXP(dg-df) - dconv
            dc2(ipha) = EXP(2*(dg-df)) - EXP(dh-df)

         ENDIF

      ENDDO

      RETURN
      END


c==============================================================

      SUBROUTINE AB_CHI2(phacur, dc1, dc2)

      DOUBLE PRECISION dc1(*), dc2(*)
      REAL phacur(*)

c  calculate the summation arrays used in the alpha matrix and beta vector.
c  For chi^2 these are particularly simple - being -phacur and 1 respectively.

c       arguments:
c       phacur  r        i: vector of T=(y_obs-y_mod)/sigma at current parameter
c       dc1     r*8      r: array used in first derivative summation
c       dc2     r*8      r: array used in second derivative summation


      INTEGER ipha

      INTEGER dgndst, dgnbne
      EXTERNAL dgndst, dgnbne

c Loop around the pha bins

      DO ipha = 1, dgnbne(dgndst())

         dc1(ipha) = -phacur(ipha)
         dc2(ipha) = 1.d0

      ENDDO

      RETURN
      END