C *********************************************
C
C Code to calculate A_UT of the proton and
C the sin(phi_S-phi)*cos(phi)-moment of A_UT for proton 
C This moment of A_UT is sensitive to GPDs H and E 
C The code also gives the unpolarized cross section
C 
C Two cases are considered:
C-------------------------
C 1) All harmonics for BH, DVCS and Interference are kept -> A_UT_Proton_Full
C 2) All harmonics for BH and only LT harmonics for DVCS and Interference
C    terms are kept (c0_DVCS, c0_I, c1_I and s1_I) -> A_UT_Proton_LT
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 January 2009
C
C *********************************************

      program AUT_Proton            
      implicit double precision (a-h,o-z)
C The notation to remind of the Trento sign convention
      dimension a_phi_Trento(20),w_phi(20)

        pi=acos(-1.0)
        amn=0.938272

C External kinematics
C Positrons: a_sign=-1; electrons: a_sign=+1

C As an example, we use a typical HERMES kinematics
      a_sign=-1.0          
      E_lepton=27.5
      write(6,65) 'Beam energy=',E_lepton
 65   format(A,2f6.2)   
      
      Qn=2.5
      write(6,65) 'Q**2=',Qn  
      xB=0.10
      write(6,65) 'x_B=',xB  
      t=-0.1
      write(6,65) 't=',t

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
   
        N_phi=20
        phi_min=0.0
        phi_max=2.0*pi

        sum_phi_full=0.0
        sum_phi_lt=0.0
        sum_cs_full=0.0

C Gaussian integration routine
        call gausab(a_phi_Trento,w_phi,phi_min,phi_max,N_phi)
        do i=1,N_phi

C Note that the notations of BMK differ from the Trento convention
C phi_Trento=pi-phi 

        phi=pi-a_phi_Trento(i)      
        call ASYMMETRY_AUT_Proton
     >(a_sign,E_lepton,xB,Qn,t,phi,aJu,aJd,
     >A_UT_Proton_Full,A_UT_Proton_LT,CS_Full)

       sum_phi_full=sum_phi_full
     >+A_UT_Proton_Full*cos(a_phi_Trento(i))*w_phi(i)

        sum_phi_lt=sum_phi_lt
     >+A_UT_Proton_LT*cos(a_phi_Trento(i))*w_phi(i)

        sum_cs_full=sum_cs_full+CS_Full*w_phi(i)

        enddo

C There is an additional minus sign since phi_S_Trento-phi_Trento=pi+phi_BMK
C and, hence, sin(phi_S_Trento-phi_Trento)=-sin(phi_BMK)
  
        A_UT_Proton_Full_sin_cos=-sum_phi_full/pi
        A_UT_Proton_LT_sin_cos=-sum_phi_lt/pi
        UNPOL_CS_Full=sum_cs_full


        write(6,65) 'Full A_UT_sin(phi_S-phi)_cos(phi)=',
     >A_UT_Proton_Full_sin_cos
      
        write(6,65) 'LT A_UT_sin(phi_S-phi)_cos(phi)=',
     >A_UT_Proton_LT_sin_cos

        write(6,65) 'UNPOL. CROSS SECTION [nb/GeV**2]=',
     >UNPOL_CS_Full

      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 (we use the Trento convention)
C aJu and aJd = total angular momentum carried by the u-quark and the d-quark

C NOTE: A_UT is sin(phi_S-phi)-moment of A_LU, which still depends on 
C the angle phi (between the leptonic and nadronic planes) 

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

        amn=0.938272
        pi=acos(-1.0)

C Total invariant energy
        s=2.0*E_lepton*amn+amn**2

        y=Qn/(2.0*amn*x*E_lepton)
        eps=4.0*amn**2*x**2/Qn

        if (y.gt.0.95) then
        write(6,*) 'y is larger 0.95'
        endif        
        
C Minimal t
        t_min=-Qn*(2.*(1.-x)*(1.-(1.+eps)**0.5)+eps)
     >/(4.*x*(1.-x)+eps)

C Kinematic K-factor

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

C Kinematic J-factor (needed for lepton propagators)
        factor_J=(1.-y-0.5*y*eps)*(1.+t/Qn)
     >-(1.-x)*(2.-y)*t/Qn

C Lepton propagators

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

        aP2=1.+t/Qn-aP1

C We keep only the leading twist harmonics
        call BH_Proton(x,Qn,y,t,c0_BH_p,c1_BH_p,c2_BH_p)

        call DVCS_Proton(x,y,c0_DVCS_p)

        call FORM_FACTOR_NUCLEON(t,F1p,F2p,F1n,F2n)
 
        call Dual_Proton(x,Qn,t,aJu,aJd,Hscript_Im_p,Hscript_Re_p,
     >Escript_Im_p,Escript_Re_p)

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

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

C Numerator: Interference

        CI_TP=1.0/(2.-x)*(x**2*F1p-(1.0-x)*t/amn**2*F2p)*Hscript_Im_p
     >+(x**2/(2.-x)*F1p+(t/(4.0*amn**2))*((2.-x)*F1p+x**2/(2.-x)*F2p))
     >*Escript_Im_p

        delta_CI_TP=t/amn**2*(F2p*Hscript_Im_p-F1p*Escript_Im_p)

        c0_I_TP=8.0*amn*(1.-y)**0.5*factor_K/Qn**0.5*(2.0-y)
     >*((2.0-y)**2/(1.0-y)*CI_TP+delta_CI_TP)

        c1_I_TP=8.0*amn*(1.0-y)**0.5/Qn**0.5*(2.-2.*y+y**2)*CI_TP

        c2_I_TP=16.0*amn*(1.-y)**0.5*factor_K/(Qn**0.5*(2.-x))
     >*(2.-y)*CI_TP*(-x)

        aInt_TP=a_sign*1.0/(x*y**3*t*aP1*aP2)     
     >*(c0_I_TP+c1_I_TP*cos(phi)+c2_I_TP*cos(2.0*phi))

        aInt_TP_LT=a_sign*1.0/(x*y**3*t*aP1*aP2)     
     >*(c0_I_TP+c1_I_TP*cos(phi))

C Numerator: DVCS

        CDVCS_TP=2./(2.-x)*2.0*(Hscript_Im_p*Escript_Re_p
     >- Hscript_Re_p*Escript_Im_p)
        
        c0_DVCS_TP=-Qn**0.5*factor_K/(amn*(1.0-y)**0.5)
     >*(2.0-2.0*y+y**2)*CDVCS_TP

        c1_DVCS_TP=-4.0*Qn**0.5*(factor_K)**2
     >/(amn*(2.0-x)*(1.-y)**0.5)*(2.0-y)*CDVCS_TP*(-x)

        aDVCS_TP=1.0/(y**2*Qn)
     >*(c0_DVCS_TP+c1_DVCS_TP*cos(phi))

        aDVCS_TP_LT=1.0/(y**2*Qn)
     >*(c0_DVCS_TP)

C Denominator: BH 
       BH_p=1.0/(x**2*y**2*t*aP1*aP2*(1.0+eps)**2)
     >*(c0_BH_p+c1_BH_p*cos(phi)+c2_BH_p*cos(2.*phi)) 
        

C Denominator: Interference
        c0_I_p=-8.0*(2.0-y)*((2.0-y)**2/(1.-y)*(factor_K)**2
     >*(F1p*Hscript_Re_p-t/(4.*amn**2)*F2p*Escript_Re_p)
     >+t/Qn*(1.-y)*(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*(2.0-2.0*y+y**2)
     >*(F1p*Hscript_Re_p-t/(4.*amn**2)*F2p*Escript_Re_p)

        c2_I_p=-16.*(factor_K)**2/(2.0-x)*(2.-y)
     >*(F1p*Hscript_Re_p-t/(4.*amn**2)*F2p*Escript_Re_p)*(-x) 

        aInt_p=a_sign*1.0/(x*y**3*t*aP1*aP2)     
     >*(c0_I_p+c1_I_p*cos(phi)+c2_I_p*cos(2.0*phi))
        
        aInt_p_LT=a_sign*1.0/(x*y**3*t*aP1*aP2)     
     >*(c0_I_p+c1_I_p*cos(phi))

C Denominator: DVCS
        DVCS_p=1.0/(y**2*Qn)
     >*(c0_DVCS_p+8.*factor_K/(2.-x)*(2.-y)*(-x)*cos(phi))
     >*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))

        DVCS_p_LT=1.0/(y**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_UT_Proton_Full=(aInt_TP+aDVCS_TP)
     >/(BH_p+aInt_p+DVCS_p) 

        A_UT_Proton_LT=(aInt_TP_LT+aDVCS_TP_LT)
     >/(BH_p+aInt_p_LT+DVCS_p_LT)

C Special conditions 
        if (A_UT_Proton_Full.gt.1.0) then
        write(6,*) 'A_UT_Proton_Full > 1 !'
        endif

        if (A_UT_Proton_LT.gt.1.0) then
        write(6,*) 'A_UT_Proton_LT > 1 !'
        endif
        
        if (BH_p+aInt_p+DVCS_p.lt.0.0) then
        write(6,*) 'FULL CROSS SECTION IS NEGATIVE !'
        endif
        

        if (BH_p+aInt_p_LT+DVCS_p_LT.lt.0.0) then
        write(6,*) 'LT CROSS SECTION IS NEGATIVE !'
        endif

C Unpolarized cross section
        alpha_em=1./137.0
        CS_Full=(alpha_em**3*x*y)/(8.0*pi*Qn*(1.0+eps)**0.5)
     >*(BH_p+aInt_p+DVCS_p)
     >*0.389379*(10.0**6)
C [Cross section] = nbarn/GeV**2

        return
        end


C BH amplitude squared
        subroutine BH_Proton(x,Qn,y,t,c0_BH_p,c1_BH_p,c2_BH_p)
        implicit double precision (a-h,o-z)

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

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

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

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            

C Free proton       
      c0_BH_p=8.*(factor_K)**2*((2.+3.*eps)*Qn/t*Comb1_p
     >+2.*x**2*Comb2_p)
     >+(2.-y)**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-0.25*y**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*(2.-y)*((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)**2*(4.*amn**2/t*Comb1_p
     >+2.0*Comb2_p)

        return
        end


        subroutine DVCS_Proton(x,y,c0_DVCS_p)
        implicit double precision (a-h,o-z)

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

        return
        end

C Proton and neutron elastic form factors
        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
      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
c        write(6,*) 'J_u and J_d out of range'
        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***********************************************************************
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