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

Author: Dag Sverre Seljebotn

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);

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
 );
 

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":

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