      real function compazw(W,E,ic)

      implicit none
      real w,e
      integer ic

C
C  Calculated W distribution, by integrating compton spectra
C
C     ic=0          all events
C     ic=1          L=0
C     ic=2          L=2
C     ic=10,11,12   as above but using 'ideal' Compton spectra with same 'z'
C     ic=20,21,22   as above but using 'ideal' Compton spectra with same laser
C     ic=30,31,32   as above but using 'ideal' Compton spectra with z=4.3
C

      real z,xmax,xmax2,ymax,ymax2,dy1,dy2
      real x1,x2,y1,y2,g1,g2,sum,p1,p2
      integer i,n,n1,n2
      data n /100/

      real compaz,compazc,compazp

      if(ic.lt.20)then
          z=compaz(0.,E,20)
      elseif(ic.lt.30)then
          z=0.018088*E
      else
          z=4.8
C          z=4.3
          endif


      xmax=z/(1.+z)
      xmax2=2.*z/(1.+2.*z)

      ymax=log(2*E*xmax/W)
      ymax2=log(2*E*xmax2/W)

      compazw=0.
      if(ymax2.le.0.)Return     

      if(ymax.le.0. .and. ic.ge.10)Return     

      if(ymax.lt.0)ymax=0

      n1=n*ymax/ymax2+1
      n2=n-n1

      dy1=ymax/n1
      dy2=(ymax2-ymax)/n2

      sum=0

      y1=0

      p1=1.
      p2=1.

      x1=0.5*W/E*exp(y1)
      x2=0.5*W/E*exp(-y1)

      If(ic.lt.10)then
        g1=compaz(x1,E,0)*compaz(x2,E,0)
      else
        g1=compazc(x1,z,0.85)*compazc(x2,z,0.85)
        endif

      if(ic.eq.1)p1=0.5+0.5*compaz(x1,E,10)*compaz(x2,E,10)
      if(ic.eq.2)p1=0.5-0.5*compaz(x1,E,10)*compaz(x2,E,10)
      if(ic.eq.11 .or. ic.eq.21 .or. ic.eq.31)p1=
     +       0.5+0.5*compazp(x1,1.0,0.85,z)*compazp(x2,1.0,0.85,z)
      if(ic.eq.12 .or. ic.eq.22 .or. ic.eq.32)p1=
     +       0.5-0.5*compazp(x1,1.0,0.85,z)*compazp(x2,1.0,0.85,z)

      do i=1,n1

        y2=y1+dy1
        if(i.eq.n1)y2=y2-0.001*dy1

        x1=0.5*W/E*exp(y2)
        x2=0.5*W/E*exp(-y2)

        If(ic.lt.10)then
          g2=compaz(x1,E,0)*compaz(x2,E,0)
        else
          g2=compazc(x1,z,0.85)*compazc(x2,z,0.85)
          endif

        if(ic.eq.1)p2=0.5+0.5*compaz(x1,E,10)*compaz(x2,E,10)
        if(ic.eq.2)p2=0.5-0.5*compaz(x1,E,10)*compaz(x2,E,10)
        if(ic.eq.11 .or. ic.eq.21 .or. ic.eq.31)p2=
     +       0.5+0.5*compazp(x1,1.0,0.85,z)*compazp(x2,1.0,0.85,z)
        if(ic.eq.12 .or. ic.eq.22 .or. ic.eq.32)p2=
     +       0.5-0.5*compazp(x1,1.0,0.85,z)*compazp(x2,1.0,0.85,z)


        sum=sum+(g1*p1+g2*p2)*dy1/2.

        y1=y1+dy1
        g1=g2
        p1=p2

        enddo

      y2=y1+0.001*dy2
      x1=0.5*W/E*exp(y2)
      x2=0.5*W/E*exp(-y2)

      If(ic.lt.10)then
        g1=compaz(x1,E,0)*compaz(x2,E,0)
      else
        g1=compazc(x1,z,0.85)*compazc(x2,z,0.85)
        endif


      if(ic.eq.1)p1=0.5+0.5*compaz(x1,E,10)*compaz(x2,E,10)
      if(ic.eq.2)p1=0.5-0.5*compaz(x1,E,10)*compaz(x2,E,10)
      if(ic.eq.11 .or. ic.eq.21 .or. ic.eq.31)p1=
     +       0.5+0.5*compazp(x1,1.0,0.85,z)*compazp(x2,1.0,0.85,z)
      if(ic.eq.12 .or. ic.eq.22 .or. ic.eq.32)p1=
     +       0.5-0.5*compazp(x1,1.0,0.85,z)*compazp(x2,1.0,0.85,z)

      do i=1,n2

        y2=y1+dy2
        if(i.eq.n2)y2=y2-0.001*dy2

        x1=0.5*W/E*exp(y2)
        x2=0.5*W/E*exp(-y2)

        If(ic.lt.10)then
          g2=compaz(x1,E,0)*compaz(x2,E,0)
        else
          g2=compazc(x1,z,0.85)*compazc(x2,z,0.85)
          endif


        if(ic.eq.1)p2=0.5+0.5*compaz(x1,E,10)*compaz(x2,E,10)
        if(ic.eq.2)p2=0.5-0.5*compaz(x1,E,10)*compaz(x2,E,10)
        if(ic.eq.11 .or. ic.eq.21 .or. ic.eq.31)p2=
     +       0.5+0.5*compazp(x1,1.0,0.85,z)*compazp(x2,1.0,0.85,z)
        if(ic.eq.12 .or. ic.eq.22 .or. ic.eq.32)p2=
     +       0.5-0.5*compazp(x1,1.0,0.85,z)*compazp(x2,1.0,0.85,z)

        sum=sum+(g1*p1+g2*p2)*dy2/2.

        y1=y1+dy2
        g1=g2
        p1=p2

        enddo

      compazw=  w/E/E * sum

      return
      end


