Many bug and warning fixes

.github/workflows/create_coverage_report.yml

@@ -43,7 +43,7 @@ jobs:
         - name: run test
           run: | 
             cd build
-            make test
+            ctest --output-on-failure
           continue-on-error: true
         - name: coverage report
           run: |

.github/workflows/create_documentation.yml

@@ -44,7 +44,7 @@ jobs:
       - name: run test
         run: | 
           cd build
-          make test
+          ctest --output-on-failure
         continue-on-error: true
       - name: coverage report
         run: |

.github/workflows/deploy_new_documentation.yml

@@ -48,7 +48,7 @@ jobs:
       - name: run test
         run: | 
           cd build
-          make test
+          ctest --output-on-failure
         continue-on-error: true
       - name: coverage report
         run: |

.github/workflows/testing_OSX.yml

@@ -54,7 +54,7 @@ jobs:
          CRPROPA_DATA_PATH: "/Users/runner/work/CRPropa3/CRPropa3/build/data"
        run: |
          cd build
-         make test
+         ctest --output-on-failure
      - name: Archive test results
        if: always()
        uses: actions/upload-artifact@v4

.github/workflows/testing_ubuntu22.yml

@@ -41,7 +41,7 @@ jobs:
       - name: Run tests
         run: |
           cd build
-          make test
+          ctest --output-on-failure
       - name: Archive test results
         if: always()
         uses: actions/upload-artifact@v4

.github/workflows/testing_ubuntu24.yml

@@ -39,7 +39,7 @@ jobs:
       - name: Run tests
         run: |
           cd build
-          make test
+          ctest --output-on-failure
       - name: Archive test results
         if: always()
         uses: actions/upload-artifact@v4
@@ -84,7 +84,7 @@ jobs:
       - name: Run tests
         run: |
           cd build
-          make test 
+          ctest --output-on-failure
       - name: Archive test results
         if: always()
         uses: actions/upload-artifact@v4

cmake/FindGooglePerfTools.cmake

@@ -7,7 +7,7 @@
 #  STACKTRACE_LIBRARIES, where to find the stacktrace library.
 #  PROFILER_LIBRARIES, where to find the profiler library.
 
-FIND_PATH(GOOGLE_PERFTOOLS_INCLUDE_DIR google/heap-profiler.h)
+FIND_PATH(GOOGLE_PERFTOOLS_INCLUDE_DIR NAMES google/heap-profiler.h gperftools/heap-profiler.h)
 
 SET(TCMALLOC_NAMES ${TCMALLOC_NAMES} tcmalloc)
 FIND_LIBRARY(TCMALLOC_LIBRARY  NAMES ${TCMALLOC_NAMES})

include/crpropa/AssocVector.h

@@ -46,12 +46,12 @@ namespace Loki
 
     namespace Private
     {
-        template <class Value, class C>
+        template <class Value, class C, class K>
         class AssocVectorCompare : public C
         {
-            typedef std::pair<typename C::first_argument_type, Value>
+            typedef K first_argument_type;
+            typedef std::pair<first_argument_type, Value>
                 Data;
-            typedef typename C::first_argument_type first_argument_type;
 
         public:
             AssocVectorCompare()
@@ -98,10 +98,10 @@ namespace Loki
     >
     class AssocVector
         : private std::vector< std::pair<K, V>, A >
-        , private Private::AssocVectorCompare<V, C>
+        , private Private::AssocVectorCompare<V, C, K>
     {
         typedef std::vector<std::pair<K, V>, A> Base;
-        typedef Private::AssocVectorCompare<V, C> MyCompare;
+        typedef Private::AssocVectorCompare<V, C, K> MyCompare;
 
     public:
         typedef K key_type;
@@ -110,20 +110,19 @@ namespace Loki
 
         typedef C key_compare;
         typedef A allocator_type;
-        typedef typename A::reference reference;
-        typedef typename A::const_reference const_reference;
+        typedef typename A::value_type& reference;
+        typedef const typename A::value_type* const_reference;
         typedef typename Base::iterator iterator;
         typedef typename Base::const_iterator const_iterator;
         typedef typename Base::size_type size_type;
         typedef typename Base::difference_type difference_type;
-        typedef typename A::pointer pointer;
-        typedef typename A::const_pointer const_pointer;
+        typedef typename A::value_type* pointer;
+        typedef const typename A::value_type* const_pointer;
         typedef typename Base::reverse_iterator reverse_iterator;
         typedef typename Base::const_reverse_iterator const_reverse_iterator;
 
         class value_compare
-            : public std::binary_function<value_type, value_type, bool>
-            , private key_compare
+            : private key_compare
         {
             friend class AssocVector;
 

include/crpropa/Random.h

@@ -92,7 +92,7 @@ namespace crpropa {
 
  Mersenne Twister random number generator -- a C++ class Random
  Based on code by Makoto Matsumoto, Takuji Nishimura, and Shawn Cokus
- Richard J. Wagner  v1.0  15 May 2003  rjwagner@writeme.com
+ Richard J. Wagner  v1.0  15 May 2003  rjwagner\@writeme.com
  */
 class Random {
 public:

include/crpropa/Vector3.h

@@ -420,7 +420,7 @@ class Vector3 {
 
 	const std::string getDescription() {
 		char buffer[256];
-		sprintf(buffer, "Vector(%.6G, %.6G, %.6G)", data[0], data[1], data[2]);
+		snprintf(buffer, 256, "Vector(%.6G, %.6G, %.6G)", data[0], data[1], data[2]);
 		return buffer;
 	}
 

include/crpropa/magneticField/turbulentField/SimpleGridTurbulence.h

@@ -30,7 +30,7 @@ class SimpleTurbulenceSpectrum : public TurbulenceSpectrum {
 	*/
 	SimpleTurbulenceSpectrum(double Brms, double lMin, double lMax,
 	                         double sIndex = 5. / 3)
-	    : TurbulenceSpectrum(Brms, lMin, lMax, constScaleBendover * lMax, sIndex, 0) {}
+	    : TurbulenceSpectrum(Brms, lMin, lMax, 1000 * lMax, sIndex, 0) {}
 	~SimpleTurbulenceSpectrum() {}
 
 	/**

include/crpropa/magneticField/turbulentField/TurbulentField.h

@@ -81,7 +81,7 @@ class TurbulenceSpectrum : public Referenced {
 
 	/**
   Computes the magnetic field coherence length
-	Obtained from the definition of \f$l_c = 1/B_{\rm rms}^2 \int_0^\infty dr\langleB(0)B^*(r)\rangle \f$
+	Obtained from the definition of \\f$l_c = 1/B_{\\rm rms}^2 \\int_0^\\infty dr\\langleB(0)B^*(r)\\rangle \\f$
 	Approximates the true value correctly as long as lBendover <= lMax/8 (~5%
   error) (for the true value the above integral should go from lMin to lMax)
 	*/

include/crpropa/module/Acceleration.h

@@ -122,13 +122,13 @@ class DirectedFlowOfScatterCenters : public StepLengthModifier {
 /// @class QuasiLinearTheory
 /// @brief Scales the steplength according to quasi linear theory.
 /// @details Following quasi-linear theory [Schlickeiser1989], the mean free
-/// path \f$\lambda\f$ of a
-///  particle with energy \f$E\f$ and charge \f$Z\f$ in a field with turbulence
-///  spectrum \f$\frac{k}{k_{\min}}^{-q}\f$ is
-///   \f$ \lambda = {\left(\frac{B}{\delta B}\right)}^2 {\left(R_G\;
-///   k_{\min}\right)}^{1-q} R_G \equiv \lambda_0 {\left( \frac{E}{1
-///   EeV}\frac{1}{Z} \right)}^{2-q} \f$
-///  where \f$R_G = \frac{E}{B Z}\f$ is the gyro-radius of the
+/// path \\f$\\lambda\\f$ of a
+///  particle with energy \\f$E\\f$ and charge \\f$Z\\f$ in a field with turbulence
+///  spectrum \\f$\\frac{k}{k_{\\min}}^{-q}\\f$ is
+///   \\f$ \\lambda = {\\left(\\frac{B}{\\delta B}\\right)}^2 {\\left(R_G\\;
+///   k_{\\min}\\right)}^{1-q} R_G \\equiv \\lambda_0 {\\left( \\frac{E}{1
+///   EeV}\\frac{1}{Z} \\right)}^{2-q} \\f$
+///  where \\f$R_G = \\frac{E}{B Z}\\f$ is the gyro-radius of the
 ///  particles.
 /// This class implements the rigidity dependent scaling factor used to modify
 /// the base step length.

libs/sophia/sophia_interface.f

@@ -2568,45 +2568,55 @@ SUBROUTINE DECPAR(LA,P0,ND,LL,P)
 
 C...generation of the masses, compute weight, if rejected try again
 240   RORD(1) = 1.D0
-      DO 260 IL1=2,ND-1
-      RSAV = RNDM(0)
-      DO 250 IL2=IL1-1,1,-1
-      IF(RSAV.LE.RORD(IL2))   GOTO 260
-250     RORD(IL2+1)=RORD(IL2)
+      DO IL1=2,ND-1
+        RSAV = RNDM(0)
+        DO IL2=IL1-1,1,-1
+          IF(RSAV.LE.RORD(IL2))   GOTO 260
+          RORD(IL2+1)=RORD(IL2)
+        ENDDO
 260     RORD(IL2+1)=RSAV
+      ENDDO
       RORD(ND) = 0.D0
       WT = 1.D0
-      DO 270 IL=ND-1,1,-1
-      PV(IL,5)=PV(IL+1,5)+P(IL,5)+(RORD(IL)-RORD(IL+1))*(PV(1,5)-PS)
-270   WT=WT*PAWT(PV(IL,5),PV(IL+1,5),P(IL,5))
+      DO IL=ND-1,1,-1
+        PV(IL,5)=PV(IL+1,5)+P(IL,5)+(RORD(IL)-RORD(IL+1))*(PV(1,5)-PS)
+        WT=WT*PAWT(PV(IL,5),PV(IL+1,5),P(IL,5))
+      ENDDO
       IF (WT.LT.RNDM(0)*WWTMAX)   GOTO 240
 
 C...Perform two particle decays in respective cm frame
-280   DO 300 IL=1,ND-1
-      PA=PAWT(PV(IL,5),PV(IL+1,5),P(IL,5))
-      UE(3)=2.D0*RNDM(0)-1.D0
-      PHI=2.D0*PI*RNDM(0)
-      UT = SQRT(1.D0-UE(3)**2)
-      UE(1) = UT*COS(PHI)
-      UE(2) = UT*SIN(PHI)
-      DO 290 J=1,3
-      P(IL,J)=PA*UE(J)
-290   PV(IL+1,J)=-PA*UE(J)
-      P(IL,4)=SQRT(PA**2+P(IL,5)**2)
-300   PV(IL+1,4)=SQRT(PA**2+PV(IL+1,5)**2)
+280   DO IL=1,ND-1
+        PA=PAWT(PV(IL,5),PV(IL+1,5),P(IL,5))
+        UE(3)=2.D0*RNDM(0)-1.D0
+        PHI=2.D0*PI*RNDM(0)
+        UT = SQRT(1.D0-UE(3)**2)
+        UE(1) = UT*COS(PHI)
+        UE(2) = UT*SIN(PHI)
+        DO J=1,3
+          P(IL,J)=PA*UE(J)
+          PV(IL+1,J)=-PA*UE(J)
+        ENDDO
+        P(IL,4)=SQRT(PA**2+P(IL,5)**2)
+        PV(IL+1,4)=SQRT(PA**2+PV(IL+1,5)**2)
+      ENDDO
 
 C...Lorentz transform decay products to lab frame
-      DO 310 J=1,4
-310   P(ND,J)=PV(ND,J)
-      DO 340 IL=ND-1,1,-1
-      DO 320 J=1,3
-320   BE(J)=PV(IL,J)/PV(IL,4)
-      GA=PV(IL,4)/PV(IL,5)
-      DO 340 I=IL,ND
-      BEP = BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
-      DO 330 J=1,3
-330   P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
-340   P(I,4)=GA*(P(I,4)+BEP)
+      DO J=1,4
+        P(ND,J)=PV(ND,J)
+      ENDDO
+      DO IL=ND-1,1,-1
+        DO J=1,3
+          BE(J)=PV(IL,J)/PV(IL,4)
+        ENDDO
+        GA=PV(IL,4)/PV(IL,5)
+        DO I=IL,ND
+          BEP = BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
+          DO J=1,3
+            P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
+          ENDDO
+          P(I,4)=GA*(P(I,4)+BEP)
+        ENDDO
+      ENDDO
 
 C...Weak decays
         IF (MAT .EQ. 1)  THEN
@@ -2621,14 +2631,17 @@ SUBROUTINE DECPAR(LA,P0,ND,LL,P)
 
 C...Boost back for rapidly moving particle
       IF (MBST .EQ. 1)   THEN
-         DO 440 J=1,3
-440      BE(J)=P0(J)/P0(4)
-         GA= P0(4)/P0(5)
-         DO 460 I=1,ND
-         BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
-         DO 450 J=1,3
-450         P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
-460         P(I,4)=GA*(P(I,4)+BEP)
+        DO J=1,3
+          BE(J)=P0(J)/P0(4)
+        ENDDO
+        GA= P0(4)/P0(5)
+        DO I=1,ND
+          BEP=BE(1)*P(I,1)+BE(2)*P(I,2)+BE(3)*P(I,3)
+          DO J=1,3
+            P(I,J)=P(I,J)+GA*(GA*BEP/(1.+GA)+P(I,4))*BE(J)
+          ENDDO
+          P(I,4)=GA*(P(I,4)+BEP)
+        ENDDO
       ENDIF
 
 C...labels for antiparticle decay
@@ -2804,8 +2817,9 @@ DOUBLE PRECISION FUNCTION PO_RNDGAM(ALAM,ETA)
       IF(RNDM(0).GT.Y**(F-1.D0)) GOTO 10
    40 IF(N.EQ.0) GOTO 70
    50 Z = 1.D0
-      DO 60 I = 1,N
-   60 Z = Z*RNDM(0)
+      DO I = 1,N
+        Z = Z*RNDM(0)
+      ENDDO
       Y = Y-LOG(Z+1.D-7)
    70 PO_RNDGAM = Y/ALAM
 

src/ProgressBar.cpp

@@ -52,7 +52,7 @@ void ProgressBar::setPosition(unsigned long position) {
 		std::string s = " - Finished at ";
 		s.append(ctime(&currentTime));
 		char fs[255];
-		std::sprintf(fs, "%c[%d;%dm Finished %c[%dm", 27, 1, 32, 27, 0);
+		std::snprintf(fs, 255, "%c[%d;%dm Finished %c[%dm", 27, 1, 32, 27, 0);
 		std::printf(stringTmpl.c_str(), fs, 100, "Needed",
 				int(tElapsed / 3600), (int(tElapsed) % 3600) / 60,
 				int(tElapsed) % 60, s.c_str());
@@ -68,7 +68,7 @@ void ProgressBar::setError() {
 	std::string s = " - Finished at ";
 	s.append(ctime(&currentTime));
 	char fs[255];
-	std::sprintf(fs, "%c[%d;%dm  ERROR   %c[%dm", 27, 1, 31, 27, 0);
+	std::snprintf(fs, 255, "%c[%d;%dm  ERROR   %c[%dm", 27, 1, 31, 27, 0);
 	std::printf(stringTmpl.c_str(), fs, _currentCount, "Needed",
 			int(tElapsed / 3600), (int(tElapsed) % 3600) / 60,
 			int(tElapsed) % 60, s.c_str());

src/Random.cpp

@@ -282,7 +282,8 @@ uint64_t Random::randInt64()
 {
 	int64_t a = randInt();
 	int64_t b = randInt();
-	return (b + a << 32);
+	// a is shifted to the left to create a proper 64 bit integer
+	return (b + (a << 32));
 }
 
 
@@ -390,11 +391,11 @@ void Random::initialize(const uint32_t seed) {
 void Random::reload() {
 	uint32_t *p = state;
 	int i;
-	for (i = N - M; i--; ++p)
+	for (i = +N - M; i--; ++p)
 		*p = twist(p[M], p[0], p[1]);
 	for (i = M; --i; ++p)
-		*p = twist(p[M - N], p[0], p[1]);
-	*p = twist(p[M - N], p[0], state[0]);
+		*p = twist(p[+M - N], p[0], p[1]);
+	*p = twist(p[+M - N], p[0], state[0]);
 
 	left = N, pNext = state;
 }

src/Variant.cpp

@@ -668,6 +668,7 @@ Variant Variant::fromString(const std::string& s, Type t) {
 			std::string s;
 			v.push_back(s);
 		}
+		return Variant(v);
 	} else {
 		std::string msg; 
 		msg += "fromString not implemented for type ";
@@ -1010,7 +1011,7 @@ Variant::operator const std::vector<Variant>&() const {
 
 #define INT_FUNCTION(to_type, fun, to)                                                                                 \
 	to Variant::fun() const {                                                                                          \
-		switch (type) {                                                                                                \
+		switch (to_type) {                                                                                                \
 			case Variant::TYPE_BOOL:                                                                                   \
 				return data._t_bool ? 1 : 0;                                                                           \
 				break;                                                                                                 \

src/magneticField/CMZField.cpp

@@ -192,7 +192,7 @@ Vector3d CMZField::getMCField(const Vector3d& pos) const {//Field in molecular c
         Bphi = - B1*exp(-z*z/Hc/Hc)*R/r;
     }
     else{
-        - B1*exp(-z*z/Hc/Hc)*R/r1*(3*r/r1- 2*r*r/r1*r1);
+        Bphi = - B1*exp(-z*z/Hc/Hc)*R/r1*(3*r/r1- 2*r*r/r1*r1);
     }
 
     b.x -= Bphi*sin(phi);

src/magneticField/MagneticField.cpp

@@ -26,7 +26,7 @@ Vector3d &PeriodicMagneticField::getExtends() {
 	return extends;
 }
 
-void PeriodicMagneticField::setExtends(const Vector3d &origin) {
+void PeriodicMagneticField::setExtends(const Vector3d &extends) {
 	this->extends = extends;
 }
 

src/magneticField/turbulentField/PlaneWaveTurbulence.cpp

@@ -59,7 +59,7 @@
 
 #ifdef ENABLE_FAST_WAVES
 #include <immintrin.h>
-#include <memory.h>
+#include <memory>
 #endif
 
 namespace crpropa {