c pole_ratio.f
      
c ratio of K+/pi+ dsigL/dt assuming pole dominance
c see eqns 3,4 of Huber 2008 paper and references listed below   

c gh 21.05.31

      implicit none

      real*8 mkg,mpg,mLg,mSg,mpig,mng,pi,hbarc
      real*8 q2g,wg,nug,pgamg,betacm,gammacm,pcmLg,ekcmLg,pcmSg,ekcmSg
      real*8 pgamcmg,nuparcmg,tminLg,tminSg,Fk,fksq
      real*8 r2_mono,Fpi_mono,r2_dip,Fpi_dip,pmono,Fpi_fit,fpisq
      real*8 lNpi,lNk,gpoleKL,gpoleKS,tt,gkLn,gkSn,gpolepi,gpinn
      real*8 dl_polepi,dl_poleKL,dl_poleKS
      real*8 ratio_KLpi,ratio_KSpi
      
      integer itt

      mkg=0.49368
      mpg=0.93827
      mLg=1.11568
      mSg=1.19264
      mpig=0.13957
      mng=0.93957
      
      pi=3.14159
      hbarc=0.197  ! (GeV-fm)

      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)

c      pcmSg  = sqrt( (wg**2 + mSg**2 - mkg**2)**2 - 4.*(wg*mSg)**2 )
c     1        /(2.*wg)
c      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
c      tminSg   = q2g-mkg**2+2.*(nuparcmg*ekcmSg-pgamcmg*pcmSg)

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 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_piNN form factor values from P. Kroll, arXiv:1602.03803
c   coupling const=13.1+/-0.2 and monopole cutoff=0.44+/-0.07    
c note the lambda_N cutoff value determined by T. Meissner, PRC 52(1995)3386,
c   is much bigger, 0.80 GeV, for similar 13.4 coupling constant.
      lNpi=0.44
      gpolepi=13.1
      
c g_kpY form factor values from Goloskokov & Kroll, EPJA 47, 112 (2011)
c they give no guidance on the cutoff
      lNk=(0.44+0.80)/2
      gpoleKL=-13.3
      gpoleKS=-3.5
      
      write(6,10)
 10   format(/,4x'Q2',7x,'W',6x'-t        K+L/pi+n',2x,'K+S/pi+n')

      do itt=1,30
         tt=-tminLg-0.05*float(itt-1)

         gpinn=gpolepi*(lNpi**2-mpig**2)/ (lNpi**2-tt)
         gkLn= gpoleKL*(lNk**2-mkg**2)  / (lNk**2-tt)
         gkSn= gpoleKS*(lNk**2-mkg**2)  / (lNk**2-tt)
      
c simplified sigL with flux factor and other common factors removed
c
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_polepi= ((gpinn**2)*fpisq)/ ((tt-mpig**2)**2)
         dl_poleKL= ((gkLn**2)*fksq)  / ((tt-mkg**2)**2)
         dl_poleKS= ((gkSn**2)*fksq)  / ((tt-mkg**2)**2)

         ratio_KLpi= dl_poleKL / dl_polepi
         ratio_KSpi= dl_poleKS / dl_polepi

         write(6,110)q2g,wg,tt,ratio_KLpi,ratio_KSpi
 110     format(3f8.3,3x,2f10.4)
      
      enddo
      end