program integrate
real*8 n, a, b, k

n = 1000
k = 1

a = 0
b = 1

call int_hopkins(a, b, k, n)

end program integrate

real*8 function f(x)
    real*8 x
    f = x
    return
end function f

real*8 function g(x)
    real*8 x
    g = x
    return
end function g

real*8 function dg(x)
    real*8 x
    dg = 0*x + 1
    return
end function

real*8 function sinc(x)
    real*8 x
    if(x < 0.000001) then
        sinc = 1
    else
        sinc = dsin(x) / x
    end if
    return
end function sinc

subroutine int_hopkins(a, b, k, nx)
    real*8 a, b, nx, k, x, px
    real*8 f, g, dg
    real*8 h, xi, sx

    complex*8 I

    h = (b - a) / nx

    px = a
    x = px + h

    I = 0
    do while (x < b)
        xi = (px + x) / 2
        sx = k * dg(xi) * (x - px)

        I = I + (x - px) * f(xi) * sinc(sx) * cdexp(CMPLX(0,1) * k * g(xi))

        px = x
        x = x + h
        end do

    write(*,*) I
    return
end subroutine int_hopkins