<SYSTEM_TASK:>
Solve the following problem using Python, implementing the functions described below, one line at a time
<END_TASK>
<USER_TASK:>
Description:
 def spharm_lm(l, m, theta, phi, normalization='4pi', kind='real', csphase=1, degrees=True):
 """ Compute the spherical harmonic function for a specific degree and order. Usage ----- ylm = spharm (l, m, theta, phi, [normalization, kind, csphase, degrees]) Returns ------- ylm : float or complex The spherical harmonic function ylm, where l and m are the spherical harmonic degree and order, respectively. Parameters l : integer The spherical harmonic degree. m : integer The spherical harmonic order. theta : float The colatitude in degrees. phi : float The longitude in degrees. normalization : str, optional, default = '4pi' '4pi', 'ortho', 'schmidt', or 'unnorm' for geodesy 4pi normalized, orthonormalized, Schmidt semi-normalized, or unnormalized spherical harmonic functions, respectively. kind : str, optional, default = 'real' 'real' or 'complex' spherical harmonic coefficients. 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 spherical harmonic functions. degrees : optional, bool, default = True If True, colat and phi are expressed in degrees. Description spharm_lm will calculate the spherical harmonic function for a specific degree l and order m, and for a given colatitude theta and longitude phi. Three parameters determine how the spherical harmonic functions are defined. normalization can be either '4pi' (default), 'ortho', 'schmidt', or 'unnorm' for 4pi normalized, orthonormalized, Schmidt semi-normalized, or unnormalized spherical harmonic functions, respectively. kind can be either 'real' or 'complex', and csphase determines whether to include or exclude (default) the Condon-Shortley phase factor. The spherical harmonic 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 spherical harmonic 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 than 0. Input value was {:s}."
            .format(repr(l))
            )

    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 kind.lower() not in ('real', 'complex'):
        raise ValueError(
            "kind must be 'real' or 'complex'. " +
            "Input value was {:s}.".format(repr(kind))
            )

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

    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

    ind = (l*(l+1))//2 + abs(m)

    if degrees is True:
        theta = _np.deg2rad(theta)
        phi = _np.deg2rad(phi)

    if kind.lower() == 'real':
        p = _legendre(l, _np.cos(theta), normalization=normalization,
                      csphase=csphase, cnorm=0, packed=True)
        if m >= 0:
            ylm = p[ind] * _np.cos(m*phi)
        else:
            ylm = p[ind] * _np.sin(abs(m)*phi)

    else:
        p = _legendre(l, _np.cos(theta), normalization=normalization,
                      csphase=csphase, cnorm=1, packed=True)
        ylm = p[ind] * (_np.cos(m*phi) + 1j * _np.sin(abs(m)*phi))  # Yl|m|

        if m < 0:
            ylm = ylm.conj()
            if _np.mod(m, 2) == 1:
                ylm = - ylm

    return ylm