.

include/sctl/fft_wrapper.hpp

@@ -364,8 +364,8 @@ template <class ValueType, class FFT_Derived> class FFT_Generic {
   }
 
   static void check_align(const Vector<ValueType>& in, const Vector<ValueType>& out) {
-    SCTL_ASSERT_MSG((((uintptr_t)& in[0]) & ((uintptr_t)(SCTL_MEM_ALIGN - 1))) == 0, "sctl::FFT: Input vector not aligned to " <<SCTL_MEM_ALIGN<<" bytes!");
-    SCTL_ASSERT_MSG((((uintptr_t)&out[0]) & ((uintptr_t)(SCTL_MEM_ALIGN - 1))) == 0, "sctl::FFT: Output vector not aligned to "<<SCTL_MEM_ALIGN<<" bytes!");
+    //SCTL_ASSERT_MSG((((uintptr_t)& in[0]) & ((uintptr_t)(SCTL_MEM_ALIGN - 1))) == 0, "sctl::FFT: Input vector not aligned to " <<SCTL_MEM_ALIGN<<" bytes!");
+    //SCTL_ASSERT_MSG((((uintptr_t)&out[0]) & ((uintptr_t)(SCTL_MEM_ALIGN - 1))) == 0, "sctl::FFT: Output vector not aligned to "<<SCTL_MEM_ALIGN<<" bytes!");
     // TODO: copy to auxiliary array if unaligned
   }
 

include/sctl/sph_harm.hpp

@@ -479,6 +479,7 @@ template <class Real> class SphericalHarmonics{
     };
     static MatrixStorage& MatrixStore(){
       static MatrixStorage storage;
+      if (!storage.Qx_.size()) storage.Resize(SCTL_SHMAXDEG);
       return storage;
     }
 };

include/sctl/sph_harm.txx

@@ -662,13 +662,13 @@ template <class Real> void SphericalHarmonics<Real>::StokesEvalSL(const Vector<R
   for (Long i = 0; i < N; i++) { // Set StokesOp
 
     Real cos_theta, csc_theta, cos_phi, sin_phi;
-    { // Set exp_phi, cos_theta, csc_theta, cos_phi, sin_phi
+    { // Set cos_theta, csc_theta, cos_phi, sin_phi
       cos_theta = cos_theta_phi[i * 2 + 0];
       csc_theta = 1 / sqrt<Real>(1 - cos_theta * cos_theta);
       cos_phi = cos(cos_theta_phi[i * 2 + 1]);
       sin_phi = sin(cos_theta_phi[i * 2 + 1]);
     }
-    Complex<Real> imag(0,1), exp_phi(cos_phi, sin_phi);
+    Complex<Real> imag(0,1), exp_phi(cos_phi, -sin_phi);
 
     for (Long m = 0; m <= p0; m++) {
       for (Long n = m; n <= p0; n++) {
@@ -706,7 +706,7 @@ template <class Real> void SphericalHarmonics<Real>::StokesEvalSL(const Vector<R
           Complex<Real> Ycsc_2 = Y_2 * csc_theta;
           if (fabs(cos_theta) == 1) {
             auto Y_csc0 = [exp_phi, cos_theta](Long n, Long m) {
-              if (m == 1) -sqrt<Real>((2*n+1)*n*(n+1)) * ((n%2==0) && (cos_theta<0) ? -1 : 1) * exp_phi;
+              if (m == 1) return -sqrt<Real>((2*n+1)*n*(n+1)) * ((n%2==0) && (cos_theta<0) ? -1 : 1) * exp_phi;
               return Complex<Real>(0, 0);
             };
             Ycsc_0 = Y_csc0(n - 1, m);
@@ -866,13 +866,13 @@ template <class Real> void SphericalHarmonics<Real>::StokesEvalDL(const Vector<R
   for (Long i = 0; i < N; i++) { // Set StokesOp
 
     Real cos_theta, csc_theta, cos_phi, sin_phi;
-    { // Set exp_phi, cos_theta, csc_theta, cos_phi, sin_phi
+    { // Set cos_theta, csc_theta, cos_phi, sin_phi
       cos_theta = cos_theta_phi[i * 2 + 0];
       csc_theta = 1 / sqrt<Real>(1 - cos_theta * cos_theta);
       cos_phi = cos(cos_theta_phi[i * 2 + 1]);
       sin_phi = sin(cos_theta_phi[i * 2 + 1]);
     }
-    Complex<Real> imag(0,1), exp_phi(cos_phi, sin_phi);
+    Complex<Real> imag(0,1), exp_phi(cos_phi, -sin_phi);
 
     for (Long m = 0; m <= p0; m++) {
       for (Long n = m; n <= p0; n++) {
@@ -910,7 +910,7 @@ template <class Real> void SphericalHarmonics<Real>::StokesEvalDL(const Vector<R
           Complex<Real> Ycsc_2 = Y_2 * csc_theta;
           if (fabs(cos_theta) == 1) {
             auto Y_csc0 = [exp_phi, cos_theta](Long n, Long m) {
-              if (m == 1) -sqrt<Real>((2*n+1)*n*(n+1)) * ((n%2==0) && (cos_theta<0) ? -1 : 1) * exp_phi;
+              if (m == 1) return -sqrt<Real>((2*n+1)*n*(n+1)) * ((n%2==0) && (cos_theta<0) ? -1 : 1) * exp_phi;
               return Complex<Real>(0, 0);
             };
             Ycsc_0 = Y_csc0(n - 1, m);
@@ -1078,7 +1078,7 @@ template <class Real> void SphericalHarmonics<Real>::StokesEvalKL(const Vector<R
       cos_phi = cos(cos_theta_phi[i * 2 + 1]);
       sin_phi = sin(cos_theta_phi[i * 2 + 1]);
     }
-    Complex<Real> imag(0,1), exp_phi(cos_phi, sin_phi);
+    Complex<Real> imag(0,1), exp_phi(cos_phi, -sin_phi);
 
     StaticArray<Real, COORD_DIM> norm0;
     { // Set norm0 <-- Q^t * norm