def legendre_lm(l, m, z, normalization='4pi', csphase=1, cnorm=0):
    """
    Compute the associated Legendre function for a specific degree and order.

    Usage
    -----
    plm = legendre_lm (l, m, z, [normalization, csphase, cnorm])

    Returns
    -------
    plm : float
        The associated Legendre functions for degree l and order m.

    Parameters
    ----------
    l : integer
        The spherical harmonic degree.
    m : integer
        The spherical harmonic order.
    z : float
        The argument of the associated Legendre functions.
    normalization : str, optional, default = '4pi'
        '4pi', 'ortho', 'schmidt', or 'unnorm' for use with geodesy 4pi
        normalized, orthonormalized, Schmidt semi-normalized, or unnormalized
        spherical harmonic functions, respectively.
    csphase : optional, integer, default = 1
        If 1 (default), the Condon-Shortley phase will be excluded. If -1, the
        Condon-Shortley phase of (-1)^m will be appended to the associated
        Legendre functions.
    cnorm : optional, integer, default = 0
        If 1, the complex normalization of the associated Legendre functions
        will be used. The default is to use the real normalization.

    Description
    -----------
    legendre_lm will calculate the associated Legendre function for a specific
    degree l and order m. The Legendre functions are used typically as a part
    of the spherical harmonic functions, and three parameters determine how
    they are defined. normalization can be either '4pi' (default), 'ortho',
    'schmidt', or 'unnorm' for use with 4pi normalized, orthonormalized,
    Schmidt semi-normalized, or unnormalized spherical harmonic functions,
    respectively. csphase determines whether to include or exclude (default)
    the Condon-Shortley phase factor. cnorm determines whether to normalize
    the Legendre functions for use with real (default) or complex spherical
    harmonic functions.

    The Legendre functions are calculated using the standard three-term
    recursion formula, and in order to prevent overflows, the scaling approach
    of Holmes and Featherstone (2002) is utilized. The resulting functions are
    accurate to about degree 2800. See Wieczorek and Meschede (2018) for exact
    definitions on how the Legendre functions are defined.

    References
    ----------
    Holmes, S. A., and W. E. Featherstone, A unified approach to the Clenshaw
    summation and the recursive computation of very high degree and order
    normalised associated Legendre functions, J. Geodesy, 76, 279-299,
    doi:10.1007/s00190-002-0216-2, 2002.

    Wieczorek, M. A., and M. Meschede. SHTools — Tools for working with
    spherical harmonics, Geochem., Geophys., Geosyst., 19, 2574-2592,
    doi:10.1029/2018GC007529, 2018.
    """
    if l < 0:
        raise ValueError(
            "The degree l must be greater or equal to 0. Input value was {:s}."
            .format(repr(l))
            )

    if m < 0:
        raise ValueError(
            "The order m must be greater or equal to 0. Input value was {:s}."
            .format(repr(m))
            )

    if m > l:
        raise ValueError(
            "The order m must be less than or equal to the degree l. " +
            "Input values were l={:s} and m={:s}.".format(repr(l), repr(m))
            )

    if normalization.lower() not in ('4pi', 'ortho', 'schmidt', 'unnorm'):
        raise ValueError(
            "The normalization must be '4pi', 'ortho', 'schmidt', " +
            "or 'unnorm'. Input value was {:s}."
            .format(repr(normalization))
            )

    if csphase != 1 and csphase != -1:
        raise ValueError(
            "csphase must be either 1 or -1. Input value was {:s}."
            .format(repr(csphase))
            )

    if cnorm != 0 and cnorm != 1:
        raise ValueError(
            "cnorm must be either 0 or 1. Input value was {:s}."
            .format(repr(cnorm))
            )

    if normalization.lower() == 'unnorm' and lmax > 85:
        _warnings.warn("Calculations using unnormalized coefficients " +
                       "are stable only for degrees less than or equal " +
                       "to 85. lmax for the coefficients will be set to " +
                       "85. Input value was {:d}.".format(lmax),
                       category=RuntimeWarning)
        lmax = 85

    if normalization == '4pi':
        p = _PlmBar(l, z, csphase=csphase, cnorm=cnorm)
    elif normalization == 'ortho':
        p = _PlmON(l, z, csphase=csphase, cnorm=cnorm)
    elif normalization == 'schmidt':
        p = _PlmSchmidt(l, z, csphase=csphase, cnorm=cnorm)
    elif normalization == 'unnorm':
        p = _PLegendreA(l, z, csphase=csphase, cnorm=cnorm)

    return p[(l*(l+1))//2+m]
Compute the associated Legendre function for a specific degree and order.

    Usage
    -----
    plm = legendre_lm (l, m, z, [normalization, csphase, cnorm])

    Returns
    -------
    plm : float
        The associated Legendre functions for degree l and order m.

    Parameters
    ----------
    l : integer
        The spherical harmonic degree.
    m : integer
        The spherical harmonic order.
    z : float
        The argument of the associated Legendre functions.
    normalization : str, optional, default = '4pi'
        '4pi', 'ortho', 'schmidt', or 'unnorm' for use with geodesy 4pi
        normalized, orthonormalized, Schmidt semi-normalized, or unnormalized
        spherical harmonic functions, respectively.
    csphase : optional, integer, default = 1
        If 1 (default), the Condon-Shortley phase will be excluded. If -1, the
        Condon-Shortley phase of (-1)^m will be appended to the associated
        Legendre functions.
    cnorm : optional, integer, default = 0
        If 1, the complex normalization of the associated Legendre functions
        will be used. The default is to use the real normalization.

    Description
    -----------
    legendre_lm will calculate the associated Legendre function for a specific
    degree l and order m. The Legendre functions are used typically as a part
    of the spherical harmonic functions, and three parameters determine how
    they are defined. normalization can be either '4pi' (default), 'ortho',
    'schmidt', or 'unnorm' for use with 4pi normalized, orthonormalized,
    Schmidt semi-normalized, or unnormalized spherical harmonic functions,
    respectively. csphase determines whether to include or exclude (default)
    the Condon-Shortley phase factor. cnorm determines whether to normalize
    the Legendre functions for use with real (default) or complex spherical
    harmonic functions.

    The Legendre functions are calculated using the standard three-term
    recursion formula, and in order to prevent overflows, the scaling approach
    of Holmes and Featherstone (2002) is utilized. The resulting functions are
    accurate to about degree 2800. See Wieczorek and Meschede (2018) for exact
    definitions on how the Legendre functions are defined.

    References
    ----------
    Holmes, S. A., and W. E. Featherstone, A unified approach to the Clenshaw
    summation and the recursive computation of very high degree and order
    normalised associated Legendre functions, J. Geodesy, 76, 279-299,
    doi:10.1007/s00190-002-0216-2, 2002.

    Wieczorek, M. A., and M. Meschede. SHTools — Tools for working with
    spherical harmonics, Geochem., Geophys., Geosyst., 19, 2574-2592,
    doi:10.1029/2018GC007529, 2018.