Aquila-Consortium/libsharp2
0
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

Change API to be able to support spin-weighted Legendre functions in the future

Browse source
This commit is contained in:
Dag Sverre Seljebotn 2017-05-22 11:02:43 +02:00
parent a93db0b1aa
commit da748f619a
4 changed files with 30 additions and 23 deletions
Download patch file Download diff file Expand all files Collapse all files

25
libsharp/sharp_legendre_table.c
View file

@ -1417,19 +1417,16 @@ New high-level wrapper
void sharp_normalized_associated_legendre_table(
ptrdiff_t m,
int spin,
ptrdiff_t lmax,
ptrdiff_t ntheta,
/* contigious 1D array of theta values to compute for,
contains ntheta values */
double *theta,
/* number of columns in out; see below. Should be >= (lmax - m). */
ptrdiff_t ncols,
/* contigiuos 2D array, in "theta-major ordering". Has `ntheta`
rows and `ncols` columns. Indexed as out[itheta * ncols + (l - m)].
If `ncols > lmax - m` then those entries are not accessed.
*/
ptrdiff_t theta_stride,
ptrdiff_t l_stride,
ptrdiff_t spin_stride,
double *out
) {
if (spin != 0) UTIL_FAIL ("sharp_normalized_associated_legendre_table: only spin=0 has been implemented so far");
Ylmgen_C ctx;
ptrdiff_t itheta, l, lmin;
@ -1444,14 +1441,14 @@ void sharp_normalized_associated_legendre_table(
Ylmgen_recalc_Ylm_sse2(&ctx);
lmin = IMIN(*ctx.firstl, lmax + 1);
for (l = m; l < lmin; ++l) {
out[itheta * ncols + l - m] = 0;
out[(itheta + 1) * ncols + l - m ] = 0;
out[itheta * theta_stride + (l - m) * l_stride + spin * spin_stride] = 0;
out[(itheta + 1) * theta_stride + (l - m) * l_stride + spin * spin_stride] = 0;
}
for (l = IMAX(lmin, m); l <= lmax; ++l) {
double v1, v2;
read_v2df(ctx.ylm_sse2[l], &v1, &v2);
out[itheta * ncols + l - m] = v1;
out[(itheta + 1) * ncols + l - m] = v2;
out[itheta * theta_stride + (l - m) * l_stride + spin * spin_stride] = v1;
out[(itheta + 1) * theta_stride + (l - m) * l_stride + spin * spin_stride] = v2;
}
}
#endif
@ -1460,10 +1457,10 @@ void sharp_normalized_associated_legendre_table(
Ylmgen_recalc_Ylm(&ctx);
lmin = IMIN(*ctx.firstl, lmax + 1);
for (l = m; l < lmin; ++l) {
out[itheta * ncols + l - m] = 0;
out[itheta * theta_stride + (l - m) * l_stride + spin * spin_stride] = 0;
}
for (l = IMAX(lmin, m); l <= lmax; ++l) {
out[itheta * ncols + l - m] = ctx.ylm[l];
out[itheta * theta_stride + (l - m) * l_stride + spin * spin_stride] = ctx.ylm[l];
}
}
Ylmgen_destroy(&ctx);

22
libsharp/sharp_legendre_table.h
View file

@ -54,27 +54,37 @@ extern "C" {
(Internally, sin(theta) is also used for part of the computation, making theta
the most convenient argument.)
\param m The m-value to compute a table for
NOTE: Support for spin-weighted Legendre functions is on the TODO-list. Only spin=0
is supported now.
\param m The m-value to compute a table for; must be >= 0
\param spin The spin parameter; pass 0 for the regular associated Legendre functions.
NOTE: This is present for future compatability, currently only 0 is supported.
\param lmax A table will be provided for l = m .. lmax
\param ntheta How many theta values to evaluate for
\param theta Contiguous 1D array of theta values
\param ncols Number of columns in the out-array; should have ncols >= (lmax - m)
\param out Contiguous 2D array that will receive the output. Each output entry
is assigned to out[itheta * ncols + (l - m)].
\param theta_stride See below
\param l_stride See below
\param spin_stride See below. "ispin" will always be 0 if spin==0, or 0 for positive spin
and 1 for the corresponding negative spin otherwise.
\param out Contiguous 3D array that will receive the output. Each output entry
is assigned to out[itheta * theta_stride + (l - m) * l_stride + ispin * spin_stride].
*/
void sharp_normalized_associated_legendre_table(
ptrdiff_t m,
int spin,
ptrdiff_t lmax,
ptrdiff_t ntheta,
/* contiguous 1D array of theta values to compute for,
contains ntheta values */
double *theta,
/* number of columns in out; see below. Should be >= (lmax - m). */
ptrdiff_t ncols,
/* contiguous 2D array, in "theta-major ordering". Has `ntheta`
rows and `ncols` columns. Indexed as out[itheta * ncols + (l - m)].
If `ncols > lmax - m` then those entries are not accessed.
*/
ptrdiff_t theta_stride,
ptrdiff_t l_stride,
ptrdiff_t spin_stride,
double *out
);

4
python/libsharp/libsharp.pxd
View file

@ -59,8 +59,8 @@ cdef extern from "sharp.h":
sharp_alm_info *alm_info, int ntrans, int flags, double *time,
unsigned long long *opcnt) nogil
void sharp_normalized_associated_legendre_table(int m, int lmax, int ntheta,
double *theta, int ncols, double *out) nogil
void sharp_normalized_associated_legendre_table(int m, int spin, int lmax, int ntheta,
double *theta, int theta_stride, int l_stride, int spin_stride, double *out) nogil
cdef extern from "sharp_geomhelpers.h":

2
python/libsharp/libsharp.pyx
View file

@ -267,5 +267,5 @@ def normalized_associated_legendre_table(int lmax, int m, theta):
if lmax < m:
raise ValueError("lmax < m")
with nogil:
sharp_normalized_associated_legendre_table(m, lmax, theta_.shape[0], &theta_[0], lmax - m + 1, &out_[0,0])
sharp_normalized_associated_legendre_table(m, 0, lmax, theta_.shape[0], &theta_[0], lmax - m + 1, 1, 1, &out_[0,0])
return out