Skip to content
Snippets Groups Projects
Commit ee483763 authored by Oleg Kalashev's avatar Oleg Kalashev
Browse files

Turbulent magnetic field parameters adjusted (as in arXiv version)

parent dc81f42b
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
src/app/crbeam/CRbeam.cpp
+ 16
− 10
View file @ ee483763
......@@ -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;
}
......
src/app/crbeam/CRbeam.h
+ 3
− 1
View file @ ee483763
......@@ -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;
......
src/lib/Deflection3D.cpp
+ 14
− 8
View file @ ee483763
......@@ -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);
......
src/lib/Deflection3D.h
+ 3
− 2
View file @ ee483763
......@@ -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;
};
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment