      Subroutine CompAZGN(Ein,Eg,Pg)
C
C Generate photon energy from CompAZ parametrization
C
C  Input:  
C  ======
C             Ein   - electron beam energy [GeV]  50.<=Abs(Ein)<=500.
C
C  Output: 
C  =======
C
C  For Ein>0  Eg    - generated photon energy
C             Pg    - photon polarization for this energy
C
C  For Ein<0  Eg    - normalization of CompAZ spectra for E=ABS(Ein)
C             Pg    - not used
C
C
      Implicit None
C
      Real Ein,Eg,Pg
C
C Tables to store binned distribution
C      
      Integer Nbin,Nbin1
      Parameter (Nbin=1000, Nbin1=Nbin+1)

      Real Ebin(Nbin1)
      
      Real Prob(Nbin)
      Real Eslope(Nbin)
    
      Real Polar(Nbin)
      Real Pslope(Nbin)

      Integer Other(Nbin)

      Real Eold,Norm
      Data Eold / 0. /

      Save Ebin,Prob,Eslope,Polar,Pslope,Other
      Save Eold,Norm

      Real Ee,z,z2,xmin,xmax,xmax2
      Integer i,j,n1,n2,n3,imin,imax
      Real x,x1,x2,dx1,dx2,dx3,dx
      Real g1,g2,p1,p2,sum
      Real pmin,pmax,r,s
      Real rtab(3)

      Real Compaz
C
C Check energy
C
      Ee=ABS(Ein)

      If(Ee.LT.50. .OR. Ee.GT.500.)then

        write(6,*)'************************************************* '
        write(6,*)' '
        write(6,*)' '
        write(6,*)'CompAZGN  energy out of range (50.-500.) : ',Ee
        write(6,*)' '
        write(6,*)' '
        write(6,*)'************************************************* '

        STOP
C      ======

        Endif

C
C Initialize for new energy
C      

      If(Ee.NE.Eold)then

        write(6,*)'CompAZGN initialized for  Ee = ',Ee

        z=compaz(0.,Ee,20)

        xmax=compaz(0.,Ee,21)
        xmax2=compaz(0.,Ee,22)
        
        z2=z*(1-xmax)
        xmin=z2/(z2+1.) * (1-xmax)

        n1=nbin*xmin/xmax2
        dx1=xmin/n1

        n2=nbin*(xmax-xmin)/xmax2
        dx2=(xmax-xmin)/n2

        n3=nbin-n1-n2
        dx3=(xmax2-xmax)/n3

        x1=0.
        Ebin(1)=x1*Ee
        dx=dx1

        x2=x1+0.001*dx
        g1=compaz(x2,Ee,0)
        p1=compaz(x2,Ee,10)

        sum=0.

        Do i=1,nbin

          x2=x1+dx
          Ebin(i+1)=x2*Ee

          if(i.eq.n1+n2)x2=x2-0.001*dx

          g2=compaz(x2,Ee,0)
          p2=compaz(x2,Ee,10)

          Prob(i)=(g1+g2)*dx/2.
          Eslope(i)=(g2-g1)/(g2+g1)

          sum=sum+(g1+g2)*dx/2.

          Polar(i)=(p1+p2)/2.
          Pslope(i)=(p2-p1)/(p2+p1)

          x1=x1+dx
          g1=g2
          p1=p2

          if(i.eq.n1)dx=dx2

          if(i.eq.n1+n2)then
            dx=dx3
            x2=x1+0.001*dx
            g1=compaz(x2,Ee,0)
            p1=compaz(x2,Ee,10)
            endif

          EndDo
C
C Normalize
C
          Norm=sum
          write(6,*)'CompAZ normalization : ',Norm

          sum=sum/nbin

          Do i=1,Nbin
            Prob(i)=Prob(i)/sum
            Other(i)=-1
            EndDo
C
C Sort table
C           
        Do i=1,Nbin
          pmin=2.
          pmax=0.
          Do j=1,Nbin
             if(Other(j).lt.0.and.Prob(j).lt.pmin)then
               imin=j
               pmin=Prob(j)
               EndIf
             if(Other(j).lt.0.and.Prob(j).ge.pmax)then
               imax=j
               pmax=Prob(j)
               EndIf
             EndDo

          other(imin)=imax
          prob(imax)=prob(imax)-(1.-prob(imin))

          EndDo

        Eold=Ee

        EndIf

       If(Ein.lt.0.)then
          Eg=Norm
          Pg=0.
          Return
          EndIf


       call Ranmar(rtab,3)
       i=nbin*rtab(1) + 1
       if(i.gt.nbin)i=nbin

       if(rtab(2).gt.prob(i))i=other(i)

       r=rtab(3)
       s=Eslope(i)

       if(s.eq.0)then
          x=r
       Else
          x=(1-s)*(1-s)+4*s*r
          x=(s-1+sqrt(x))/2./s
          EndIf

       Eg=Ebin(i)+x*(Ebin(i+1)-Ebin(i))

       Pg=Polar(i)*(1+(2*x-1.)*Pslope(i))

       Return
       End


