Example to calculate photon-ALP oscillations from NGC 1275 in a cavity field
This notebook demonstrates how to calculate the photon-ALP transition probability for NGC 1275, the central AGN of the Perseus cluster. The assumed B-field environments are the same as in https://arxiv.org/abs/1908.03084, and include the cluster field and the magnetic field of the Milky Way. Electron density model parameters taken from https://iopscience.iop.org/article/10.1086/374923/fulltext/57155.text.html.
Note that in order to obtain the field used in 1908.03084, you have to set \(\theta\) to 225° instead of 45° (adding 180°). This effectively inverts the sign of all field components. While this change of \(\theta\) and the sign is irrelevant for the conversion probability, it changes the field along the line of sight and inverts the sign of the RM, as well.
If you haven’t installed gammaALPs already, you can do so using pip. Just uncomment the line below:
[1]:
#!pip install gammaALPs
[2]:
from gammaALPs.core import Source, ALP, ModuleList
from gammaALPs.base import environs, transfer
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.gridspec as gridspec
from ebltable.tau_from_model import OptDepth
[3]:
%matplotlib inline
Set the ALP
Initialize an ALP object, that stores the ALP mass \(m\) (in neV) and the coupling \(g\) (in \(10^{-11}\mathrm{GeV}^{-1}\)).
[4]:
m, g = 2, 100
alp = ALP(m,g)
Set the source
Set the source properties (redshift and sky coordinates) in the Source containier
[5]:
ngc1275 = Source(z = 0.017559, ra = '03h19m48.1s', dec = '+41d30m42s')
Init the module list
Initialize the list of transfer modules that will store the different magnetic field environments.
Energies are supplied in GeV as numpy.ndarray
[6]:
logEMeV = np.linspace(0.46709351, 6.53207135, num=183, endpoint=True)
EGeV = np.power(10, logEMeV - 3.)
Now initialize the initial photon polarization. Since we are dealing with a gamma-ray source, no ALPs are initially present in the beam (third diagonal element is zero). The polarization density matrix is normalized such that its trace is equal to one, \(\mathrm{Tr}(\rho_\mathrm{in}) = 1\).
[7]:
pin = np.diag((1.,1.,0.)) * 0.5
Add modules:
Now we add propagation modules for the cluster, the EBL, and the Galactic magnetic field; run the module. Repeat for only the ICM field.
[8]:
alp = ALP(m, g)
ml = ModuleList(alp, ngc1275, pin = pin, EGeV = EGeV)
ml.add_propagation("ICMStructured",
0, # position of module counted from the source.
B0 = 8.3, # rms of B field
R = 93,
theta = 225,
theta_rad=False,
pa=147,
pa_rad=False,
n0 = 3.9e-2, # normalization of electron density
n2 = 4.05e-3, # second normalization of electron density, see Churazov et al. 2003, Eq. 4
r_abell = 500., # extension of the cluster
r_core = 80., # electron density parameter, see Churazov et al. 2003, Eq. 4
r_core2 = 280., # electron density parameter, see Churazov et al. 2003, Eq. 4
beta = 1.2, # electron density parameter, see Churazov et al. 2003, Eq. 4
beta2= 0.58, # electron density parameter, see Churazov et al. 2003, Eq. 4
eta = 0 # scaling of B-field with electron denstiy
)
ml.add_propagation("EBL",1, eblmodel = 'dominguez') # EBL attenuation comes second, after beam has left cluster
ml.add_propagation("GMF",2, model = 'pshirkov', model_sym = 'BSS') # finally, the beam enters the Milky Way Field
px, py, pa = ml.run()
# again with only structured field
ml2 = ModuleList(alp, ngc1275, pin = pin, EGeV = EGeV)
ml2.add_propagation("ICMStructured",
0, # position of module counted from the source.
B0 = 8.3, # rms of B field
R = 93,
theta = 225,
theta_rad=False,
pa=147,
pa_rad=False,
n0 = 3.9e-2, # normalization of electron density
n2 = 4.05e-3, # second normalization of electron density, see Churazov et al. 2003, Eq. 4
r_abell = 500., # extension of the cluster
r_core = 80., # electron density parameter, see Churazov et al. 2003, Eq. 4
r_core2 = 280., # electron density parameter, see Churazov et al. 2003, Eq. 4
beta = 1.2, # electron density parameter, see Churazov et al. 2003, Eq. 4
beta2= 0.58, # electron density parameter, see Churazov et al. 2003, Eq. 4
eta = 0 # scaling of B-field with electron denstiy
)
px2, py2, pa2 = ml2.run()
environs.py: 578 --- INFO: Using inputted chi
environs.py:1196 --- INFO: Using inputted chi
core.py: 658 --- INFO: Running Module 0: <class 'gammaALPs.base.environs.MixICMStructured'>
core.py: 658 --- INFO: Running Module 2: <class 'gammaALPs.base.environs.MixGMF'>
/Users/mey/Python/gammaALPs/gammaALPs/base/transfer.py:799: UserWarning: Not all values of linear polarization are real values!
warnings.warn("Not all values of linear polarization are real values!")
/Users/mey/Python/gammaALPs/gammaALPs/base/transfer.py:802: UserWarning: Not all values of circular polarization are real values!
warnings.warn("Not all values of circular polarization are real values!")
environs.py: 578 --- INFO: Using inputted chi
core.py: 658 --- INFO: Running Module 0: <class 'gammaALPs.base.environs.MixICMStructured'>
Plot the output
[9]:
plt.plot(EGeV, px[0] + py[0], label='ICM + EBL + GMF')
plt.plot(EGeV, px2[0] + py2[0], label='ICM')
plt.xscale('log')
plt.yscale('linear')
plt.ylim(0.4, 1.05)
plt.xlim(3e-2, 1e2)
plt.xlabel('Energy (GeV)')
plt.ylabel('Photon survival probability')
plt.legend(loc = 0, fontsize = 'medium')
plt.annotate(r'$m_a = {0:.1f}\,\mathrm{{neV}}, g_{{a\gamma}} = {1:.1f} \times 10^{{-11}}\,\mathrm{{GeV}}^{{-1}}$'.format(ml.alp.m,ml.alp.g),
xy = (0.95,0.1), size = 'x-large', xycoords = 'axes fraction', ha = 'right')
plt.subplots_adjust(left = 0.2)
Plot the B fields
Plot r vs B field components of structured field.
[10]:
r = ml.modules[0].r
b_r = ml.modules[0]._Bfield_model.b_r
b_theta = ml.modules[0]._Bfield_model.b_theta
b_phi = ml.modules[0]._Bfield_model.b_phi
plt.plot(r, b_r, label='$B_r$')
plt.plot(r, b_theta, label=r'$B_\theta$')
plt.plot(r, b_phi, label=r'$B_\phi$')
plt.xlabel(r'$r [\SI{}{\kilo\parsec}]$')
plt.ylabel(r'$B [\SI{}{\micro\gauss}]$')
plt.legend()
[10]:
<matplotlib.legend.Legend at 0x13eeaab50>
Error in callback <function _draw_all_if_interactive at 0x1396e8af0> (for post_execute):
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/pyplot.py:197, in _draw_all_if_interactive()
195 def _draw_all_if_interactive() -> None:
196 if matplotlib.is_interactive():
--> 197 draw_all()
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/_pylab_helpers.py:132, in Gcf.draw_all(cls, force)
130 for manager in cls.get_all_fig_managers():
131 if force or manager.canvas.figure.stale:
--> 132 manager.canvas.draw_idle()
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/backend_bases.py:1893, in FigureCanvasBase.draw_idle(self, *args, **kwargs)
1891 if not self._is_idle_drawing:
1892 with self._idle_draw_cntx():
-> 1893 self.draw(*args, **kwargs)
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/backends/backend_agg.py:388, in FigureCanvasAgg.draw(self)
385 # Acquire a lock on the shared font cache.
386 with (self.toolbar._wait_cursor_for_draw_cm() if self.toolbar
387 else nullcontext()):
--> 388 self.figure.draw(self.renderer)
389 # A GUI class may be need to update a window using this draw, so
390 # don't forget to call the superclass.
391 super().draw()
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/artist.py:95, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
93 @wraps(draw)
94 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 95 result = draw(artist, renderer, *args, **kwargs)
96 if renderer._rasterizing:
97 renderer.stop_rasterizing()
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
69 if artist.get_agg_filter() is not None:
70 renderer.start_filter()
---> 72 return draw(artist, renderer)
73 finally:
74 if artist.get_agg_filter() is not None:
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/figure.py:3154, in Figure.draw(self, renderer)
3151 # ValueError can occur when resizing a window.
3153 self.patch.draw(renderer)
-> 3154 mimage._draw_list_compositing_images(
3155 renderer, self, artists, self.suppressComposite)
3157 for sfig in self.subfigs:
3158 sfig.draw(renderer)
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/image.py:132, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
130 if not_composite or not has_images:
131 for a in artists:
--> 132 a.draw(renderer)
133 else:
134 # Composite any adjacent images together
135 image_group = []
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
69 if artist.get_agg_filter() is not None:
70 renderer.start_filter()
---> 72 return draw(artist, renderer)
73 finally:
74 if artist.get_agg_filter() is not None:
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/axes/_base.py:3070, in _AxesBase.draw(self, renderer)
3067 if artists_rasterized:
3068 _draw_rasterized(self.figure, artists_rasterized, renderer)
-> 3070 mimage._draw_list_compositing_images(
3071 renderer, self, artists, self.figure.suppressComposite)
3073 renderer.close_group('axes')
3074 self.stale = False
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/image.py:132, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
130 if not_composite or not has_images:
131 for a in artists:
--> 132 a.draw(renderer)
133 else:
134 # Composite any adjacent images together
135 image_group = []
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
69 if artist.get_agg_filter() is not None:
70 renderer.start_filter()
---> 72 return draw(artist, renderer)
73 finally:
74 if artist.get_agg_filter() is not None:
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/axis.py:1395, in Axis.draw(self, renderer, *args, **kwargs)
1393 # Shift label away from axes to avoid overlapping ticklabels.
1394 self._update_label_position(renderer)
-> 1395 self.label.draw(renderer)
1397 self._update_offset_text_position(tlb1, tlb2)
1398 self.offsetText.set_text(self.major.formatter.get_offset())
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
69 if artist.get_agg_filter() is not None:
70 renderer.start_filter()
---> 72 return draw(artist, renderer)
73 finally:
74 if artist.get_agg_filter() is not None:
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/text.py:748, in Text.draw(self, renderer)
745 renderer.open_group('text', self.get_gid())
747 with self._cm_set(text=self._get_wrapped_text()):
--> 748 bbox, info, descent = self._get_layout(renderer)
749 trans = self.get_transform()
751 # don't use self.get_position here, which refers to text
752 # position in Text:
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/text.py:381, in Text._get_layout(self, renderer)
379 clean_line, ismath = self._preprocess_math(line)
380 if clean_line:
--> 381 w, h, d = _get_text_metrics_with_cache(
382 renderer, clean_line, self._fontproperties,
383 ismath=ismath, dpi=self.figure.dpi)
384 else:
385 w = h = d = 0
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/text.py:69, in _get_text_metrics_with_cache(renderer, text, fontprop, ismath, dpi)
66 """Call ``renderer.get_text_width_height_descent``, caching the results."""
67 # Cached based on a copy of fontprop so that later in-place mutations of
68 # the passed-in argument do not mess up the cache.
---> 69 return _get_text_metrics_with_cache_impl(
70 weakref.ref(renderer), text, fontprop.copy(), ismath, dpi)
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/text.py:77, in _get_text_metrics_with_cache_impl(renderer_ref, text, fontprop, ismath, dpi)
73 @functools.lru_cache(4096)
74 def _get_text_metrics_with_cache_impl(
75 renderer_ref, text, fontprop, ismath, dpi):
76 # dpi is unused, but participates in cache invalidation (via the renderer).
---> 77 return renderer_ref().get_text_width_height_descent(text, fontprop, ismath)
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/backends/backend_agg.py:217, in RendererAgg.get_text_width_height_descent(self, s, prop, ismath)
213 return super().get_text_width_height_descent(s, prop, ismath)
215 if ismath:
216 ox, oy, width, height, descent, font_image = \
--> 217 self.mathtext_parser.parse(s, self.dpi, prop)
218 return width, height, descent
220 font = self._prepare_font(prop)
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/mathtext.py:79, in MathTextParser.parse(self, s, dpi, prop, antialiased)
77 prop = prop.copy() if prop is not None else None
78 antialiased = mpl._val_or_rc(antialiased, 'text.antialiased')
---> 79 return self._parse_cached(s, dpi, prop, antialiased)
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/mathtext.py:100, in MathTextParser._parse_cached(self, s, dpi, prop, antialiased)
97 if self._parser is None: # Cache the parser globally.
98 self.__class__._parser = _mathtext.Parser()
--> 100 box = self._parser.parse(s, fontset, fontsize, dpi)
101 output = _mathtext.ship(box)
102 if self._output_type == "vector":
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/_mathtext.py:2165, in Parser.parse(self, s, fonts_object, fontsize, dpi)
2162 result = self._expression.parseString(s)
2163 except ParseBaseException as err:
2164 # explain becomes a plain method on pyparsing 3 (err.explain(0)).
-> 2165 raise ValueError("\n" + ParseException.explain(err, 0)) from None
2166 self._state_stack = []
2167 self._in_subscript_or_superscript = False
ValueError:
r [\SI{}{\kilo\parsec}]
^
ParseFatalException: Unknown symbol: \SI, found '\' (at char 3), (line:1, col:4)
---------------------------------------------------------------------------
ValueError Traceback (most recent call last)
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/IPython/core/formatters.py:340, in BaseFormatter.__call__(self, obj)
338 pass
339 else:
--> 340 return printer(obj)
341 # Finally look for special method names
342 method = get_real_method(obj, self.print_method)
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/IPython/core/pylabtools.py:152, in print_figure(fig, fmt, bbox_inches, base64, **kwargs)
149 from matplotlib.backend_bases import FigureCanvasBase
150 FigureCanvasBase(fig)
--> 152 fig.canvas.print_figure(bytes_io, **kw)
153 data = bytes_io.getvalue()
154 if fmt == 'svg':
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/backend_bases.py:2158, in FigureCanvasBase.print_figure(self, filename, dpi, facecolor, edgecolor, orientation, format, bbox_inches, pad_inches, bbox_extra_artists, backend, **kwargs)
2155 # we do this instead of `self.figure.draw_without_rendering`
2156 # so that we can inject the orientation
2157 with getattr(renderer, "_draw_disabled", nullcontext)():
-> 2158 self.figure.draw(renderer)
2159 if bbox_inches:
2160 if bbox_inches == "tight":
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/artist.py:95, in _finalize_rasterization.<locals>.draw_wrapper(artist, renderer, *args, **kwargs)
93 @wraps(draw)
94 def draw_wrapper(artist, renderer, *args, **kwargs):
---> 95 result = draw(artist, renderer, *args, **kwargs)
96 if renderer._rasterizing:
97 renderer.stop_rasterizing()
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
69 if artist.get_agg_filter() is not None:
70 renderer.start_filter()
---> 72 return draw(artist, renderer)
73 finally:
74 if artist.get_agg_filter() is not None:
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/figure.py:3154, in Figure.draw(self, renderer)
3151 # ValueError can occur when resizing a window.
3153 self.patch.draw(renderer)
-> 3154 mimage._draw_list_compositing_images(
3155 renderer, self, artists, self.suppressComposite)
3157 for sfig in self.subfigs:
3158 sfig.draw(renderer)
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/image.py:132, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
130 if not_composite or not has_images:
131 for a in artists:
--> 132 a.draw(renderer)
133 else:
134 # Composite any adjacent images together
135 image_group = []
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
69 if artist.get_agg_filter() is not None:
70 renderer.start_filter()
---> 72 return draw(artist, renderer)
73 finally:
74 if artist.get_agg_filter() is not None:
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/axes/_base.py:3070, in _AxesBase.draw(self, renderer)
3067 if artists_rasterized:
3068 _draw_rasterized(self.figure, artists_rasterized, renderer)
-> 3070 mimage._draw_list_compositing_images(
3071 renderer, self, artists, self.figure.suppressComposite)
3073 renderer.close_group('axes')
3074 self.stale = False
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/image.py:132, in _draw_list_compositing_images(renderer, parent, artists, suppress_composite)
130 if not_composite or not has_images:
131 for a in artists:
--> 132 a.draw(renderer)
133 else:
134 # Composite any adjacent images together
135 image_group = []
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
69 if artist.get_agg_filter() is not None:
70 renderer.start_filter()
---> 72 return draw(artist, renderer)
73 finally:
74 if artist.get_agg_filter() is not None:
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/axis.py:1395, in Axis.draw(self, renderer, *args, **kwargs)
1393 # Shift label away from axes to avoid overlapping ticklabels.
1394 self._update_label_position(renderer)
-> 1395 self.label.draw(renderer)
1397 self._update_offset_text_position(tlb1, tlb2)
1398 self.offsetText.set_text(self.major.formatter.get_offset())
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/artist.py:72, in allow_rasterization.<locals>.draw_wrapper(artist, renderer)
69 if artist.get_agg_filter() is not None:
70 renderer.start_filter()
---> 72 return draw(artist, renderer)
73 finally:
74 if artist.get_agg_filter() is not None:
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/text.py:748, in Text.draw(self, renderer)
745 renderer.open_group('text', self.get_gid())
747 with self._cm_set(text=self._get_wrapped_text()):
--> 748 bbox, info, descent = self._get_layout(renderer)
749 trans = self.get_transform()
751 # don't use self.get_position here, which refers to text
752 # position in Text:
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/text.py:381, in Text._get_layout(self, renderer)
379 clean_line, ismath = self._preprocess_math(line)
380 if clean_line:
--> 381 w, h, d = _get_text_metrics_with_cache(
382 renderer, clean_line, self._fontproperties,
383 ismath=ismath, dpi=self.figure.dpi)
384 else:
385 w = h = d = 0
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/text.py:69, in _get_text_metrics_with_cache(renderer, text, fontprop, ismath, dpi)
66 """Call ``renderer.get_text_width_height_descent``, caching the results."""
67 # Cached based on a copy of fontprop so that later in-place mutations of
68 # the passed-in argument do not mess up the cache.
---> 69 return _get_text_metrics_with_cache_impl(
70 weakref.ref(renderer), text, fontprop.copy(), ismath, dpi)
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/text.py:77, in _get_text_metrics_with_cache_impl(renderer_ref, text, fontprop, ismath, dpi)
73 @functools.lru_cache(4096)
74 def _get_text_metrics_with_cache_impl(
75 renderer_ref, text, fontprop, ismath, dpi):
76 # dpi is unused, but participates in cache invalidation (via the renderer).
---> 77 return renderer_ref().get_text_width_height_descent(text, fontprop, ismath)
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/backends/backend_agg.py:217, in RendererAgg.get_text_width_height_descent(self, s, prop, ismath)
213 return super().get_text_width_height_descent(s, prop, ismath)
215 if ismath:
216 ox, oy, width, height, descent, font_image = \
--> 217 self.mathtext_parser.parse(s, self.dpi, prop)
218 return width, height, descent
220 font = self._prepare_font(prop)
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/mathtext.py:79, in MathTextParser.parse(self, s, dpi, prop, antialiased)
77 prop = prop.copy() if prop is not None else None
78 antialiased = mpl._val_or_rc(antialiased, 'text.antialiased')
---> 79 return self._parse_cached(s, dpi, prop, antialiased)
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/mathtext.py:100, in MathTextParser._parse_cached(self, s, dpi, prop, antialiased)
97 if self._parser is None: # Cache the parser globally.
98 self.__class__._parser = _mathtext.Parser()
--> 100 box = self._parser.parse(s, fontset, fontsize, dpi)
101 output = _mathtext.ship(box)
102 if self._output_type == "vector":
File ~/Python/miniconda3/envs/gammaALPs-0.4/lib/python3.9/site-packages/matplotlib/_mathtext.py:2165, in Parser.parse(self, s, fonts_object, fontsize, dpi)
2162 result = self._expression.parseString(s)
2163 except ParseBaseException as err:
2164 # explain becomes a plain method on pyparsing 3 (err.explain(0)).
-> 2165 raise ValueError("\n" + ParseException.explain(err, 0)) from None
2166 self._state_stack = []
2167 self._in_subscript_or_superscript = False
ValueError:
r [\SI{}{\kilo\parsec}]
^
ParseFatalException: Unknown symbol: \SI, found '\' (at char 3), (line:1, col:4)
<Figure size 640x480 with 1 Axes>
We plot the electron density. Further, nel is used for calculating the RM at a later point.
[11]:
nel = ml.modules[0]._nelicm(r) # in cm^-3
plt.plot(r, nel)
plt.yscale('log')
plt.xscale('log')
plt.ylabel('$n_\mathrm{el}$ (cm$^{-3}$)')
plt.xlabel('$r$ (kpc)')
[11]:
Text(0.5, 0, '$r$ (kpc)')
We can further plot a field line to illustrate that the field really is structured, opposed to a gaussian turbulence field.
For this, we use the equation determining a single field line \(\newcommand{\vert}{|}\newcommand{\vect}{\mathbf}\vect{x}_{i+1} = \vect{x}_{i} + \frac{\vect{B}(\vect{x}_i)}{\left \vert \vect{B}(\vect{x}_i) \right \vert} \Delta s\) for a small \(\Delta s\).
[12]:
import numpy as np
import matplotlib as mpl
from scipy.spatial.transform import Rotation as R
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.collections import LineCollection
from matplotlib.colors import ListedColormap, BoundaryNorm
Functions and constants
In the following block, we define the neccessary functions (field components etc.), some constants and transformations of coordinates and vector fields to be used.
[13]:
def field_line(x_0, steps, ds):
x = np.zeros((steps, 3))
x[0] = x_0
for i in range(steps-1):
x[i+1] = integrate(x[i], ds)
return x
def integrate(x_prev, ds):
'''Returns next point along field line, given some starting point x_prev and a small increment ds'''
r, theta, phi = cart_to_sphere(x_prev[0], x_prev[1], x_prev[2])
B_sph = B_sphere(r, theta, phi)
B_cart = np.matmul(trafo(r, theta, phi), B_sph.transpose()).transpose()
norm = np.linalg.norm(B_cart)
x_new = x_prev + B_cart / norm * ds
return x_new
def cart_to_sphere(x, y, z):
'''Coordinate transformation from cartesian to spherical.'''
r = np.sqrt(x**2 + y**2 + z**2)
theta = np.arccos(z / r)
phi = np.arctan2(y, x)
return r, theta, phi
def sphere_to_cart(r, theta, phi):
'''Coordinate transformation from spherical to cartesian'''
x = r * np.sin(theta) * np.cos(phi)
y = r * np.sin(theta) * np.sin(phi)
z = r * np.cos(theta)
return x, y, z
def trafo(r, theta, phi):
'''Returns trafo used for vector fields, i.e. Field A cartesian = S * Field A spherical. S is returned'''
mat = np.array([[np.sin(theta) * np.cos(phi), np.cos(theta) * np.cos(phi), -np.sin(phi)],
[np.sin(theta) * np.sin(phi), np.cos(theta) * np.sin(phi), np.cos(phi)],
[np.cos(theta), -np.sin(theta), 0]])
return mat
'''Constants used for field definitions'''
alpha = 5.7634591968
c = 1
F_0 = c * (alpha * np.cos(alpha) - np.sin(alpha)) * alpha**2
def b_r(r, theta):
val = 2 * np.cos(theta) * f(r) / r**2
return val
def b_theta(r, theta):
val = - np.sin(theta) * f_prime(r) / r
return val
def b_phi(r, theta):
val = alpha * np.sin(theta) * f(r) / r
return val
def f(r):
return c * (alpha * np.cos(alpha * r) - np.sin(alpha * r) / r) \
- F_0 * r**2 / alpha**2
def f_prime(r):
return c * ( - alpha**2 * np.sin(alpha * r) \
- alpha * np.cos(alpha * r) / r \
+ np.sin(alpha * r) / r**2) \
- 2 * F_0 * r / alpha**2
def B_sphere(r, theta, phi):
return np.array([[b_r(r, theta), b_theta(r, theta), b_phi(r, theta)]])
def norm(r):
return np.sqrt(r[:, :, :, 0]**2)
Plotting the field line
We plot an “interactive” field line, in the sense that you can rotate it around in space!
Because pyplot lacks a proper function for coloured 3D plots, we use instead a scatter with some colormap.
[14]:
%matplotlib notebook
# so that we can have an interactive plot!
fig = plt.figure()
ax = fig.gca(projection='3d')
ax.set_xlabel('$x$')
ax.set_ylabel('$y$')
ax.set_zlabel('$z$')
# start a field line at r=0.2, theta=pi / 2 (on the equator), phi=0
# choose some stepsize ds=0.001 and integrate 100000 steps
ds = 0.001
steps = 100000
#if you want to plot multiple field lines, change this for loop
for phi in np.linspace(0, 2 * np.pi, num=1):
x_0 = sphere_to_cart(0.2, np.pi / 2, phi)
x = field_line(x_0, steps, ds)
#print(x)
#plot only every 10th point to save memory and time
line = ax.scatter(x[::10, 0], x[::10, 1], x[::10, 2],
c=np.linalg.norm(x[::10], axis=-1),
cmap=plt.get_cmap('viridis'),
s=1, marker='.')
cbar = fig.colorbar(line, ax=ax, label='$r$')
ax.set_xticks((-0.5, 0, 0.5))
ax.set_yticks((-0.5, 0, 0.5))
ax.set_zticks((-0.4, 0, 0.4))
ax.view_init(elev=25, azim=-107)
---------------------------------------------------------------------------
TypeError Traceback (most recent call last)
Cell In[14], line 5
2 # so that we can have an interactive plot!
4 fig = plt.figure()
----> 5 ax = fig.gca(projection='3d')
6 ax.set_xlabel('$x$')
7 ax.set_ylabel('$y$')
TypeError: gca() got an unexpected keyword argument 'projection'
Calculate rotation measure
Faraday RM can be calculated via a method of the structured field model, arguments areelectron density evaluated at r = structured_field.r in cm^-3
[ ]:
ml.modules[0].Bfield_model.rotation_measure(nel) # in rad * m^-2
Taylor et al. (2006) found RM values between 6500 and 7500 rad m^-2
[ ]: