C *********************************************
C
C Code to calculate A_LU for Xe-131 and the proton 
C 
C We keep all harmonic for BH amplitude and 
C twist-2 harmonics for DVCS and Interference: 
C c0_DVCS, c0_I, c1_I and s1_I
C 
C Using the dual parameterization of GPDs, 
* V. Guzey and T. Teckentrup, Phys. Rev. D74, 054027 (2006) [hep-ph/0607099]
* V. Guzey, T. Teckentrup, e-Print: arXiv:0810.3899


C Note: used the non-factorized Regge-motivated t-dependence 
C
C V. Guzey, Theory Group, Jefferson Lab, vguzey@jlab.org
C Jan. 2009
C
C ************************************************      

      program ALU_on_Xe131            
      implicit double precision (a-h,o-z)

       dimension t_array(15) 

        data t_array /-0.005,-0.0075,-0.01,-0.02,-0.03,-0.04,-0.05,
     >-0.06,-0.07,-0.08,-0.09,-0.1,-0.2,-0.3,-0.4 /

        pi=acos(-1.0)

C External kinematics
C Positrons: a_sign=-1; electrons: a_sign=+1
      a_sign=-1.0
     
      E_lepton=27.5 
      write(6,65) 'Beam energy=',E_lepton
 65   format(A,2f6.3)   
      Qn=1.7
      write(6,65) 'Q**2=',Qn  
      xB=0.065
      write(6,65) 'x_B=',xB  
c      t=-0.1  
c      write(6,65) 't=',t  
        phi=pi/2.0
      write(6,65) 'phi=',phi  

C J_u and J_d
C There are 49 allowed combinations
C (J_u,J_d)=(0.6,1),(0.6,0.6),(0.6,0.2),(0.6, 0),(0.6,-0.2),(0.6,-0.6),(0.6.-1)
C           (0.4,1),(0.4,0.6),(0.4,0.2),(0.4, 0),(0.4,-0.2),(0.4,-0.6),(0.4.-1)
C           (0.2,1),(0.2,0.6),(0.2,0.2),(0.2, 0),(0.2,-0.2),(0.2,-0.6),(0.2.-1)
C             ..................................................................
C           (-0.6,1),(-0.6,0.6),(-0.6,0.2),(-0.6, 0),(-0.6,-0.2),(-0.6,-0.6),(-0.6.-1)
               
        aJu=0.0
        aJd=0.0

C        phi=pi-phi_Trento

        Nt=15
        do it=1,Nt
        t=t_array(it)

        call ASYMMETRY_ALU_Xe131
     >(a_sign,E_lepton,xB,Qn,t,phi,aJu,aJd,A_LU_Full,
     >A_LU_CohEnr,A_LU_Coh,A_LU_Incoh,A_LU_Proton)


        write(6,44)t,A_LU_Full/A_LU_Proton,A_LU_Coh/A_LU_Proton,
     >A_LU_CohEnr/A_LU_Proton,A_LU_Incoh/A_LU_Proton
 44     format(10f10.4)
        
        enddo

      stop
      end

C Main subroutine:
C a_sign = the sign in front of the interf. term (a_sign=1 electrons; a_sign=-1 positrons)
C E_lepton = lepton energy
C x = x_Bjorken
C Qn= Q**2 [GeV**2]
C t = 4-momentum transfer
C phi =angle between the lepton and production planes in the BMK notation
C aJu and aJd = total angular momentum carried by the u-quark and the d-quark

C Use the results of  
C A.V. Belitsky, D. Mueller, A. Kirchner, Nucl. Phys. B629, 323 (2002) [hep-ph/0112108]
        subroutine ASYMMETRY_ALU_Xe131
     >(a_sign,E_lepton,x,Qn,t,phi,aJu,aJd,A_LU_Full,
     >A_LU_CohEnr,A_LU_Coh,A_LU_Incoh,A_LU_Proton)
        implicit double precision (a-h,o-z)

        amn=0.938272
        pi=acos(-1.0)

C Kr-84
        A=131.0
        Z=54.0

        y=Qn/(2.0*x*amn*E_lepton)
        y_p=y

        if (y.gt.0.95) then
        write(6,*) 'y is larger 0.95'
        endif

        eps=4.0*amn**2*x**2/Qn
                
C Minimal t
        t_min=-Qn*(2.*(1.-x/A)*(1.-(1.+eps)**0.5)+eps)
     >/(4.*x/A*(1.-x/A)+eps)

        t_min_p=-Qn*(2.*(1.-x)*(1.-(1.+eps)**0.5)+eps)
     >/(4.*x*(1.-x)+eps)

        factor_K=(-t/Qn*(1.-x/A)*(1.-y-0.25*y**2*eps)
     >*(1.-t_min/t)*((1.+eps)**0.5+(4.*x/A*(1.-x/A)+eps)
     >/(4.*(1.-x/A))*(t-t_min)/Qn))**0.5

        factor_K_p=(-t/Qn*(1.-x)*(1.-y_p-0.25*y_p**2*eps)
     >*(1.-t_min_p/t)*((1.+eps)**0.5+(4.*x*(1.-x)+eps)
     >/(4.*(1.-x))*(t-t_min_p)/Qn))**0.5


        factor_J=(1.-y-0.5*y*eps)*(1.+t/Qn)
     >-(1.-x/A)*(2.-y)*t/Qn

        factor_J_p=(1.-y_p-0.5*y_p*eps)*(1.+t/Qn)
     >-(1.-x)*(2.-y_p)*t/Qn

C Lepton propagators

        aP1=-(factor_J+2.0*factor_K*cos(phi))
     >/(y*(1.+eps))

        aP2=1.+t/Qn-aP1

        aP1_p=-(factor_J_p+2.0*factor_K_p*cos(phi))
     >/(y_p*(1.+eps))

        aP2_p=1.+t/Qn-aP1_p

        call BH(x,Qn,y,y_p,t,c0_BH,c1_BH,c2_BH,
     >c0_BH_p,c1_BH_p,c2_BH_p,c0_BH_n,c1_BH_n,c2_BH_n)

        call INTERF(x,Qn,y,y_p,t,c0_I,c1_I)

        call DVCS(x,y,y_p,c0_DVCS,c0_DVCS_p)

C Note the argument shift A/(A-1)
C  V.~Guzey and M.~Strikman,``DVCS on spinless nuclear targets in impulse approximation,''  Phys.\ Rev.\  C {\bf 68}, 015204 (2003) [arXiv:hep-ph/0301216]

        call FF(131.0/130.0*t,F_A_CohEnr)
        call FF(t,F_A)       

        call FORM_FACTOR_NUCLEON(t,F1p,F2p,F1n,F2n)

        xi_N=x/(2.-x)
        xi_A=x/(2.*A-x)

        call Dual_Proton(x,Qn,t,aJu,aJd,Hscript_Im_p,
     >Hscript_Re_p,Escript_Im_p,Escript_Re_p)

        call Dual_Neutron(x,Qn,t,aJu,aJd,
     >Hscript_Im_n,Hscript_Re_n,
     >Escript_Im_n,Escript_Re_n)

C IMPORTANT!!! 
C Guzey & Teckentrup sign convention for CFF differs by an overall sign
C from BMK

C Since we use the BMK expressions for the DVCS asymmetry, we need to reverse the sign
        Hscript_Im_p=-Hscript_Im_p
        Hscript_Re_p=-Hscript_Re_p
        Escript_Im_p=-Escript_Im_p
        Escript_Re_p=-Escript_Re_p

        Hscript_Im_n=-Hscript_Im_n
        Hscript_Re_n=-Hscript_Re_n
        Escript_Im_n=-Escript_Im_n
        Escript_Re_n=-Escript_Re_n


C Building A_LU
C Note: a_sign=1 corresponds to electrons
C       a_sign=-1 corresponds to positrons


C Coherent contrib. to the numerator
        DeltaI_CohEnr=a_sign/(x*y**3*t*aP1*aP2)
     >*8.0*factor_K*y*(2.0-y)*Z*F_A_CohEnr*(2.0-x/A)/(2.0-x)
     >*((Z-1.0)*(Hscript_Im_p+t/(4.0*amn**2)*Escript_Im_p)
     >+(A-Z)*(Hscript_Im_n+t/(4.0*amn**2)*Escript_Im_n))*F_A_CohEnr
     >*sin(phi)

        DeltaI_Coh=a_sign/(x*y**3*t*aP1*aP2)
     >*8.0*factor_K*y*(2.0-y)*Z*F_A*(2.0-x/A)/(2.0-x)
     >*(Z*(Hscript_Im_p+t/(4.0*amn**2)*Escript_Im_p)
     >+(A-Z)*(Hscript_Im_n+t/(4.0*amn**2)*Escript_Im_n))*F_A
     >*sin(phi)


C Incoherent contrib. to the numerator
        DeltaI_Incoh=a_sign/(x*y_p**3*t*aP1_p*aP2_p)
     >*8.*factor_K_p*y_p*(2.-y_p)*sin(phi)
     >*(Z*(F1p*Hscript_Im_p-t/(4.*amn**2)*F2p*Escript_Im_p)
     >+(A-Z)*(F1n*Hscript_Im_n-t/(4.*amn**2)*F2n*Escript_Im_n))

C Free proton
        DeltaI_p=a_sign/(x*y_p**3*t*aP1_p*aP2_p)
     >*8.*factor_K_p*y_p*(2.-y_p)*sin(phi)
     >*(F1p*Hscript_Im_p-t/(4.*amn**2)*F2p*Escript_Im_p)

C BH Coherent
       BH_CohEnr=1.0/(x**2*y**2*t*aP1*aP2*(1.+eps)**2)
     >*Z*(Z-1.)
     >*(c0_BH+c1_BH*cos(phi)+c2_BH*cos(2.*phi))
     >*(F_A_CohEnr)**2  

        BH_Coh=1.0/(x**2*y**2*t*aP1*aP2*(1.+eps)**2)
     >*Z**2*(c0_BH+c1_BH*cos(phi)+c2_BH*cos(2.*phi))
     >*(F_A)**2  

C BH Incoherent        
        BH_Incoh=1.0/(x**2*y_p**2*t
     >*aP1_p*aP2_p*(1.+eps)**2)
     >*((c0_BH_p+c1_BH_p*cos(phi)+c2_BH_p*cos(2.*phi))*Z
     >+(A-Z)*(c0_BH_n+c1_BH_n*cos(phi)
     >+c2_BH_n*cos(2.*phi)))

C Free proton
        BH_p=1.0/(x**2*y_p**2*t*aP1_p*aP2_p*(1.+eps)**2)
     >*(c0_BH_p+c1_BH_p*cos(phi)+c2_BH_p*cos(2.*phi))


C Interf. denominator, Coherent
        aInt_CohEnr=a_sign/(x*y**3*t*aP1*aP2)
     >*Z*F_A_CohEnr*(2.0-x/A)/(2.0-x)
     >*(c0_I+c1_I*cos(phi))
     >*((Z-1.0)*(Hscript_Re_p+t/(4.0*amn**2)*Escript_Im_p)
     >+(A-Z)*(Hscript_Re_n+t/(4.0*amn**2)*Escript_Im_n))
     >*F_A_CohEnr

       aInt_Coh=a_sign/(x*y**3*t*aP1*aP2)
     >*Z*F_A*(2.0-x/A)/(2.0-x)
     >*(c0_I+c1_I*cos(phi))
     >*(Z*(Hscript_Re_p+t/(4.0*amn**2)*Escript_Im_p)
     >+(A-Z)*(Hscript_Re_n+t/(4.0*amn**2)*Escript_Im_n))*F_A 
    
C Interf. denom., Incoherent
C Proton harmonics
        c0_I_p=-8.0*(2.0-y_p)*((2.0-y_p)**2/(1.-y_p)
     >*(factor_K_p)**2
     >*(F1p*Hscript_Re_p-t/(4.*amn**2)*F2p*Escript_Re_p)
     >+t/Qn*(1.-y_p)*(2.-x)
     >*((F1p*Hscript_Re_p-t/(4.*amn**2)*F2p*Escript_Re_p)
     >-x**2/(2.0-x)**2*(F1p+F2p)*(Hscript_Re_p+Escript_Re_p)))

        c1_I_p=-8.0*factor_K_p*(2.0-2.0*y_p+y_p**2)
     >*(F1p*Hscript_Re_p-t/(4.*amn**2)*F2p*Escript_Re_p)


C Neutron harmonics
        c0_I_n=-8.0*(2.0-y_p)*((2.0-y_p)**2/(1.-y_p)
     >*(factor_K_p)**2
     >*(F1n*Hscript_Re_n-t/(4.*amn**2)*F2n*Escript_Re_n)
     >+t/Qn*(1.-y_p)*(2.-x)
     >*((F1n*Hscript_Re_n-t/(4.*amn**2)*F2n*Escript_Re_n)
     >-x**2/(2.0-x)**2*(F1n+F2n)*(Hscript_Re_n+Escript_Re_n)))

        c1_I_n=-8.0*factor_K_p*(2.0-2.0*y_p+y_p**2)
     >*(F1n*Hscript_Re_n-t/(4.*amn**2)*F2n*Escript_Re_n)


        aInt_Incoh=a_sign/(x*y_p**3*t*aP1_p*aP2_p)
     >*((c0_I_p+c1_I_p*cos(phi))*Z
     >+(c0_I_n+c1_I_n*cos(phi))*(A-Z))
    
C Free proton
        aInt_p=a_sign/(x*y_p**3*t*aP1_p*aP2_p)
     >*(c0_I_p+c1_I_p*cos(phi))

C DVCS Coh.

        DVCS_CohEnr=1.0/(y**2*Qn)
     >*(2.-x/A)**2/(2.0-x)**2
     >*c0_DVCS
     >*(Z*(Z-1.0)*((Hscript_Im_p+t/(4.0*amn**2)*Escript_Im_p)**2
     >+(Hscript_Re_p+t/(4.0*amn**2)*Escript_Re_p)**2)
     >+(A-Z)*(A-Z-1.0)*((Hscript_Im_n+t/(4.0*amn**2)*Escript_Im_n)**2
     >+(Hscript_Re_n+t/(4.0*amn**2)*Escript_Re_n)**2)
     >+2.0*Z*(A-Z)*((Hscript_Im_p+t/(4.0*amn**2)*Escript_Im_p)
     >*(Hscript_Im_n+t/(4.0*amn**2)*Escript_Im_n)
     >+(Hscript_Re_p+t/(4.0*amn**2)*Escript_Re_p)
     >*(Hscript_Re_n+t/(4.0*amn**2)*Escript_Re_n)))
     >*(F_A_CohEnr)**2


        DVCS_Coh=1.0/(y**2*Qn)
     >*(2.-x/A)**2/(2.0-x)**2
     >*c0_DVCS
     >*((Z*(Hscript_Im_p+t/(4.0*amn**2)*Escript_Im_p)
     >+(A-Z)*(Hscript_Im_n+t/(4.0*amn**2)*Escript_Im_n))**2
     >+(Z*(Hscript_Re_p+t/(4.0*amn**2)*Escript_Re_p)
     >+(A-Z)*(Hscript_Re_n+t/(4.0*amn**2)*Escript_Re_n))**2)
     >*(F_A)**2


C DVCS Incoh.
        DVCS_Incoh=1.0/(y_p**2*Qn)
     >*c0_DVCS_p
     >*(Z*(1.0/(2.0-x)**2*(4.*(1.-x)*(Hscript_Im_p**2+Hscript_Re_p**2)
     >-x**2*2.0*(Hscript_Im_p*Escript_Im_p+Hscript_Re_p*Escript_Re_p)
     >-(x**2+(2.0-x)**2*t/(4.0*amn**2))
     >*(Escript_Im_p**2+Escript_Re_p**2)))
     >+(A-Z)*(1.0/(2.0-x)**2*(4.*(1.-x)
     >*(Hscript_Im_n**2+Hscript_Re_n**2)
     >-x**2*2.0*(Hscript_Im_n*Escript_Im_n+Hscript_Re_n*Escript_Re_n)
     >-(x**2+(2.0-x)**2*t/(4.0*amn**2))
     >*(Escript_Im_n**2+Escript_Re_n**2))))

C Free proton
        DVCS_p=1.0/(y_p**2*Qn)
     >*c0_DVCS_p
     >*(1.0/(2.0-x)**2*(4.*(1.-x)*(Hscript_Im_p**2+Hscript_Re_p**2)
     >-x**2*2.0*(Hscript_Im_p*Escript_Im_p+Hscript_Re_p*Escript_Re_p)
     >-(x**2+(2.0-x)**2*t/(4.0*amn**2))
     >*(Escript_Im_p**2+Escript_Re_p**2)))
     

        A_LU_Full=(DeltaI_CohEnr+DeltaI_Incoh)
     >/
     >(BH_CohEnr+BH_Incoh+aInt_CohEnr+aInt_Incoh
     >+DVCS_CohEnr+DVCS_Incoh)

C Coherent-enriched
       A_LU_CohEnr=(DeltaI_CohEnr)
     >/
     >(BH_CohEnr+aInt_CohEnr+DVCS_CohEnr)

C Purely Coherent
       A_LU_Coh=(DeltaI_Coh)
     >/
     >(BH_Coh+aInt_Coh+DVCS_Coh)  

c Incoherent
       A_LU_Incoh=(DeltaI_Incoh)
     >/
     >(BH_Incoh+aInt_Incoh+DVCS_Incoh) 


C Only proton

       A_LU_Proton=DeltaI_p
     >/(BH_p+aInt_p+DVCS_p)


        return
        end


C BH amplitude squared: both nuclear and proton and neutron harmonics
        subroutine BH(x,Qn,y,y_p,t,c0_BH,c1_BH,c2_BH,
     >c0_BH_p,c1_BH_p,c2_BH_p,c0_BH_n,c1_BH_n,c2_BH_n)
        implicit double precision (a-h,o-z)

        amn=0.938272
        eps=4.0*amn**2*x**2/Qn

        A=131.0
        Z=54.0
        
        t_min=-Qn*(2.*(1.-x/A)*(1.-(1.+eps)**0.5)+eps)
     >/(4.*x/A*(1.-x/A)+eps)

        t_min_p=-Qn*(2.*(1.-x)*(1.-(1.+eps)**0.5)+eps)
     >/(4.*x*(1.-x)+eps)

        factor_K=(-t/Qn*(1.-x/A)*(1.-y-0.25*y**2*eps)
     >*(1.-t_min/t)*((1.+eps)**0.5+(4.*x/A*(1.-x/A)+eps)
     >/(4.*(1.-x/A))*(t-t_min)/Qn))**0.5

        factor_K_p=(-t/Qn*(1.-x)*(1.-y_p-0.25*y_p**2*eps)
     >*(1.-t_min_p/t)*((1.+eps)**0.5+(4.*x*(1.-x)+eps)
     >/(4.*(1.-x))*(t-t_min_p)/Qn))**0.5

C Spin-0 nucleus
C No nuclear FF in this expressions. The ff is added later
      c0_BH=((2.-y)**2+y**2*(1.+eps)**2)
     >*(4.*x**2*amn**2/t+4.*(1.-x/A)+(4.*x/A+eps)*t/Qn)
     >+32.*x**2*amn**2/t*(factor_K)**2
     >+2.*eps*(4.*(1.-y)*(3.+2.*eps)+y**2*(2.-eps**2))
     >-4.*(x/A)**2*(2.0-y)**2*(2.+eps)*t/Qn    
     
      c1_BH=-8.*factor_K*(2.-y)*(2.*x/A+eps-4.*x**2*amn**2/t)

      c2_BH=32.*(factor_K)**2*x**2*amn**2/t 


C Proton and neutron form factors
      call FORM_FACTOR_NUCLEON(t,F1p,F2p,F1n,F2n)

      Comb1_p=(F1p)**2-t/(4.*amn**2)*(F2p)**2
      Comb2_p=(F1p+F2p)**2            
      Comb1_n=(F1n)**2-t/(4.*amn**2)*(F2n)**2
      Comb2_n=(F1n+F2n)**2 

C Free proton       
      c0_BH_p=8.*(factor_K_p)**2*((2.+3.*eps)*Qn/t*Comb1_p
     >+2.*(x)**2*Comb2_p)
     >+(2.-y_p)**2*((2.+eps)*(4.*x**2*amn**2/t*(1.+t/Qn)**2
     >+4.0*(1.-x)*(1.+x*t/Qn))*Comb1_p
     >+4.*(x)**2*(x+(1.-x+0.5*eps)*(1.-t/Qn)**2
     >-x*(1.-2.0*x)*t**2/Qn**2)*Comb2_p)
     >+8.*(1.+eps)*(1.-y_p-0.25*y_p**2*eps)
     >*(2.0*eps*(1.-t/(4.*amn**2))*Comb1_p
     >-x**2*(1.-t/Qn)**2*Comb2_p)

      c1_BH_p=8.*factor_K_p*(2.-y_p)*((4.*x**2*amn**2/t
     >-2.*x-eps)*Comb1_p
     >+2.*x**2*(1.-(1.-2.*x)*t/Qn)*Comb2_p)

      c2_BH_p=8.*x**2*(factor_K_p)**2*(4.*(amn)**2/t*Comb1_p
     >+2.0*Comb2_p)

C Free neutron
      c0_BH_n=8.*(factor_K_p)**2*((2.+3.*eps)*Qn/t*Comb1_n
     >+2.*(x)**2*Comb2_n)
     >+(2.-y_p)**2*((2.+eps)*(4.*x**2*amn**2/t*(1.+t/Qn)**2
     >+4.0*(1.-x)*(1.+x*t/Qn))*Comb1_n
     >+4.*(x)**2*(x+(1.-x+0.5*eps)*(1.-t/Qn)**2
     >-x*(1.-2.0*x)*t**2/Qn**2)*Comb2_n)
     >+8.*(1.+eps)*(1.-y_p-0.25*y_p**2*eps)
     >*(2.0*eps*(1.-t/(4.*amn**2))*Comb1_n
     >-x**2*(1.-t/Qn)**2*Comb2_n)

      c1_BH_n=8.*factor_K_p*(2.-y_p)*((4.*x**2*amn**2/t
     >-2.*x-eps)*Comb1_n
     >+2.*x**2*(1.-(1.-2.*x)*t/Qn)*Comb2_n)

      c2_BH_n=8.*x**2*(factor_K_p)**2*(4.*amn**2/t*Comb1_n
     >+2.0*Comb2_n)

        return
        end

C Interference: only nuclear harmonics
        subroutine INTERF(x,Qn,y,y_p,t,c0_I,c1_I)
        implicit double precision (a-h,o-z)

        amn=0.938272
        eps=4.0*amn**2*x**2/Qn

        A=131.0
        Z=54.0     

        t_min=-Qn*(2.*(1.-x/A)*(1.-(1.+eps)**0.5)+eps)
     >/(4.*x/A*(1.-x/A)+eps)

        factor_K=(-t/Qn*(1.-x/A)*(1.-y-0.25*y**2*eps)
     >*(1.-t_min/t)*((1.+eps)**0.5+(4.*x/A*(1.-x/A)+eps)
     >/(4.*(1.-x/A))*(t-t_min)/Qn))**0.5

        t_min_p=-Qn*(2.*(1.-x)*(1.-(1.+eps)**0.5)+eps)
     >/(4.*x*(1.-x)+eps)

        factor_K_p=(-t/Qn*(1.-x)*(1.-y_p-0.25*y_p**2*eps)
     >*(1.-t_min_p/t)*((1.+eps)**0.5+(4.*x*(1.-x)+eps)
     >/(4.*(1.-x))*(t-t_min_p)/Qn))**0.5

C Nucleus
      c0_I=-8.*t/Qn*(2.-y)*((2.0-x/A)*(1.-y)
     >-(1.-x/A)*(2.-y)**2*(1.-t_min/t))

      c1_I=-8.*factor_K*(2.0-2.0*y+y**2)  

        return
        end

        subroutine DVCS(x,y,y_p,c0_DVCS,c0_DVCS_p)
        implicit double precision (a-h,o-z)

        A=131.0
        Z=54.0

        c0_DVCS=2.0*(2.0-2.0*y+y**2)

        c0_DVCS_p=2.0*(2.0-2.0*y_p+y_p**2)

        return
        end


C Charge FF of Xe-131
        subroutine FF(t,F_A)
        implicit double precision (a-h,o-z)
        dimension r(100),wr(100)
        pi=acos(-1.0)

        sum=0.0
        N=100
        r_min=0.0
        r_max=10.0
        call gausab(r,wr,r_min,r_max,N)
        do i=1,N
        sum=sum+sin((Abs(t))**0.5/0.197327*r(i))/
     >((Abs(t))**0.5/0.197327)
     >*rho(r(i))*r(i)*wr(i)
        enddo

        F_A=4.*pi*sum

        return
        end

        double precision function rho(rr)
        implicit double precision (a-h,o-z)

C Xe-131
        w=0.3749       
        z=2.6776
        c=5.3376
        rho0=0.00112617

        rho=rho0*(1.0+w*rr**2/c**2)/(1.0+Exp((rr**2-c**2)/z**2))

        return
        end


        subroutine FORM_FACTOR_NUCLEON(t,F1p,F2p,F1n,F2n)
        implicit double precision (a-h,o-z)

        amn=0.938272
        ammp=2.793
        ammn=-1.913

        F1p=1./(1.0-t/0.71)**2*1./(1.0-t/(4.*amn**2))
     >*(1.-t/(4.*amn**2)*2.79)

        F2p=1./(1.0-t/0.71)**2*1./(1.0-t/(4.*amn**2))
     >*(ammp-1.)

        F1n=1./(1.0-t/0.71)**2*1./(1.0-t/(4.*amn**2))
     >*(-t/(4.*amn**2)*ammn)

        F2n=1./(1.0-t/0.71)**2*1./(1.0-t/(4.*amn**2))*ammn
             
        return
        end
        
C Main subroutine for CFFs of GPDs H and E of the Proton
      subroutine Dual_Proton(xB,Q2,t,aJu,aJd,
     >Hscript_Im_p,Hscript_Re_p,Escript_Im_p,Escript_Re_p)
      implicit double precision (a-h,o-z)

      Parameter (N_E_data_files=49)
      Character FileName1(N_E_data_files)*40,FileName2*40

* Arrays of variables from the read data files
      dimension xB_array(27),Q2_array(19),t_array(13),
     >Hscript_Re_array(13,27,19),Hscript_Im_array(13,27,19)

* Temporary arrays used for 3D interpolation
      dimension Q2tmp_Im(19),Q2tmp_Re(19),xBtmp_Im(27),
     >xBtmp_Re(27),ttmp_Im(13),ttmp_Re(13)

* The same for the GPD E
      dimension Q2E_array(10),Escript_Re_array(13,27,10),
     >Escript_Im_array(13,27,10)

      dimension Q2Etmp_Im(10),Q2Etmp_Re(10),xBEtmp_Im(27),
     >xBEtmp_Re(27),tEtmp_Im(13),tEtmp_Re(13)


      data FileName1 /
     >'Escript_dual_guidal_Ju06_Jd1_2009.dat',
     >'Escript_dual_guidal_Ju06_Jd06_2009.dat',
     >'Escript_dual_guidal_Ju06_Jd02_2009.dat',
     >'Escript_dual_guidal_Ju06_Jd00_2009.dat',
     >'Escript_dual_guidal_Ju06_Jdm02_2009.dat',
     >'Escript_dual_guidal_Ju06_Jdm06_2009.dat',
     >'Escript_dual_guidal_Ju06_Jdm1_2009.dat', 
     >'Escript_dual_guidal_Ju04_Jd1_2009.dat',
     >'Escript_dual_guidal_Ju04_Jd06_2009.dat',
     >'Escript_dual_guidal_Ju04_Jd02_2009.dat',
     >'Escript_dual_guidal_Ju04_Jd00_2009.dat',
     >'Escript_dual_guidal_Ju04_Jdm02_2009.dat',
     >'Escript_dual_guidal_Ju04_Jdm06_2009.dat',
     >'Escript_dual_guidal_Ju04_Jdm1_2009.dat',  
     >'Escript_dual_guidal_Ju02_Jd1_2009.dat',
     >'Escript_dual_guidal_Ju02_Jd06_2009.dat',
     >'Escript_dual_guidal_Ju02_Jd02_2009.dat',
     >'Escript_dual_guidal_Ju02_Jd00_2009.dat',
     >'Escript_dual_guidal_Ju02_Jdm02_2009.dat',
     >'Escript_dual_guidal_Ju02_Jdm06_2009.dat',
     >'Escript_dual_guidal_Ju02_Jdm1_2009.dat',
     >'Escript_dual_guidal_Ju00_Jd1_2009.dat',
     >'Escript_dual_guidal_Ju00_Jd06_2009.dat',
     >'Escript_dual_guidal_Ju00_Jd02_2009.dat',
     >'Escript_dual_guidal_Ju00_Jd00_2009.dat',
     >'Escript_dual_guidal_Ju00_Jdm02_2009.dat',
     >'Escript_dual_guidal_Ju00_Jdm06_2009.dat',
     >'Escript_dual_guidal_Ju00_Jdm1_2009.dat',  
     >'Escript_dual_guidal_Jum02_Jd1_2009.dat',
     >'Escript_dual_guidal_Jum02_Jd06_2009.dat',
     >'Escript_dual_guidal_Jum02_Jd02_2009.dat',
     >'Escript_dual_guidal_Jum02_Jd00_2009.dat',
     >'Escript_dual_guidal_Jum02_Jdm02_2009.dat',
     >'Escript_dual_guidal_Jum02_Jdm06_2009.dat',
     >'Escript_dual_guidal_Jum02_Jdm1_2009.dat',
     >'Escript_dual_guidal_Jum04_Jd1_2009.dat',
     >'Escript_dual_guidal_Jum04_Jd06_2009.dat',
     >'Escript_dual_guidal_Jum04_Jd02_2009.dat',
     >'Escript_dual_guidal_Jum04_Jd00_2009.dat',
     >'Escript_dual_guidal_Jum04_Jdm02_2009.dat',
     >'Escript_dual_guidal_Jum04_Jdm06_2009.dat',
     >'Escript_dual_guidal_Jum04_Jdm1_2009.dat',   
     >'Escript_dual_guidal_Jum06_Jd1_2009.dat',
     >'Escript_dual_guidal_Jum06_Jd06_2009.dat',
     >'Escript_dual_guidal_Jum06_Jd02_2009.dat',
     >'Escript_dual_guidal_Jum06_Jd00_2009.dat',
     >'Escript_dual_guidal_Jum06_Jdm02_2009.dat',
     >'Escript_dual_guidal_Jum06_Jdm06_2009.dat',
     >'Escript_dual_guidal_Jum06_Jdm1_2009.dat'/  

      common /Dual_Data/xB_array,Q2_array,t_array,
     >Hscript_Re_array,Hscript_Im_array,Q2E_array,
     >Escript_Re_array,Escript_Im_array

C Reading the data file for CFF H
      open(1,file='Hscript_dual_guidal_2009.dat',
     >status='unknown')

      Nt=13
      Nx=27
      Nq=19     

      do it=1,Nt
      do ix=1,Nx
      do iq=1,Nq

      read(1,*)xB_array(ix),Q2_array(iq),t_array(it),
     >Hscript_Im_array(it,ix,iq),Hscript_Re_array(it,ix,iq)

      enddo
      enddo
      enddo

      close(1)

C Determining the appropriate data file for the given (J_u,J_d) 
C Go through 49 possibilities
      if (aJu.eq.0.6.and.aJd.eq.1.0) then
      FileName2=FileName1(1)
      else
      if (aJu.eq.0.6.and.aJd.eq.0.6) then
      FileName2=FileName1(2)
      else
      if (aJu.eq.0.6.and.aJd.eq.0.2) then
      FileName2=FileName1(3)
      else
      if (aJu.eq.0.6.and.aJd.eq.0.0) then
      FileName2=FileName1(4)
      else   
      if (aJu.eq.0.6.and.aJd.eq.-0.2) then
      FileName2=FileName1(5)
      else 
      if (aJu.eq.0.6.and.aJd.eq.-0.6) then
      FileName2=FileName1(6)
      else
      if (aJu.eq.0.6.and.aJd.eq.-1.0) then
      FileName2=FileName1(7)
      else
      if (aJu.eq.0.4.and.aJd.eq.1.0) then
      FileName2=FileName1(8)
      else
      if (aJu.eq.0.4.and.aJd.eq.0.6) then
      FileName2=FileName1(9)
      else
      if (aJu.eq.0.4.and.aJd.eq.0.2) then
      FileName2=FileName1(10)
      else
      if (aJu.eq.0.4.and.aJd.eq.0.0) then
      FileName2=FileName1(11)
      else   
      if (aJu.eq.0.4.and.aJd.eq.-0.2) then
      FileName2=FileName1(12)
      else 
      if (aJu.eq.0.4.and.aJd.eq.-0.6) then
      FileName2=FileName1(13)
      else
      if (aJu.eq.0.4.and.aJd.eq.-1.0) then
      FileName2=FileName1(14)
      else          
      if (aJu.eq.0.2.and.aJd.eq.1.0) then
      FileName2=FileName1(15)
      else
      if (aJu.eq.0.2.and.aJd.eq.0.6) then
      FileName2=FileName1(16)
      else
      if (aJu.eq.0.2.and.aJd.eq.0.2) then
      FileName2=FileName1(17)
      else
      if (aJu.eq.0.2.and.aJd.eq.0.0) then
      FileName2=FileName1(18)
      else   
      if (aJu.eq.0.2.and.aJd.eq.-0.2) then
      FileName2=FileName1(19)
      else 
      if (aJu.eq.0.2.and.aJd.eq.-0.6) then
      FileName2=FileName1(20)
      else
      if (aJu.eq.0.2.and.aJd.eq.-1.0) then
      FileName2=FileName1(21)
      else 
      if (aJu.eq.0.0.and.aJd.eq.1.0) then
      FileName2=FileName1(22)
      else
      if (aJu.eq.0.0.and.aJd.eq.0.6) then
      FileName2=FileName1(23)
      else
      if (aJu.eq.0.0.and.aJd.eq.0.2) then
      FileName2=FileName1(24)
      else
      if (aJu.eq.0.0.and.aJd.eq.0.0) then
      FileName2=FileName1(25)
      else   
      if (aJu.eq.0.0.and.aJd.eq.-0.2) then
      FileName2=FileName1(26)
      else 
      if (aJu.eq.0.0.and.aJd.eq.-0.6) then
      FileName2=FileName1(27)
      else
      if (aJu.eq.0.0.and.aJd.eq.-1.0) then
      FileName2=FileName1(28)
      else
      if (aJu.eq.-0.2.and.aJd.eq.1.0) then
      FileName2=FileName1(29)
      else
      if (aJu.eq.-0.2.and.aJd.eq.0.6) then
      FileName2=FileName1(30)
      else
      if (aJu.eq.-0.2.and.aJd.eq.0.2) then
      FileName2=FileName1(31)
      else
      if (aJu.eq.-0.2.and.aJd.eq.0.0) then
      FileName2=FileName1(32)
      else   
      if (aJu.eq.-0.2.and.aJd.eq.-0.2) then
      FileName2=FileName1(33)
      else 
      if (aJu.eq.-0.2.and.aJd.eq.-0.6) then
      FileName2=FileName1(34)
      else
      if (aJu.eq.-0.2.and.aJd.eq.-1.0) then
      FileName2=FileName1(35)
      else 
      if (aJu.eq.-0.4.and.aJd.eq.1.0) then
      FileName2=FileName1(36)
      else
      if (aJu.eq.-0.4.and.aJd.eq.0.6) then
      FileName2=FileName1(37)
      else
      if (aJu.eq.-0.4.and.aJd.eq.0.2) then
      FileName2=FileName1(38)
      else
      if (aJu.eq.-0.4.and.aJd.eq.0.0) then
      FileName2=FileName1(39)
      else   
      if (aJu.eq.-0.4.and.aJd.eq.-0.2) then
      FileName2=FileName1(40)
      else 
      if (aJu.eq.-0.4.and.aJd.eq.-0.6) then
      FileName2=FileName1(41)
      else
      if (aJu.eq.-0.4.and.aJd.eq.-1.0) then
      FileName2=FileName1(42)
      else 
       if (aJu.eq.-0.6.and.aJd.eq.1.0) then
      FileName2=FileName1(43)
      else
      if (aJu.eq.-0.6.and.aJd.eq.0.6) then
      FileName2=FileName1(44)
      else
      if (aJu.eq.-0.6.and.aJd.eq.0.2) then
      FileName2=FileName1(45)
      else
      if (aJu.eq.-0.6.and.aJd.eq.0.0) then
      FileName2=FileName1(46)
      else   
      if (aJu.eq.-0.6.and.aJd.eq.-0.2) then
      FileName2=FileName1(47)
      else 
      if (aJu.eq.-0.6.and.aJd.eq.-0.6) then
      FileName2=FileName1(48)
      else
      if (aJu.eq.-0.6.and.aJd.eq.-1.0) then
      FileName2=FileName1(49)
      else 
      write(6,*) 'Ju and Jd out or range'  

        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif


      open(2,file=FileName2,status='unknown')

      Nt=13
      Nx=27
      NqE=10     

      do it=1,Nt
      do ix=1,Nx
      do iq=1,NqE

      read(2,*)a,Q2E_array(iq),a,
     >Escript_Im_array(it,ix,iq),Escript_Re_array(it,ix,iq)

      enddo
      enddo
      enddo

      close(2)

             
* 3D interpolation done separately in all three variables
      do it=1,Nt
      do ix=1,Nx

      do iq=1,Nq
       
      Q2tmp_Im(iq)=Hscript_Im_array(it,ix,iq)
      Q2tmp_Re(iq)=Hscript_Re_array(it,ix,iq)   
      enddo

      do iq=1,NqE
       
      Q2Etmp_Im(iq)=Escript_Im_array(it,ix,iq)
      Q2Etmp_Re(iq)=Escript_Re_array(it,ix,iq)   
      enddo

      call linear(Nq,Q2_array,Q2tmp_Im,Q2,xBtmp_Im(ix)) 
      call linear(Nq,Q2_array,Q2tmp_Re,Q2,xBtmp_Re(ix))   

      call linear(NqE,Q2E_array,Q2Etmp_Im,Q2,xBEtmp_Im(ix)) 
      call linear(NqE,Q2E_array,Q2Etmp_Re,Q2,xBEtmp_Re(ix)) 

      enddo

      call linear(Nx,xB_array,xBtmp_Im,xB,ttmp_Im(it))  
      call linear(Nx,xB_array,xBtmp_Re,xB,ttmp_Re(it)) 

      call linear(Nx,xB_array,xBEtmp_Im,xB,tEtmp_Im(it))  
      call linear(Nx,xB_array,xBEtmp_Re,xB,tEtmp_Re(it)) 

      enddo

      call linear(Nt,t_array,ttmp_Im,t,Hscript_Im_p) 
      call linear(Nt,t_array,ttmp_Re,t,Hscript_Re_p)  

      call linear(Nt,t_array,tEtmp_Im,t,Escript_Im_p) 
      call linear(Nt,t_array,tEtmp_Re,t,Escript_Re_p)

      return
      end

* 1D linear interpolation
      subroutine linear(n,xa,ya,x,y)      
      implicit double precision (a-h,o-z)
      integer n,ns,i

      dimension xa(n),ya(n)

C just in case there is some rounding error
      y=0.0
      
C Later added by T. Teckentrup      
      ns=0
      
      dif=abs(x-xa(1))

      do i=1,n-1

C Finding the closest index i
      dift=abs(x-xa(i)) 
      if (dift.le.dif) then
      ns=i
      dif=dift
      endif
      enddo     

C x could be greater or smaller than xa(i)
      if (x.gt.xa(ns)) then   

      y1=ya(ns)
      y2=ya(ns+1)
      x1=xa(ns)
      x2=xa(ns+1) 

      y=((y1-y2)*x+(y2*x1-y1*x2))/(x1-x2)

      else

      y1=ya(ns-1)
      y2=ya(ns)
      x1=xa(ns-1)
      x2=xa(ns) 

      y=((y1-y2)*x+(y2*x1-y1*x2))/(x1-x2)   

      endif

      return
      end

C Main subroutine for CFFs of GPDs H and E of the Neutron
      subroutine Dual_Neutron(xB,Q2,t,aJu,aJd,
     >Hscript_Im_n,Hscript_Re_n,Escript_Im_n,Escript_Re_n)
      implicit double precision (a-h,o-z)

      Parameter (N_E_data_files=49)
      Character FileName1_n(N_E_data_files)*46,
     >FileName2_n*46

* Arrays of variables from the read data files
      dimension xB_array_n(27),Q2_array_n(19),t_array_n(13),
     >Hscript_Re_array_n(13,27,19),Hscript_Im_array_n(13,27,19)

* Temporary arrays used for 3D interpolation
      dimension Q2tmp_Im_n(19),Q2tmp_Re_n(19),xBtmp_Im_n(27),
     >xBtmp_Re_n(27),ttmp_Im_n(13),ttmp_Re_n(13)

* The same for the GPD E
      dimension Q2E_array_n(10),Escript_Re_array_n(13,27,10),
     >Escript_Im_array_n(13,27,10)

      dimension Q2Etmp_Im_n(10),Q2Etmp_Re_n(10),xBEtmp_Im_n(27),
     >xBEtmp_Re_n(27),tEtmp_Im_n(13),tEtmp_Re_n(13)

      data FileName1_n /
     >'Escript_dual_guidal_Ju06_Jd1_neutron_2009.dat',
     >'Escript_dual_guidal_Ju06_Jd06_neutron_2009.dat',
     >'Escript_dual_guidal_Ju06_Jd02_neutron_2009.dat',
     >'Escript_dual_guidal_Ju06_Jd00_neutron_2009.dat',
     >'Escript_dual_guidal_Ju06_Jdm02_neutron_2009.dat',
     >'Escript_dual_guidal_Ju06_Jdm06_neutron_2009.dat',
     >'Escript_dual_guidal_Ju06_Jdm1_neutron_2009.dat', 
     >'Escript_dual_guidal_Ju04_Jd1_neutron_2009.dat',
     >'Escript_dual_guidal_Ju04_Jd06_neutron_2009.dat',
     >'Escript_dual_guidal_Ju04_Jd02_neutron_2009.dat',
     >'Escript_dual_guidal_Ju04_Jd00_neutron_2009.dat',
     >'Escript_dual_guidal_Ju04_Jdm02_neutron_2009.dat',
     >'Escript_dual_guidal_Ju04_Jdm06_neutron_2009.dat',
     >'Escript_dual_guidal_Ju04_Jdm1_neutron_2009.dat',  
     >'Escript_dual_guidal_Ju02_Jd1_neutron_2009.dat',
     >'Escript_dual_guidal_Ju02_Jd06_neutron_2009.dat',
     >'Escript_dual_guidal_Ju02_Jd02_neutron_2009.dat',
     >'Escript_dual_guidal_Ju02_Jd00_neutron_2009.dat',
     >'Escript_dual_guidal_Ju02_Jdm02_neutron_2009.dat',
     >'Escript_dual_guidal_Ju02_Jdm06_neutron_2009.dat',
     >'Escript_dual_guidal_Ju02_Jdm1_neutron_2009.dat',
     >'Escript_dual_guidal_Ju00_Jd1_neutron_2009.dat',
     >'Escript_dual_guidal_Ju00_Jd06_neutron_2009.dat',
     >'Escript_dual_guidal_Ju00_Jd02_neutron_2009.dat',
     >'Escript_dual_guidal_Ju00_Jd00_neutron_2009.dat',
     >'Escript_dual_guidal_Ju00_Jdm02_neutron_2009.dat',
     >'Escript_dual_guidal_Ju00_Jdm06_neutron_2009.dat',
     >'Escript_dual_guidal_Ju00_Jdm1_neutron_2009.dat',  
     >'Escript_dual_guidal_Jum02_Jd1_neutron_2009.dat',
     >'Escript_dual_guidal_Jum02_Jd06_neutron_2009.dat',
     >'Escript_dual_guidal_Jum02_Jd02_neutron_2009.dat',
     >'Escript_dual_guidal_Jum02_Jd00_neutron_2009.dat',
     >'Escript_dual_guidal_Jum02_Jdm02_neutron_2009.dat',
     >'Escript_dual_guidal_Jum02_Jdm06_neutron_2009.dat',
     >'Escript_dual_guidal_Jum02_Jdm1_neutron_2009.dat',
     >'Escript_dual_guidal_Jum04_Jd1_neutron_2009.dat',
     >'Escript_dual_guidal_Jum04_Jd06_neutron_2009.dat',
     >'Escript_dual_guidal_Jum04_Jd02_neutron_2009.dat',
     >'Escript_dual_guidal_Jum04_Jd00_neutron_2009.dat',
     >'Escript_dual_guidal_Jum04_Jdm02_neutron_2009.dat',
     >'Escript_dual_guidal_Jum04_Jdm06_neutron_2009.dat',
     >'Escript_dual_guidal_Jum04_Jdm1_neutron_2009.dat',   
     >'Escript_dual_guidal_Jum06_Jd1_neutron_2009.dat',
     >'Escript_dual_guidal_Jum06_Jd06_neutron_2009.dat',
     >'Escript_dual_guidal_Jum06_Jd02_neutron_2009.dat',
     >'Escript_dual_guidal_Jum06_Jd00_neutron_2009.dat',
     >'Escript_dual_guidal_Jum06_Jdm02_neutron_2009.dat',
     >'Escript_dual_guidal_Jum06_Jdm06_neutron_2009.dat',
     >'Escript_dual_guidal_Jum06_Jdm1_neutron_2009.dat'/  

      common /Dual_Data/xB_array_n,Q2_array_n,t_array_n,
     >Hscript_Re_array_n,Hscript_Im_array_n,Q2E_array_n,
     >Escript_Re_array_n,Escript_Im_array_n

C Reading the data file for CFF H
      open(11,file='Hscript_dual_guidal_neutron_2009.dat',
     >status='unknown')

      Nt=13
      Nx=27
      Nq=19     

      do it=1,Nt
      do ix=1,Nx
      do iq=1,Nq

      read(11,*)xB_array_n(ix),Q2_array_n(iq),t_array_n(it),
     >Hscript_Im_array_n(it,ix,iq),Hscript_Re_array_n(it,ix,iq)

      enddo
      enddo
      enddo

      close(11)

C Determining the appropriate data file for the given (J_u,J_d) 
C Go through 49 possibilities
      if (aJu.eq.0.6.and.aJd.eq.1.0) then
      FileName2_n=FileName1_n(1)
      else
      if (aJu.eq.0.6.and.aJd.eq.0.6) then
      FileName2_n=FileName1_n(2)
      else
      if (aJu.eq.0.6.and.aJd.eq.0.2) then
      FileName2_n=FileName1_n(3)
      else
      if (aJu.eq.0.6.and.aJd.eq.0.0) then
      FileName2_n=FileName1_n(4)
      else   
      if (aJu.eq.0.6.and.aJd.eq.-0.2) then
      FileName2_n=FileName1_n(5)
      else 
      if (aJu.eq.0.6.and.aJd.eq.-0.6) then
      FileName2_n=FileName1_n(6)
      else
      if (aJu.eq.0.6.and.aJd.eq.-1.0) then
      FileName2_n=FileName1_n(7)
      else
      if (aJu.eq.0.4.and.aJd.eq.1.0) then
      FileName2_n=FileName1_n(8)
      else
      if (aJu.eq.0.4.and.aJd.eq.0.6) then
      FileName2_n=FileName1_n(9)
      else
      if (aJu.eq.0.4.and.aJd.eq.0.2) then
      FileName2_n=FileName1_n(10)
      else
      if (aJu.eq.0.4.and.aJd.eq.0.0) then
      FileName2_n=FileName1_n(11)
      else   
      if (aJu.eq.0.4.and.aJd.eq.-0.2) then
      FileName2_n=FileName1_n(12)
      else 
      if (aJu.eq.0.4.and.aJd.eq.-0.6) then
      FileName2_n=FileName1_n(13)
      else
      if (aJu.eq.0.4.and.aJd.eq.-1.0) then
      FileName2_n=FileName1_n(14)
      else          
      if (aJu.eq.0.2.and.aJd.eq.1.0) then
      FileName2_n=FileName1_n(15)
      else
      if (aJu.eq.0.2.and.aJd.eq.0.6) then
      FileName2_n=FileName1_n(16)
      else
      if (aJu.eq.0.2.and.aJd.eq.0.2) then
      FileName2_n=FileName1_n(17)
      else
      if (aJu.eq.0.2.and.aJd.eq.0.0) then
      FileName2_n=FileName1_n(18)
      else   
      if (aJu.eq.0.2.and.aJd.eq.-0.2) then
      FileName2_n=FileName1_n(19)
      else 
      if (aJu.eq.0.2.and.aJd.eq.-0.6) then
      FileName2_n=FileName1_n(20)
      else
      if (aJu.eq.0.2.and.aJd.eq.-1.0) then
      FileName2_n=FileName1_n(21)
      else 
      if (aJu.eq.0.0.and.aJd.eq.1.0) then
      FileName2_n=FileName1_n(22)
      else
      if (aJu.eq.0.0.and.aJd.eq.0.6) then
      FileName2_n=FileName1_n(23)
      else
      if (aJu.eq.0.0.and.aJd.eq.0.2) then
      FileName2_n=FileName1_n(24)
      else
      if (aJu.eq.0.0.and.aJd.eq.0.0) then
      FileName2_n=FileName1_n(25)
      else   
      if (aJu.eq.0.0.and.aJd.eq.-0.2) then
      FileName2_n=FileName1_n(26)
      else 
      if (aJu.eq.0.0.and.aJd.eq.-0.6) then
      FileName2_n=FileName1_n(27)
      else
      if (aJu.eq.0.0.and.aJd.eq.-1.0) then
      FileName2_n=FileName1_n(28)
      else
      if (aJu.eq.-0.2.and.aJd.eq.1.0) then
      FileName2_n=FileName1_n(29)
      else
      if (aJu.eq.-0.2.and.aJd.eq.0.6) then
      FileName2_n=FileName1_n(30)
      else
      if (aJu.eq.-0.2.and.aJd.eq.0.2) then
      FileName2_n=FileName1_n(31)
      else
      if (aJu.eq.-0.2.and.aJd.eq.0.0) then
      FileName2_n=FileName1_n(32)
      else   
      if (aJu.eq.-0.2.and.aJd.eq.-0.2) then
      FileName2_n=FileName1_n(33)
      else 
      if (aJu.eq.-0.2.and.aJd.eq.-0.6) then
      FileName2_n=FileName1_n(34)
      else
      if (aJu.eq.-0.2.and.aJd.eq.-1.0) then
      FileName2_n=FileName1_n(35)
      else 
      if (aJu.eq.-0.4.and.aJd.eq.1.0) then
      FileName2_n=FileName1_n(36)
      else
      if (aJu.eq.-0.4.and.aJd.eq.0.6) then
      FileName2_n=FileName1_n(37)
      else
      if (aJu.eq.-0.4.and.aJd.eq.0.2) then
      FileName2_n=FileName1_n(38)
      else
      if (aJu.eq.-0.4.and.aJd.eq.0.0) then
      FileName2_n=FileName1_n(39)
      else   
      if (aJu.eq.-0.4.and.aJd.eq.-0.2) then
      FileName2_n=FileName1_n(40)
      else 
      if (aJu.eq.-0.4.and.aJd.eq.-0.6) then
      FileName2_n=FileName1_n(41)
      else
      if (aJu.eq.-0.4.and.aJd.eq.-1.0) then
      FileName2_n=FileName1_n(42)
      else 
       if (aJu.eq.-0.6.and.aJd.eq.1.0) then
      FileName2_n=FileName1_n(43)
      else
      if (aJu.eq.-0.6.and.aJd.eq.0.6) then
      FileName2_n=FileName1_n(44)
      else
      if (aJu.eq.-0.6.and.aJd.eq.0.2) then
      FileName2_n=FileName1_n(45)
      else
      if (aJu.eq.-0.6.and.aJd.eq.0.0) then
      FileName2_n=FileName1_n(46)
      else   
      if (aJu.eq.-0.6.and.aJd.eq.-0.2) then
      FileName2_n=FileName1_n(47)
      else 
      if (aJu.eq.-0.6.and.aJd.eq.-0.6) then
      FileName2_n=FileName1_n(48)
      else
      if (aJu.eq.-0.6.and.aJd.eq.-1.0) then
      FileName2_n=FileName1_n(49)
      else 
      write(6,*) 'Ju and Jd out or range'  

        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif
        endif


      open(12,file=FileName2_n,status='unknown')

      Nt=13
      Nx=27
      NqE=10     

      do it=1,Nt
      do ix=1,Nx
      do iq=1,NqE

      read(12,*)a,Q2E_array_n(iq),a,
     >Escript_Im_array_n(it,ix,iq),Escript_Re_array_n(it,ix,iq)

      enddo
      enddo
      enddo

      close(12)

             
* 3D interpolation done separately in all three variables
      do it=1,Nt
      do ix=1,Nx

      do iq=1,Nq
       
      Q2tmp_Im_n(iq)=Hscript_Im_array_n(it,ix,iq)
      Q2tmp_Re_n(iq)=Hscript_Re_array_n(it,ix,iq)   
      enddo

      do iq=1,NqE
       
      Q2Etmp_Im_n(iq)=Escript_Im_array_n(it,ix,iq)
      Q2Etmp_Re_n(iq)=Escript_Re_array_n(it,ix,iq)   
      enddo

      call linear(Nq,Q2_array_n,Q2tmp_Im_n,Q2,xBtmp_Im_n(ix)) 
      call linear(Nq,Q2_array_n,Q2tmp_Re_n,Q2,xBtmp_Re_n(ix))   

      call linear(NqE,Q2E_array_n,Q2Etmp_Im_n,Q2,xBEtmp_Im_n(ix)) 
      call linear(NqE,Q2E_array_n,Q2Etmp_Re_n,Q2,xBEtmp_Re_n(ix)) 

      enddo

      call linear(Nx,xB_array_n,xBtmp_Im_n,xB,ttmp_Im_n(it))  
      call linear(Nx,xB_array_n,xBtmp_Re_n,xB,ttmp_Re_n(it)) 

      call linear(Nx,xB_array_n,xBEtmp_Im_n,xB,tEtmp_Im_n(it))  
      call linear(Nx,xB_array_n,xBEtmp_Re_n,xB,tEtmp_Re_n(it)) 

      enddo

      call linear(Nt,t_array_n,ttmp_Im_n,t,Hscript_Im_n) 
      call linear(Nt,t_array_n,ttmp_Re_n,t,Hscript_Re_n)  

      call linear(Nt,t_array_n,tEtmp_Im_n,t,Escript_Im_n) 
      call linear(Nt,t_array_n,tEtmp_Re_n,t,Escript_Re_n)

      return
      end

C***********************************************************************
     
C Integration

      SUBROUTINE GAUSAB( X, W, A, B, NPT )

C=======================================================================
C                               /b
C     Gaussian quadratures for  | f(x)dx.   NPT must be even & <100.
C                              /a
C
C=======================================================================

      implicit double precision (a-h,o-z)
      dimension X(1), W(1)
      
    
      COMMON / GAUSCM / XX(100), WW(100)
C-----------------------------------------------------------------------
      NPT2 = NPT/2

      IF( A .EQ. 0. .AND. B .GT. .5E20 ) THEN
      CALL GAUS(NPT,1)
      DO 5 N = 1, NPT
      X(N) = XX(N)
      W(N) = WW(N)
5     CONTINUE
      RETURN
      END IF

      CALL GAUS(NPT,0)

      C1 = (B-A)/2.
      C2 = (B+A)/2.

      DO 10 N = 1, NPT2
      M = NPT2 + 1 - N
      X(N) = - XX(M)*C1 + C2
      W(N) =   WW(M)*C1
10    CONTINUE

      DO 20 N = NPT2+1, NPT
      M = N - NPT2
      X(N) = XX(M)*C1 + C2
      W(N) = WW(M)*C1
20    CONTINUE

      RETURN

      END
C=======================================================================
        SUBROUTINE GAUS(NN,ITYPE)
C
C       STANDARD SETTINGS ARE: IALF=IBTA=0,0<NN<100,
C               ITYPE=0 : -1 TO +1 POINTS
C               ITYPE=1 :  0 TO INFINITY
C
      implicit double precision (a-h,o-z)

      DATA IALF, IBTA / 0, 0 /
C***********************************************************************
C
C     GAUJAC COMPUTES GAUSS - JACOBI INTEGRATION WEIGHTS AND NODES
C
C***********************************************************************
C
      DIMENSION X(100),A(100),B(100),C(100)
      DIMENSION XX(100), WW(100)
      DIMENSION G(100)
C
        COMMON/GAUSCM/XX,WW
C
      PI=3.141592653589793238D0
C
      ISUM=IALF+IBTA
      IDIF=IBTA-IALF
      FNUM=ISUM*IDIF
      ALF=IALF
      BTA=IBTA
C
C     CALCULATE COEFFICIENTS REQUIRED BY JACOBI
C
      DO 20 I=1,100
      I2=I+I
      IM=I-1
      DEN=(ISUM+I2)*(ISUM+I2-2)
      IF(ISUM.EQ.0) GO TO 67
      B(I)=FNUM/DEN
      GO TO 68
67    B(I)=0.D0
68    FNUMP=4*IM*(IALF+IM)*(IBTA+IM)*(ISUM+IM)
      DENP=(ISUM+I2-1)*(ISUM+I2-2)**2*(ISUM+I2-3)
      IF(((ISUM.EQ.0).OR.(ISUM.EQ.1)).AND.(I.EQ.1)) GO TO 69
      C(I)=FNUMP/DENP
      GO TO 20
69    C(I)=0.D0
20    CONTINUE

C-----------------------------------------------------------------------
C     For outputting quadratures
C-----------------------------------------------------------------------
C      TYPE 34
C34    FORMAT(1X,'FILENAME FOR INTERMEDIATE RESULTS =? (A10)')
C      ACCEPT 35,FNAMX
C35    FORMAT(A10)
C      OPEN(UNIT=3,FILE=FNAMX,STATUS='NEW')
C-----------------------------------------------------------------------

      CALL JACOBI(NN,X,G,ALF,BTA,B,C,EPS,CSX,CSA,TSX,TSA)

C-----------------------------------------------------------------------
C     For outputting quadratures
C-----------------------------------------------------------------------
C      TYPE 70
C70    FORMAT(1X,'OUTPUT FILENAME ? (A10)')
C      ACCEPT 71,FNAM
C71    FORMAT(A10)
C      OPEN(UNIT=7,FILE=FNAM,STATUS='NEW')
C-----------------------------------------------------------------------

C     OUTPUT DATA STATEMENT OF WEIGHTS AND NODES
C
      IF(ITYPE.EQ.0) GO TO 130
      DO 100 I=1,NN
      II=NN+1-I
C     XX(I)=(1.0D+00+X(II))/(1.0D+00-X(II))
C     WW(I)=2.0D+00*G(II)/((1.0D+00-X(II))**(IALF+2)*
C    1     (1.0D+00+X(II))**IBTA)
      ZZ=PI*(1.0D+00+X(II))/4.0D+00
      XX(I)=SIN(ZZ)/COS(ZZ)
      WW(I)=G(II)*PI*0.25D+00/(COS(0.25D+00*PI*(1.0D+00
     1      +X(II)))**2)

C-----------------------------------------------------------------------
C      WRITE(7,120)I,I,XX(I),WW(I)
120   FORMAT(6X,6HDATAX(,I2,6H,M),W(,I2,4H,M)/,D22.15,1H,,
     1      D22.15,1H/)
C-----------------------------------------------------------------------

100   CONTINUE
        RETURN
 130  CONTINUE
      NI=NN/2
      DO 140 I=1,NI
      II=NI+1-I
      XX(I)=X(II)
      WW(I)=G(II)

C-----------------------------------------------------------------------
C      WRITE(7,120)I,I,XX(I),WW(I)
C-----------------------------------------------------------------------

140   CONTINUE
      RETURN
      END


      SUBROUTINE JACOBI(NN,X,A,ALF,BTA,B,C,EPS,CSX,CSA,TSX,TSA)
c      IMPLICIT REAL*8 (A-H,O-Z)
        implicit double precision (a-h,o-z)
C


C***********************************************************************
C
C     CALCULATES THE ZEROS X(I) OF THE NN-TH ORDER JACOBI POLYNOMIAL
C     PN(ALF,BTA) FOR THE SEGMENT (-1,1)
C     THE LARGEST ZERO WILL BE STORED IN X(1).
C     ALSO CALCULATES THE CORRESPONDING COEFFICIENTS A(I) OF THE
C     NN-TH ORDER GAUSS-JACOBI QUADRATURE FORMULA OF DEGREE 2*NN-1.
C     THIS SUBROUTINE MUST BE GIVEN THE COEFFICIENTS:
C
C                (ALF+BTA)(BTA-ALF)
C        B(N) = -----------------------------
C               (ALF+BTA+2N)(ALF+BTA+2N-2)
C
C                4(N-1)(ALF+N-1)(BTA+N-1)(ALF+BTA+N-1)
C        C(N) = ---------------------------------------------
C               (ALF+BTA+2N-1)(ALF+BTA+2N-2)**2(ALF+BTA+2N-3)
C
C
C     IN THE RECURSION RELATION
C
C        P(N) = (X-B(N))*P(N-1)-C(N)*P(N-2)
C
C     FOR ALL N LESS THAN OR EQUAL TO THE HIGHEST DEGREE NN.
C
C        CSX = CALC SUM X(I)      TSX = TRUE SUM X(I)
C        CSA = CALC SUM A(I)      TSA = TRUE SUM A(I)
C
C***********************************************************************
C
      DIMENSION X(100),A(100),B(100),C(100)
      FN=NN
      CSX=0.D0
      CSA=0.D0
      BETA=EXP(FLGAMA(ALF+1.D0)+FLGAMA(BTA+1.D0)-FLGAMA(ALF+BTA+2.D0))
      CC=2.**(ALF+BTA+1.)*BETA
      TSX=FN*(BTA-ALF)/(ALF+BTA+2.*FN)
      TSA=CC
      DO 1 J=2,NN
1     CC=CC*C(J)
      DO 12 I=1,NN
      IF(I-1)12,2,3
C        LARGEST ZERO
2     AN=ALF/FN
      BN=BTA/FN
      R1=(1.D0+ALF)*(2.78/(4.+FN*FN)+.768*AN/FN)
      R2=1.D0+1.48*AN+.96*BN+.452*AN*AN+.83*AN*BN
      XT=1.D0-R1/R2
      GO TO 11
3     IF(I-2)12,4,5
C        SECOND ZERO
4     R1=(4.1+ALF)/((1.+ALF)*(1.+.156*ALF))
      R2=1.D0+.06*(FN-8.)*(1.D0+.12*ALF)/FN
      R3=1.D0+.012*BTA*(1.D0+.25*ABS(ALF))/FN
      RATIO=R1*R2*R3
      XT=XT-RATIO*(1.D0-XT)
      GO TO 11
5     IF(I-3) 12,6,7
C        THIRD ZERO
6     R1=(1.67+.28*ALF)/(1.D0+.37*ALF)
      R2=1.D0+.22*(FN-8.D0)/FN
      R3=1.D0+8.*BTA/((6.28D0+BTA)*FN*FN)
      RATIO=R1*R2*R3
      XT=XT-RATIO*(X(1)-XT)
      GO TO 11
7     IF(NN-I-1)10,9,8
C        MIDDLE ZEROS
8     XT=3.*X(I-1)-3.*X(I-2)+X(I-3)
      GO TO 11
C        SECOND LAST ZERO
9     R1=(1.D0+.235*BTA)/(.766D0+.119*BTA)
      R2=1./(1.D0+.639*(FN-4.D0)/(1.D0+.71*(FN-4.D0)))
      R3=1./(1.D0+20.*ALF/((7.5D0+ALF)*FN*FN))
      RATIO=R1*R2*R3
      XT=XT+RATIO*(XT-X(I-2))
      GO TO 11
C        LAST ZERO
10    R1=(1.D0+.37*BTA)/(1.67D0+.28*BTA)
      R2=1./(1.D0+.22*(FN-8.D0)/FN)
      R3=1./(1.D0+8.*ALF/((6.28D0+ALF)*FN*FN))
      RATIO=R1*R2*R3
      XT=XT+RATIO*(XT-X(I-2))
C
11    CALL GSROOT(XT,NN,ALF,BTA,DPN,PN1,B,C,EPS)
      X(I)=XT
      A(I)=CC/(DPN*PN1)
C-----------------------------------------------------------------------
C      WRITE(3,20)ALF,BTA,NN,I,XT,A(I)
C-----------------------------------------------------------------------

      CSX=CSX+XT
      CSA=CSA+A(I)
12    CONTINUE

C-----------------------------------------------------------------------
C      WRITE(3,22)CSX,CSA,TSX,TSA
20    FORMAT(1X,2F6.2,2I3,2(1X,D26.18),1X,(1X,D26.18))
22    FORMAT(1H0,/,5X,'CSX =',D25.18,5X,'CSA =',D25.18,/,
     1 5X,'TSX =',D25.18,5X,'TSA =',D25.18)
C-----------------------------------------------------------------------

      RETURN
      END

      SUBROUTINE GSROOT(X,NN,ALF,BTA,DPN,PN1,B,C,EPS)
c      IMPLICIT REAL*8 (A-H,O-Z)
        implicit double precision (a-h,o-z)


C***********************************************************************
C
C     IMPROVES THE APPROXIMATE ROOT X
C     IN ADDITION WE ALSO OBTAIN:
C      DPN = DERIVATIVE OF P(N) AT X
C      PN1 = VALUE OF P(N-1) AT X
C
C***********************************************************************
C
      DIMENSION B(100),C(100)
      ITER=0
1     ITER=ITER+1
      CALL RECR(P,DP,PN1,X,NN,ALF,BTA,B,C)
      D=P/DP
      X=X-D
      IF(ABS(D)-EPS)3,3,2
2     IF(ITER-10)1,3,3
3     DPN=DP
      RETURN
      END

      SUBROUTINE RECR(PN,DPN,PN1,X,NN,ALF,BTA,B,C)
c      IMPLICIT REAL*8 (A-H,O-Z)
        implicit double precision (a-h,o-z)


      DIMENSION B(100),C(100)
      P1=1.D0
      P=X+(ALF-BTA)/(ALF+BTA+2.D0)
      DP1=0.D0
      DP=1.D0
      DO 1 J=2,NN
      Q=(X-B(J))*P-C(J)*P1
      DQ=(X-B(J))*DP+P-C(J)*DP1
      P1=P
      P=Q
      DP1=DP
      DP=DQ
1     CONTINUE
      PN=P
      DPN=DP
      PN1=P1
      RETURN
      END

       double precision FUNCTION FLGAMA(W)
c      IMPLICIT REAL*8 (A-H,O-Z)
        implicit double precision (a-h,o-z)

C
C***********************************************************************
C
C     CALCULATES LOG(BASE E)GAMMA(W) FOR W REAL AND GAMMA(W) POSITIVE.
C     USES STIRLING'S APPROXIMATION
C     ACCURATE TO ABOUT 12 SIGNIFICANT PLACES
C
C***********************************************************************
C
      PI=3.141592653589793238D0
      X=W
      M=0
      FK=-1.D0
      IF(X-.5D0)1,1,2
C        W LESS EQ .5
1     M=1
      XPI=X*PI
      X=1.D0-X
2     FK=FK+1.D0
      IF(X+FK-6.D0)2,2,3
3     Z=X+FK
      ZZ=Z*Z
C        LOG GAMMA(Z), Z GREATER 6.
      Y=(Z-.5)*LOG(Z)-Z+.9189385332047+(((((((-10861851./ZZ
     1 +2356200.)/ZZ-704820.)/ZZ+309400.)/ZZ-218790.)/ZZ
     2 +291720.)/ZZ-1021020.)/ZZ+30630600.)/Z/367567200.
      IF(FK)6,6,4
4     IK=FK
      DO 5 I=1,IK
      FK=FK-1.
5     Y=Y-LOG(X+FK)
6     IF(M)7,11,7
7     P=PI/SIN(XPI)
      IF(P)8,8,10
8     WRITE(6,9)W
9     FORMAT(2X,6HGAMMA(D11.4,13H) IS NEGATIVE)
      Y=0.
      GO TO 11
10    Y=LOG(P)-Y
11    FLGAMA=Y
      RETURN
      END
c========

c****** Interpolation
c  N: number of data points, (X(1),Y(1)), ....., (X(N),Y(N))
c  At the point xx, the interpolated value, yy, is obtained
c-----------------------------------------------------------
c  call spline(N,X,Y,B,C,D) (only once this procedure is necessary !)
c   
c  yy=seval(N,xx,X,Y,B,C,D)
c-----------------------------------------------------------
c
c***************************************************************************
C ---------------------------------------------------------------------
      SUBROUTINE SPLINE(N,X,Y,B,C,D)
C ---------------------------------------------------------------------
c***************************************************************************
C CALCULATE THE COEFFICIENTS B,C,D IN A CUBIC SPLINE INTERPOLATION.
C INTERPOLATION SUBROUTINES ARE TAKEN FROM
C G.E. FORSYTHE, M.A. MALCOLM AND C.B. MOLER,
C COMPUTER METHODS FOR MATHEMATICAL COMPUTATIONS (PRENTICE-HALL, 1977).
c      IMPLICIT REAL*8(A-H,O-Z)
      implicit double precision (A-H,O-Z)
      DIMENSION X(5000),Y(5000),B(5000),C(5000),D(5000)
      NM1=N-1
      IF(N.LT.2) RETURN
      IF(N.LT.3) GO TO 250
      D(1)=X(2)-X(1)
      C(2)=(Y(2)-Y(1))/D(1)
      DO 210 I=2,NM1
        D(I)=X(I+1)-X(I)
        B(I)=2.0D0*(D(I-1)+D(I))
        C(I+1)=(Y(I+1)-Y(I))/D(I)
        C(I)=C(I+1)-C(I)
 210             CONTINUE
      B(1)=-D(1)
      B(N)=-D(N-1)
      C(1)=0.0D0
      C(N)=0.0D0
      IF(N.EQ.3) GO TO 215
      C(1)=C(3)/(X(4)-X(2))-C(2)/(X(3)-X(1))
      C(N)=C(N-1)/(X(N)-X(N-2))-C(N-2)/(X(N-1)-X(N-3))
      C(1)=C(1)*D(1)**2.0D0/(X(4)-X(1))
      C(N)=-C(N)*D(N-1)**2.0D0/(X(N)-X(N-3))
 215       CONTINUE
      DO 220 I=2,N
        T=D(I-1)/B(I-1)
        B(I)=B(I)-T*D(I-1)
        C(I)=C(I)-T*C(I-1)
 220             CONTINUE
      C(N)=C(N)/B(N)
      DO 230 IB=1,NM1
        I=N-IB
        C(I)=(C(I)-D(I)*C(I+1))/B(I)
 230             CONTINUE
      B(N)=(Y(N)-Y(NM1))/D(NM1)+D(NM1)*(C(NM1)+2.0D0*C(N))
      DO 240 I=1,NM1
        B(I)=(Y(I+1)-Y(I))/D(I)-D(I)*(C(I+1)+2.0D0*C(I))
        D(I)=(C(I+1)-C(I))/D(I)
        C(I)=3.0D0*C(I)
 240             CONTINUE
      C(N)=3.0D0*C(N)
      D(N)=D(N-1)
      RETURN
 250       CONTINUE
      B(1)=(Y(2)-Y(1))/(X(2)-X(1))
      C(1)=0.0D0
      D(1)=0.0D0
      B(2)=B(1)
      C(2)=0.0D0
      D(2)=0.0D0
      RETURN
      END
c
c***************************************************************************
C ---------------------------------------------------------------------
c      REAL*8 FUNCTION SEVAL(N,XX,X,Y,B,C,D)
       double precision FUNCTION SEVAL(N,XX,X,Y,B,C,D)
C ---------------------------------------------------------------------
c***************************************************************************
C CALCULATE THE DISTRIBUTION AT XX BY CUBIC SPLINE INTERPOLATION.
c      IMPLICIT REAL*8(A-H,O-Z)
       IMPLICIT double precision (A-H,O-Z)
      DIMENSION X(5000),Y(5000),B(5000),C(5000),D(5000)
      DATA I/1/
      IF(I.GE.N) I=1
      IF(XX.LT.X(I)) GO TO 310
      IF(XX.LE.X(I+1)) GO TO 330
 310       CONTINUE
      I=1
      J=N+1
 320       CONTINUE
      K=(I+J)/2
      IF(XX.LT.X(K)) J=K
      IF(XX.GE.X(K)) I=K
      IF(J.GT.I+1) GO TO 320
 330       CONTINUE
      DX=XX-X(I)
      SEVAL=Y(I)+DX*(B(I)+DX*(C(I)+DX*D(I)))
      RETURN
      END