如何将（r，球谐）空间中表示的数据插值到规则的笛卡尔网格？

Question

我正在尝试将球谐函数插值到三次笛卡尔网格。

我的球形伪光谱模拟的输出数据在

Nr

和

rMin

之间具有

rMax

径向级别，每个级别都包含一组经度和纬度的有限阶球谐函数。球谐函数通过三角截断映射到包含

Ni

纬度和

Nj

经度的物理球面网格。

域名如下：

• 径向水平仪：

  rMin <= r(k) <= rMax

  ，带索引

  1 <= k <= Nr

• 球谐函数（三角截断，没有变换造成的混叠）：

  • Nm = (Nj-1)/3

  • 0 <= m <= Nm

  • m <= l <= Nm

  • nlm == (nm+1)*(nm+2)/2

    （

    l

    、

    m

    组合总数）

数据数组：

• 光谱形式：

  complex*16, dimension( 1:nlm, 1:Nr ) :: foo_spectral

• 笛卡尔形式：

  real*8, dimension( 1:Nx, 1:Ny, 1:Nz ) :: foo_cartesian

我正在寻找一种准确而有效的方法，将数据从其光谱表示形式插值到边长为

2*rMax

的立方笛卡尔网格中，以便球形域完全适合内部。然而，我只想在球体内插值 ：对于对应于 r<rMin

 或 

rMax<r

 的点，立方网格应该具有 

OUTSIDE_DOMAIN

 值。 

目前，我必须将数据从其光谱表示（球面谐波：

foo(Nr,nlm)

）转换为物理表示（球面网格：

foo(Nr,Ni,Nj)

），然后使用 IDL 中的 

QHULL 例程从物理进行插值，球形网格到物理立方网格 (foo(Nx,Ny,Nz)

)（注意 

Nx==Ny==Nz

 代表立方网格）。 

我的数据大小大于我现有的代码（用 IDL 编写）可以处理的大小，并且对于我的目的而言，无需转换为球形空间。我想要一种更直接的独立方法，例如不依赖于 IDL。

对于如何做到这一点有什么想法吗？我愿意使用开源库，但如果不必这样做就好了。

提前致谢！

Answer 1

我强烈建议为此使用库；球谐变换很难高效且准确地完成，并且您的第一次尝试不太可能像现有例程一样好。

同事评价很高的一个库是

SHTns，它既可以为您进行综合（逆变换），又可以在任意点进行插值（对于任何给定的 shell）。它有 Fortran 绑定。您仍然必须自己处理多个径向壳，以一种或另一种方式（可能通过做您现在正在做的事情 - 将所有内容转换为球形网格，然后使用标准插值方法进入立方网格），并且虽然这有点棘手，但它比球谐变换部分要简单得多。