如何在Healpy中获取地图的θ、φ二阶导数？

Great question! Let's walk through how to get those second-order derivatives (including mixed θφ, θθ, and φφ) in Healpy. First, a quick note: Healpy doesn't have a direct built-in function that matches the der2 optional output from Fortran Healpix's alm2map_der—but we have two reliable workarounds to achieve this.

Option 1: Use First-Order Derivatives + Numerical Differentiation #

Healpy's hp.alm2map_der1 function gives you the first-order θ and φ derivatives of a map. We can build on these to compute second-order derivatives numerically:

• θθ second derivative: Take the θ-direction numerical derivative of the first-order θ derivative map
• φφ second derivative: Take the φ-direction numerical derivative of the first-order φ derivative map
• θφ mixed derivative: Take the φ-direction derivative of the first-order θ derivative map (or vice versa—they should be theoretically equivalent, with minor numerical differences)

Here's a concrete code example:

import healpy as hp
import numpy as np

# Set up a test scenario
nside = 64
lmax = 3 * nside - 1
alm = hp.sphtfunc.synalm(np.ones(lmax + 1), lmax=lmax)  # Random test alm

# Get first-order derivatives from Healpy
original_map, dtheta_map, dphi_map = hp.alm2map_der1(alm, nside=nside, lmax=lmax)

# Calculate θθ second derivative
theta, _ = hp.pix2ang(nside, np.arange(hp.nside2npix(nside)))
avg_dtheta = np.mean(np.diff(theta))  # Approximate θ step (HEALPix pixels aren't perfectly uniform)
dtheta2_map = np.gradient(dtheta_map, avg_dtheta, axis=0)

# Calculate φφ second derivative
pixels_per_ring = hp.nside2npix(nside) // (2 * nside)
avg_dphi = 2 * np.pi / pixels_per_ring  # φ step per ring
dphi2_map = np.gradient(dphi_map, avg_dphi, axis=1)

# Calculate θφ mixed derivative
dthetadphi_map = np.gradient(dtheta_map, avg_dphi, axis=1)

Keep in mind: Since HEALPix pixels don't form a regular grid, this is an approximation—accuracy improves with higher nside.

Option 2: Compute Directly from Spherical Harmonic Coefficients #

For more precise results, we can use the analytic properties of spherical harmonics to derive second-order derivative coefficients from the original alm, then map those to a healpy map.

• φφ second derivative: The second φ-derivative of a spherical harmonic (Y_{lm}) is simply (-m^2 Y_{lm}). So we can modify the alm directly:

  # Compute alm for φφ derivative
  alm_phi2 = np.copy(alm)
  l, m = hp.sphtfunc.Alm.getlm(lmax)
  alm_phi2 *= -m ** 2
  # Map to healpy map
  dphi2_map = hp.alm2map(alm_phi2, nside=nside, lmax=lmax)
  

• θθ and θφ derivatives: These require recursive relations for the θ-derivative of spherical harmonics. You can reference the Fortran Healpix source code for alm2map_der to translate the der2 calculation logic into Python—this will give you exact matches to the Fortran output.

Bonus: Wrap the Fortran Function Directly #

If you need perfect parity with the Fortran alm2map_der der2 output, you can use f2py to wrap the original Fortran subroutine into a Python-callable module. This skips any approximation and uses the exact code you're familiar with.

内容的提问来源于stack exchange，提问作者johnhenry