
      SUBROUTINE dopile(frame, npiled, g0, alpha, psffrac, nregions, 
     &                  ne, ear, ifl, photar, outar, tmpar, conv, pfrac)

      INTEGER ne, ifl, npiled, nregions
      REAL ear(0:ne), photar(ne), frame, g0, alpha, psffrac
      DOUBLE PRECISION outar(ne), tmpar(ne), conv(ne), pfrac(0:npiled)

c Subroutine to do the actual pile-up calculation. Split out from the
c main routine to make dynamic memory easier to handle. 
c
c Arguments :
c     frame      R       i: Frame time (in seconds)
c     npiled     I       i: Max. number of photons to pile up
c     g0         R       i: Grade 0 fraction
c     alpha      R       i: Index for parametrization of grade morphing
c     psffrac    R       i: Fraction of the flux to pile-up
c     nregions   I       i: Number of independent regions
c     ne         I       i: Number of energy bins
c     ifl        I       i: Dataset number (not used)
c     photar     R     i/r: Predicted counts/sec in each energy range
c     outar      D       w: Workspace array
c     tmpar      D       w: Workspace array
c     conv       D       w: Workspace array
c     pfrac      D       w: Workspace array

      DOUBLE PRECISION factor, pnorm, piledflux, sum, incounts
      DOUBLE PRECISION expfac, fract
      INTEGER ie, je, ipile, ioff

      CHARACTER wrtstr*255

c This is the fraction of the input counts that are used to calculate the 
c pile-up

      fract = psffrac/nregions

c This is an offset we will require in the convolution in case ear(0) is
c non-zero. Note that we assume the energy array is linear and starts at
c a multiple of the energy bin size.

      ioff = -ear(0)/(ear(1)-ear(0))

c Set the output array and the temporary arrays for the single photon case
c and find the sum of the input array

      pnorm = 0.d0
      DO ie = 1, ne
         outar(ie) = DBLE(photar(ie)*fract)
         tmpar(ie) = outar(ie)
         pnorm = pnorm + outar(ie)
      ENDDO

      incounts = pnorm * frame
      WRITE(wrtstr,*) 'Pileup: input model count rate/frame/region = ', 
     &           incounts
      CALL xwrite(wrtstr, 15)

c Calculate the exponential factor. 

      expfac = EXP(-1 * frame * pnorm / g0)

      pfrac(0) = 0.0d0
      pfrac(1) = pnorm * frame * expfac
      sum = pfrac(1)

c Loop round pile-ups

      factor = 1.0d0

      DO ipile = 2, npiled

c Do the convolution (should probably do an fft here but use BFI for now).

         DO ie = 1, ne
            conv(ie) = 0.
            DO je = 1, ie+ioff
               conv(ie) = conv(ie) + DBLE(photar(je)*fract)
     &                               *tmpar(ie-je+1+ioff)
            ENDDO
         ENDDO

c Add the latest iteration of the temporary array into the summation

         factor = factor * DBLE(alpha * frame/float(ipile))
         pfrac(ipile) = 0.d0
         DO ie = 1, ne
            tmpar(ie) = conv(ie)
            piledflux = tmpar(ie)*factor
            outar(ie) = outar(ie) + piledflux
            pfrac(ipile) = pfrac(ipile) + piledflux
         ENDDO

         pfrac(ipile) = pfrac(ipile) * frame * expfac
         sum = sum + pfrac(ipile)

      ENDDO

c Now multiply by the exponential term and add in the unpiled fraction

      DO ie = 1, ne
         piledflux = piledflux + outar(ie)*expfac
         photar(ie) = nregions*SNGL(outar(ie)*expfac) 
     &                + (1.-psffrac)*photar(ie)
      ENDDO

      WRITE(wrtstr,*) 'Pileup: output model count rate/frame/region = ', 
     &                piledflux * frame
      CALL xwrite(wrtstr, 15)

c Write out diagnostic information if required. The second number is the
c probability of seeing that number of photons piled up in a frame and the
c third number is the probability of an observed event being due to that
c number of piled up photons.

      IF ( sum .GT. 0.0 ) THEN
         WRITE(wrtstr,'(i4,3(1x,1pg12.6))') 0, expfac, pfrac(0), 
     &                                      pfrac(0)/sum
         CALL xwrite(wrtstr, 25)
         pnorm = expfac
         DO ipile = 1, npiled
            pnorm = pnorm * incounts / ipile
            WRITE(wrtstr,'(i4,3(1x,1pg12.6))') ipile, pnorm, 
     &                           pfrac(ipile), pfrac(ipile)/sum
            CALL xwrite(wrtstr, 25)
         ENDDO
      ENDIF

      RETURN
      END


