lsum Function

  function lsum ( lmin, order )
  !! Function that calculates the sum of higher order terms, \(c_n(\ell)\),
  !! provided by [[ldep]] from \(\ell_{\mathrm{min}}\) to \(\infty\).

    integer(ip), intent(in) :: lmin
    !! The value of \(\ell\) at which to start the sum,
    !! \(\ell_{\mathrm{min}}\).
    integer(ip), intent(in) :: order
    !! The order of the higher order term, \(n\).
    real(wp) :: lsum
    !! The return value is \[\sum_{\ell=\ell_{\mathrm{min}}}^{\infty} c_n(\ell) =
    !!   \begin{cases}
    !!     \frac{\ell_{\mathrm{min}}}{4\ell_{\mathrm{min}}^2-1}, & n=2 \\
    !!     \frac{\ell_{\mathrm{min}}}{48\ell_{\mathrm{min}}^4-
    !!       120\ell_{\mathrm{min}}^2+27}, & n=3 \\
    !!     \frac{\ell_{\mathrm{min}}}{5 (64\ell_{\mathrm{min}}^6-
    !!       560\ell_{\mathrm{min}}^4+1036\ell_{\mathrm{min}}^2-225)}, & n=4
    !!   \end{cases} \].
    real(wp) :: lm

    lm = real(lmin,wp)
    select case (order)
    case (2)
      lsum = lm/(4.0_wp*lm**2-1.0_wp)
    case (4)
      lsum = lm/(48.0_wp*lm**4-120.0_wp*lm**2+27.0_wp)
    case (6)
      lsum = lm/(5.0_wp &
                 *(64.0_wp*lm**6-560.0_wp*lm**4+1036.0_wp*lm**2-225.0_wp))
    case default
      stop "error: lsum called with unsupported order"
    end select

  end function lsum