ugly as all hell, but adjointness seems to be fine

interpol_ng/demo.py

@@ -4,7 +4,7 @@ import pysharp
 import time
 import matplotlib.pyplot as plt
 
-np.random.seed(20909)
+np.random.seed(2099)
 
 def nalm(lmax, mmax):
     return ((mmax+1)*(mmax+2))//2 + (mmax+1)*(lmax-mmax)
@@ -24,24 +24,26 @@ def compress_alm(alm,lmax):
     res[lmax+2::2] = np.sqrt(2)*alm[lmax+1:].imag
     return res
 
-
 def myalmdot(a1,a2,lmax,mmax,spin):
     return np.vdot(compress_alm(a1,lmax),compress_alm(a2,lmax))
 def mydot(a1,a2,spin):
     return np.vdot(theta_extend(a1,spin),theta_extend(a2,spin))
 
-lmax=512
-mmax=lmax
-kmax=12
+lmax=30
+kmax=13
 
 
 # get random sky a_lm
 # the a_lm arrays follow the same conventions as those in healpy
-slmT = random_alm(lmax, mmax)
+slmT = random_alm(lmax, lmax)
 
 # build beam a_lm
 blmT = random_alm(lmax, kmax)
-ptg = np.array([[1.2,0.01,1.],[1.4,0.7,2.]])
+ptg=np.random.uniform(0.,1.,3*1000000).reshape(1000000,3)
+ptg[:,0]*=np.pi
+ptg[:,1]*=2*np.pi
+ptg[:,2]*=2*np.pi
+#ptg = np.array([[0.129,0.01,1.],[3.1,0.7,2.]])
 foo = interpol_ng.PyInterpolator(slmT,blmT,lmax, kmax, epsilon=1e-6, nthreads=1)
 bar=foo.interpol(ptg)
 print(foo.Nphi(),foo.Nphi0())

interpol_ng/interpol_ng.cc

@@ -336,7 +336,7 @@ arr.apply([](T &v){v=0.;});
       auto ainfo = sharp_make_triangular_alm_info(lmax,lmax,1);
 
       // move stuff from border regions onto the main grid
-      for (size_t i=0; i<ntheta+2*supp; ++i)
+      for (size_t i=0; i<cube.shape(0); ++i)
         for (size_t j=0; j<supp; ++j)
           for (size_t k=0; k<cube.shape(2); ++k)
             {
@@ -352,7 +352,20 @@ arr.apply([](T &v){v=0.;});
             cube.v(supp+1+i,j+supp,k) += fct*cube(supp-1-i,j2+supp,k);
             cube.v(supp+ntheta-2-i, j+supp,k) += fct*cube(supp+ntheta+i,j2+supp,k);
             }
-
+for (size_t k=0; k<cube.shape(2); ++k)
+{
+double fct = (((k+1)/2)&1) ? -1 : 1;
+for (size_t j=0,j2=nphi/2; j<nphi/2; ++j,++j2)
+  {
+  if (j2>=nphi) j2-=nphi;
+  double tval = (cube(supp,j+supp,k) + fct*cube(supp,j2+supp,k));
+  cube.v(supp,j+supp,k) = tval;
+  cube.v(supp,j2+supp,k) = fct*tval;
+  tval = (cube(supp+ntheta-1,j+supp,k) + fct*cube(supp+ntheta-1,j2+supp,k));
+  cube.v(supp+ntheta-1,j+supp,k) = tval;
+  cube.v(supp+ntheta-1,j2+supp,k) = fct*tval;
+  }
+}
       vector<double>lnorm(lmax+1);
       for (size_t i=0; i<=lmax; ++i)
         lnorm[i]=sqrt(4*pi/(2*i+1.));
@@ -362,8 +375,8 @@ arr.apply([](T &v){v=0.;});
       decorrect(m1,0);
    for (size_t j=0; j<nphi0; ++j)
      {
-     m1.v(0,j)*=0.5;;
-     m1.v(ntheta0-1,j)*=0.5;;
+     m1.v(0,j)*=0.5;
+     m1.v(ntheta0-1,j)*=0.5;
      }
       sharp_alm2map_adjoint(a1.Alms().vdata(), m1.data(), *ginfo, *ainfo, 0, nthreads);
       for (size_t m=0; m<=lmax; ++m)
@@ -379,13 +392,13 @@ arr.apply([](T &v){v=0.;});
         decorrect(m2,k);
    for (size_t j=0; j<nphi0; ++j)
      {
-     m1.v(0,j)*=0.5;;
-     m1.v(ntheta0-1,j)*=0.5;;
+     m1.v(0,j)*=0.5;
+     m1.v(ntheta0-1,j)*=0.5;
      }
    for (size_t j=0; j<nphi0; ++j)
      {
-     m2.v(0,j)*=0.5;;
-     m2.v(ntheta0-1,j)*=0.5;;
+     m2.v(0,j)*=0.5;
+     m2.v(ntheta0-1,j)*=0.5;
      }
 
         sharp_alm2map_spin_adjoint(k, a1.Alms().vdata(), a2.Alms().vdata(), m1.data(),