module spherical_harmonics
    implicit none
    integer, parameter :: dp = kind(1.0d0)
    ! Use numerical constant for PI for maximum compatibility
    real(dp), parameter :: PI = 3.14159265358979323846264338327950288_dp

contains

    ! Calculates the spherical harmonics Y_lm(theta, phi)
    ! l: orbital angular momentum quantum number
    ! m: magnetic quantum number
    ! theta: polar angle (radians)
    ! phi: azimuthal angle (radians)
    function ylm(l, m, theta, phi) result(res)
        integer, intent(in) :: l, m
        real(dp), intent(in) :: theta, phi
        complex(kind=dp) :: res

        real(dp) :: p_val
        real(dp) :: cos_theta

        if (abs(m) > l) then
            res = cmplx(0.0_dp, 0.0_dp, dp)
            return
        end if

        cos_theta = cos(theta)

        ! Calculate normalized associated Legendre polynomials N_lm * P_lm(cos_theta)
        p_val = normalized_legendre_p(l, abs(m), cos_theta)

        if (m >= 0) then
            res = p_val * cmplx(cos(m*phi), sin(m*phi), dp)
        else
            ! Use the property Y_{l,-m} = (-1)^m * conjg(Y_lm)
            if (mod(m, 2) == 0) then
                res = p_val * cmplx(cos(m*phi), sin(m*phi), dp)
            else
                res = -p_val * cmplx(cos(m*phi), sin(m*phi), dp)
            end if
        end if

    end function ylm

    ! Calculate associated Legendre polynomials with normalization factors using recurrence relations
    ! N_lm = sqrt((2l+1)/4pi * (l-m)!/(l+m)!)
    function normalized_legendre_p(l, m, x) result(res)
        integer, intent(in) :: l, m
        real(dp), intent(in) :: x
        real(dp) :: res

        real(dp) :: p_mm, p_mp1m, p_ll
        real(dp) :: somx2
        integer :: i

        ! Assume m >= 0
        somx2 = sqrt((1.0_dp - x)*(1.0_dp + x))

        ! Calculation of P_mm
        p_mm = 1.0_dp / sqrt(4.0_dp * PI)
        do i = 1, m
            p_mm = -p_mm * somx2 * sqrt(real(2*i+1, dp) / real(2*i, dp))
        end do

        if (l == m) then
            res = p_mm
            return
        end if

        ! Calculation of P_{m+1, m}
        p_mp1m = x * sqrt(real(2*m+3, dp)) * p_mm

        if (l == m + 1) then
            res = p_mp1m
            return
        end if

        ! Calculation of P_ll using recurrence relations
        ! Normalized form of (l-m)P_l = x(2l-1)P_{l-1} - (l+m-1)P_{l-2}
        do i = m + 2, l
            p_ll = (x * sqrt(real(4*i**2 - 1, kind=dp)) * p_mp1m &
                   - sqrt(real((2*i+1)*(i+m-1)*(i-m-1), kind=dp) / real(2*i-3, kind=dp)) * p_mm) &
                   / sqrt(real(i**2 - m**2, kind=dp))
            p_mm = p_mp1m
            p_mp1m = p_ll
        end do

        res = p_mp1m

    end function normalized_legendre_p

end module spherical_harmonics
!! Usage example
!program test_ylm
!    use spherical_harmonics
!    implicit none
!    integer :: l, m
!    real(dp) :: theta, phi
!    complex(kind=dp) :: y
!
!    l = 5
!    m = 1
!    theta = PI / 4.0_dp  ! 45 degrees
!    phi = PI / 3.0_dp    ! 60 degrees
!
!    y = ylm(l, m, theta, phi)
!
!    print *, "Testing Y_lm for l=", l, " m=", m
!    print *, "theta=", theta, " phi=", phi
!    print *, "Result Y_lm:", y
!    print *, "Magnitude:", abs(y)
!
!end program test_ylm