c dsigL_pole.f
      
c calculation of K+ pole contribution to dsigL/dt
c see eqns 3,4 of Huber 2008 paper and references listed below   

c gh 16.07.07

      implicit none

      real*8 mkg,mpg,mlg,msg,pi,hbarc,eesq
      real*8 q2g,wg,nug,pgamg,betacm,gammacm,pcmLg,ekcmLg,pcmSg,ekcmSg
      real*8 pgamcmg,nuparcmg,tminLg,tminSg,Fk,fksq
      real*8 lN,gpoleL,gpoleS,tt,gkLn,gkSn,nn,dl_poleL,dl_poleS
      real*8 fact_q,fact_w,fact_t,sigl_L,sigl_S
      
      integer itt

      mkg=0.49368
      mpg=0.93827
      mLg=1.11568
      mSg=1.19264
      
      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

      pcmLg  = sqrt( (wg**2 + mLg**2 - mkg**2)**2 - 4.*(wg*mLg)**2 )
     1        /(2.*wg)
      ekcmLg= sqrt(pcmLg**2 + mkg**2)

      pcmSg  = sqrt( (wg**2 + mSg**2 - mkg**2)**2 - 4.*(wg*mSg)**2 )
     1        /(2.*wg)
      ekcmSg= sqrt(pcmSg**2 + mkg**2)

! calculate t_min
      pgamcmg = (pgamg-betacm*nug)*gammacm
      nuparcmg= nug/gammacm-betacm*pgamcmg                     !p.59 of notes

      tminLg   = q2g-mkg**2+2.*(nuparcmg*ekcmLg-pgamcmg*pcmLg) !p.105
      tminSg   = q2g-mkg**2+2.*(nuparcmg*ekcmSg-pgamcmg*pcmSg)

c fit to experimental Fk
c use monopole formula from Goloskokov & Kroll, EPJA 47, 112 (2011)
c essentially Fpi multiplied by 0.9 to account for flavor-symmetry breaking
      Fk=0.9/(1.+q2g/0.462)
      fksq=Fk**2

c g_kpY form factor values from Goloskokov & Kroll, EPJA 47, 112 (2011)
      lN=0.44
      gpoleL=13.3
      gpoleS=-3.5
      
      write(6,10)
 10   format(/,4x'Q2',7x,'W',6x'-t       dsigL_pole',2x,
     1         'dsigL  %dsigL_pole   dsigS_pole',2x,
     2         'dsigS  %dsigS_pole')

      do itt=1,30

         tt=-tminLg-0.02*float(itt-1)
         gkLn=gpoleL*(lN**2-mkg**2)/(lN**2-tt)
         gkSn=gpoleS*(lN**2-mkg**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
c according to Favart et al, aXiv:1511.04535 the only changes from pion pole 
c should be the meson mass, gkYn and lN
         dl_poleL=(4.*hbarc*eesq*(gkLn**2)*q2g*fksq*(-tt))/
     1        (nn*(tt-mkg**2)**2)
         dl_poleS=(4.*hbarc*eesq*(gkSn**2)*q2g*fksq*(-tt))/
     1        (nn*(tt-mkg**2)**2)

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

c-----------------------------------------------------------------------
c factorized p(e,e'K+)Lambda model described in Tanja Horn Dec 2007 report 
c model is based on Q^2, W dependence of Bebek and Brauel data, using 
c updated Koltenuk fits to improve epsilon, etc. dependence

c Lambda
         fact_q = 4.50/(q2g+2.67)**2
         fact_w = 4.1959*pcmLg/(wg*(wg**2-mpg**2))
         fact_t = 1.10*exp(-2.1*(tt-tminLg))
         sigl_L = fact_q*fact_t*fact_w*(1./3.)

c Sigma
         fact_q = 0.1587/(q2g+0.785)**2
         fact_w = 1.3176*pcmSg/(wg*(wg**2-mpg**2))
         fact_t = 1.0*exp(-1.0*(tt-tminSg))
         sigl_S = fact_q*fact_t*fact_w*(1./3.)

         write(6,110)q2g,wg,tt,dl_poleL,sigl_L,dl_poleL/sigl_L*100,
     1     dl_poleS,sigl_S,dl_poleS/sigl_S*100
 110     format(3f8.3,3x,3f10.4,3x,3f10.4)
      
      enddo
      end