diff --git a/src/app/crbeam/CRbeam.cpp b/src/app/crbeam/CRbeam.cpp
index d30120b564768dfb52eb50aa957e7127796502e7..e9b528b3ad46258f78ebe11845d68a179b6b2dd3 100644
--- a/src/app/crbeam/CRbeam.cpp
+++ b/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,7 +150,8 @@ 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" },
- { PTauPrint, CmdInfo::FLAG_NULL, "-ptau", "--print-tau", NULL, "print tau" },
+ { 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.) {
- if (fRandomizeEGMF) {
- Deflection3D* defl = new Deflection3D();
+ 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(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
diff --git a/src/app/crbeam/CRbeam.h b/src/app/crbeam/CRbeam.h
index 6affaf2d5680c6309d6c527d7426265072d0a4a9..9c94ccc12e0c133f4253cc10cad6bfed9845cb9f 100644
--- a/src/app/crbeam/CRbeam.h
+++ b/src/app/crbeam/CRbeam.h
@@ -95,6 +95,7 @@ private:
double fNmodesEGMF;
bool fRandomizeEGMF;
bool fTurbulentEGMF;
+ uint fDeflectionAccuracy;
bool fLogging;
bool fTrajectoryLogging;
int fLogThread;
diff --git a/src/lib/Deflection3D.cpp b/src/lib/Deflection3D.cpp
index eaec1474c7aacd85b4eb953361abd01eb5283dcf..3febb8daf09a8917140a10cbc05051cf6b986b38 100644
--- a/src/lib/Deflection3D.cpp
+++ b/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();
diff --git a/src/lib/Deflection3D.h b/src/lib/Deflection3D.h
index 7a2938dfffa2673e79b17caa29f7d2de86108295..dd1dcfa35c8cf0c826088870fc85f547b5f0df54 100644
--- a/src/lib/Deflection3D.h
+++ b/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;