c ------------------------------------------------------------------------
c
c  ** NOTE - FOR DETECTORS WITH LIMITED RESPONSES, THE ARRAY OF
c     ENERGIES SHOULD BE EXTENDED WITH THE 'ENERGIES' COMMAND. 
c
c   FOR EXAMPLE, WITH RXTE:
c   >    energies 0.2 50 500 lin 
c   
c
c   This routine is designed to give a physically motivated
c   prescription to model Comptonization of an input seed spectrum.
c   Photon number is preserved in the output.  
c
c     USER INPUT:
c
c     param(1) - Photon index of the Comptonized component in the spectrum. 
c     param(2) - Normalization of the spectrum. 
c     param(3) - Flag - >0 for up-scattering only (SIMPL-1), <=0 for both
c                 up-scattering and down-scattering (SIMPL-2).
c------------------------------------------------------------------------
c     XSPEC INPUT
c     
c     2 Vectors:
c     -----------------
c     1) Energy
c     2) Flux (Energy)
c
c
c     EAR(0:NE) - array of observed energy bins in keV
c     NE        - # observed energy bins
c     PARAM     - list of parameters
c     PHOTAR() - an array of photon flux in each bin of energy PHOTAR(NE)
c     PHOTER() - formal output for errors
c
c
c
c
c  EXAMPLE:
c      model phabs(simpl(diskbb))
c
c==========================================================================

      subroutine simpl(ear,ne,param,IFL,photar,photer) 

      integer i,j,k
      double precision norm,gamma,gnormUP,gnormDN,gamma1,gamma2
      real ear(0:ne),param(3),photar(ne),photer(ne)
      real tmparr(ne),enavgam1(ne),enavgam2(ne)
      real engam1(0:ne),engam2(0:ne)


c--------------------------
      INCLUDE 'xspec.inc'
c--------------------------

      
      norm=param(2)
      gamma=param(1)
      switch=param(3)
      
      
c     WORKAROUND to avoid pole at gamma=1
      IF (gamma .eq. 1) then 
         gamma=1.001
      ENDIF

      gamma1=gamma-1.
      gamma2=gamma+2.
      
c     Initialize arrays
      
      do 10 i=0,ne
         engam1(i)=ear(i)**(-gamma1)
         if (i.ne.0) then
            tmparr(i)=0.
            enavgam1(i)=(0.5*(ear(i-1)+ear(i)))**gamma1
         endif
 10   continue
      
      
      IF (switch .gt. 0) THEN   !UP-SCATTERING ONLY
c     Loop over energy bins
         do 20 i=1,ne
            
c     Calculate scattered contribution to the same bin
            tmparr(i)=tmparr(i)+photar(i)*(1.-enavgam1(i)*engam1(i))
            
c     Loop over all higher bins and calculate scattered contributions
            do 20 j=i+1,ne
               tmparr(j)=tmparr(j)+enavgam1(i)*(engam1(j-1)-engam1(j))
     $              *photar(i)
 20         continue
            
            
         ELSE IF (switch .le. 0) THEN !UP-SCATTERING & DOWN-SCATTERING 
            
            gnormUP=(gamma+2.)/(1.+2.*gamma)
            gnormDN=(gamma-1.)/(1.+2.*gamma)
     
            do 30 k=0,ne
               engam2(k)=ear(k)**(gamma2)
               if (k .ne. 0) then
                  enavgam2(k)=(0.5*(ear(k-1)+ear(k)))**(-gamma2)
               endif
 30         continue
            
c     Loop over energy bins
            do 40 i=1,ne
c     Calculate scattered contribution to the same bin
               tmparr(i)=tmparr(i)+photar(i)*((1.-enavgam1(i)*engam1(i))
     $              *gnormUP+gnormDN*(1.-engam2(i)*enavgam2(i) ))
c     Loop over all bins and calculate scattered contributions
           do 40 j=1,ne
              if (j .lt. i) then
                 tmparr(j)=tmparr(j)+enavgam2(i)*(engam2(j)-engam2(j-1))
     $                *photar(i)*gnormDN
              else if (j .gt. i) then
                 tmparr(j)=tmparr(j)+enavgam1(i)*(engam1(j-1)-engam1(j))
     $                *photar(i)*gnormUP
              endif
 40        continue
               
        ENDIF
                  
C     ----------------------------------------
C     OUTPUT TO CHECK THE SCRIPT - uncomment to use
C     ----------------------------------------
C      open(unit=90,file='sop.out')      
C      do i=1,300
C         write(90,*) i,ear(i),norm,gamma,tmparr(i),photar(i)
C      enddo
C      CLOSE(90)
C     ----------------------------------------
      do i = 1, ne
         photar(i) = (1.-norm)*photar(i)+norm*tmparr(i)
      enddo

      RETURN
      
      END