program call_pb
c
      implicit none
c
      real*8 e0,theta,ep ! units of GeV,radians,GeV
      real*8 xsecqe,xsec_inel,sigmot,sigma
      real*8 w,q2,thr,enu,f2,R,f1,wsq
      real*8 asym,xn_aren
      real*8 xbj,eps,p1,p2,p3,p4,p5,mp,mp2
      integer i
c 
      mp=.93827
      mp2=mp*mp
       open(2,file='deut_asym_kin.dat')
      do i=1,100
         read(2,*,end=999) q2,xbj,eps,p1,p2,p3,p4,p5,w,theta  ! units GeV,GeV2,deg
         thr=theta/180.*3.14159
         e0=5.7525
         enu=(w*w+q2-mp2)/2./mp
         ep=e0-enu
         call QUASIY8(1.d0,2.0d0,E0,EP,theta,xsecqe)
         xsecqe=xsecqe*1.0d-03 ! convert to nb/sr/MeV
         wsq=w*w
         call f1f2in06(1.d0,2.d0,Q2,wsq,f1,F2,R)
	 xsec_inel=sigmot(E0*1000.,thr)*f2*(1.d0/enu+
     &     2.d0*enu*(1.d0+q2/enu**2)*(tan(thr/2.d0))**2/(q2*(1d0+r))) ! 6/03
         xsec_inel=xsec_inel*1.0d+4  ! nb/sr/MeV
         write(*,'(5(f7.4,1x),6e15.5)') e0,ep,theta,q2,w,xsecqe,xsec_inel,xsecqe+xsec_inel,f1,f2,r
      enddo
c
 999  end
      
*sigmot
      real*8 function sigmot(e,thr)
      implicit none
      real*8 e, thr, alph, hbarc
      alph=1.d0/137.03604d0
      hbarc=197.3286d0
      sigmot=(alph*hbarc*cos(thr/2.d0)/2.d0/e/sin(thr/2.d0)**2)**2
c  fm**2/sr
      return
      end