Convert values of spin-weighted function on a grid to mode weights

    Parameters
    ----------
    map : array_like, complex, shape (..., Ntheta, Nphi)
        Values of the spin-weighted function on grid points of the sphere.  This array may have more
        than two dimensions, where initial dimensions may represent different times, for example, or
        separate functions on the sphere.  The final two dimensions should give the values of the
        function, in the order described below in the 'Notes' section.
    s : int or array, int, shape (...)
        Spin weight of the function.  If `amp` is multidimensional and this is an array, its
        dimensions must match the first dimensions of `map`, and the different values are the spin
        weights of the different functions represented by those dimensions.  Otherwise, if `map` is
        multidimensional and `s` is a single integer, all functions are assumed to have the same
        spin weight.
    lmax : int
        The largest `ell` value present in the input array.


    Returns
    -------
    salm : ndarray, complex, shape (..., (lmax+1)**2)
        Mode weights of the spin-weighted function.  This array is shaped like the input `map` array,
        but has one less dimension.  The final dimension describes the values of the mode weights on
        the corresponding sphere, as described below in the 'Notes' section.


    See also
    --------
    spinsfast.map2salm : Roughly the inverse of this function.


    Notes
    -----

    The input data represent the values on this grid of spherical coordinates:

        np.array([[map(theta, phi)
                   for phi in np.linspace(0.0, 2*np.pi, num=2*lmax+1, endpoint=False)]
                  for theta in np.linspace(0.0, np.pi, num=2*lmax+1, endpoint=True)])

    The input is converted to a contiguous complex numpy array if necessary.

    The output `salm` data are given in increasing order of `ell` value, always starting with
    (ell, m) = (0, 0) even if `s` is nonzero, proceeding to (1, -1), (1, 0), (1, 1), etc.
    Explicitly, the ordering matches this:

        [map_lm(ell, m) for ell in range(lmax+1) for m in range(-ell, ell+1)]

    Note that `map2salm` and `salm2map` are not true inverses of each other for several reasons.
    First, modes with `ell < |s|` should always be zero; they are simply assumed to be zero on input
    to `salm2map`.  It is possible to define a `map` function that violates this assumption -- for
    example, having a nonzero average value over the sphere, if the function has nonzero spin `s`,
    this is impossible.  Also, it is possible to define a map of a function with so much angular
    dependence that it cannot be captured with the given `lmax` value.  For example, a discontinuous
    function will never be perfectly resolved.


    Example
    -------
    >>> s = -2
    >>> lmax = 8
    >>> theta_phi = np.array([[[theta, phi]
                               for phi in np.linspace(0.0, 2*np.pi, num=2*lmax+1, endpoint=False)]
                              for theta in np.linspace(0.0, np.pi, num=2*lmax+1, endpoint=True)])
    >>> map = np.array([[np.sqrt(3/(8*np.pi)) * np.sin(tp[0]) for tp in _] for _ in theta_phi])
    >>> salm = spinsfast.map2salm(map, s, lmax)

def map2salm(map, s, lmax):
    """Convert values of spin-weighted function on a grid to mode weights

    Parameters
    ----------
    map : array_like, complex, shape (..., Ntheta, Nphi)
        Values of the spin-weighted function on grid points of the sphere.  This array may have more
        than two dimensions, where initial dimensions may represent different times, for example, or
        separate functions on the sphere.  The final two dimensions should give the values of the
        function, in the order described below in the 'Notes' section.
    s : int or array, int, shape (...)
        Spin weight of the function.  If `amp` is multidimensional and this is an array, its
        dimensions must match the first dimensions of `map`, and the different values are the spin
        weights of the different functions represented by those dimensions.  Otherwise, if `map` is
        multidimensional and `s` is a single integer, all functions are assumed to have the same
        spin weight.
    lmax : int
        The largest `ell` value present in the input array.


    Returns
    -------
    salm : ndarray, complex, shape (..., (lmax+1)**2)
        Mode weights of the spin-weighted function.  This array is shaped like the input `map` array,
        but has one less dimension.  The final dimension describes the values of the mode weights on
        the corresponding sphere, as described below in the 'Notes' section.


    See also
    --------
    spinsfast.map2salm : Roughly the inverse of this function.


    Notes
    -----

    The input data represent the values on this grid of spherical coordinates:

        np.array([[map(theta, phi)
                   for phi in np.linspace(0.0, 2*np.pi, num=2*lmax+1, endpoint=False)]
                  for theta in np.linspace(0.0, np.pi, num=2*lmax+1, endpoint=True)])

    The input is converted to a contiguous complex numpy array if necessary.

    The output `salm` data are given in increasing order of `ell` value, always starting with
    (ell, m) = (0, 0) even if `s` is nonzero, proceeding to (1, -1), (1, 0), (1, 1), etc.
    Explicitly, the ordering matches this:

        [map_lm(ell, m) for ell in range(lmax+1) for m in range(-ell, ell+1)]

    Note that `map2salm` and `salm2map` are not true inverses of each other for several reasons.
    First, modes with `ell < |s|` should always be zero; they are simply assumed to be zero on input
    to `salm2map`.  It is possible to define a `map` function that violates this assumption -- for
    example, having a nonzero average value over the sphere, if the function has nonzero spin `s`,
    this is impossible.  Also, it is possible to define a map of a function with so much angular
    dependence that it cannot be captured with the given `lmax` value.  For example, a discontinuous
    function will never be perfectly resolved.


    Example
    -------
    >>> s = -2
    >>> lmax = 8
    >>> theta_phi = np.array([[[theta, phi]
                               for phi in np.linspace(0.0, 2*np.pi, num=2*lmax+1, endpoint=False)]
                              for theta in np.linspace(0.0, np.pi, num=2*lmax+1, endpoint=True)])
    >>> map = np.array([[np.sqrt(3/(8*np.pi)) * np.sin(tp[0]) for tp in _] for _ in theta_phi])
    >>> salm = spinsfast.map2salm(map, s, lmax)

    """
    import numpy as np
    map = np.ascontiguousarray(map, dtype=np.complex128)
    salm = np.empty(map.shape[:-2]+(N_lm(lmax),), dtype=np.complex128)
    if map.ndim>2:
        s = np.ascontiguousarray(s, dtype=np.intc)
        if s.ndim != map.ndim-2 or np.product(s.shape) != np.product(map.shape[:-2]):
            s = s*np.ones(map.shape[:-2], dtype=np.intc)
        _multi_map2salm(map, salm, s, lmax)
    else:
        _map2salm(map, salm, s, lmax)
    return salm