c dsigL_pole.f
      
c calculation of pion pole contribution to dsigL/dt
c see eqns 3,4 of Huber 2008 paper      

c gh 16.06.12

      implicit none

      real*8 mpig,mpg,mng,pi,hbarc,eesq
      real*8 q2g,wg,nug,pgamg,betacm,gammacm,pcmg,epicmg
      real*8 pgamcmg,nuparcmg,tming
      real*8 r2_mono,Fpi_mono,r2_dip,Fpi_dip,pmono,Fpi_fit,fpisq
      real*8 lN,gpole,tt,gpinn,nn,dl_pole
      real*8 polefactor,fpi,q2fpi2,wfactor,sigl
      
      integer itt

      mpig=0.13957
      mpg=0.93827
      mng=0.93957
      
      pi=3.14159
      hbarc=0.197  ! (GeV-fm)
      eesq=4.*pi*hbarc/137.

      write(6,100)
 100  format(' Enter Q^2, W in GeV')
      read(5,*)q2g,wg

      nug = (wg**2 + q2g - mpg**2)/(2.*mpg)
      pgamg = sqrt( nug**2 + q2g)

c find speed of virtual photon+proton c.m. frame
      betacm  = pgamg/(nug+mpg)
      gammacm = (nug+mpg)/wg

      pcmg  = sqrt( (wg**2 + mng**2 - mpig**2)**2 - 4.*(wg*mng)**2 )
     1        /(2.*wg)
      epicmg= sqrt(pcmg**2 + mpig**2)

! calculate t_min
      pgamcmg = (pgamg-betacm*nug)*gammacm
      nuparcmg= nug/gammacm-betacm*pgamcmg                     !p.59 of notes
      tming   = q2g-mpig**2+2.*(nuparcmg*epicmg-pgamcmg*pcmg)  !p.105

c fit to experimental Fpi
c monoopole+dipole fit eqns 8-10 of Huber 2008 paper

c Amendolia fits, units fm^2
      r2_mono=0.431
      Fpi_mono=1./(1.+(r2_mono*q2g)/(6*hbarc**2))

      r2_dip=0.411
      Fpi_dip=1./(1.+(r2_dip*q2g)/(12*hbarc**2))**2
      
      pmono=0.85  ! best fit is 85% monopole + 15% dipole
      Fpi_fit=pmono*Fpi_mono+(1-pmono)*Fpi_dip
      fpisq=Fpi_fit**2

c g_piNN (monopole) form factor from T. Meissner, PRC 52(95)3386.
c      lN=0.800
c      gpole=13.4
c g_piNN (monopole) form factor used in the Bonn potential
c      lN=1.30
c g_piNN form factor (central values) in Favart et al, arXiv:1511.04535
c      lN=0.50
c      gpole=13.3
c g_piNN form factor values from P. Kroll, arXiv:1602.03803
      lN=0.44
      gpole=13.1
      
      write(6,10)
 10   format(/,8x'Q2',11x,'W',10x'-t      dsigL_pole',5x,
     1         'dsigL  %dsigL_pole')

      do itt=1,30

         tt=-tming-0.02*float(itt-1)
         gpinn=gpole*(lN**2-mpig**2)/(lN**2-tt)
      
c NN (flux factor according to Vanderhaeghen)
         nn=32.*pi*(wg**2-mpg**2)*
     1    sqrt( (wg**2-mpg**2)**2+q2g**2+2.*q2g*(wg**2+mpg**2) )

c BTM calculation of Actor, Korner, Bender, Il.Nuo.Cim 24A(74),369
         dl_pole=(4.*hbarc*eesq*(gpinn**2)*q2g*fpisq*(-tt))/
     1        (nn*(tt-mpig**2)**2)

c convert from (fm/GeV)^2to mub/(GeV^2)
         dl_pole=dl_pole*(10**4)

c-----------------------------------------------------------------------
c Tanja's sig2006 global fit, including higher Q2 
         
c pion pole term
         polefactor=abs(tt)/((abs(tt)+mpig**2)**2)
         
         fpi=1./(1.+1.77*q2g+0.12*q2g**2)
         q2fpi2=q2g*fpi**2
         sigl =19.8*polefactor*q2fpi2*exp((-3.66*(abs(tt))))

c W rescaling according to 1/(w**2-mp**2)**2
         wfactor=8.539/(wg**2-mpg**2)**2
         
         sigl=sigl*wfactor

         write(6,110)q2g,wg,tt,dl_pole,sigl,dl_pole/sigl*100
 110     format(6f12.4)
      
      enddo
      end