      Real Function Compaz(x,E,icase)
C
      Implicit none
C
      Real x,E
      Integer icase
C
C Input parameters:
C -----------------
C
C     x               Fraction of beam energy carried by photon  ( 0 < x < 1 )
C
C     E               Electron beam energy [GeV]
C
C     Icase           Choice of output distribution:
C
C     Icase = 0       Full spectra - sum of three contributions
C             1       Main contribution from direct Compton process
C             2       Contribution from scattering with 2 photons
C             3       Contribution from scattering on secondary electrons
C
C            10       Polarisation of the full spectra
C            11       Polarisation of the 1st component (direct Compton)
C            12       Polarisation of the 2nd component (two photon)
C            13       Polarisation of the 3rd component (secondary electron)
C
C            20       Effective z parameter
C            21       maximum x value for direct Compton
C            22       maximum x value for scattering with 2 photons
C
C========================================================================
C
C Parameters of the distribution, partly fitted to the spectra from V.Telnov
C
      Real par(14)
      Data par /
     +   18.08800    ,                              ! Laser parameter
     +    0.06594662 ,   0.7060851   ,              ! Nonlinear corrections
     +    0.7115863  ,  -0.6776124   ,  0.000000 ,  ! Total normalization
     +    0.5540019  ,  37.38912     ,              ! 2nd component fraction
     +    0.7257064  ,   1.517959    ,              ! 3rd component fraction
     +    3.315865   ,   1.387423    ,              ! Corrections
     +    0.85       ,   1.00        /              ! Polarizations
C
C========================================================================
C
      Real Ee,z0,ksi,z
      Real Tot,Tot2,Tot3
      Real Br2,gpow,pole,poll,pol,pole2
C
      Real xmax,xmax2,g
      Real gflux1,gflux2,gflux3,gflux
      Real Pol1,Pol2,Pol3,Ptot
      Real g1,g2,dg,p1,p2,dp
      Real ymin,ymin2,ymax,y1,y2,dy
C
      Logical DoPolar
C
      integer iy,ny
      Parameter (ny=100)
C
      Real compazc,compazp
C
C=========================================
C
C Check input flag
C
      If(.NOT.((Icase.GE.0  .and. Icase.LE.3 ) .OR.
     +         (Icase.GE.10 .and. Icase.LE.13) .OR.
     +         (Icase.GE.20 .and. Icase.LE.22))     )Then
         write(6,*)' COMPAZ: Invalid distribution flag : ',icase
         Compaz=0.
         Return
         Endif
C
C Check input energy
C       
      If(E.LT.50. .OR. E.GT.500.)then
        write(6,*)' COMPAZ : Beam energy outside range (50.-500.): ',E
        Compaz=0.
        Return
        Endif
C
C Polarisation calculation flag
C
      DoPolar =  (Icase.GE.10 .and. Icase.LE.13) 
C
C=========================================
C
C  Set parameter values for given energy
C  -------------------------------------
C
C Beam energy in TeV
C
      Ee = e/1000.
C
C Nominal value of z = 4*Ee*El/me^2
C
      z0=par(1)*ee
C
C correction nonlinear effects
C
      ksi=par(2)+par(3)*ee
C
C effective z
C
      z=z0/(1+ksi)
C
C Maximum x values for different contributions
C
      xmax=z/(1.+z)
      xmax2=2.*z/(1.+2.*z)
C
      If(Icase.EQ.20)Then
         Compaz=z
         Return
      ElseIf(Icase.EQ.21)Then
         Compaz=xmax
         Return
      ElseIf(Icase.EQ.22)Then
         Compaz=xmax2
         Return
         Endif
C
C-------------------------------------------------------------
C
C normalization of main compton component
C
      Tot=par(4)+ee*par(5)+ee*ee*par(6)
C
C normalization of additional contributions
C
      Tot2=Tot*par(7)*(1-exp(-ksi*ksi*par(8)))
      Tot3=Tot*(par(9)+ee*par(10))
C
C Corrections to Compton spectra
C
      Br2=par(11)
      gpow=par(12)
C
C Electron/laser polarisation
C
      Pole=par(13)
      Poll=par(14)
C
      Pol=Pole*Poll
C
C==========================================
C
C Leading Compton term
C --------------------
C
      gflux1=0.
      pol1=0.

      If(x.gt.0. .and. x.le.xmax .and. Tot .gt. 0.)then

        gflux1=compazc(x,z,Pol)
        g=(z/x - z -1)
        if(g.gt.0.)then
          g=sqrt(g)
        else
          g=0.
          endif

        gflux1=gflux1*exp( -Br2 * g*g/8.)

        If(DoPolar) pol1=compazp(x,Poll,Pole,z)

        EndIf
C
      If(Icase.eq.1)then
         compaz=Tot*gflux1
         return
      ElseIf(Icase.eq.11)then
         compaz=pol1
         return
         endif
C
C------------------------------------------------------------------
C 
C Scattering on 2 photons
C -----------------------
C
      gflux2=0.
      pol2=0.
      
      If(x.gt.0. .and. x.lt.xmax2 .and. Tot2 .gt. 0.)then

        gflux2=compazc(x,2.*z,Pol)
        g= 2.*z/x - 2.*z -1
        if(g.gt.0.)then
          g=sqrt(g)
        else
          g=0.
          endif

        gflux2=gflux2*exp( -Br2 * g*g/8.)* g**gpow

        If(DoPolar .and. x.gt.xmax)then
           Pol2=compazp(x,poll,Pole,2.*z)
        elseIf(DoPolar)then
           Pol2=compazp(xmax,poll,Pole,2.*z)
           Endif
C
        EndIf
C 
      If(Icase.eq.2)then
         compaz=Tot2*gflux2
         return
      ElseIf(Icase.eq.12)then
         compaz=pol2
         return
         endif
C
C----------------------------------------------------------------
C
C Electron rescattering
C ---------------------
C
      gflux3=0.
      pol3=0.
      
      If(x.gt.0. .and. x.le.xmax .and. Tot3 .gt. 0.)then

        ymin=0.5*x*(1+sqrt(4./x/z + 1.))

        ymin2=1./(z+1)
        if(ymin2.gt.ymin)ymin=ymin2

        ymax=1.0

        If(ymax.gt.ymin)then

          Pole2= 0.

          y1=ymin
          y2=y1
C
C To avoid numerical problems
C
          y1=y1*1.0001

          g1= log(y1*z+1.)/log(z+1.)/y1  * 
     +      compazc(x/y1,y1*z,pole2*Poll)*compazc(1-y1,z,Pol)

          If(DoPolar) p1=compazp(x/y1,poll,Pole2,y1*z)
C
          dy=(ymax-ymin)/ny

          Do iy=1,ny       
          
            y2=y2+dy
            g2=log(y2*z+1.)/log(z+1.)/y2  * 
     +        compazc(x/y2,y2*z,pole2*poll)*compazc(1-y2,z,Pol)
            
            dg = (g1+g2)*dy/2.
            dp = (g1*p1+g2*p2)*dy/2.

            gflux3=gflux3 + dg

            If(DoPolar)then
               p2=compazp(x/y2,poll,Pole2,y2*z)
               dp = (g1*p1+g2*p2)*dy/2.

               Pol3=Pol3 + dp

               p1=p2
               EndIf

            y1=y2
            g1=g2

            enddo

          EndIf

        EndIf
C
      If(Icase.eq.3)then
         compaz=Tot3*gflux3
         return
      ElseIf(Icase.eq.13 .and. gflux3.gt.0.)then
         compaz=pol3/gflux3
         return
      ElseIf(Icase.eq.13)then
         compaz=0.
         return
         endif
C
C----------------------------------------------------------------
C
       Gflux=Tot*gflux1 + Tot2*gflux2 + Tot3*gflux3
C
       If(Icase.EQ.0)then
          compaz=Gflux
       Elseif(Gflux.GT.0. .and. Icase.eq.10)then
          Ptot=Tot*gflux1*pol1 + Tot2*gflux2*pol2 + Tot3*pol3
          compaz=Ptot/Gflux
       Else
          compaz=0.
          Endif 
C
       Return
       End
C
C
C*******************************************************
C
      Real Function Compazc(x,z,Pol)
C
C 'Ideal' spectrum of compton backscattered photons
C
      Implicit none
      Real x,z,Pol
C
      Real xmax,r,y,gflux,norm,g,corr
C
      xmax=z/(1.+z)
C
      If(x.lt.0. .or. x.gt.xmax)then
        Compazc=0.
        Return
        EndIf
C
C Flux normalisation (corrected!  formula in paper is wrong)
C
      Norm = (z*z*z + 18*z*z + 32*z + 16)/(2*z*(z+1)*(z+1)) +
     +       (1 - 4/z - 8/z/z)*log(1+z) +
     +   Pol*( 2  +  z*z/(2*(z+1)*(z+1)) - (1+2/z)*log(1+z) )
C
C
C Additional variables, for easier calculation
C
         r = x/(1-x)
         y = 1 - x + 1/(1-x)
C
C Photon flux, as given in hep-ph/9508397
C
         gflux= y + 4*r/z*(r/z-1) + Pol*r*(2-x)*(2*r/z-1)
C
         Compazc=gflux/norm
C
         Return
         End

C---------------------------------------------------------
C
      Real Function Compazp(x,Pl,Pb,z)
C
C Calculate spectrum of compton backscattered photons
C
      Implicit none
      Real x,Pl,Pb,z
C
      Real xmax,r,y,gflux,gpol,Pol
C
      xmax=z/(1.+z)
C
      If(x.lt.0. .or. x.gt.xmax)then
        Compazp=0.
        Return
        EndIf
C
      Pol = Pl*Pb
      r = x/(1-x)
      y = 1 - x + 1/(1-x)
C
C Photon polarised flux
C
         gpol=Pb*r*(1+(1-x)*(2*r/z - 1)**2) + Pl*y*(2*r/z-1)
C
C Photon flux
C
         gflux= y + 4*r/z*(r/z-1) + Pol*r*(2-x)*(2*r/z-1)
C
         Compazp=gpol/gflux
C
         Return
         End
