How can I interpolate data represented in (r,spherical harmonic) space to a regular cartesian grid (F90)?

Question

I'm trying to interpolate spherical harmonics to a cubic, Cartesian grid.

The output data of my spherical, pseudo-spectral simulation has Nr radial levels between rMin and rMax, each containing a set of finite-order spherical harmonics for longitude and latitude. The spherical harmonics are mapped to a physical spherical grid containing Ni latitudes and Nj longitudes via a triangular truncation.

The domain is as follows:

• Radial levels: rMin <= r(k) <= rMax, with indexing 1 <= k <= Nr
• Spherical harmonics (triangular truncation, without aliasing from transform):
  • Nm = (Nj-1)/3
  • 0 <= m <= Nm
  • m <= l <= Nm
  • nlm == (nm+1)*(nm+2)/2 (the total number of l,m combinations)

Data arrays:

• Spectral form: complex*16, dimension( 1:nlm, 1:Nr ) :: foo_spectral
• Cartesian form: real*8, dimension( 1:Nx, 1:Ny, 1:Nz ) :: foo_cartesian

I'm looking for an accurate and efficient way to interpolate the data from its spectral representation to a cubic Cartesian grid with edge-length 2*rMax, such that the spherical domain fits perfectly inside. I only want to interpolate within the sphere, however: for points corresponding to r<rMin or rMax<r, the cubic grid should have OUTSIDE_DOMAIN values.

Currently, I have to transform the data from its spectral representation (spherical harmonics: foo(Nr,nlm)) to a physical representation (spherical grid: foo(Nr,Ni,Nj)), and then use a QHULL routine in IDL to interpolate from the physical, spherical grid to the physical, cubic grid (foo(Nx,Ny,Nz)) (note that Nx==Ny==Nz for a cubic grid).

The size of my data is larger than my existing code (written in IDL) can handle, and converting to spherical space is unnecessary for my purposes. I'd like a more direct method that is stand-alone -- not dependent on IDL, for instance.

Any thoughts about how this could be done? I'm willing to use open-source libraries, but it would be nice to not have to.

Thanks in advance!

Answer 1

I would strongly recommend using libraries for this; the spherical harmonic transform is hard to do efficiently and accurately and it's unlikely your first attempts will be anything like as good as existing routines.

One library a colleague thinks quite highly of is SHTns, which will both do the synthesis (inverse transform) for you and do the interpolation (for any given shell) at an arbitrary point. It has fortran bindings. You'd still have to handle the multiple radial shells yourself, one way or another (probably by doing what you're doing now - transform everything onto a spherical grid, and then use standard interpolation methods to get onto a cubic grid), and while that is a little tricky to do right it's much more straightforward than the spherical harmonic transform part.