Merge pull request #197 from EXP-code/covarianceBugfix

Covariance bugfix

Author: michael-petersen

CMakeLists.txt

@@ -2,7 +2,7 @@ cmake_minimum_required(VERSION 3.25) # Needed for CUDA, MPI, and CTest features
 
 project(
   EXP
-  VERSION "7.9.1"
+  VERSION "7.9.2"
   HOMEPAGE_URL https://github.com/EXP-code/EXP
   LANGUAGES C CXX Fortran)
 

doc/exp.cfg

@@ -48,7 +48,7 @@ PROJECT_NAME           = EXP
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         = 7.9.1
+PROJECT_NUMBER         = 7.9.2
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a

doc/exp.cfg.breathe

@@ -48,7 +48,7 @@ PROJECT_NAME           = EXP
 # could be handy for archiving the generated documentation or if some version
 # control system is used.
 
-PROJECT_NUMBER         = 7.9.1
+PROJECT_NUMBER         = 7.9.2
 
 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a

exputil/EmpCylSL.cc

@@ -4102,8 +4102,20 @@ void EmpCylSL::accumulate(double r, double z, double phi, double mass,
     }
 
     if (compute and covar) {
-      Eigen::VectorXcd vec = std::complex<double>(mcos, msin) *
-	vc[id].row(mm).transpose() * norm;
+      int size = vc[id].row(mm).size();
+      if (size != NORDER) {
+        throw GenericError("size of vectors must match", __FILE__, __LINE__, 1039, false);
+      }
+      Eigen::VectorXcd vec(size);
+      
+      Eigen::VectorXd vC = vc[id].row(mm).transpose() * norm;
+      Eigen::VectorXd vS = vs[id].row(mm).transpose() * norm;
+      
+      // For m = 0, the azimuthal dependence vanishes; set vS to zero to ensure a purely real result
+      if (mm==0) vS.setZero();
+      
+      vec.real() = vC*mcos + vS*msin;
+      vec.imag() = vC*msin - vS*mcos;
 
       VC[id][whch][mm] += mass * vec;
       MV[id][whch][mm] += mass * (vec * vec.adjoint());