The angular coordinates on a spherical shell are the co-latitude $\theta$ and the longitude $\phi$ .

The spatial discretization is defined by a call to shtns_init; or to shtns_set_grid or shtns_set_grid_auto after shtns_create.

A field A (or one of its component in case of a vector field) is discretized on an ordered grid, which consists of

NPHI equally spaced nodes in longitude, spanning the range of angle between 0 (included) and $2\pi$ /MRES (excluded),
NLAT gauss nodes or equally spaced nodes (depending on shtns_type) in latitude, spanning angles from 0 to $\pi$ .

There are constraints on the sizes NLAT and NPHI. The sampling theorem requires NPHI > 2*MMAX, with MMAX the highest Fourier mode to be resolved. Similarly NLAT > LMAX is required for a Gauss-Legendre quadrature.

When dealing with non-linear equation, or more generally when analysing data that is not band-limited to LMAX and/or MMAX, anti-aliasing constraints exist. For a given non-linear order of the equation N (typically 2), the anti-aliasing conditions are :

NPHI > (N+1)*MMAX
NLAT > (N+1)*LMAX/2 for Gauss-Legendre quadrature
NLAT > N*LMAX for our DCT quadrature

Furthermore, FFT is more efficient for certain values of NPHI and/or NLAT. If you do not require specific NLAT or NPHI, you can let SHTns choose them for you, by using shtns_set_grid_auto with nphi and/or nlat set to zero.

Identifying a grid point with its indices it and ip in latitudinal and longitudinal direction respectively (indices start at 0), one has :

$\theta$ = acos(ct[it]) where ct = shtns_info::ct is initialized by the call to shtns_init
$\phi$ = ip * 2 $\pi$ /(NPHI*MRES) and PHI_RAD and PHI_DEG can be used to get the phi angle corresponding to ip.
How to access $A(\theta,\phi)$ depends on the data layout, but it is always a contiguous array of double precision floating point values (64 bit)

Currently, three data layouts are supported. Contiguous longitudes, Contiguous latitudes, and Native layout.

Contiguous longitudes

In this layout, increasing longitudes are stored next to each other for each latitude. That is $A(\theta,\phi)$ = A[it*NPHI + ip] in C or A(ip,it) in Fortran. Use SHT_PHI_CONTIGUOUS to instruct shtns_init to use this spatial data layout :

shtns_init ( sht_gauss | SHT_PHI_CONTIGUOUS, ... )

will tell shtns to precompute everything for a gauss grid, and spatial data stored with longitude varying fastest.

Contiguous latitudes

In this layout, increasing latitude are stored next to each other for each longitude. That is $A(\theta,\phi)$ = A[ip*NLAT + it] in C or A(it,ip) in Fortran. Use SHT_THETA_CONTIGUOUS to instruct shtns_init to use this spatial data layout.

Native layout

The native way of storing spatial field data (which will help you achieve the best performance with SHTns) is the same as the Contiguous latitudes layout, except that it requires you to allocate slightly more memory for a field than the NLAT*NPHI double values. Namely NLAT*(NPHI/2+1)*2 doubles are required (instead of NLAT*NPHI) to be able tu use the in-place FFT of FFTW.

In Fortran this means you will allocate data as if the phi direction had (NPHI/2+1)*2 points instead fo NPHI. The additional space, is located at the end of the useful data, and you don't need to worry about it.

To instruct shtns_init that your spatial data has been set up using this layout, use SHT_NATIVE_LAYOUT.

Note: One must be careful when allocating spatial data, as the Fourier transform requires some additional space.
Precisely, the number of required double is NSPAT_ALLOC = NLAT*(NPHI/2+1)*2. The first NLAT*NPHI values are the spatial data.; In addition, spatial data must be allocated using the fftw_malloc function, for example using :
A = fftw_malloc( NSPAT_ALLOC * sizeof(double) );; It is not recommended to use the native layout with Fortran or Python.