SUBROUTINE calcml(phamod, phaobs, phacor, phabkg, region, brgion,
     &                  areasc, breasc, likely)

      IMPLICIT NONE

      REAL phamod(*), phaobs(*), phacor(*), phabkg(*)
      REAL region(*), brgion(*), areasc(*), breasc(*)
      REAL likely

c      calcml      kaa 9 January 1989
c            XSPEC subroutine to calculate likelihood for a particular
c            measurement array
c!!  try and put in ML formula from Cash - this only works if
c    errors are only counting statistics on source
c    The pathological case of model(E)=0 is handled by setting
c    log(model(E))=-32.

c  6/23/95 use normalization suggested by John Castor to give a
c          goodness-of-fit

c      phamod      r            i: model values
c      phaobs      r            i: observed count rate
c      phacor      r            i: correction on observed count rate
c      phabkg      r            i: background count rate
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      likely      r            r: value of likelihood


      DOUBLE PRECISION dconv, dsum, dobs, dmod, ldobs, ldmod
      DOUBLE PRECISION dbconv, dbobs, dtconv, dbfact, dfact
      INTEGER ipha, idset
      CHARACTER wrtstr*255
      LOGICAL qwarn

      DOUBLE PRECISION lprior, wcntrb
      REAL dgtime, dgbtme
      INTEGER dgndst, dgnbnb, dgnbne
      LOGICAL fgqcrf
      EXTERNAL dgndst, dgnbnb, dgnbne, dgtime, dgbtme
      EXTERNAL fgqcrf, lprior, wcntrb

      dsum = 0.d0

      DO idset = 1, dgndst()

         IF ( fgqcrf(idset) ) THEN
            WRITE(wrtstr,'(a,i4)') 
     &        'Correction data being ignored for dataset ', idset
            CALL xwrite(wrtstr, 10)
         ENDIF

c Check whether any of the source datapoints are negative - if so we should
c not be running cstat so warn the user

         qwarn = .FALSE.
         DO ipha = dgnbnb(idset), dgnbne(idset)
            IF ( phaobs(ipha) .LT. 0. ) qwarn = .TRUE.
         ENDDO
         IF ( qwarn ) THEN
            WRITE(wrtstr,'(a,i4,a,a)') 
     &        'Dataset ', idset, ' has bins with negative counts - ',
     &        'cstat will not work correctly'
            CALL xwrite(wrtstr, 10)
         ENDIF

         dfact = DBLE(dgtime(idset))

c First calculate the case of no background

         IF ( dgbtme(idset) .EQ. 0. ) THEN

            DO ipha = dgnbnb(idset), dgnbne(idset)

               dconv = dfact * DBLE(areasc(ipha))

               IF (phaobs(ipha) .LE. 0. ) THEN

                  dsum = dsum + dconv*DBLE(phamod(ipha))

               ELSE

                  dobs = dconv * DBLE(phaobs(ipha))
                  ldobs = log(dobs)

c Put a floor under the model of 1.e-5 counts/bin - should be well below
c the systematics limit

                  dmod = MAX(dconv * DBLE(phamod(ipha)), 1.d-5)
                  ldmod = log(dmod)

                  dsum = dsum + dmod - dobs + dobs*(ldobs-ldmod)

               ENDIF

            ENDDO

c Now if there is a background file defined

         ELSE

c Check for negative bins in the background file

            qwarn = .FALSE.
            DO ipha = dgnbnb(idset), dgnbne(idset)
               IF ( phabkg(ipha) .LT. 0. ) qwarn = .TRUE.
            ENDDO
            IF ( qwarn ) THEN
               WRITE(wrtstr,'(a,i4,a,a)') 
     &           'Dataset ', idset, ' has background bins with ',
     &           'negative counts - cstat will not work correctly'
               CALL xwrite(wrtstr, 10)
            ENDIF

            dbfact = DBLE(dgbtme(idset))

c Loop round channels in this dataset

            DO ipha = dgnbnb(idset), dgnbne(idset)

               dconv  = dfact*DBLE(areasc(ipha))
               dbconv = dbfact*DBLE(breasc(ipha))
     &                        *DBLE(brgion(ipha))/DBLE(region(ipha))
               dtconv = dconv + dbconv

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

               dsum = dsum + wcntrb(dobs, dbobs, dmod, dconv, 
     &                              dbconv, dtconv)

c end loop over channels

            ENDDO

         ENDIF

c end loop over datasets

      ENDDO

c Include the Bayesian prior term (the log of the prior probability)

      likely = 2 * SNGL(dsum) - 2 * lprior()

      RETURN
      END

c ****************************************************************************

      FUNCTION wcntrb(dobs, dbobs, dmod, dconv, dbconv, dtconv)

      DOUBLE PRECISION wcntrb, dobs, dbobs, dmod, dconv, dbconv, dtconv

c Routine to calculate the contribution to W/2 from a given channel

      DOUBLE PRECISION da, db, dc, dq, df

c Handle the special cases - first if the source and background are zero

      IF ( dobs .EQ. 0 .AND. dbobs .EQ. 0 ) THEN

         wcntrb = dconv*dmod

c if source is non-zero but background is zero

      ELSEIF ( dobs .NE. 0 .AND. dbobs .EQ. 0 ) THEN

         IF ( dmod .LE. dobs/dtconv ) THEN

            wcntrb = - dmod*dbconv - dobs*log(dconv/dtconv)

         ELSE

            wcntrb = dmod*dconv - dobs*log(dmod*dconv)
     &                            - dobs*(1-log(dobs))

         ENDIF

c if source is zero and background is non-zero

      ELSEIF ( dobs .EQ. 0 .AND. dbobs .NE. 0 ) THEN

         wcntrb = dmod*dconv - dbobs*log(dbconv/dtconv)

c if the model is zero but the source and background are non-zero

      ELSEIF ( dmod .EQ. 0 .AND. dobs .NE. 0 .AND. dbobs .NE. 0 ) THEN

         wcntrb = dobs*log(dobs) + dbobs*log(dbobs) 
     &            - dobs*log(dconv*(dobs+dbobs)/dtconv)
     &            - dbobs*log(dbconv*(dobs+dbobs)/dtconv)

c now the standard solution with source, background, and model non-zero

      ELSE

c method from Numerical Recipes for solving a quadratic to avoid round-off
c error in the case when b^2 >> ac.

         da = dtconv
         db = dtconv*dmod - dobs - dbobs
         dc = -dbobs*dmod
         dq = -0.5*(db + SIGN(1.0d0,db)*SQRT( db*db - 4*da*dc ))
         df = dq/dtconv
         IF ( df .LT. 0.d0 ) df = dc/dq

         wcntrb = dconv*dmod + dtconv*df - dobs*(1-log(dobs)) 
     &                       - dbobs*(1-log(dbobs))

         IF ( (dconv*dmod+dconv*df) .GT. 0.d0 ) THEN
            wcntrb = wcntrb - dobs * log(dconv*dmod+dconv*df)
         ELSE
            wcntrb = wcntrb - dobs * (-32)
         ENDIF

         IF ( (dbconv*df) .GT. 0.d0 ) THEN
            wcntrb = wcntrb - dbobs * log(dbconv*df)
         ELSE
            wcntrb = wcntrb - dbobs * (-32)
         ENDIF

      ENDIF

      END