correlation length test added, deflection accuracy param added

src/app/crbeam/CRbeam.cpp

@@ -78,7 +78,7 @@ enum ClineParams
     PNparticles,PBatchSize,
     PRedshift,PFilamentZ,PEnergy,PMinEnergy,PPowerLaw,PMinSourceEnergy,PPrimary,PNoEMcascade,PSourceEvolM,
 	PBackground,PBackgroundMult,PExtraBackground,PExtraBackgroundPhysical,PEBLCut,PMonoCmb,PFixedCmbT,PExtDeltaZconst,PExtPowerLow,PExtDeltaZexp,
-	POutput,POverwriteOutput,PThinning,PRescalePPP,PEGMF,PEGMFLmin,PEGMFLmax,PEGMFNmodes,PRandomEGMF,PTurbulentEGMF,PTauPrint,POutputSuffix
+	POutput,POverwriteOutput,PThinning,PRescalePPP,PEGMF,PEGMFLmin,PEGMFLmax,PEGMFNmodes,PRandomEGMF,PTurbulentEGMF,PDeflectionAccuracy,PTauPrint,POutputSuffix
 };
 #define xstr(s) str(s)
 #define str(s) #s
@@ -150,6 +150,7 @@ static CmdInfo commands[] = {
 		{ PEGMFNmodes,    CmdInfo::FLAG_ARGUMENT, "-mfnm",   "--nmodesEGMF",   "300",    "Number of modes in magnetic field spectrum" },
 		{ PRandomEGMF,    CmdInfo::FLAG_NULL,     "-mfr",   "--randomEGMF",   NULL,   "randomize EGMF (use different EGMF configurations for different initial particles)" },
 		{ PTurbulentEGMF,    CmdInfo::FLAG_NULL,     "-mft",   "--turbulentEGMF",   NULL,   "use turbulent EGMF with Kolmogorov spectrum" },
+        { PDeflectionAccuracy,    CmdInfo::FLAG_ARGUMENT, "-mfac",   "--accEGMF",   "10",    "Deflection simulation accuracy factor" },
 		{ PTauPrint,    CmdInfo::FLAG_NULL,     "-ptau",   "--print-tau",   NULL,   "print tau" },
         { POutputSuffix, CmdInfo::FLAG_ARGUMENT, "-os", "--output-suffix", "", "Output dir name suffix (appended to autogenerated names)"},
 		{ 0 },
@@ -224,7 +225,7 @@ CRbeam::CRbeam(int argc, char** argv):
 		fNmodesEGMF(0),
 		fRandomizeEGMF(false),
 		fTurbulentEGMF(false),
-		//fMaxDeflection(0.5/180.*M_PI),
+        fDeflectionAccuracy(10),
 		fLogging(false),
         fTrajectoryLogging(false),
 		fLogThread(-1),
@@ -258,6 +259,7 @@ CRbeam::CRbeam(int argc, char** argv):
 	fNmodesEGMF = cmd(PEGMFNmodes);
 	fRandomizeEGMF = cmd(PRandomEGMF);
 	fTurbulentEGMF = cmd(PTurbulentEGMF);
+	fDeflectionAccuracy = (int)cmd(PDeflectionAccuracy);
 	//fMaxDeflection = cmd(PMaxDeflection);
 	//fMaxDeflection *= (M_PI/180.);
 
@@ -497,15 +499,21 @@ int CRbeam::run()
 	int Bslot = -1;
     SmartPtr<MagneticField> mf;
 	if(fEGMF>0.) {
+        std::cerr << "fDeflectionAccuracy " << fDeflectionAccuracy << std::endl;
+	    if (fTurbulentEGMF){
+            CosmoTime t0;
+            auto l_cor = TurbulentMF::MeanCorLength(rand, fLminEGMF, fLmaxEGMF, fEGMF, fNmodesEGMF, t0);
+            std::cerr << "l_cor estimate: " << l_cor << " Mpc" << std::endl;
+	    }
         if (fRandomizeEGMF) {
-			Deflection3D* defl = new Deflection3D();
+			Deflection3D* defl = new Deflection3D(fDeflectionAccuracy);
 			Bslot = defl->MFSlot();
 			pe.AddInteraction(defl);
 		}
 		else {
 			mf = fTurbulentEGMF ? ((MagneticField*) new TurbulentMF(rand, fLminEGMF, fLmaxEGMF, fEGMF, fNmodesEGMF)) :
 								((MagneticField*) new MonochromaticMF(fLmaxEGMF, fEGMF, Beta, Alpha, Theta, Phi));
-			pe.AddInteraction(new Deflection3D(mf));
+			pe.AddInteraction(new Deflection3D(mf, fDeflectionAccuracy));
 		}
 	}
 	//// generating initial state

src/app/crbeam/CRbeam.h

@@ -95,6 +95,7 @@ private:
 	double			fNmodesEGMF;
 	bool 			fRandomizeEGMF;
 	bool 			fTurbulentEGMF;
+	uint            fDeflectionAccuracy;
 	bool            fLogging;
 	bool 			fTrajectoryLogging;
 	int             fLogThread;

src/lib/Deflection3D.cpp

@@ -27,7 +27,6 @@
 #include "Deflection3D.h"
 #include "Test.h"
 
-
 namespace Interactions {
     using namespace mcray;
     using namespace Utils;
@@ -213,6 +212,33 @@ namespace Interactions {
         return 0;
     }
 
+    double MagneticField::CorrelationLength(std::vector<double> x, std::vector<double> aDirection, const CosmoTime &aTime,
+                                            double aMaxValMpc) {
+        double step = 0.1*MinVariabilityScale(aTime);
+        int max_steps = (int)(aMaxValMpc*units.Mpc/step + 0.5);
+        if (max_steps == 0)
+            return step/units.Mpc;
+        auto norm = sqrt(aDirection[0]*aDirection[0] + aDirection[1]*aDirection[1] + aDirection[2]*aDirection[2]);
+        for(int dim_i=0; dim_i<3; dim_i++)
+            aDirection[dim_i] *= (step/norm);
+
+        auto iniB = std::vector<double>(3);
+        auto curB = std::vector<double>(3);
+        GetValueGauss(x.data(), aTime, iniB);
+        for(int step_i=0; step_i<max_steps; step_i++){
+            for(int dim_i=0; dim_i<3; dim_i++)
+                x[dim_i] += aDirection[dim_i];
+            GetValueGauss(x.data(), aTime, curB);
+            double sc_product = 0.;
+            for(int dim_i=0; dim_i<3; dim_i++)
+                sc_product += curB[dim_i]*iniB[dim_i];
+            if (sc_product < 0)
+                return step*(step_i+1)/units.Mpc;
+        }
+        return aMaxValMpc;
+    }
+
+
     int MonochromaticMF::UnitTest() {
         std::vector<double> B(3);
         CosmoTime t(0);
@@ -294,6 +320,18 @@ namespace Interactions {
         }
     }
 
+    double TurbulentMF::MeanCorLength(Randomizer &aRandomizer, double aLmin_Mpc, double aLmax_Mpc, double aVariance_Gauss,
+                                      int aNmodes, const CosmoTime &aTime, int aNsamples, double aMaxValMpc){
+        double l_cor_sum = 0.;
+        auto x = std::vector<double>({0.,0.,0.});
+        auto dx = std::vector<double>({0.,0.,1.});
+        for(int i=0; i<aNsamples; i++){
+            TurbulentMF B(aRandomizer, aLmin_Mpc, aLmax_Mpc, aVariance_Gauss, aNmodes);
+            l_cor_sum += B.CorrelationLength(x, dx, aTime, aMaxValMpc);
+        }
+        return l_cor_sum/aNsamples;
+    }
+
     int TurbulentMF::UnitTest() {
         if(!cosmology.IsInitialized())
             cosmology.Init();

src/lib/Deflection3D.h

@@ -44,6 +44,10 @@ namespace Interactions {
         virtual double MinVariabilityScale(const CosmoTime &aTime) const = 0;
 
         int Print(double lambdaMpc, std::ostream& aOutput=std::cout);
+
+        double CorrelationLength(std::vector<double> x, std::vector<double> aDirection, const CosmoTime &aTime,
+                                         double aMaxValMpc=10000);
+
     };
 
     class Deflection3D : public DeflectionInteraction {
@@ -110,6 +114,8 @@ namespace Interactions {
         virtual void GetValueGauss(const double *x, const CosmoTime &aTime, std::vector<double> &outValue) const;
         virtual double MinVariabilityScale(const CosmoTime &aTime) const;
         static int UnitTest();
+        static double MeanCorLength(Randomizer &aRandomizer, double aLmin_Mpc, double aLmax_Mpc, double aVariance_Gauss,
+                                    int aNmodes, const CosmoTime &aTime, int aNsamples=1000, double aMaxValMpc=10000);
     private:
         double fLmin;
         double fLmax;