diff --git a/src/app/crbeam/CRbeam.cpp b/src/app/crbeam/CRbeam.cpp
index 8987c5042e7937e177b3e4783bffa98d342e1b19..b10e19d746b777653352ae3e822efde3a1dca14c 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,PEGMFL,PRandomEGMF,PTurbulentEGMF,PTauPrint,POutputSuffix
+ POutput,POverwriteOutput,PThinning,PRescalePPP,PEGMF,PEGMFLmin,PEGMFLmax,PEGMFNmodes,PRandomEGMF,PTurbulentEGMF,PTauPrint,POutputSuffix
};
#define xstr(s) str(s)
#define str(s) #s
@@ -145,7 +145,9 @@ static CmdInfo commands[] = {
{ PRescalePPP, CmdInfo::FLAG_ARGUMENT, "-pt", "--ppp-thinning", "10", "PPP Rescale Coefficient (double > 0, actual number of secondaries per interaction)" },
{ PEGMF, CmdInfo::FLAG_ARGUMENT, "-mf", "--EGMF", "1e-15", "Extragalactic magnetic field in Gauss" },
- { PEGMFL, CmdInfo::FLAG_ARGUMENT, "-mfl", "--lEGMF", "1", "Extragalactic magnetic field correlation length in Mpc at z=0" },
+ { PEGMFLmin, CmdInfo::FLAG_ARGUMENT, "-mflmin", "--lminEGMF", "0.05", "Extragalactic magnetic field minimum scale in Mpc at z=0" },
+ { PEGMFLmax, CmdInfo::FLAG_ARGUMENT, "-mflmax", "--lmaxEGMF", "5.0", "Extragalactic magnetic field maximum scale in Mpc at z=0" },
+ { 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" },
@@ -217,7 +219,9 @@ CRbeam::CRbeam(int argc, char** argv):
fOverwriteOutput(false),
fFixedCmb(false),
fEGMF(0.),
- fLcorEGMF(1.),
+ fLminEGMF(0.),
+ fLmaxEGMF(0.),
+ fNmodesEGMF(0),
fRandomizeEGMF(false),
fTurbulentEGMF(false),
//fMaxDeflection(0.5/180.*M_PI),
@@ -249,7 +253,9 @@ CRbeam::CRbeam(int argc, char** argv):
fPPPrescaleCoef = cmd(PRescalePPP);
fEGMF = cmd(PEGMF);
- fLcorEGMF = cmd(PEGMFL);
+ fLminEGMF = cmd(PEGMFLmin);
+ fLmaxEGMF = cmd(PEGMFLmax);
+ fNmodesEGMF = cmd(PEGMFNmodes);
fRandomizeEGMF = cmd(PRandomEGMF);
fTurbulentEGMF = cmd(PTurbulentEGMF);
//fMaxDeflection = cmd(PMaxDeflection);
@@ -316,8 +322,8 @@ int CRbeam::run()
if(fTurbulentEGMF)
outputDir += "Turb";
outputDir += ToString(fEGMF);
- outputDir += "_Lc";
- outputDir += ToString(fLcorEGMF);
+ outputDir += "_Lmax";
+ outputDir += ToString(fLmaxEGMF);
}
else
{
@@ -495,8 +501,8 @@ int CRbeam::run()
pe.AddInteraction(defl);
}
else {
- mf = fTurbulentEGMF ? ((MagneticField*) new TurbulentMF(rand, fLcorEGMF, fEGMF)) :
- ((MagneticField*) new MonochromaticMF(fLcorEGMF, fEGMF, Beta, Alpha, Theta, Phi));
+ mf = fTurbulentEGMF ? ((MagneticField*) new TurbulentMF(rand, fLminEGMF, fLmaxEGMF, fEGMF, fNmodesEGMF)) :
+ ((MagneticField*) new MonochromaticMF(fLmaxEGMF, fEGMF, Beta, Alpha, Theta, Phi));
pe.AddInteraction(new Deflection3D(mf));
}
}
@@ -587,8 +593,8 @@ int CRbeam::run()
if(fRandomizeEGMF)
{
- MagneticField* mf = fTurbulentEGMF ? ((MagneticField*) new TurbulentMF(rand, fLcorEGMF, fEGMF)) :
- ((MagneticField*) new MonochromaticMF(rand, fLcorEGMF, fEGMF));
+ MagneticField* mf = fTurbulentEGMF ? ((MagneticField*) new TurbulentMF(rand, fLminEGMF, fLmaxEGMF, fEGMF, fNmodesEGMF)) :
+ ((MagneticField*) new MonochromaticMF(rand, fLmaxEGMF, fEGMF));
fields[i%fBatchSize] = mf; //store till the end of current batch
particle.interactionData[Bslot] = mf;
}
diff --git a/src/app/crbeam/CRbeam.h b/src/app/crbeam/CRbeam.h
index a905ac2b328b9ade35e5839f8b2e084879302693..6affaf2d5680c6309d6c527d7426265072d0a4a9 100644
--- a/src/app/crbeam/CRbeam.h
+++ b/src/app/crbeam/CRbeam.h
@@ -90,7 +90,9 @@ private:
bool fOverwriteOutput;
bool fFixedCmb;
double fEGMF;
- double fLcorEGMF;
+ double fLminEGMF;
+ double fLmaxEGMF;
+ double fNmodesEGMF;
bool fRandomizeEGMF;
bool fTurbulentEGMF;
bool fLogging;
diff --git a/src/lib/Deflection3D.cpp b/src/lib/Deflection3D.cpp
index ef30f02ee871ff30fd48045601efb9233a473feb..eaec1474c7aacd85b4eb953361abd01eb5283dcf 100644
--- a/src/lib/Deflection3D.cpp
+++ b/src/lib/Deflection3D.cpp
@@ -260,28 +260,34 @@ namespace Interactions {
for(std::vector<SafePtr<MonochromaticMF> >::const_iterator pw = fWaves.begin(); pw!=fWaves.end(); pw++){
(*pw)->GetValueGauss(x, aTime, wave);
for(uint i=0; i<3; i++)
- outValue[i] += (wave[i]);
+ outValue[i] += std::sqrt(2)*(wave[i]);
}
}
double TurbulentMF::MinVariabilityScale(const CosmoTime &aTime) const {
- return fLc*0.1;
+ return fLmax*0.02;
}
- TurbulentMF::TurbulentMF(Randomizer &aRandomizer, double aLcor_Mpc, double aVariance_Gauss, double aMultK):
- fLc(aLcor_Mpc*units.Mpc)
+ TurbulentMF::TurbulentMF(Randomizer &aRandomizer, double aLmin_Mpc, double aLmax_Mpc, double aVariance_Gauss, int aNmodes):
+ fLmin(aLmin_Mpc*units.Mpc),
+ fLmax(aLmax_Mpc*units.Mpc)
{
- ASSERT(aMultK>1.);
+ ASSERT(aNmodes>1.);
+ ASSERT(aLmax_Mpc>aLmin_Mpc);
const double gamma = 11./3.;
std::vector<double> norms;
std::vector<double> Ks;
double sumNorm = 0;
- for(double k = 1./fLc/16.; k*fLc < 4e4; k*=aMultK){//here we limit the range of waves by condition normK>~1e-3
+ double aMultK = pow(aLmax_Mpc/aLmin_Mpc, 1./(aNmodes-1.));
+ double k = 2*M_PI/fLmax;
+ // std::cout << aLmin_Mpc << '\t' << aLmax_Mpc << '\t' << aNmodes << '\t' << aMultK << '\n';
+ for(uint i=0; i<aNmodes; i++){
Ks.push_back(k);
- double normK = k*k*k/(1.+pow(k*fLc,gamma));
+ double normK = (k*fLmax)*(k*fLmax)*(k*fLmax)*pow(k*fLmax,-gamma);
norms.push_back(normK);
sumNorm += normK;
fWaves.push_back(0);
+ k *= aMultK;
}
for(uint i=0; i<norms.size(); i++){
fWaves[i] = new MonochromaticMF(aRandomizer, 2.*M_PI/Ks[i]/units.Mpc, aVariance_Gauss*sqrt(norms[i]/sumNorm));
@@ -297,7 +303,7 @@ namespace Interactions {
Randomizer r;
double lambdaMpc = 1;
double Bgauss = 1e-15;
- SmartPtr<MagneticField> mf = new TurbulentMF(r, lambdaMpc, Bgauss);
+ SmartPtr<MagneticField> mf = new TurbulentMF(r, lambdaMpc/20, lambdaMpc*5, Bgauss, 100);
Deflection3D defl(mf);
std::ofstream outB("TurbulentMF_B.txt");
mf->Print(lambdaMpc, outB);
diff --git a/src/lib/Deflection3D.h b/src/lib/Deflection3D.h
index 97feefcaf251169aa672839ea927b604f0dc5be2..7a2938dfffa2673e79b17caa29f7d2de86108295 100644
--- a/src/lib/Deflection3D.h
+++ b/src/lib/Deflection3D.h
@@ -106,12 +106,13 @@ namespace Interactions {
class TurbulentMF : public MagneticField{
public:
- TurbulentMF(Randomizer &aRandomizer, double aLcor_Mpc, double aVariance_Gauss, double aMultK=2./* steps in K */);
+ TurbulentMF(Randomizer &aRandomizer, double aLmin_Mpc, double aLmax_Mpc, double aVariance_Gauss, int aNmodes);
virtual void GetValueGauss(const double *x, const CosmoTime &aTime, std::vector<double> &outValue) const;
virtual double MinVariabilityScale(const CosmoTime &aTime) const;
static int UnitTest();
private:
- double fLc;
+ double fLmin;
+ double fLmax;
std::vector<SafePtr<MonochromaticMF> > fWaves;
};