lukasmerten
commented
10 months ago Hi @AmadoraMillalen I finally found some time to take a look at your issue.
The ToroidalHaloField has default settings that are very likely not working for Galactic scale simulations:
ToroidalHaloField(double b0 = 1., double z0 = 1., double z1 = 1., double r0 = 1.) {
setParameters(b0, z0, z1, r0);
}If you rerun your simulation use some values for b0, z0, z1, and r0 that are better fitting positions on Galactic scales (kpc). Since I am not the author of that magnetic field model, I cannot suggest any concrete values, but the default values just lead a magnetic field strength of |B|=0.
The LogarithmicSpiralField actually had a bug in the part that models the logarithmic decay with Galactic height. This lead to a very large magentic field strength outside of the Galactic plane a probably trapped all particles during propagation.
This should be fixed with #459. However, please check the settings of the default parameters carefully. To me the magnetic field strength in the Galactic plane (z=0) seems to be very large.
Hi, I was doing magnification plots using: sim.add(PropagationCK(BField(), 1e-4, 0.1 parsec, 10 parsec))
With JF12Field, PT11Field, PT11Field, CMZField, TF17Field I got deflections, but with ToroidalHaloField and LogarithmicSpiralField I don't get any deflection,
... for BField, ax in zip(BFieldList, axes): energy = E * EeV
sim = ModuleList() sim.add(PropagationCK(BField(), 1e-4, 0.1 parsec, 10 parsec)) # def using Runge-Kutta obs = Observer() obs.add(ObserverSurface( Sphere(Vector3d(0), 20* kpc) )) ...
pid = - nucleusId(A,Z) position = Vector3d(-8.5, 0, 0) * kpc # Earth
lb_final = np.empty((npix, 2)) print('simulating {} particles'.format(npix))
for i in range(npix): v = healpy.pix2vec(nside,i) direction = Vector3d(v[0], v[1], v[2]) p = ParticleState(pid, energy, position, direction) c = CandidateRefPtr(Candidate(p)) cv.push_back(c)
sim.run(cv) for j in range(len(cv)): c = cv[j] d = c.current.getDirection() lb_final[j , 0] = d.getPhi() lb_final[j , 1] = d.getTheta() l_comp = [] b_comp = [] for i in range(npix): l_comp.append( np.nan_to_num(lb_final[ i, 0])) b_comp.append( np.nan_to_num(np.pi/2 - lb_final[ i , 1])) lb_final = SkyCoord(l_comp,b_comp,frame='galactic',unit='radian').galactic
pix_final = hpx.skycoord_to_healpix(lb_final) extgal = healpy.ma(np.zeros(npix)) for i in pix_final: if extgal[i] == 0: extgal[i] = 1 else: extgal[i] +=1 if (extgal == 0).any(): extgal[extgal == 0] = 0.02 if counter==4: title = f'{type(BField()).name}ASS' else: title = f'{type(BField()).name}' plt.axes(ax)
map_out = healpy.sphtfunc.smoothing(np.log10(extgal), fwhm = 0.017, iter = 1)
healpy.mollview(map_out, coord='G', title=title ,hold=True, cmap= 'jet' )
counter+=1