program hydro
c this program calculates the production of antihydrogen in 
c collisions of antinuclei with heavy targets (can also be used
c calculate K-shell capture in HI collisions).
c UNITS: hbar = c = m_e = 1
c   electron  and positron energies and momenta in units of m_e
c   e^2=1/137
c   cross sections in units of picobarns
c C.A. Bertulani version  1/October/1997
c-------------------------------------------------------------
      implicit real*8(a-h,o-z)
      parameter(n=101,nk=101)
      dimension sig(n,nk),efree(n),theta(nk),vec2(n),vec1(n),
	&          phot(n,nk,n),vec3(nk),vec4(n),vec5(nk)
      real*8 kz,kz2,kt,kt2
      open(12,file='hydro.out',status='unknown')
c inputs from screen
      print*,' enter value of ZP, ZT and GAMMA'
      read(*,*)zp,zt,gamma
c constants    
      pi=acos(-1.)
      e2=1./137.
      beta=sqrt(gamma*gamma-1.)/gamma
      a0=1./zp/e2
      cross=1.49d15
      factor=1./pi*(zt/beta)**2*e2 
      kz=1./beta/gamma**2
      kz2=kz*kz
c integration limits
      emax=25.
      emin=1.000001
      tmax=pi
      tmin=0.
      qtmax=5.
      qtmin=0.1/gamma**2
c integration step size
      de=(emax-emin)/n
      dt=(tmax-tmin)/nk
      dqt=(qtmax-qtmin)/n
c electron polar angle and electron energy
      do i=1,nk
        theta(i)=(i-1)*dt+tmin
      enddo
      do i=1,n
        efree(i)=(i-1)*de+emin
      enddo
c ********** start calculation of (d Sigma/d E_e d Omega_e)
c ******************************
c loop over electron energy
      do i=1,n
        omega=efree(i)+1.
        omgv=omega/gamma/beta
        qz=omega/beta
        qz2=qz*qz 
        pfree=sqrt(efree(i)*efree(i)-1.)
        em1=efree(i)-1.
        p2=pfree*pfree
c loop over electron polar angle
        do k=1,nk
          pz=pfree*cos(theta(k))
          pt=pfree*sin(theta(k))
c loop over virtual photon momentum 
          do jq=1,n
            qt=(jq-1)*dqt+qtmin
            qt2=qt*qt
            q2=qt2+qz2
            q2g2=qt2+omgv*omgv
            q2p2=q2-p2
            kt=qt/omega
            kt2=kt*kt
            aux=32.*e2/a0**5*pfree*omega/qt2
	      qmpz=qz-pz
            v1=p2+q2-2*pz*qz
            v2=2.*pt*qt
            v3=pz*qz-p2
	      v4=v2
c ** integral over virtual photon polarization angle was
c    performed using MAPLE V  ** 
c transverse (virtual) photo production of anti-H
            phott=aux*cross*kt2*f1(v1,v2,v3,v4,q2p2,efree(i),qt,pt)
c longitudinal (virtual) photo production of anti-H
            photz=aux*cross*kz2*f2(v1,v2,v3,v4,q2p2,efree(i),pz,qmpz)
c interference 
            photzt=aux*cross*2*kz/omega*
     #      f3(v1,v2,q2p2,efree(i),pz,qmpz,pt,qt)
c here is the total (virtual) photo cross section
            sum=phott+photz+photzt
	      phot(i,k,jq)=sum
            vec2(jq)=2.*qt*qt2*sum/q2g2**2
          enddo          
c end loop over virtual photon momentum: calculate its integral
          call arsimp(nk,vec2,dqt,sum)
          sig(i,k)=sin(theta(k))*factor*sum/omega           
        enddo
c end loop over electron polar angle: sig(i,k) is (d Sigma/d E_e d Theta_e)
      enddo
c end loop over electron energy
c******************************
c ****** integration of (d Sigma/d E_e d Theta_e) over Theta_e 
      do i=1,n
        do k=1,nk
          vec3(k)=sig(i,k)
        enddo
        call arsimp(nk,vec3,dt,sum)
        vec4(i)=sum
      enddo
c ****** integration of (virtual) photo cross section over Theta_e 
      do i=1,n,10
	  write(12,*)'photo cross sections: (E_e)  x  (qt)  x  (sigma)'
        do iq=1,n 
          qt=(iq-1)*dqt+qtmin
          do k=1,nk
            vec5(k)=sin(theta(k))*phot(i,k,iq)
c	if(iq.eq.1.and.i.eq.1)write(12,*)cos(theta(k)),phot(i,k,iq)
          enddo
          call arsimp(nk,vec5,dt,sum)
c cross section by real photons
          pfree=sqrt(efree(i)*efree(i)-1.)
	    realph=cross*4.*pi*e2/a0**5*pfree/(efree(i)+1.)**4
     &            *(efree(i)**2+2./3.*efree(i)
     &            +4./3.-(efree(i)+2.)/pfree*dlog(efree(i)+pfree)) 	         
          write(12,10)efree(i),qt,sum,realph
10        format(1x,4(g10.3,1x))
	  enddo
      enddo
      write(12,*)' '
c output of (d Sigma/d E_e)
     	write(12,*)'heavy ion cross sections'
      print*,' gamma=', gamma
      print*,' here follows (d Sigma/d E_e)' 
      write(12,*)' gamma=', gamma
      write(12,*)' here follows (d Sigma/d E_e)' 
      do j=1,n
        print*,efree(j),vec4(j)
        write(12,*)efree(j),vec4(j)
      enddo 
c integration of d Sigma/d E_e  over E_e 
      call arsimp(n,vec4,de,sum)
      res=sum
c equivalent photon approximation
      t1=gamma*gamma
      t2 = t1-1
      do i=1,n
      t4 = (efree(i)+1)**2
      t9 = efree(i)**2
      t13 = dlog(4*t2/t4+1)-4*t2/(4*t1+2*efree(i)-3+t9)
      pfree=sqrt(efree(i)*efree(i)-1.)
	  realph=cross*4.*pi*e2/a0**5*pfree/(efree(i)+1.)**4
     &          *(efree(i)**2+2./3.*efree(i)
     &          +4./3.-(efree(i)+2.)/pfree*dlog(efree(i)+pfree)) 	         
	  vec1(i)=t13*factor/(efree(i)+1.)*realph
	enddo
      call arsimp(n,vec1,de,sum)
      epa=sum
      sbrod=zt**2*1.42*(dlog(gamma*gamma+1)-gamma*gamma/(gamma*gamma+1))
c output of total cross section + compare with Brodsky's formula
      print*,' Cross section in picobarns + EPA results:'
      print*,res,epa,sbrod 
      write(12,*)' Cross section in picobarns + EPA results:'
      write(12,*)res,epa,sbrod 
      end 
c
c=====================================================================
c     
      real*8 function f1(v1,v2,v3,v4,q2p2,e,qt,pt)
      implicit real*8(a-h,o-z)
      t2 = v1**2
      t3 = t2**2
      t4 = t2*v1
      t6 = v2**2
      t7 = t6*t2
      t8 = q2p2**2
      t9 = t2*t8
      t12 = t6**2
      t13 = t6*t8
      t16 = 1/t8
      t18 = v1-v2
      t19 = v1+v2
      t22 = 1/sqrt(t18*t19)
      t23 = t18**2
      t25 = t22/t23
      t26 = t19**2
      t27 = 1/t26
      t30 = qt**2
      t36 = pt**2
      t42 = 1/q2p2
      t52 = e*q2p2
      t56 = e**2
      t57 = t3*v1
      t66 = 1/t26/t19
      t68 = 1/t23/t18
      t69 = t66*t68
      t77 = v3*q2p2
      t83 = t6*v2
      t92 = e*v3
      t95 = q2p2*t6
      t105 = e*v4
      t106 = t4*v2
      t109 = v1*t83
      t118 = -2*v3*t4*t8-3*v3*v1*t13-2*t77*t3+t77*t7+t77*t12-4*v4*v2*t9-
     #v4*t83*t8-3*v4*t4*v2*q2p2+3*v4*v1*t83*q2p2+2*t92*q2p2*t3-t92*t95*t
     #2-t92*q2p2*t12+2*t92*t57-4*t92*t4*t6+2*t92*v1*t12+3*t105*t106*q2p2
     #-3*t105*t109*q2p2+2*t105*v2*t3-4*t105*t83*t2+2*t105*t12*v2
      t134 = qt*e
      t137 = pt*e
      f1 = (e-1)*(0.3141593E1*(2*t3+4*q2p2*t4-4*t7+2*t9-4*v1*t6*q2p2+2
     #*t12+t13)*t16*t25*t27/4-0.3141593E1*(2*t30*t2+t30*t6-6*pt*qt*v1*v2
     #+t36*t2+2*t36*t6)*t27*t25*t42)+(e+1)*0.3141593E1*(2*t4*t8+3*v1*t8*
     #t6-4*t52*t3+2*t52*t7+2*t52*t12+2*t56*t57-4*t56*t4*t6+2*t56*v1*t12)
     #*t69*t22*t16+0.3141593E1*t118*t66*t68*t22*t16+2*pt*0.3141593E1*(4*
     #qt*v2*q2p2*t2+qt*t83*q2p2-pt*t4*q2p2-4*pt*v1*t95-3*t134*t106+3*t13
     #4*t109+t137*t3+t137*t7-2*t137*t12)*t69*t22*t42
      return
      end		 
c
c=====================================================================
c     
      real*8 function f2(v1,v2,v3,v4,q2p2,e,pz,qmpz)
      implicit real*8(a-h,o-z)
      t2 = v1**2
      t3 = t2**2
      t4 = t2*v1
      t6 = v2**2
      t7 = t6*t2
      t8 = q2p2**2
      t9 = t2*t8
      t12 = t6**2
      t13 = t6*t8
      t16 = 1/t8
      t18 = v1-v2
      t19 = v1+v2
      t22 = 1/sqrt(t18*t19)
      t23 = t18**2
      t24 = 1/t23
      t25 = t22*t24
      t26 = t19**2
      t27 = 1/t26
      t30 = qmpz**2
      t31 = 2*t2+t6
      t46 = e*q2p2
      t50 = e**2
      t51 = t3*v1
      t60 = 1/t26/t19
      t62 = 1/t23/t18
      t71 = v3*q2p2
      t77 = t6*v2
      t86 = e*v3
      t99 = e*v4
      t112 = -2*v3*t4*t8-3*v3*v1*t13-2*t71*t3+t71*t7+t71*t12-4*v4*v2*t9-
     #v4*t77*t8-3*v4*t4*v2*q2p2+3*v4*v1*t77*q2p2+2*t86*q2p2*t3-t86*q2p2*
     #t6*t2-t86*q2p2*t12+2*t86*t51-4*t86*t4*t6+2*t86*v1*t12+3*t99*t4*v2*
     #q2p2-3*t99*v1*t77*q2p2+2*t99*v2*t3-4*t99*t77*t2+2*t99*t12*v2
      f2 = (e-1)*(0.3141593E1*(2*t3+4*q2p2*t4-4*t7+2*t9-4*v1*t6*q2p2+2
     #*t12+t13)*t16*t25*t27/4-t30*t31*0.3141593E1*t27*t24*t22/q2p2)+(e+1
     #)*0.3141593E1*(2*t4*t8+3*v1*t8*t6-4*t46*t3+2*t46*t7+2*t46*t12+2*t5
     #0*t51-4*t50*t4*t6+2*t50*v1*t12)*t60*t62*t22*t16+0.3141593E1*t112*t
     #60*t62*t22*t16+2*pz*qmpz*t31*0.3141593E1*t27*t25
      return
      end		 
c
c=====================================================================
c     
      real*8 function f3(v1,v2,q2p2,e,pz,qmpz,pt,qt)
      implicit real*8(a-h,o-z)
      t2 = v1**2
      t3 = v2**2
      t7 = v1+v2
      t8 = t7**2
      t11 = v1-v2
      t12 = t11**2
      t18 = 1/sqrt(t2-t3)
      t19 = 1/t8/t7/t12/t11*t18
      t21 = 1/q2p2
      t22 = e*t21
      t25 = 1/t8
      t26 = 1/t12
      t27 = t25*t26
      t31 = (4*t2+t3)*v2*0.3141593E1*t19-3*t22*v1*v2*t27*t18*0.3141593E1
      t34 = qt**2
      t40 = 2*t2+t3
      t45 = 1/sqrt(t11*t7)
      t54 = e-1
      t56 = t21*qmpz
      f3 = pt*qt*qmpz*t31+pz*t34*((2*t2+3*t3)*v1*0.3141593E1*t19-t22*t4
     #0*0.3141593E1*t25*t26*t45)-pz*qt*pt*t31-t54*t34*t56*t40*0.3141593E
     #1*t27*t45+t54*qt*pt*t56*t31
      return
      end		 
c
c=====================================================================
c Integration of a real function by Simpson's rule. N (must be odd) is
c the number of mesh points, A is a vector containing the integrand 
c values on the mesh and H is the step size.
C********
       subroutine arsimp(n,a,h,sum)
       implicit real*8 (a-h,o-z)
       dimension a(n)
       l=n
       s=a(1)-a(l)
       do i=2,l,2
         s=s+4.*a(i)+2.*a(i+1)
       enddo
       sum=h*s/3.
       return
       end