subroutine comp(E,nu)
      real E,nu
      include 'comp.inc'
      Ebeam=E
      Pbeam=sqrt(Ebeam**2-me**2)
      Ki=nu
c     Ki=2*2.4e-9
c     Ki=3.43e-9
      kin=(Ebeam+Pbeam)*Ki/me
      end

      real function Elab(th)
      real th
      include 'comp.inc'
      thbeam=th*pi/180
      kout=kin/(1+(kin/me)*(1-cos(thbeam)))
      Kf=(Ebeam-Pbeam*cos(thbeam))*kout/me
      thlab=atan2(kout*sin(thbeam),Kf)
      Elab=Kf
      end

      real function tlab(th)
      real th
      include 'comp.inc'
      thbeam=th*pi/180
      kout=kin/(1+(kin/me)*(1-cos(thbeam)))
      Kf=(Ebeam-Pbeam*cos(thbeam))*kout/me
      thlab=atan2(kout*sin(thbeam),Kf)
      tlab=thlab*Ebeam/me
      end

      real function tbeam(E)
      real E
      include 'comp.inc'
      real Z
      Z=1-E/(Pbeam-E)*(me/kin)
      if (abs(Z).lt.1) then
        thbeam=acos(Z)
      else if (Z.gt.1) then
        thbeam=0
      else
        thbeam=pi
      endif
      kout=kin/(1+(kin/me)*(1-cos(thbeam)))
      Kf=(Ebeam-Pbeam*cos(thbeam))*kout/me
      thlab=atan2(kout*sin(thbeam),Kf)
      tbeam=thbeam*180/pi
      end 

      real function dsigdE(E)
      real E
      include 'comp.inc'
      if (E.gt.Elab(180.)) then
        dsigdE=0
        return
      endif
      call tbeam(E) 
      dsigdE=2*pi*Compton(kin,thbeam)
     &      /abs((Ebeam-Pbeam*cos(thbeam))*(kout/me)**2
     &          -(kout*Pbeam/me))
      end

      real function dNgdE(E)
      real E
      include 'comp.inc'
      real P,G,rB,rC,I,L
      integer N
      if (E.gt.Elab(180.)) then
        dNgdE=0
        return
      endif
      call tbeam(E) 
      P=1000    !laser power in Watts (peak times duty factor)
      G=1       !cavity gain factor
      tau=2e-12 !laser pulse length (s) times crossing factor
      rC=10e-6  !neck radius of cavity beam in m
      rB=10e-6  !electron beam radius in m
      I=1.0e-6  !electron beam current in A
      L=0.0450  !effective length of cavity (m)
      N=1       !number of passes through beam
c
c this is the formula for the c.w. cavity case
c     dNgdE=N*P*G/(Ki*1.6e+7)*dsigdE(E)
c    +      /(2*pi*(rB**2+rC**2))
c    +      *(I/1.6e-2)
c    +      *(2*L/3e8)
c
c this is the formula for the pulsed cavity case
      dNgdE=N*P*G/(Ki*1.6e+7)*dsigdE(E)
     +      /(2*pi*(rB**2+rC**2))
     +      *(I/1.6e-2)
     +      *tau
      end

*DECK,Compton.
*CMZ :          09/11/97  22.57.28  by  Rick Jones
*-- Author :  Richard Jones
      Real Function Compton(kin,theta)
      Implicit none
      Real kin,theta
 
C   Purpose: This function calculates the cross section for Compton scattering
C      to first order in alpha.  The incoming electron is taken to be at rest
C      and the incoming photon directed along the +z axis, the respective
C      polarizations being specified by initial-state spin-density matrices.
C      The final state is specified by the scattering angle theta of the
C      outgoing photon, taken to lie in the x-z plane, and the spin-density
C      matrices of the final-state particles. The meaning of a "final-state
C      spin-density matrix" supplied as an input is the specification of the
C      way the sum over final-state helicities is done to calculate the cross
C      section.  The cross section is differential in
C                     d(Omega) = d(cos(theta)) d(phi).
C      The formalism for this program is taken from Bjorken and Drell
C      "Relativistic Quantum Mechanics", p.127-130.
C
C   Explicit Inputs:  kin = initial photon momentum (GeV/c)
C                       in lab frame where electron is initially at rest,
C                       nominal direction is along z (GeV/c)
C                     theta = photon scattering angle in lab frame (radians)
C   Explicit Outputs: Compton = differential cross section in lab angles
C                       d(sigma)/d(Omega) expressed in microbarns
C
C   Implicit Inputs:  elDMi(2,2) = incoming electron spin-density matrix, where
C                         ( 1  0 )      ( 0  0 )        ( 1/2  0  )
C                         ( 0  0 ) ,    ( 0  1 ) , and  (  0  1/2 )
C                       means polarization along +z axis, -z axis, and
C                       unpolarized, respectively.
C                     phDMi(2,2) = incoming photon spin-density matrix, with the
C                       same definition as for the electron above.
C                     elDMf(2,2) = final electron spin-density matrix, where
C                         ( 1  0 )      ( 0  0 )        ( 1  0 )
C                         ( 0  0 ) ,    ( 0  1 ) , and  ( 0  1 )
C                       pick out the positive, negative, and both helicity
C                       components in the final state.
C                     phDMf(2,2) = final photon spin-density matrix, with
C                       same definition as for the final-state electron above.
C   Implicit Outputs: /Comp/ common, for details see below.
C========================================================================
 
C   Global Declarations
C-------------------------
 
C Common blocks for storing detailed results
      Common /Ckin/ Pfi(5),kout(4)
      Common /Cspin/ elDMi(2,2),phDMi(2,2),elDMf(2,2),phDMf(2,2)
      Real Pfi,kout
      Complex elDMi,phDMi,elDMf,phDMf
 
C Common blocks for storing standard matrices
      Complex i_,_i
      Parameter (i_=(0.,1.),_i=(0.,-1.))
      Common /Cdirac/ gamma(4,4,4),J(4,4,3),K(4,4,3)
      Complex gamma,J,K
      Data gamma / 0, 0, 0, 1,
     &             0, 0, 1, 0,
     &             0,-1, 0, 0,
     &            -1, 0, 0, 0,
     &                         0, 0, 0,_i,
     &                         0, 0, i_,0,
     &                         0, i_,0, 0,
     &                        _i, 0, 0, 0,
     &                                     0, 0, 1, 0,
     &                                     0, 0, 0,-1,
     &                                    -1, 0, 0, 0,
     &                                     0, 1, 0, 0,
     &                                                 1, 0, 0, 0,
     &                                                 0, 1, 0, 0,
     &                                                 0, 0,-1, 0,
     &                                                 0, 0, 0,-1   /
      Data J / 0, 1, 0, 0,
     &         1, 0, 0, 0,
     &         0, 0, 0, 1,
     &         0, 0, 1, 0,
     &                     0,_i, 0, 0,
     &                     i_,0, 0, 0,
     &                     0, 0, 0,_i,
     &                     0, 0, i_,0,
     &                                 1, 0, 0, 0,
     &                                 0,-1, 0, 0,
     &                                 0, 0, 1, 0,
     &                                 0, 0, 0,-1   /
      Data K / 0, 0, 0, i_,
     &         0, 0, i_,0,
     &         0, i_,0, 0,
     &         i_,0, 0, 0,
     &                     0, 0, 0, 1,
     &                     0, 0,-1, 0,
     &                     0, 1, 0, 0,
     &                    -1, 0, 0, 0,
     &                                 0, 0, i_,0,
     &                                 0, 0, 0,_i,
     &                                 i_,0, 0, 0,
     &                                 0,_i, 0, 0   /
 
C   Local Declarations
C------------------------
      Real me
      Parameter (me=511e-6)
      Real alpha,hbarc
      Parameter (alpha=7.3e-3,hbarc=0.197)
      Double Precision pi
      Parameter (pi=3.1415926535898)
 
      Complex epsin(4,2),epsout(4,2)
      Real ome,eta
      Complex U(4,4),V(4,4)
      Complex W(4,4,2,2)
      Complex trace,sum
      Real seta,ceta
      Real come,some
      Integer mu,nu
      Integer si,sii,sf,sff
 
C   Executable Code
C----------------------
c solve the elastic scattering kinematics
      kout(4)=kin/(1+(kin/me)*(1-cos(theta)))
      kout(3)=kout(4)*cos(theta)
      kout(2)=0
      kout(1)=kout(4)*sin(theta)
      Pfi(1)=-kout(1)
      Pfi(2)=-kout(2)
      Pfi(3)=-kout(3)+kin
      Pfi(5)=sqrt(Pfi(1)**2+Pfi(2)**2+Pfi(3)**2)
      Pfi(4)=sqrt(me**2+Pfi(5)**2)
 
c          N O T E S   O N   D E N S I T Y   M A T R I C E S
c-------------------------------------------------------------------------
c This program calculates the scattering amplitude with the y axis normal
c to the scattering plane.  If one is looking at the inclusive photon
c spectrum then he wants to average over all final electron directions.
c This average is a NOOP for the unpolarized or circularly polarized cases,
c but as soon as a second plane is introduced to the problem, either by
c plane polarization of the incoming or outgoing photon, then an averaging
c procedure becomes necessary.  The general form for a photon density
c matrix is:
c                    rho = 1/2 * (1  +  P . sigma)
c where
c 1 is the unit 2x2 matrix, P is a 3-vector representing the photon
c polarization and sigma are the Pauli matrices.  The values for P are
c       circular +/- :   P = +/- ( 0, 0, 1 )
c       plane at phi :   P = ( cos(2phi),  sin(2phi),  0 )
c All observables for plane polarization cancel out to the unpolarized
c case unless both initial and final photons are plane-polarized.  In
c that case the result for parallel (+) and perpendicular (-) planes is
c
c            W(+/-) = S(1/2,1) +/- S(x,x) +/- S(y,y)
c
c This means that the plane polarization A = (W(+)-W(-))/(W(+)+W(-))
c is given by
c                 A = ( S(x,x) + S(y,y) ) / ( 2 S(1/2,1) )
c
c I have used S(1/2,1) to represent the unpolarized cross section.

c initialize the density matrices (if not already loaded)                 !ZZZ
      trace=elDMi(1,1)+elDMi(2,2)
      If (abs(trace-1).gt.1e-4) then
        elDMi(1,1)=0.5
        elDMi(1,2)=0.0
        elDMi(2,1)=0.0
        elDMi(2,2)=0.5
      endIf
      trace=phDMi(1,1)+phDMi(2,2)
      If (abs(trace-1).gt.1e-4) then
        phDMi(1,1)=0.5
        phDMi(1,2)=0.0
        phDMi(2,1)=0.0
        phDMi(2,2)=0.5
      endIf
      trace=elDMf(1,1)+elDMf(2,2)
      If (abs(trace-1).gt.1e-4) then                                      !ZZZ
        elDMf(1,1)=1.0
        elDMf(1,2)=0.0
        elDMf(2,1)=0.0
        elDMf(2,2)=1.0
      endIf
      trace=phDMf(1,1)+phDMf(2,2)
      If (abs(trace-1).gt.1e-4) then
        phDMf(1,1)=1.0
        phDMf(1,2)=0.0
        phDMf(2,1)=0.0
        phDMf(2,2)=1.0
      endIf
 
c calculate final-state parameters:
c   eta = rapidity of final-state electron in lab
c   ome = angle from +z axis to Pfi (within the (z,-x) half-plane)
c   epsin(3,i) = incoming polarization vector for helicities + (i=1),- (i=2).
c   epsout(3,i) = outgoing polarization vector for helicities + (1), - (2).
      eta=Pfi(5)/Pfi(4)
      eta=log(eta+sqrt(eta**2+1))  ! poor man's arcsinh function
      ome=pi-theta
      seta=sinh(eta/2)
      ceta=cosh(eta/2)
      some=sin(ome/2)
      come=cos(ome/2)
      epsin(1,1)=1/sqrt(2.)
      epsin(2,1)=i_/sqrt(2.)
      epsin(3,1)=0
      epsin(4,1)=0
      epsin(1,2)=1/sqrt(2.)
      epsin(2,2)=-i_/sqrt(2.)
      epsin(3,2)=0
      epsin(4,2)=0
      epsout(1,1)=conjg(cos(theta)*epsin(1,1)+sin(theta)*epsin(3,1))
      epsout(2,1)=conjg(epsin(2,1))
      epsout(3,1)=conjg(cos(theta)*epsin(3,1)-sin(theta)*epsin(1,1))
      epsout(4,1)=0
      epsout(1,2)=conjg(cos(theta)*epsin(1,2)+sin(theta)*epsin(3,2))
      epsout(2,2)=conjg(epsin(2,2))
      epsout(3,2)=conjg(cos(theta)*epsin(3,2)-sin(theta)*epsin(1,2))
      epsout(4,2)=0
 
c calculate the W matrices:
c   W(mu,nu,si,sf) = exp(-i_*ome*Jy/2)
c                  * exp(-i_*eta*(sin(ome)*Kx/2-cos(ome)*Kz/2)
c                  *((epsout!(sf) epsin!(si) kin!)/kin
c                   +(epsin!(si) epsout!(sf) kout!)/kout))
c   where a! means a-slash (dotted with gamma matrices).
      Do si=1,2
        Do sf=1,2
          Do mu=1,4
            Do nu=1,4
              U(mu,nu)=gamma(mu,nu,4)-gamma(mu,nu,3)
              V(mu,nu)=(kout(4)*gamma(mu,nu,4)
     &                 -kout(1)*gamma(mu,nu,1)
     &                 -kout(2)*gamma(mu,nu,2)
     &                 -kout(3)*gamma(mu,nu,3))/kout(4)
            endDo
          endDo
          Call slasher(epsin(1,si),U)
          Call slasher(epsout(1,sf),U)
          Call slasher(epsout(1,sf),V)
          Call slasher(epsin(1,si),V)
          Do mu=1,4
            Do nu=1,4
              W(mu,nu,si,sf)=U(mu,nu)+V(mu,nu)
            endDo
          endDo
          Do mu=1,4
            Do nu=1,4
              U(mu,nu)=-i_*seta*(sin(ome)*K(mu,nu,1)
     &                          -cos(ome)*K(mu,nu,3))
              V(mu,nu)=-i_*some*J(mu,nu,2)
            endDo
            U(mu,mu)=ceta+U(mu,mu)
            V(mu,mu)=come+V(mu,mu)
          endDo
          call times(U,W(1,1,si,sf))
          call times(V,W(1,1,si,sf))
        endDo
      endDo
 
c now do the traces (upper components only)
      sum=0
      Do si=1,2
        Do sii=1,2
          Do sf=1,2
            Do sff=1,2
              Do mu=1,2
                Do nu=1,2
                  U(mu,nu)=elDMi(1,mu)*conjg(W(1,nu,sii,sff))
     &                    +elDMi(2,mu)*conjg(W(2,nu,sii,sff))
                  V(mu,nu)=elDMf(1,mu)*W(nu,1,si,sf)
     &                    +elDMf(2,mu)*W(nu,2,si,sf)
                endDo
              endDo
              trace=V(1,1)*U(1,1)+V(1,2)*U(2,1)
     &             +V(2,1)*U(1,2)+V(2,2)*U(2,2)
              sum=sum+phDMf(sf,sff)*phDMi(si,sii)*trace
            endDo
          endDo
        endDo
      endDo
 
c return with the differential cross section (in ub)
      Compton=sum*(hbarc*alpha*kout(4)/(2*me*kin))**2
      Compton=Compton*1e+4
      End
*DECK,slasher.
*CMZ :          09/11/97  22.57.28  by  Rick Jones
*-- Author :  Richard Jones
      Subroutine slasher(vector,matrix)
      Complex vector(4),matrix(4,4)
      Common /Cdirac/ gamma(4,4,4),J(4,4,3),K(4,4,3)
      Complex gamma,J,K
      Complex work(4,4)
      Integer mu,nu
      Do mu=1,4
        Do nu=1,4
          work(mu,nu)=vector(4)*gamma(mu,nu,4)
     &               -vector(1)*gamma(mu,nu,1)
     &               -vector(2)*gamma(mu,nu,2)
     &               -vector(3)*gamma(mu,nu,3)
        endDo
      endDo
      Call times(work,matrix)
      End
*DECK,times.
*CMZ :          09/11/97  22.57.28  by  Rick Jones
*-- Author :  Richard Jones
      Subroutine times(factor,product)
      Complex factor(4,4),product(4,4)
      Complex work(4,4)
      Integer mu,nu
      Do mu=1,4
        Do nu=1,4
          work(mu,nu)=factor(1,nu)*product(mu,1)
     &               +factor(2,nu)*product(mu,2)
     &               +factor(3,nu)*product(mu,3)
     &               +factor(4,nu)*product(mu,4)
        endDo
      endDo
      Do mu=1,4
        Do nu=1,4
          product(mu,nu)=work(mu,nu)
        endDo
      endDo
      End