program Pi0Phot
! List of subroutines
!   DiffCross(Eg,t, DSG)
!       return the differential cross section from F amplitudes
!   Ampl_F(s,t, Fvec)
!       return the F amplitudes
!   Cos2T(Eg,cost,mass, t,coef,t0,t1)
!       return cos for given t
!   T2Cos(Eg,t,mass, cos,coef)
!       return t for given cos

    implicit double precision (a-h,o-z)
    double complex                  :: II = (0.d0,1.0d0)
    double precision, dimension(4)  :: mass
    double precision, dimension(6)  :: paralp, parcou
    integer                         :: infile=9, outfile=10
    open (infile,  FILE = 'PARAM.TXT' , FORM = 'FORMATTED',STATUS = 'OLD')
    open (outfile, FILE = 'output.txt' , FORM = 'FORMATTED',STATUS = 'UNKNOWN')

    PI=3.141592653589793238462643383279502884197
    xmg   = 0.0     ! photon mass
    xmpi = 0.138    ! pion mass
    xmp  = 0.938    ! proton mass
    mass = (/ xmg,xmp,xmpi,xmp /)

    ! parameters: trajectories and couplings
    ! the (linear) trajectories are
    !   vector pole : alp1 = paralp(1) + t*paralp(2)
    !   vector cut  : alp2 = paralp(3) + t*paralp(4)
    !   axial  pole : alp3 = paralp(5) + t*paralp(6)
    ! the couplings are (p=pole and c=cut)
    !   parcoup = (/ g1p, g4p, g1c, g4c, g2p, g2c /)
    paralp = (/ 0.442457d0, 1.09991d0, 0.46181d0, 0.16654d0, -0.193332d0, 1.0213d0  /)
    parcou = (/ 3.8873d0, -10.1403d0, -1.76187d0, -3.58089d0, -8.887d0, 0.0d0/)


    read(infile,*) Eg, ist, vmin, vmax, dv
    write(*,'(a20,f6.3,a6)') 'Beam energy Eg = ', Eg, ' GeV'
    Print*, ''
    Print*, 'Units: t (GeV^2),   Dsig/Dt (micro barn/GeV^2),  Dsig/DOmega (micro barn)  '
    Print*, ''
    Print*, '    t         cos        Dsig/Dt        Dsig/DOmega  '

    if (ist==1) then
    ! the X variable is t
        ! CHECK: if the tmin and tmax are outside [t1,t0]
        call Cos2T(Eg,1.0d0,mass, t,coef,t0,t1)
        if (vmin<t1) then
            vmin = t1
        endif
        if (vmax>t0) then
            vmax = t0
        endif
        ! end of CHECK
        nb = nint((vmax-vmin)/dv) + 1
        t = vmax
        do i=1,nb
            call DiffCross(Eg, t, dsgdt)
            call T2Cos(Eg,t,mass, cost,coef)
            write(*,'(2f10.3, 2e15.3)') t,cost, dsgdt, dsgdt/coef
            write(outfile,'(2f10.3, 2e15.3)') t,cost, dsgdt, dsgdt/coef
            t = t - dv
        end do
    elseif (ist==0) then
    ! the X variable is cost
        nb = nint((vmax-vmin)/dv) + 1
        cost = vmin
        do i=1,nb
            call Cos2T(Eg,cost,mass, t,coef,t0,t1)
            call DiffCross(Eg, t, dsgdt)
            write(*,'(2f10.3, 2e15.3)') t,cost, dsgdt, dsgdt/coef
            write(outfile,'(2f10.3, 2e15.3)') t,cost, dsgdt, dsgdt/coef
            cost = cost + dv
        end do
    endif

    ! make the plots
    call SYSTEM('/opt/local/bin/gnuplot gnuplot.txt')

    close(infile)
    close(outfile)

contains

subroutine DiffCross(Eg,t, DSG)
! return the differential cross section
! Input:
!   Eg      beam energy in the lab frame
!   t       the madelstam variable
! Output:
!   DSG     differential cross section
    implicit double precision (a-h,o-z)
    double complex tc,kt,pt,cost, f01p,f11p,f01m,f11m, res1, Fvec
    dimension Fvec(4)

    s = xmp**2 + 2*xmp*Eg
    som = 2*xmp**2 + xmpi**2
    tc = complex(t,0)
    kt = (t - xmpi**2)/(2*sqrt(tc))
    pt = sqrt(tc-4*xmp**2)/2.0
    cost = (2*s + t - som)/(4*kt*pt)
    sint2 = 1- abs(cost)**2

    call Ampl_F(s,t, Fvec)

    f01p = kt*sqrt(tc)/2.0*Fvec(1)
    f11p = kt/2.0*Fvec(3);
    f01m = pt*kt*Fvec(2);
    f11m = -pt*kt*sqrt(tc)*Fvec(4)

    res1 = 2*(abs(f01p)**2  + abs(f01m)**2)*abs(sint2)
    res1 = res1  + abs(f11p - f11m)**2*(1 - cost)**2 + abs(f11p + f11m)**2*(1 + cost)**2
    DSG = real(res1)/(64*PI*xmp**2*Eg**2)*389.4

end subroutine

subroutine Ampl_F(s,t, Fvec)
! rerurn the F amplitudes
! Input:
!   s       the madelstam variable
!   t       the madelstam variable
! Output:
!   Fvec    vector with the four F amplitudes
!           in term of CGLN A amplitudes
!           F1 =    - A1 + 2 MP A2
!           F2 =      A1 +    t A2
!           F3 = 2 MP A1 -    t A4
!           F4 =      A3
    implicit double precision (a-h,o-z)
    double complex   :: Fvec
    dimension Fvec(4)
    double complex F1,F2,F3,F4, xi1,xi2,xi3

    g1p = parcou(1)
    g4p = parcou(2)
    g1c = parcou(3)
    g4c = parcou(4)
    g2p = parcou(5)
    g2c = parcou(6)

    alp1 = paralp(1) + paralp(2)*t      ! vector pole trajectory
    alp2 = paralp(3) + paralp(4)*t      ! vector cut trajectory
    alp3 = paralp(5) + paralp(6)*t      ! axial pole trajectory

    xi1 = gamma(1 - alp1)*( 1 - exp(-II*PI*alp1) )/2.d0*s**(alp1-1)
    xi2 = gamma(1 - alp2)*( 1 - exp(-II*PI*alp2) )/2.d0*s**(alp2-1)/log(s)
    xi3 = gamma(1 - alp3)*( 1 - exp(-II*PI*alp3) )/2.d0*s**(alp3-1)

    F1 = -(t*g1p - 2*xmp*g4p)*xi1 - (t*g1c - 2*xmp*g4c)*xi2
    F2 = t*g2p*xi3
    F3 = t*(2*xmp*g1p - g4p)*xi1 + t*(2*xmp*g1c - g4c)*xi2
    F4 = 0

    Fvec = (/ F1, F2, F3, F4 /)
end subroutine





subroutine Cos2T(Eg,cost,mass, t,coef,t0,t1)
! return cos for given t
! Input:
!   Eg      beam energy in the lab frame
!   cost    cosine of the scattering angle in the center of mass frame
!   mass    vector with the four masses for beam + target --> meson  + recoil
!           mass(1) = beam
!           mass(2) = target
!           mass(3) = meson
!           mass(4) = recoil
! Output:
!   t       the madelstam variable
!   coef    the Jacobian coef= PI/(qi qf)
!   t0      t1 < t < t0
!   t1      t0 and t1 are negatif |t0|<|t1|

    implicit double precision (a-h,o-z)
    double precision, dimension(4):: mass
    xma = mass(1)
    xmb = mass(2)
    xmc = mass(3)
    xmd = mass(4)

    s = xma**2 + xmb**2 + 2*xmb*Eg
    qi = sqrt( (s - (xma + xmb)**2)*(s - (xma - xmb)**2)/(4.0*s) )
    qf = sqrt( (s - (xmc + xmd)**2)*(s - (xmc - xmd)**2)/(4.0*s) )
    som = xma**2 + xmb**2 + xmc**2 + xmd**2
    t = 2.0*qi*qf*cost + som/2.0 - s/2.0 - (xma**2 - xmb**2)*(xmc**2 - xmd**2)/(2.0*s)
    coef = PI/(qi*qf)
    t0 = (xma**2 - xmb**2 - xmc**2 + xmd**2)**2/4./s - (qi - qf)**2
    t1 = (xma**2 - xmb**2 - xmc**2 + xmd**2)**2/4./s - (qi + qf)**2
end subroutine

subroutine T2Cos(Eg,t,mass, cos,coef)
! return t for given cos
! Input:
!   Eg      beam energy in the lab frame
!   t       the madelstam variable
!   mass    vector with the four masses for beam + target --> meson  + recoil
!           mass(1) = beam
!           mass(2) = target
!           mass(3) = meson
!           mass(4) = recoil
! Output:
!   cost    cosine of the scattering angle in the center of mass frame
!   coef    the Jacobian coef= PI/(qi qf)

    implicit double precision (a-h,o-z)
    double precision, dimension(4):: mass
    xma = mass(1)
    xmb = mass(2)
    xmc = mass(3)
    xmd = mass(4)

    s = xma**2 + xmb**2 + 2*xmb*Eg
    qi = sqrt( (s - (xma + xmb)**2)*(s - (xma - xmb)**2)/(4.0*s) )
    qf = sqrt( (s - (xmc + xmd)**2)*(s - (xmc - xmd)**2)/(4.0*s) )
    som = xma**2 + xmb**2 + xmc**2 + xmd**2
    cos = (t + (xma**2 - xmb**2)*(xmc**2 - xmd**2)/(2.0*s) - som/2.0 + s/2.0 )/(2.0*qi*qf)
    coef = PI/(qi*qf)

end subroutine


end program Pi0Phot
