SUBROUTINE xsgaul(ear, ne, param, ifl, photar, photer)

      INTEGER ne, ifl
      REAL ear(0:ne), param(2), photar(ne), photer(ne)

C---
C XSPEC model subroutine
C Simple gaussian line profile
C---
C see ADDMOD for parameter descriptions
C number of model parameters:2
C       1       EL      line energy (in energy units, e.g. keV)
C       2       W	line width (sigma) (in energy units)
C intrinsic energy range: none
C       N.B. the line may have significant flux at unphysical (i.e. negative)
C            energies when EL<~W.
C algorithm:
C       n(E) = (1./(sqrt(2pi*W*W)) * exp(-0.5*((E-EL)/W)**2)) dE
C       N.B. the norm is equal to the total integrated number of counts in
C            the line.
C       If W is less than or equal to zero, the line is treated as a delta
C            function.
C       N.B. when the energy spacing is much greater than W and the stepsize
C            for EL then the partial derivative determinations can lead
C            to fit error conditions.  The solution is to increase the
C            stepsize for EL
C---
C 17 Dec 1983 - rashafer
C 1.1: An attempt to use simple linear interpolation for the
C      response to better handle very narrow lines.  This
C      is not strictly correct, and effectively the lines will have
C      widths of order of the ear spacing at that point.
C 4/19/98  kaa  Reverted to putting whole line in a single bin if the
C               width is small. This means for lines with width < ear
C               spacing the energy can only be determined to the ear
C               spacing. This change is necessary to determine line widths
C               correctly using high resolution spectroscopy.
C---

      REAL crtsig
      PARAMETER (crtsig=6.0)
      REAL el, w, emin, emax, winv, alow, erf, ahi
      REAL emidie, emidje, delta
      INTEGER ie, ielow, je

c Initialize the model array

      DO ie = 1, ne
         photar(ie) = 0.
      ENDDO

c Calculate the line

      CALL gaussline(ear, ne, param(1), param(2), 1.0, crtsig, 
     &               .FALSE., photar)

      RETURN
      END

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

      SUBROUTINE gaussline(ear, ne, ecenter, width, norm, crtsig, 
     &                     qspeedy, photar)

      IMPLICIT NONE

      INTEGER ne

      REAL ear(0:ne), photar(ne)
      REAL ecenter, width, norm, crtsig

      LOGICAL qspeedy

c Add the input gaussian line to the photar array. This is broken out from
c the xsgaul function to make it available to other models which may need
c to calculate gaussian line shapes. If qspeedy is true then uses a tabulation
c otherwise does a proper calculation of the erf function.

      REAL emin, emax, winv, alow, ahi
      REAL emidie, emidje, delta
      INTEGER ie, ielow, je, ien

      DOUBLE PRECISION energyfrac
      REAL erf
      INTEGER bisearch
      EXTERNAL erf, bisearch

c first trap case of zero width line

      IF ( width .LE. 0. ) THEN
         ie = bisearch(ne, ear, ecenter)
         IF ( ie .GT. 0 ) photar(ie) = photar(ie) + norm
         RETURN
      ENDIF

c Set the line energy, width and the energy range over which the
c line will be calculated (essentially +/- crtsig sigma)

      width = max(0., width)
      emin = ecenter - crtsig*width
      emax = ecenter + crtsig*width

c If the line calculation energies lie outside the response energy range
c then give up

      IF ( (emax.LT.ear(0)) .OR. (emin.GE.ear(ne)) ) RETURN


c This is the speedy method using tabulated gaussian data

      IF ( qspeedy ) THEN

c Find the bin with the line center. If line center is below first bin then
c don't calculate lower part of the line. If line center is above first bin
c then just calculate part of line within energy range.

         IF ( ecenter .LT. ear(0) ) THEN
            ien = 0
         ELSEIF ( ecenter .GT. ear(ne) ) THEN
            ien = ne
         ELSE
            ien = bisearch(ne, ear, ecenter)
         ENDIF

c Do the low energy part of the line

         ielow = ien
         alow = 0
         DO WHILE ( ielow .GE. 1 )
            ahi = energyfrac(ecenter, width, ear(ielow-1))
            IF ( ahi .GE. 0. ) THEN
                photar(ielow) = photar(ielow) +
     &                   (ahi-alow)*0.5*norm
            ELSE

c Too many sigma away so stop now and add the rest of the Gaussian into this
c bin. Not strictly correct but shouldn't matter and ensures that the total
c flux is preserved.

               photar(ielow) = photar(ielow) +
     &                      (1.-alow)*0.5*norm
               ielow = 1
            ENDIF
            alow = ahi
            ielow = ielow - 1
         ENDDO

c Find the bin with the line center. If line center is above the last bin then
c don't calculate upper part of the line. If line center is below first bin
c then just calculate part of line within energy range.

         IF ( ecenter .LT. ear(0) ) THEN
            ien = 1
         ELSEIF ( ecenter .GT. ear(ne) ) THEN
            ien = ne + 1
         ELSE
            ien = bisearch(ne, ear, ecenter)
         ENDIF

c Do the high energy part of the line

         ielow = ien
         alow = 0
         DO WHILE ( ielow .LE. ne )
            ahi = energyfrac(ecenter, width, ear(ielow))
            IF ( ahi .GE. 0. ) THEN
               photar(ielow) = photar(ielow) +
     &                 (ahi-alow)*0.5*norm
            ELSE

c Too many sigma away so stop now and add the rest of the Gaussian into this
c bin. Not strictly correct but shouldn't matter and ensures that the total
c flux is preserved.

               photar(ielow) = photar(ielow) +
     &                  (1.-alow)*0.5*norm
               ielow = ne
            ENDIF
            alow = ahi
            ielow = ielow + 1
         ENDDO

c and this is the slower method using erf calculation.

      ELSE

         winv = 1./width/sqrt(2.)

c If there is part of the gaussian below the response energy range then
c store that in alow.

         ielow = 0
         IF (emin.LE.ear(0)) THEN
            alow = 0.5*erf(winv*(ear(0)-ecenter))
            ielow = -1
            ie = 1
         ELSE
            alow = -0.5
         ENDIF

c find the start of the line

         IF (ielow.EQ.0) ielow = bisearch(ne, ear, emin)
         IF (ielow .NE. -1) ie = ielow

c integrate the gaussian for the energy ranges covering the line

         DO WHILE (ear(ie).LT.emax)
            ahi = 0.5*erf(winv*(ear(ie)-ecenter))
            photar(ie) = photar(ie) + norm*(ahi - alow)
            alow = ahi
            ie = ie + 1
            IF (ie.GT.ne) RETURN
         ENDDO

c if all of the model was in a single bin

         IF (ie.EQ.ielow) THEN
            photar(ie) = photar(ie) + norm
         ELSE
            photar(ie) = photar(ie) + norm*(0.5 - alow)
         ENDIF

      ENDIF

      RETURN
      END

c ****************************************************************
      FUNCTION energyfrac(linecentre, linewidth, energy)

      IMPLICIT NONE

      DOUBLE PRECISION energyfrac
      REAL linecentre, linewidth, energy

c Function to return amount of energy emitted between line centre and energy

      DOUBLE PRECISION sigma
      INTEGER minindex
      DOUBLE PRECISION remainder

      INTEGER GAUSSMAXSTEP
      DOUBLE PRECISION GAUSSMAX, GAUSSSTEP
      PARAMETER (GAUSSMAXSTEP=30)
      PARAMETER (GAUSSMAX=3)
      PARAMETER (GAUSSSTEP=GAUSSMAX*1.0/GAUSSMAXSTEP)

      DOUBLE PRECISION intGaus(0:GAUSSMAXSTEP)

c Integrated Gaussian probability distribution (x=0->3 in steps of 0.1)

      DATA intGaus /0, 0.07965567455405775, 0.1585194188782062,
     & 0.23582284437790513, 0.3108434832206486, 0.38292492254802624, 
     & 0.4514937644998527, 0.5160726955538542, 0.5762892028332065,
     & 0.6318797493064812, 0.6826894921370861, 0.7286678781072344,
     & 0.7698606595565836, 0.8063990308287794, 0.8384866815324581,
     & 0.8663855974622838, 0.8904014166008838, 0.910869074482914,
     & 0.9281393617741482, 0.9425668803679965, 0.9544997361036416,
     & 0.9642711588743671, 0.9721931049730028, 0.9785517799566483,
     & 0.9836049281508079, 0.9875806693484477, 0.9906776239525625,
     & 0.9930660523939188, 0.9948897393391443, 0.9962683733992319,
     & 0.9973002039367398/

c Compute how many sigmas away we are
      sigma = abs( (energy-linecentre)/linewidth )

c Now interpolate from the table
      minindex = INT(sigma/GAUSSSTEP)

c If we're past the edge of the tabulated data return -1.
      IF ( minindex .GE. GAUSSMAXSTEP ) THEN
         energyfrac = -1
         RETURN
      ENDIF

      remainder = (sigma - DBLE(minindex)*GAUSSSTEP) * (1./GAUSSSTEP)

c Do the interpolation
      energyfrac = (1.-remainder)*intGaus(minindex) +
     &             remainder*intGaus(minindex+1)

      RETURN
      END

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

      FUNCTION bisearch(ne, ear, energy)

      INTEGER bisearch, ne
      REAL ear(0:ne), energy

c Function to do a bisection search for which element of the energy array
c the input energy lies in.

      INTEGER low, high


      IF ( energy .LT. ear(0) .OR. energy .GT. ear(ne) ) THEN
         bisearch = -1
         RETURN
      ENDIF

      low = 0
      high = ne

      DO WHILE ( (high-low) .GT. 1 )

         bisearch = (low + high) / 2
         IF ( energy .GT. ear(bisearch-1) ) THEN
            low = bisearch
         ELSE
            high = bisearch
         ENDIF

      ENDDO

      bisearch = low

      RETURN
      END