From da748f619a864f0889ad100230696b8912dc4f7c Mon Sep 17 00:00:00 2001 From: Dag Sverre Seljebotn Date: Mon, 22 May 2017 11:02:43 +0200 Subject: [PATCH] Change API to be able to support spin-weighted Legendre functions in the future --- libsharp/sharp_legendre_table.c | 25 +++++++++++-------------- libsharp/sharp_legendre_table.h | 22 ++++++++++++++++------ python/libsharp/libsharp.pxd | 4 ++-- python/libsharp/libsharp.pyx | 2 +- 4 files changed, 30 insertions(+), 23 deletions(-) diff --git a/libsharp/sharp_legendre_table.c b/libsharp/sharp_legendre_table.c index 5420899..7fef085 100644 --- a/libsharp/sharp_legendre_table.c +++ b/libsharp/sharp_legendre_table.c @@ -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); diff --git a/libsharp/sharp_legendre_table.h b/libsharp/sharp_legendre_table.h index ccabf52..de4cdd6 100644 --- a/libsharp/sharp_legendre_table.h +++ b/libsharp/sharp_legendre_table.h @@ -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 ); diff --git a/python/libsharp/libsharp.pxd b/python/libsharp/libsharp.pxd index bdc296b..fb20f5c 100644 --- a/python/libsharp/libsharp.pxd +++ b/python/libsharp/libsharp.pxd @@ -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": diff --git a/python/libsharp/libsharp.pyx b/python/libsharp/libsharp.pyx index d4bbeb2..e1c39bb 100644 --- a/python/libsharp/libsharp.pyx +++ b/python/libsharp/libsharp.pyx @@ -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