      program lomon_fit
c
c E. L. Lomon in PRC66 045501 (2002)
c
c  Models GKex(01) , GKex(01-) , GKex(02L) , GKex(02S)
c
c add GKex(05) from 
c
c
      implicit none

      integer*4 i,im
      real*8 grhorat(5),kapparho(5),gomegarat(5),kappaomega(5)
      real*8 kappaomegap(5),gomegaprat(5)
      real*8 gphirat(5),kappaphi(5),muphi(5)
      real*8 lambda1(5),lambdad(5),lambda2(5)
      real*8 lambdaqcd(5),n(5)
      real*8 kappanu,kappas,mrho,momega,mphi,momegap,mrhop
      real*8 q2,qt2
      real*8 f1rho,f2rho,f1d,f2d,f1iv,f2iv
      real*8 f1omega,f2omega,f1phi,f2phi
      real*8 f1is,f2is
      real*8 t1,t2,t3
      real*8 f1p,f2p,gep,gmp,tau,gd,mp,mu,ratio,gen,gmn,f1n,f2n,mun
c
      data grhorat  / 0.0636, 0.0598, 0.0608, 0.0401,0.007208 /
      data kapparho /-0.4175,-15.9227,5.3038, 6.8190,12.0 /
      data gomegarat/ 0.7918, 0.6981, 0.6896, 0.6739,0.7021 /
      data kappaomega / 5.1109,1.9333,-2.8585,0.8762,0.4027 /
      data gomegaprat / 0.0,   0.0,  0.2346, 0.2552,0.164  /
      data kappaomegap / 0.0,  0.0, 18.2284, 1.4916,-2.973/ 
      data gphirat    /-0.3011,-0.5270,-0.1852,-0.1676,-.1711/
      data kappaphi   /13.4385,2.3241,13.0037,7.0172,0.01 /
      data muphi      / 1.1915,1.5113,0.6848, 0.8544,0.2/
      data lambda1    / 0.9660,1.1276,0.9441,0.9407,.93088 /
      data lambdad    / 1.3406,1.8598,1.2350,1.211,1.181 /
      data lambda2    / 2.1382,1.2255,2.8268,2.7891,2.6115 /
      data lambdaqcd  / 0.1163, 0.1315, 0.150, 0.150,0.150 /
      data n          / 1.0,    0.8709, 1.0,  1.0 ,1.0   /
c parameter fixed for all fits
      kappanu=3.706
      kappas=-0.12
      momega=0.784
      mphi=1.019
      mrho=0.776
      mrhop=1.45
      momegap=1.419
      mp=0.9382796
      mu=2.79
      mun=-1.91
c
      write(*,*) ' Which model (1-5)?'
      write(*,*) '1=GKex(01)  2=GKex(01-)  3=GKex(02L)  
     1  4=GKex(02S)  4=GKex(05) '
      read(*,*) im

      do i = 1,3000
         q2=(i)*0.01
         tau=q2/(4.0*mp**2)
         qt2=q2*log((lambdad(im)**2+q2)/lambdaqcd(im)**2)/
     &        log(lambdad(im)**2/lambdaqcd(im)**2)
         f1rho=(lambda1(im)**2/(lambda1(im)**2+qt2))*
     &        (lambda2(im)**2/(lambda2(im)**2+qt2))
         f1omega=f1rho
         f1phi=f1rho*(q2/(lambda1(im)**2+q2))**1.5
         f1d=(lambdad(im)**2/(lambdad(im)**2+qt2))*
     &        (lambda2(im)**2/(lambda2(im)**2+qt2))
         f2rho=(lambda1(im)**2/(lambda1(im)**2+qt2))*
     &        (lambda2(im)**2/(lambda2(im)**2+qt2))**2
         f2omega=f2rho
         f2phi=f2rho*((lambda1(im)**2*(q2+muphi(im)**2))/
     &        (muphi(im)**2*(lambda1(im)**2+q2)))**1.5
         f2d=(lambdad(im)**2/(lambdad(im)**2+qt2))*
     &        (lambda2(im)**2/(lambda2(im)**2+qt2))**2
        t1=n(im)/2.*((1.0317+0.0875/((1+q2/0.3176)**2))/(1+q2/0.5496))
         t2=grhorat(im)*mrhop**2/(mrhop**2+q2)*f1rho
         t3=(1-1.1192*n(im)/2.-grhorat(im))*f1d
c         if(abs(q2-5.6).le.0.05)write(*,*)q2,t1,t2,t3,t1+t2+t3
         f1iv=t1*f1rho+t2+t3         
         t1=n(im)/2.*((5.7824+0.3907/((1+q2/0.1422)))/(1+q2/0.5362))
         t2=kapparho(im)*grhorat(im)*mrhop**2/(mrhop**2+q2)*f2rho
c         if(q2.ge.9.9)write(*,*)kappanu,6.1731*n,f2d
         t3=(kappanu-6.1731*n(im)/2.-kapparho(im)*grhorat(im))*f2d
c         if(abs(q2-5.6).le.0.05)write(*,*)q2,t1,t2,t3,t1+t2+t3
         f2iv=t1*f2rho+t2+t3
         t1=(gomegarat(im)*(momega**2/(momega**2+q2))+
     &        gomegaprat(im)*(momegap**2/(momegap**2+q2)))*f1omega
         t2=gphirat(im)*(mphi**2/(mphi**2+q2))*f1phi
         t3=(1-gomegarat(im)-gomegaprat(im))*f1d
c         if(abs(q2-5.6).le.0.05)write(*,*)q2,t1,t2,t3,t1+t2+t3
         f1is=t1+t2+t3
         t1=(kappaomega(im)*gomegarat(im)*(momega**2/(momega**2+q2))+
     &        kappaomegap(im)*gomegaprat(im)*(momegap**2/(momegap**2+
     &        q2)))*f2omega
         t2=kappaphi(im)*gphirat(im)*(mphi**2/(mphi**2+q2))*f2phi
         t3=(kappas-kappaomega(im)*gomegarat(im)-
     &     kappaomegap(im)*gomegaprat(im)-kappaphi(im)*gphirat(im))*f2d
c         if(abs(q2-5.6).le.0.05)write(*,*)q2,t1,t2,t3,t1+t2+t3
         f2is=t1+t2+t3
         f1p=(f1is+f1iv)/2.0
         f2p=(f2is+f2iv)/2.0
         f1n=(f1is-f1iv)/2.0
         f2n=(f2is-f2iv)/2.0
         gep=f1p-tau*f2p
         gmp=f1p+f2p
         gen=f1n-tau*f2n
         gmn=f1n+f2n
         gd=1.0/(1.0+q2/0.71)**2
         ratio=mu*gep/gmp
         write(41,'(6f15.5)')
     >   q2,gmp/(mu*gd),ratio,gmn/(mun*gd),gen,gmp/gmn
      enddo


      stop
      end

