phase

module global_vars
    integer :: N = 100, l = 2
    real :: R = 10.0, V_0 = 2, a = 5.0, lk = 0.1, rk = 10.0, dk = 0.05, k
    save
end module

function V(x)
    use global_vars
    real :: V
    real, intent(in) :: x
    if (x < a) then
        V = -V_0
    else
        V = 0.0
    end if
    return
end function

recursive function jl(x, l) result(ans)
    real, intent(in) :: x
    real :: ans
    integer, intent(in) :: l
    if (l == 0) then
        ans = sin(x) / x
    else if (l == 1) then
        ans = sin(x) / x**2 - cos(x) / x
    else
        ans = (2 * (l - 1) + 1) * jl(x, l - 1) / x  - jl(x, l - 2)
    end if
end function

function jl_hat(x, l)
    real, intent(in) :: x
    integer, intent(in) :: l
    real :: jl, jl_hat
    jl_hat = jl(x, l) * x
    return
end function

recursive function nl(x, l) result(ans)
    real, intent(in) :: x
    real :: ans
    integer, intent(in) :: l
    if (l == 0) then
        ans = -cos(x) / x
    else if (l == 1) then
        ans = -cos(x) / x**2 - sin(x) / x
    else
        ans = (2 * (l - 1) + 1) * nl(x, l - 1) / x  - nl(x, l - 2)
    end if
end function

function nl_hat(x, l)
    real, intent(in) :: x
    integer, intent(in) :: l
    real :: nl, nl_hat
    nl_hat = nl(x, l) * x
end function

function dfac(x)
    integer, intent(in) :: x
    integer :: i, dfac
    dfac = mod(x, 2)
    i = dfac + 2
    do while (i <= x)
        dfac = dfac * i
        i = i + 2
    end do
    return
end function

function F(x, delta)
    use global_vars
    real :: nl, jl, jl_hat, nl_hat, F, V
    real, intent(in) :: x, delta
    F = -V(x) * (cos(delta) * jl_hat(k * x, l) - sin(delta) * nl_hat(k * x, l))**2 / k
    return
end function

subroutine RungeKutta(h, N, Y)
    implicit none
    integer :: i, N
    real :: k, F, h, x, S(4), Y(N)
    do i = 1, N - 1
        x = h * i
        S(1) =       F(x          , Y(i)                 )
        S(2) = 2.0 * F(x + h / 2.0, Y(i) + h * S(1) / 2.0)
        S(3) = 2.0 * F(x + h / 2.0, Y(i) + h * S(2) / 2.0)
        S(4) =       F(x + h      , Y(i) + h * S(3)      )
        Y(i + 1) = Y(i) + h * sum(S) / 6.0
    end do
end subroutine

program phase
    use global_vars
    implicit none
    integer :: i, out = 42, dfac
    real :: V, h
    real, allocatable :: DELTA(:)
    allocate(DELTA(N))
    h = R / N
    k = lk
    open(unit = out, file = "output.txt")
    write(out, *) "k               delta"
    do while (k <= rk)
        DELTA(1) = -V(0.0) * (k * h)**(2 * l + 1) * h**2 / ((2 * l + 3) * dfac(2*l + 1)**2)
        call RungeKutta(h, N, DELTA)
        write(out, *) k, DELTA(N)
        k = k + dk
    end do
    close(out)
end program