c $Id: calfak.f,v 1.11 2004/07/01 12:34:37 kaa Exp $

      SUBROUTINE CALFAK(phaobs, phabkg, varobs, varbkg, phamod, areasc)

      IMPLICIT NONE

      REAL phaobs(*), phabkg(*), varobs(*), varbkg(*)
      REAL phamod(*), areasc(*)

C---
C XSPEC subroutine to calculate fake data, similar to the old
C FAKEIT subroutine (but now restructured).  N.B.  This assumes
C that both data AND background are in the Gaussian regime.
C Because of the way INDATA calculates errors, negative counts
C are set to 0.
C---
C PHAOBS    I/O  On input initialized to zero or the expected
C                contribution of the correction file no. count rate
C                On output is the faked count rate.
C PHABKG    I    Input with the expected background contribution to
C                the count rate
C VAROBS    I/O  On input initialized to zero. On output is zero if Poisson
C                statistics are used and equal to the sigma if Gaussian are used.
C VARBKG    I    Input with errors on expected background contribution. This
C                array is only used with Gaussian statistics.
C PHAMOD    I    The current model count rate.
C AREASC    I    The AREASCAL values
C---
C 26 april 1985 - rashafer
C  7-Dec-1988 - Data and background calculated as Poisson deviates - kaa
C 17-Nov-1993 - Uses new improved random number generator. Also adds the
C               background to the model prediction then uses that as mean
C               for random number calculation rather than getting two separate
C               random numbers and adding them as before.
C 21-Jul-1997   Calculate RATE when counting statistics is not set  -han
C---
      REAL expos, grasp, gaunum
      INTEGER ifile, ipha
      LOGICAL qbgpoiss

      REAL dgtime, dgbtme
      INTEGER dgndst, dgnbnb, dgnbne
      LOGICAL dgqfst, dgdtps
      EXTERNAL dgndst, dgnbnb, dgnbne, dgqfst, dgtime, dgbtme, dgdtps

c Loop round the data files

      DO ifile = 1, dgndst()

c Get information for this file

         expos = dgtime(ifile)

c Calculate the model predicted counts

         DO ipha = dgnbnb(ifile), dgnbne(ifile)
            grasp = expos * areasc(ipha) 
            phaobs(ipha) = grasp * (phaobs(ipha) + phamod(ipha))
         ENDDO

c If the background is Poisson then add in background to get mean and 
c randomize if required

         IF ( dgbtme(ifile) .EQ. 0 .OR. dgdtps(ifile, 'b') ) THEN

            IF ( dgbtme(ifile) .NE. 0. ) THEN
               DO ipha = dgnbnb(ifile), dgnbne(ifile)
                  grasp = expos * areasc(ipha) 
                  phaobs(ipha) = phaobs(ipha) + grasp * phabkg(ipha)
               ENDDO
            ENDIF

c If the randomization is required then get Poisson numbers based
c on these means

            IF ( dgqfst() ) THEN
               CALL gtpoir( phaobs(dgnbnb(ifile)),
     &                     (dgnbne(ifile)-dgnbnb(ifile)+1) )
            ENDIF

c Convert back to a count rate
 
            DO ipha = dgnbnb(ifile), dgnbne(ifile)
               IF ( areasc(ipha) .NE. 0. ) THEN
                  grasp = expos * areasc(ipha)
                  phaobs(ipha) = phaobs(ipha) / grasp
               ENDIF
            ENDDO

c Else if the background is non-Poisson then do the source and background
c separately and sum them

         ELSEIF ( .NOT.dgdtps(ifile, 'b') ) THEN

c If the randomization is required then get Poisson numbers based
c on source means

            IF ( dgqfst() ) THEN
               CALL gtpoir( phaobs(dgnbnb(ifile)),
     &                     (dgnbne(ifile)-dgnbnb(ifile)+1) )
            ENDIF

c Set the variances

            DO ipha = dgnbnb(ifile), dgnbne(ifile)
               varobs(ipha) = phaobs(ipha)
            ENDDO

c Now add on the background. If randomization is required then we use Gaussian
c statistics (note that gtgaus draws a number from G(0,1))

            DO ipha = dgnbnb(ifile), dgnbne(ifile)
               gaunum = 0.0
               IF ( dgqfst() ) CALL gtgaus(gaunum, 1)
               grasp = expos * areasc(ipha) 
               phaobs(ipha) = phaobs(ipha) + 
     &                 grasp*(phabkg(ipha)+gaunum*sqrt(varbkg(ipha)))
               varobs(ipha) = varobs(ipha) + grasp*grasp*varbkg(ipha)
            ENDDO

c Convert back to a count rate
 
            DO ipha = dgnbnb(ifile), dgnbne(ifile)
               IF ( areasc(ipha) .NE. 0. ) THEN
                  grasp = expos * areasc(ipha)
                  phaobs(ipha) = phaobs(ipha) / grasp
                  varobs(ipha) = varobs(ipha) / grasp / grasp 
               ENDIF
            ENDDO

         ENDIF

c End loop round the files

      ENDDO

      RETURN
      END