dissipation

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

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

function F(x)
    use global_vars
    implicit none
    real :: F, V
    real, intent(in) :: x
    F = V(x) + l * (l + 1) / x**2 - k**2
    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) / x * jl(x, l - 1) - 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

subroutine Simpson(h, f, N, ans)
    implicit none
    integer :: i, N
    real :: h
    complex :: f(N), ans
    ans = 0.0
    do i = 2, N - 1, 2
        ans = ans + f(i - 1) + 4.0 * f(i) + f(i + 1)
    end do
    ans = ans * h / 3.0
end subroutine

program scattering
    use global_vars
    implicit none
    integer :: i, j, out = 42
    real :: F, V, h, jl_hat,  pi = 3.1415926
    complex :: Al, delta
    integer, allocatable :: ipiv(:)
    complex, allocatable :: X(:),  mA(:, :)
    allocate(X(N))
    allocate(ipiv(N))
    allocate(mA(N, N))
    h = R / N
    k = lk
    open(unit = out, file = "output.txt")
    do while (k <= rk)
        !mA = -d^2/dr^2 + F(r) - k^2
        mA(:, :) = 0.0
        do i = 1, N
            if(i > 1) mA(i, i - 1) = 12.0 / h**2  - F(h * (i - 1))
            mA(i, i) = -24.0 / h**2 - 10.0 * F(h * i)
            if(i < N) mA(i, i + 1) = 12.0 / h**2 - F(h * (i + 1))
        end do
        mA(N, N) = mA(N, N) + exp((0.0, 1.0) * k * h) * (12.0 / h**2 - F(h * (N + 1)))
        !X = -V(r)jl(kr)
        do i = 1, N
            X(i) = 10.0 * V(h * i) * jl_hat(k * i * h, l)
            if(i > 1) X(i) = X(i) + V(h * (i - 1)) * jl_hat(k * (i - 1) * h, l)
            if(i < n) X(i) = X(i) + V(h * (i + 1)) * jl_hat(k * (i + 1) * h, l)
        end do
        !X = chi_l
        call cgesv(N, 1, mA, N, ipiv, X, N, i)
        !X = jl(kr)V(r)f_l(r, k) = jl(kr)V(r)(jl(kr) + chi_l(r, k))
        do i =1,n
            X(i) = jl_hat(k * h * i, l) * V(h * i) * (jl_hat(k * h * i, l) + X(i))
        end do
        call Simpson(h, X, N, Al)
        Al = -1.0 * Al / k**2
        delta = log(1 + (0.0, 2.0) * k * Al) / (0.0, 2.0)
        write(out, *) k, real(delta)
        k = k + dk
    end do
    close(out)
end program