separability_matrix reports pseudocylindrical, Bonne and HEALPix projections as separable #20578

Description

separability_matrix (and is_separable) say that the two outputs of the pseudocylindrical projections, the Bonne projection and HEALPix each depend on only one input. They do not: on these projections the meridians are curved, so the projected x depends on latitude. Any code that uses the separability matrix, for example a gwcs axis-correlation matrix built from a Pix2Sky_* model, therefore decides incorrectly that the pixel x axis is independent of latitude.

Expected behavior

For Sanson-Flamsteed, x = φ cos θ. Analytically, ∂x/∂θ = -φ sin θ · π/180, which is -0.1309 at (φ, θ) = (15°, 30°). So output 0 depends on both inputs, and the matrix should be [[True, True], [False, True]] (is_separable → [False, True]). The model's own values confirm this: x is 15, 12.99 and 7.5 for the same longitude at latitudes 0, 30 and 60.

How to Reproduce

from astropy.modeling import models, separability_matrix, is_separable

m = models.Sky2Pix_SansonFlamsteed()   # x = lon * cos(lat), y = lat
print(separability_matrix(m))
print(is_separable(m))

# x changes when only the latitude changes
print(m(15, 0)[0], m(15, 30)[0], m(15, 60)[0])

[[ True False]
 [False  True]]
[ True  True]
15.0 12.99038105676658 7.500000000000002

Versions

import astropy
astropy.system_info()

# Paste the result here