Magnetic field setting by Dolag field

[wilson-2015]

I am trying to plot the cosmic ray energy spectrum with magnetic field by Dolag but I am not getting the expected flat spectrum.

I am simulating within the energy range of 3 EeV to 150 EeV without interactions. The spectrum has strong modulations at low energy but straightens towards higher energy as shown in the plot attached Un1.pdf.

I don't understand what I am Missing out in the magnetic field setting kindly have a look at the section of my script below and advice please:

---

filename_bfield = ('/common/workspace/sim/bfield/dolag_B_54-186Mpc_440b.raw')
gridOrigin = Vector3d(54, 54, 54)
boxOrigin gridSize = 440
size = 186Mpc
b_factor = 1.0
randomSeed = 42 turbSpectrum = SimpleTurbulenceSpectrum(8 nG, 60 kpc, 800 kpc, 5./3) gridprops = GridProperties(gridOrigin, 440, 30 * kpc) vgrid=Grid3f(gridprops) loadGrid(vgrid, filename_bfield, b_factor) BField = SimpleGridTurbulence(turbSpectrum, gridprops, randomSeed)

simulation setup

sim = ModuleList() sim.add(PropagationCK(BField, 1e-2, 10 kpc, 100 kpc))

[Laruku25]

This is an issue that's preventing me from finish my paper too. Nobody knows what parameters to plug in for B factor, origin and even the size might not be 186*Mpc.

In dolag's paper it says it's about 115Mpc. Then you have CRPROPAs paper using 132Mpc. The file says 54-186Mpc. Does that mean there's a high resolution of the field up to 54Mpc? And 186*Mpc is total size?

The B factor is another issue which I have resolved, the correct value is B=1e-1. Then I can match the maximum field strength and the minium field strength outlined by dolag's paper, the maximum is ~1uG which is in clusters.

Now even if I get the correct B_factor there's another issue:

You don't get the same deflections that Dolag does. A 50*EeV proton for example should onl be bent by 5 degrees according to his paper and only if it goes through clusters but Look at the graph attached which has the following settings. Note that I've used the backtrack method to get the trajectory.

B=1e-1 (this replicates the min and max field strengths of Dolag's paper) origin=(0.5size,0.5size,0.5*size) basically the middle of the box. If you put the origin at 0,0,0, then you will start at the corner the box which is wrong.

The labels are in Tesla, you can change them to uG (to match Dolag's Paper) multiplying by 10e4(Tesla to Gauss)*10e6(uG).

That's definitely not the straight line path Dolag was preaching (note that this graph is in 2D so you can't see the impact of the clusters well in this image, if you change the B_factor to 1e-2, then the impact of the clusters on angles are more obvious in a 2D graph but they are still not 5 degrees they will look more or less the same as this graph, I will attach different graphs with different B factors later).

By the way why are you using a turbulent field on top of the dolag field, maybe I am missing something. Please let us know anything you may know about this Dolag's file.

I think our best would be to use the Dolag's quimby file but unfortunately the download link hasn't been working for a while now. I do know people have published stuff using that file though so it must be working well!

Also note that the origin at 0,0,0 doesn't match the actual position of us as an observer, i tried your 54*Mpc, and the simulation still starts outside of clusters like this graph. :(

This is the trajectory that you get with your B_factor of 1.0. It's for a 50*EeV.

This is the trajectory with Bfactor=1.0 and a 150*EeV proton.

Anyways I think your particles are propagating all over the place so at low energies the deflections are super strong so there's a higher chance of detection at the observer you choose, at higher energies the deflections are still all over the place so they will not flatten, you will still detect quite a lot because they still make a lot of turns though not as much as low energy particles.

Can you tell me more about your simulation? Why does the spectrum have to flatten? I have a lot to learn :(.

[wilson-2015]

Hi Laruku, Thanks for your detailed explanation and sorry for a late reply. Starting with your last question, why the spectrum have to flatten?, this is because I am throwing particles with a flat spectral index of -1 which I will re-weight later to a higher spectral index. And just like the particles detected at the source, I would expect particles detected at the observer to posses a flat spectrum for unweighted events and a power law spectrum for the re-weighted events irrespective of the energy. However this seems to be a factor of the magnetic field in use and i would like to try Dolag's Quimby file as well.

I don't think I have enough knowledge about Dolag field that anyone can count on. However, I am simulating cosmic ray particles from specific astrophysical sources from the universe and my interest is how these candidates spread around the source. The Dolag field seems to be too strong and deflects particles all over the place that is why I have turbulent field as well. Nevertheless, I am still not so sure if this is the right approach.

Finally you have suggested that the origin=(0.5size, 0.5size, 0.5*size) and that the labels are in Tesla. Which labels are you referring to here? I thought size is already given in Mpc.

Regards, Wilson.

[Laruku25]

Hey Wilson,

by labels I mean the numbers you see on the colorbars are in Tesla. The right field strenght (im pretty sure of it) is 1uG or 1e10 T at the red spots (clusters) just like the dolag paper from there everything will calibrate to the real values using the B factor.

I did not know the turbulent field will make the dolag field weaker I might try that. What Im trying now is to reduce the Bfactor to 1e-3 using this factor I can get a 5 degree deflection for a 40*EeV proton when it goes through a cluster just like Dolag's paper.

This means the field strenghts will be severely reduced and it may not make sense but well I'm finding that more useful than to have the Bfactor to 1e-1 or 1, it makes finding any correlation to sources impossible!

Oh by the way, by size i mean size=186*Mpc. I believed this is the length of the edge of the simulation box so setting the origin to x,y,z=186/2,186/2,186/2 Mpc will position your observer at the middle/origin of the box.

Now checking out CRPROPA's paper again it seems that the Dolag Field was changed to match the strengths of Miniati. The Dolag field file kept the original structure of Dolag which follows the distribution of matter bla bla bla but the strengths especially around clusters were changed to encompass bigger regions (bigger red spots bigger area with 1uG ).

Anyways, thanks for your explanation I understand why you want a flat spectrum. Why don't you try using 132Mpc (this is in crpropas paper) and set B=1e-3 and origin at s=132/2,132/2,132/2 Mpc and take out the turbulent field, how long does your simulation take? Can you easily experiment like that or is it too computionally intensive?

One thing that would really help me out is to know whether or not the origin here in the graphs I showed you represents our position in the local universe, I'm not an expert in astronomy so I wouldn't know I would guess the earth should be somewhere in the red spots. Once I have that information I will just run the simulation with b=1e-3 size=186/2 or 132/2 but thue position of earth in that map is important for me.

[wilson-2015]

Hi Laruku,

Thanks again for your detailed and helpful explanation. Yes, today I used your values, set B=1e-3 and commended out turbulent field. This gave me a flat spectrum that I was expecting. However, the magnetic field seems so reduced that the arrival direction of events are not visible and out of expectation. I tried B=1e-2, and seems to produce good results for my case.

I also wanted to visualize the Dolag magnetic field on an histogram. In the manual I used this syntax:

dumpGridToTxt(Bfield.getGrid(), 'myfield.txt'),

however the txt file generated is so large that it can not be used to plot the histogram by interactive python. In case you have some knowledge on how this is done I will be glad to know.

Thanks and regards, Wilson.

[Laruku25]

glad its working for you, i think it may not be fartched to use different bfactors since we only know the strength at galaxy clusters and cores and not really outside so deflections may not be as large as many suggest.

By visualizing a magnetic field on a histogram, what do you mean? Can you give me an example perhaps from another paper of the kind of graph you want? I can make the code for you np (its probably done anyways maybe i just have to make a few tweaks or maybe not depending on what you want).

Also based on the graph I posted, do you think that origin represents our galaxies location in the local unvierse or dolags map? Do you have an astrnomer friend you can ask maybe :P

[wilson-2015]

thanks once more. I will look around in case I find someone with more knowledge on your plots I will let you know.

The histogram I wanted looks like this dwl.pdf. The problem is the large file that I get when I dump grid to text. If there is a different way of plotting this I will really appreciate to know.

[Laruku25]

Try to get the values into an array using B.getField(Vector3d(a,b,c)) and make those vectors increase in small intervals as detailed as you want it to be, then use that array for your histogram.

The graphs I showed you represent the entire size/shape of dolag's field by the way. If it's 132Mpc for you then make the vectors change as low as 1kpc if you want.

lemme know if that helps you out.

[wilson-2015]

Hey Laruku!,

You will pardon me for my poor programming skills here. The way I understood is that starting from the centre of the box (0,0,0), I increment the vectors in small sizes (preferably 1kpc) up to 75 Mpc which is the full size of the box.

The quickest code I am trying think about doesn't work maybe some inputs will be helpful:

BField = MagneticFieldGrid(vgrid) a=b=c = np.linspace(0, 75, 0.001) Bfield = BField.getField(Vector3d(a,b,c) * Mpc) / nG with open('Bfield.txt', 'w') as outfile: for fname in Bfield: with open(fname) as infile: for line in infile: outfile.write(line)

[JonPaulLundquist]

I too appear unable to reproduce the results of arXiv:astro-ph/0410419. All of the important information is either missing or contradictory as already discussed. On another github issue thread it was suggested to use a b_factor = 109 nG. This gives an average field of 5x10^-6 microgauss which from figures like figure 12 seems to be about correct. Yet, I can't seem to find any object corresponding to Virgo, Perseus, or Coma with fields on the order of 1 microgauss. It seems that the conclusion is the b_factor is somewhere between 0.0001 to 1. The size is 132 Mpc to 184 Mpc. And the center is either at Vector3d(0,0,0) or Vector3d(size/2,size/2,size/2).

How are you making those magnetic field plots? Just B.getField(Vector3d(a,b,c))? My results so far don't really look like that.

[wilson-2015]

Hi! I hope the average field of 5x10^-6 microgauss, you are referring to Brms and not average mean because this will be too strong. With b_factor in that range I too get an average field of around 1 picogauss which is in agreement with the results from cumulative filling factors.

I could not plot the attached plot with B.getField(Vector3d(a,b,c)) an error is raised when I load in the values of a, b, c: """TypeError: Wrong number or type of arguments for overloaded function 'new_Vector3d'""

The solution that I am thinking about is how to simply convert a .raw file to a readable format.

[JonPaulLundquist]

A larger b_factor results in a larger average magnetic field and I'm using a smaller value. My initialization is:

obsPosition = Vector3d(0,0,0) gridSize = 440 size = 186*Mpc b_factor = 10*9 nG spacing = size/(gridSize) vgrid = Grid3f( obsPosition, gridSize, spacing ) loadGrid(vgrid, filename_bfield, b_factor) Bfield = MagneticFieldGrid(vgrid)

Your suggestion of b_factor = 10^-2 results in: B_ave = meanFieldStrength(vgrid) / nG = 0.4 nG = 5x10^-4 microgauss B_rms = rmsFieldStrength(vgrid) / nG = 118 nG = 0.118 microgauss

This is too large according to the voids shown in Figures 11 and 12 in the Dolag et al. paper which have ~10^-6 microgauss.

I have been using a smaller b_factor = 1 Gauss = 10^9 * nG = 10^-4 (100 times smaller than yours) that results in similarly 100 times smaller values:

B_ave = 5x10^-6 microgauss B_rms = 1.2x10^-3 microgauss.

[Laruku25]

Then I can match the maximum field strength and the minium field strength outlined by dolag's paper, the maximum is ~1uG which is in clusters.

Hye John, if you want to match Dolag's max and min field strength you can use B=0.1 , then you'll get the values of the first figure I posted.

10^-10 in clusters virgo persa etc which in gauss is 10^-10*10^4 (Tesla to Gauss conversion) = 10^-6 = 1micG

and you do the same for the lowest values. I've gotten all the field values in an array and then use the python max and min function and matched the min and max field strengths in clusters and voids.

The problem with this is that the deflections are too strong even for a proton, I believe this is because the file they have given us to download here does not represent the true field strength distribution of dolag, what i mean by this is that they have made the yellow/red spots wider and bigger causing bigger deflections and may have increased the field strengths around clusters (keeping the maximum and minimum). They did this because they believed that Dolag's field estimates were wrong and followed Miniati's advice there. In summary, that field represents a mix of Dolag and Miniati's fields. They've kept Dolag's matter distribution but they've added Miniatis field strengths estimates in it.

I've got this from a CRPROPA paper and what's what I understood from it.

@wilson-2015

Bro I've been busy with some stuff give me a few days and I'll help you with your histogram code.

Forgot to tell you guys that the quimby file apparently has the exact field distribution and strength of Dolag, unfortunately, the link no longer works. But I bet some people have it stored somewhere I do too, it's somwhere on an old laptop that no longer works so I need to extract all the files in the hard-disk(if it works) to get it.

[JonPaulLundquist]

Hello Laruku25, So the reference to https://arxiv.org/pdf/astro-ph/0410419.pdf for the Dolag field on the CRPropa website https://crpropa.github.io/CRPropa3/pages/howto_cite_crpropa.html is not correct? I've sent an email to Dolag and one of the CRPropa guys for information about this file.

A b_factor = 0.1 results in an average field of 5 nG which would indeed result in much larger deflections than Dolag et al (2005). The minimum field should be ~10^-8.5 microgauss and the average ~10^-6 according to figures 11 and 12 and the first paragraph of page 26. Just picking some random points I'm getting ~10^-3 microgauss.

You state you have all the field values in an array. Is that how you're making your figures? Did you fill that array using the getField() CRPropa function? Both the Bfield and vgrid of the code above aren't python arrays. How did you choose your points?

[Laruku25]

No, the reference isnt necessarily wrong. It does follow the structure and shape of the paper but the strengths have been slightly changed while keeping the minimum and maximum values intact. At least that's I believe. So naturally you will get different different average values for the fields so don't expect to match anything on dolag's paper except maybe the minium and max fields (since I think clusters and void values are much more easily measured/estimated than in other places).

The quimby dolag file I believe hasn't done any modification so it stays true to the paper.

I got all the values using the getField put them in array and chose the points to be separated by small intervals, the smaller the interval the more time it takes to load up the figure.

I did the same thing for the miniati file and was able to reproduce the exact same graph on miniati papers.

https://indico.cern.ch/event/504078/contributions/2329750/attachments/1353692/2044782/VanVliet_UHECR2016.pdf

The miniati paper did a graph like mine but they did not use the file given here they used a different seed to load up a different field, their magnetic fields don't follow the local universe like dolag does but the CRPROPA guys believe they have the right strengths so they used that and super imposed it on Dolag's file who does follow the distribution of matter in the local universe.

Well it's what I believe, only they know whats up.

[Laruku25]

https://forge.physik.rwth-aachen.de/public/quimby/mhd/ th quimby file was available through this link, maybe you can email or call the guys in charge of the repository, im sure they have it somewhere. https://git.rwth-aachen.de/3pia/forge/quimby

[wilson-2015]

Helo Laruku,

I hope you will get some time soon and look at my issue, I am really eager I want to know whether that histogram can be plotted with getField approach. Otherwise, I have been hoping to get an alternative way of reading a .raw file but I have just not succeeded yet.

[JonPaulLundquist]

Coming back to this - I believe this is how you can load the magnetic field magnitude into numpy vectors so you can analyze it with ease.

from crpropa import Mpc, nG, Vector3d, Grid3f, loadGrid, MagneticFieldGrid, meanFieldStrength, rmsFieldStrength, dumpGridToTxt
from numpy import arange, diff, zeros, fromstring, round, where, std, logical_and, log10
from numpy.linalg import norm
import fileinput
import matplotlib.pyplot as plt
from astropy.coordinates import SkyCoord, cartesian_to_spherical

filename_bfield = 'dolag_B_54-186Mpc_440b.raw'
gridSize = 440
size = 186*Mpc
b_factor = 10**9 * nG
spacing = size/(gridSize)
obsPosition = Vector3d(0,0,0)                       #this makes defining sources around the observer easier
boxOrigin = Vector3d(-0.5*size,-0.5*size,-0.5*size) #this makes defining sources around the observer easier

vgrid = Grid3f(boxOrigin, gridSize, spacing )
loadGrid(vgrid, filename_bfield, b_factor)
Bfield = MagneticFieldGrid(vgrid)
aveField = meanFieldStrength(vgrid) / (1000*nG)
rmsField = rmsFieldStrength(vgrid) / (1000*nG)

dumpGridToTxt(vgrid, 'DolagField.txt')

steps = arange(-(size/2)/Mpc,(size/2)/Mpc,(size/Mpc)/gridSize)
steps = steps + min(diff(steps))/2

x = zeros((steps.size**3), dtype='float16')
y = zeros((steps.size**3), dtype='float16')
z = zeros((steps.size**3), dtype='float16')
i = 0
for j in range(steps.size):
    for k in range(steps.size):
        for m in range(steps.size):
            x[i] = steps[j]
            y[i] = steps[k]
            z[i] = steps[m]
            i = i + 1

D = norm([x, y, z], axis = 0) #Distances from coordinates

Bmag = zeros(steps.size**3, dtype='float16')
i = 0
for line in fileinput.input(['myfield.txt']):
    Bmag[i] = norm(fromstring(line, sep=' '))/(1000*nG)
    i = i + 1

print(round(mean(Bmag)*10)/10==round(aveField*10)/10) #Prints TRUE
print(round(std(Bmag)*10)/10==round(rmsField*10)/10) #Prints TRUE

I believe this gives you the positions and magnitudes in microGauss without using an insane amount of memory. The average field is 4.7x10^-6 microGauss consistent with Dolag et al. The maximum field is 8.4 microGauss also consistent with Dolag et al.

Questions remain:

ind = where(Bmag>1) returns three positions with fields over 1 microGauss at 78 Mpc, 90 Mpc, and 132 Mpc. So if those are supposed to correspond to Virgo, Perseus, and Coma the distances are way off...

ind = logical_and(D>53.5, D<54.1) Around the distance to Virgo results in a max(Bmag[ind]) = 0.0004 microGauss which is too small...

Plotting like this:

supergalactic = cartesian_to_spherical(x[ind],y[ind],z[ind])
supergalactic = SkyCoord(sgl=supergalactic[2], sgb=supergalactic[1], frame='supergalactic', unit='rad')

fig = plt.figure(figsize=(8,6))
ax = fig.add_subplot(111, projection="aitoff")
s = ax.scatter(supergalactic.sgl.wrap_at('180d').rad, supergalactic.sgb.rad,c=np.log10(Bmag[ind])-6,s=0.15)
ax.set_xticklabels(['','240','','300','','0','','60','','120',''], fontweight='semibold')
ax.set_yticklabels(['','-60','','-30','','','','30','','60',''], fontweight='semibold')
ax.grid(True)
cbar = fig.colorbar(s, orientation='vertical')
cbar.set_label(label='log$_{10}$(Gauss)', weight='bold')
plt.savefig('MapDolagField_54Mpc.png')

Looks like this:

[wilson-2015]

Thanks for the great code. I have followed it and printed TRUE, TRUE meaning that both the aveField and rmsField as per Bmag are equivalent to that of the vgrid. I now expected that plotting Bmag in a histogram should give me a nice plot like this new.pdf but this is not what I am obtaining. Is there something maybe I am missing out? Off course when I print(Bmag) give a vector with zeros.

Questions remains:-maybe Laruku at one point will comment and shed more light on this. Which data have you used for "ind" to generate the sky map?

[JonPaulLundquist]

The indices were picking out a distance. ind = logical_and(D>53.5, D<54.1). This was the wrong distance, not sure what I was thinking... I believe the b_factor may also be in the name of the file: B = 54 Gauss.

filename_bfield = 'dolag_B_54-186Mpc_440b.raw'
gridSize = 440
size = 186*Mpc
b_factor = 54*10**9 * nG
spacing = size/(gridSize)
boxOrigin = Vector3d(-0.5*size,-0.5*size,-0.5*size) #this makes defining sources around the observer easier

vgrid = Grid3f(boxOrigin, gridSize, spacing )
loadGrid(vgrid, filename_bfield, b_factor)
Bfield = MagneticFieldGrid(vgrid)

dumpGridToTxt(vgrid, 'DolagField.txt')

Bmag = zeros(steps.size**3, dtype='float32')
i = 0
for line in fileinput.input(['DolagField.txt']):
    Bmag[i] = norm(fromstring(line, sep=' '))/(1000*nG)
    i = i + 1

This gives an average Bmag = 2.5x10^-4 microGauss which appears consistent with your linked histogram but my minimum is smaller and maximum larger. Doing the histogram like so:

plt.hist(np.log10(Bmag*10**-6), bins=np.arange(-12,-6,0.001),density=True, log=True)

Looks like this:

Very similar but has more values at the high and low ranges. I'm guessing they sampled using the CRpropa getField() function. What is the paper you link to?

Changing from float16 to float32 will get rid of zeros in Bmag if your b_factor is small. This b_factor = 56 Gauss results in a maximum field of 23.7 microGauss roughly at the location of the Virgo Supercluster at a distance of 27.62 Mpc. That's about 5 Mpc off from Dolag et al. Figure 26 but the correct magnitude.... The plot looks like this:

I can find nothing that seems to correspond to Perseus. I've tried assuming supergalactic, galactic, and equatorial coordinates and flipping my "steps" vector.

[Laruku25]

Helo Laruku,

I hope you will get some time soon and look at my issue, I am really eager I want to know whether that histogram can be plotted with getField approach. Otherwise, I have been hoping to get an alternative way of reading a .raw file but I have just not succeeded yet.

Hey Wilson Im sorry I cant get back to you. Im about to travel to europe and kind of scared tbh getting my papers ready and this covid stuff too. But I see that you've gotten the help you needed from John. In case you need to graph fields here's the code I used. I am not too versed in programming either I only know python and only basic calculations:

` from crpropa import* import numpy as np import math import matplotlib.pyplot as plt import matplotlib.cm as cm from matplotlib.ticker import MaxNLocator from matplotlib.colors import BoundaryNorm from matplotlib.ticker import LogFormatter from matplotlib.colors import LogNorm

def loadDolagField(): file="Dolag.raw" size=132Mpc gridOrigin=Vector3d(0.5size,0.5size,0.5size) gridSize=440 b_factor=1e-1 spacing=size/(gridSize) vgrid=VectorGrid(gridOrigin,gridSize,spacing) loadGrid(vgrid,file,b_factor) EGMF=MagneticFieldGrid(vgrid) return EGMF

def loadMiniatiField(): file="Miniati.raw" size=50*Mpc/.67

gridOrigin=Vector3d(0,0,0)

gridOrigin=Vector3d(0.5*size,0.5*size,0.5*size)
gridSize=256
b_factor=1e-4
spacing=size/(gridSize)
vgrid=VectorGrid(gridOrigin,gridSize,spacing)
loadGrid(vgrid,file,b_factor)
EGMF=MagneticFieldGrid(vgrid)
return EGMF

def getValues(x,y,z,B): z=z R=B.getField(Vector3d(x,y,z)) R=R.getX()2+R.getY()2+R.getZ()**2 R=math.sqrt(R) return R

def getXYZ(min,max,interval,z): M=loadDolagField() x=np.arange(min,max,interval) y=np.arange(min,max,interval) x_center=0.5(x[:-1]+x[1:]) y_center=0.5(y[:-1]+y[1:]) X,Y=np.meshgrid(x_center,y_center) Z = [[None for p in range(len(x)-1) ] for r in range(len(x)-1) ] for i,val in enumerate(Z): for j,bol in enumerate(Z): Z[i][j]=getValues(X[i][j],Y[i][j],z,M) Z=np.array(Z) x=x/Mpc y=y/Mpc return x,y,Z

def GenerateTracks(particle,energy,direction,EGField,distance): position=Vector3d(0,0,0) if particle=='Proton': particle=-nucleusId(1,1) if particle=='Helium': particle=-nucleusId(2,1) if particle=='Oxygen': particle=-nucleusId(16,8) if particle=='Iron': particle=-nucleusId(56,26) if particle=='Silicon': particle=-nucleusId(28,14) if EGField=='Dolag': EGMF=loadDolagField() if EGField=='Miniati': EGMF=loadMiniatiField() p=ParticleState(particle,energy,position,direction) sim=ModuleList() sim.add(PropagationCK(EGMF,1e-4,60kpc,800kpc)) print('ok') obs=Observer() obs.add(ObserverLargeSphere(Vector3d(0),float(distance*Mpc))) sim.add(obs) output = TextOutput(EGField+'.txt', Output.Trajectory3D) sim.add(output) c=Candidate(p) sim.run(c) print c

GenerateTracks('Proton',40EeV,Vector3d(1,1,0),'Dolag',30) x,y,Z=getXYZ(-30Mpc,30Mpc,.1Mpc,132/.2*Mpc) data=np.genfromtxt('Dolag.txt',names=True) a,b = data['X'],data['Y'] fig=plt.figure(figsize=(12,5)) cmap = plt.get_cmap('rainbow') plot=plt.pcolormesh(x,y,Z,cmap=cmap,norm=LogNorm()) plt.colorbar(plot) plt.plot(a,b) plt.show() '

@JonPaulLundquist

You may want to read this paper. https://arxiv.org/pdf/1603.07142.pdf . Section 4.2 talks about the extragalactic fields and it may be talking about this file in particular.

I will get back to this stuff later but thanks for your awesome code, it's definitely going to help us get more acquianted with astrophysical graphs in python.

[JonPaulLundquist]

I don't believe the field in https://arxiv.org/pdf/1603.07142.pdf is the same file. That covers (132 Mpc)^3 volume (not 186 Mpc) and is a grid of 256^3 points (not 440).

filename_bfield = 'dolag_B_54-186Mpc_440b.raw'
gridSize = 256
size = 132*Mpc
b_factor = 54*10**9 * nG
spacing = size/(gridSize)
obsPosition = Vector3d(0,0,0)                       #this makes defining sources around the observer easier
boxOrigin = Vector3d(-0.5*size,-0.5*size,-0.5*size) #this makes defining sources around the observer easier

vgrid = Grid3f(boxOrigin, gridSize, spacing )
loadGrid(vgrid, filename_bfield, b_factor)

Results in "RuntimeError: loadGrid: file and grid size do not match" because the file has 440^3 grid points.

Also, I notice that you're defining the corner of your field box to be gridOrigin=Vector3d(0.5size,0.5size,0.5size) but your observer is at position=Vector3d(0,0,0). This is the opposite of the recommendation here: https://crpropa.github.io/CRPropa3/pages/example_notebooks/extragalactic_fields/MHD_models.v4.html and means that particles won't be travelling through the field the whole time. Personally, I think my obsPosition and boxOrigin makes more sense. You also probably want to add sim.add( PeriodicBox( boxOrigin, boxSize ) ).

[JonPaulLundquist]

I just started a new report on this since this seems to be getting overlooked. https://github.com/CRPropa/CRPropa3/issues/379

I think that:

gridSize = 440
size = 132*Mpc
b_factor = 100*10**9 * nG

And the file name and the reference on the CRPropa site should be changed.

[wilson-2015]

Hello Paul!, First thanks for simplifying the code for the histogram and sorry for a late reply. I was still going through and now happy that everything works well.

Its great that you have opened a new ticket on Dolag field, actually a lot needs clarification. As you mentioned earlier, I made some adjustment in my script on boxOrigin and obsPosition to reflect parameters in your code.

Bro Laruku thanks for your code yesterday. I will go through and in case of any issue I will get in touch.

Thanks.

[Laruku25]

I don't believe the field in https://arxiv.org/pdf/1603.07142.pdf is the same file. That covers (132 Mpc)^3 volume (not 186 Mpc) and is a grid of 256^3 points (not 440).

filename_bfield = 'dolag_B_54-186Mpc_440b.raw'
gridSize = 256
size = 132*Mpc
b_factor = 54*10**9 * nG
spacing = size/(gridSize)
obsPosition = Vector3d(0,0,0)                       #this makes defining sources around the observer easier
boxOrigin = Vector3d(-0.5*size,-0.5*size,-0.5*size) #this makes defining sources around the observer easier

vgrid = Grid3f(boxOrigin, gridSize, spacing )
loadGrid(vgrid, filename_bfield, b_factor)

Results in "RuntimeError: loadGrid: file and grid size do not match" because the file has 440^3 grid points.

Also, I notice that you're defining the corner of your field box to be gridOrigin=Vector3d(0.5size,0.5size,0.5size) but your observer is at position=Vector3d(0,0,0). This is the opposite of the recommendation here: https://crpropa.github.io/CRPropa3/pages/example_notebooks/extragalactic_fields/MHD_models.v4.html and means that particles won't be travelling through the field the whole time. Personally, I think my obsPosition and boxOrigin makes more sense. You also probably want to add sim.add( PeriodicBox( boxOrigin, boxSize ) ).

Hey John I have doubts myself yes the 256 doesn't work only 440b works but i believe the magnetic field structure and shapes are from that paper, if they were from dolag's paper, you'd get the small deflection angles with heavy nuclei and and protons of 40*EeV the only way to get those is to get the B_field factor you've posted or something lower like that but if you put that B_factor you don't get the minimum and maximum field strengths, do note that CRPROPA uses Tesla instead of Gauss too. I've graphed the dolag's field shape and a trajectory and you can see that with an observer position at 0,0,0 as long as you center the box in the origin, i get the trajcetory to start in the middle.

Also about the repeating box function, i think crpropa does that by default, i will use this same code to span a bigger region and you'll see it automatically repeats the structure.

I've seen that you've wrote a new post i will read it later.

You're making good progress John, this is helping me a lot.

[JonPaulLundquist]

I've graphed the dolag's field shape and a trajectory and you can see that with an observer position at 0,0,0 as long as you center the box in the origin, i get the trajcetory to start in the middle.

The center of the box is not the gridOrigin/boxOrigin parameter though. It's apparently the corner. See the CRPropa example: https://crpropa.github.io/CRPropa3/pages/example_notebooks/extragalactic_fields/MHD_models.v4.html Your posted code has gridOrigin>obsPosition in x, y, and z. So I don't see how the observer would be within the box...

If you want to simulate a trajectory to the Milky Way through the local universe with this simulated magnetic field you need to do something like the given example:

gridOrigin = Vector3d(0,0,0)             ## origin of the 3D data, preferably at boxOrigin
vgrid = Grid3f( gridOrigin, gridSize, spacing )

obsPosition = Vector3d(0.5*size,0.5*size,0.5*size) # position of observer, MW is in center of constrained simulations
obs.add( ObserverSmallSphere( obsPosition, obsSize ) )

That's straight from the CRPropa website. gridOrigin < obsPosition. Or something like what I'm doing:

gridOrigin = Vector3d(-0.5*size, -0.5*size, -0.5*size)             ## origin of the 3D data, preferably at boxOrigin
vgrid = Grid3f( gridOrigin, gridSize, spacing )

obsPosition = Vector3d(0, 0, 0) # position of observer, MW is in center of constrained simulations
obs.add( ObserverSmallSphere( obsPosition, obsSize ) )

Also, gridOrigin < obsPosition. This way I can simply define my isotopically emitting sources by some distance from (0,0,0) and they will be within the local universe field:

direction = Vector3d()
direction.setRThetaPhi(D, theta, phi)
s.add(SourcePosition(direction))

I don't think my version of CRPropa automatically uses a periodic box as it's in the website example. I also thought it was done automatically and I wasted a lot of time not using it. Now my observer efficiency has increased dramatically due to more magnetic deflection.

As for the field magnitude - looking at all grid points directly is more accurate instead of using B.getField() (which uses interpolation so your values will always be less then the local maximums unless you're really lucky and exactly select the right point). Since "54" in the file name evidently doesn't have any relevance I think something like b_factor = 10*10*9 nG works best to approximately match Dolag. You get local maximum fields on the order of ~5 microG for both Virgo and Perseus which is roughly what Dolag has in his paper. The maximum you get is 85 microG at (194.84902946, -25.51684215) galactic coords. at 64 Mpc which is pretty large. Not sure what that is supposed to correspond to. The average is ~5x10^-5 microGauss.

[JonPaulLundquist]

I know this thread has been dead awhile but I just made a cumulative fraction plot similar to Dolag et. al middle figure 1 (https://iopscience.iop.org/article/10.1088/1475-7516/2005/01/009/pdf) with the b_factor = 10*10*9 nG. It looks very similar but there's a lower fraction of locations with B < ~10^-6 Gauss. Their plots aren't the whole field though.

[lukasmerten]

Thanks @wilson-2015, @Laruku25, and @JonPaulLundquist for your lively discussions during the last weeks. Unfortunately, this issue has reached a point were to many things are mixed that a clear resolution of each of the points seems not possible for now. Therefore it will be closed. However, we would like to encourage you to continue your discussions via email.

Please note that the CRPropa issues are not meant for discussion of post-processing questions how to plot the magnetic field?) or interpretation of physics results. That’s not because these questions are not relevant or interesting but simply because the group of developers is too small to provide any support in this direction.

Any of the provided tabulated magnetic fields are more or less external to CRPropa. We try to provide them to our best knowledge, but we cannot take any responsibility for any of the physical parameters of these models.

Please note: What was used in the CRPropa 3 paper was a model combining the Miniati magnetic field--density relation with the Dolag distribution, as explained in the paper. This model was meant to be illustrative and has nothing to do with a preference for one magnetic field or another, which is actually not well understood.