|
4 | 4 | from typing import TypeAlias, TypeVar, cast |
5 | 5 |
|
6 | 6 | import jax |
7 | | -from jax import lax |
8 | | -from jax.scipy.special import sph_harm |
9 | 7 | from jaxtyping import ArrayLike, Shaped |
10 | 8 |
|
11 | 9 | import quaxed.numpy as jnp |
12 | 10 |
|
13 | 11 | import galax._custom_types as gt |
| 12 | +from galax.potential._src.builtin.multipole import compute_Ylm |
14 | 13 |
|
15 | 14 | BatchableIntSz0: TypeAlias = Shaped[gt.IntSz0, "*#batch"] |
16 | 15 |
|
@@ -53,16 +52,6 @@ def psi_of_r(r: T) -> T: |
53 | 52 | # ============================================================================= |
54 | 53 |
|
55 | 54 |
|
56 | | -@partial(jax.numpy.vectorize, signature="(n),(),()->(n)", excluded=(3,)) |
57 | | -@partial(jax.jit, static_argnames=("m_max",)) # TODO: should l,m be static? |
58 | | -def _real_Ylm(theta: gt.SzN, l: gt.IntSz0, m: gt.IntSz0, m_max: int) -> gt.SzN: |
59 | | - # TODO: sph_harm only supports scalars, even though it returns an array! |
60 | | - theta = jnp.atleast_1d(theta) |
61 | | - return sph_harm( |
62 | | - m, jnp.atleast_1d(l), theta=jnp.zeros_like(theta), phi=theta, n_max=m_max |
63 | | - ).real |
64 | | - |
65 | | - |
66 | 55 | @partial(jax.jit, static_argnames=("m_max",)) |
67 | 56 | def real_Ylm( |
68 | 57 | theta: gt.SzAny, l: BatchableIntSz0, m: BatchableIntSz0, m_max: int = 100 |
@@ -90,4 +79,5 @@ def real_Ylm( |
90 | 79 | ``(n),(),()->(n)``. |
91 | 80 | """ |
92 | 81 | # TODO: raise an error if m > m_max |
93 | | - return _real_Ylm(theta, l, m, m_max) |
| 82 | + _real_Ylm, _ = compute_Ylm(l, m, theta=theta, phi=jnp.zeros_like(theta), l_max=m_max) |
| 83 | + return _real_Ylm |
0 commit comments