C Copy of Tito's file to compute coeffs for
c physics asymmetry, but this version for hydrogen
c (so much simpler). 
cc Jul04: EJB, set up for G0 kinematics rather than for SAMPLE

c It computes the coefficients a0-a3 with the asymmetry
c written as
c        A = a0 + a1*GEs + a2*GMs + a3*GAe
c I checked that the coefficients agree reasonably well with
c the one Raphael produced a few years ago.
c In this formulation, GAe(=GAe(T=1)) contains no s-quarks, and delta_s is
c buried in a0. So to get the "no-s-quark" asymmetry one has
c to set GAs=0 in the code, and then one has to put in an appropriate
c value of GAe. I left commented out how I think that GAe should
c be computed.
c The nucleon electromagnetic ffs are dipole with Galster for GEn.


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c     Input parameters: Beam energy (GeV)
c                       GAe 
c                       GMs
c
c      RT1 = -0.227              ! Zhu et al (MS-bar)
c
ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
      program aphys_lh2_g0

      Implicit NONE
      real energy, GAe, GMs, GEs, Aphystot, q2
      integer i

      energy = 3.03

      GMs = 0.
      GEs = 0.
      GAe = 0.


c gratuitous nested subroutine leftover from old code.
        call get_aphystot(energy, GEs, GMs, GAe, Aphystot)

c      write(*,*)q2, Aphystot


      End

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

      subroutine get_aphystot(energy, GEs, GMs, GAe, Aphystot)
      Implicit NONE
      real energy, GAe, GMs, GEs
      real d1, d2, d3
      real q2
      real Ah0, Ah1, Ah2, Ah3
      real Aphystot
      real c1, c2, c3

      integer i

      do q2=0.1,1.0,0.05

c 1=coeff in front of GEs
c 2  GMs
c 3  GAe

      d1=1.
      d2=0.
      d3=0.
      call get_Ap(energy, q2, d1, d2, d3, Ah1)
      d1=0.
      d2=1.
      d3=0.
      call get_Ap(energy, q2, d1, d2, d3, Ah2)
      d1=0.
      d2=0.
      d3=1.
      call get_Ap(energy, q2, d1, d2, d3, Ah3)
      d1=0.
      d2=0.
      d3=0.
      call get_Ap(energy, q2, d1, d2, d3, Ah0)

      
c      write(*,110)q2, ah0, ah1, ah2, ah3
      c1 = ah1 - ah0
      c2 = ah2 - ah0
      c3 = ah3 - ah0
      write(*,110) q2, ah0, c1, c2, c3


      enddo

 110  format(5(1x,f8.3))
 1212 format(1x,5(1x,0pf9.4))

      return
      End

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c
c     elastic scattering
c
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc


      subroutine get_Ap(energy,Q2, GEs, GMs, GAe, Ap)
      Implicit NONE

      real energy, GAe, GMs, GEs, Ap
      real Q2, mp, ee, ep, th
      real sth2, thetae, sth
      real aa, ae, am, aden
      real epsilon, epsilonp, tau


      real GEp, GMp, GEn, GMn,  GEZp, GMZp, GEZn, GMZn
      real GAz, GAs, GAZp, GAZn
      real GEp0, GMp0, GEn0, GMn0
      real GA0
      real RVp, RVn, RVs, RA0, GAs0, RA1
      real MV2, MA2

      real GF, alpha
      real s2w, hbarc2
      parameter (mp=938.272E-3)
      parameter (GF=1.166E-5)
      parameter (alpha=7.29735E-3)
      parameter (GEp0=1.0E0, GMp0=2.79E0, GEn0=0.0E0, GMn0=-1.91E0)
      parameter (GA0=1.2670E0)
      parameter (GAs = -0.1)
c      parameter (GAs = 0.)

      parameter (RVp=-0.054e0, RVn=-0.0143E0, RVs =0.0E0) 
      parameter (RA0=0.03E0)
      parameter (RA1=-0.227)
      parameter (MV2=0.711e0, MA2 = 1.0E0)
      parameter (s2w=0.23113E0, hbarc2 = .389379E6)
      real pi
      parameter (pi=3.14159265358979E0)
      integer i

c Energy in GeV

         ee = energy

         ep = ee - q2/2/mp
         sth2=q2/4/ee/ep
         sth=sqrt(sth2)
         th = 2*asin(sth)

         thetae=th*180/pi

         tau = Q2/4./mp**2.
         epsilon = 1./(1.+2.*(1.+tau)*tan(th/2.)**2.)
         epsilonp = sqrt((1.-epsilon**2.)*tau*(1.+tau))

c dipole form factors and (old) Galster for GEn
c
         GEp = GEp0/(1.+Q2/MV2)**2.
         GMp = GMp0*GEp
         GEn = -GMn0*GEp*tau/(1.+5.6*tau) !Galster parameterization
         GMn = GMn0*GEp

         GEZp = (1.-4.*s2w)*(1.+RVp)*GEp-(1.+RVn)*GEn-(1.+RVs)*GEs
         GMZp = (1.-4.*s2w)*(1.+RVp)*GMp-(1.+RVn)*GMn-(1.+RVs)*GMs
         GEZn = (1.-4.*s2w)*(1.+RVp)*GEn-(1.+RVn)*GEp-(1.+RVs)*GEs
         GMZn = (1.-4.*s2w)*(1.+RVp)*GMn-(1.+RVn)*GMp-(1.+RVs)*GMs 

         GAz = GA0/(1.+Q2/MA2)**2.
c         GAe = -GAz*(1+RA1) 

         GAZp = GAe + RA0 + GAs

         ae = epsilon*GEp*GEZp
         am = tau*GMp*GMZp
         aa = -(1.-4.*s2w)*epsilonp*GMp*GAZp
         aden = epsilon*GEp**2.+tau*GMp**2.

         Ap = -GF*Q2/4./pi/alpha/sqrt(2.)*((ae + am + aa)/aden)*1.e6


      return
      end