GalMag - A Python tool for computing realistic galactic magnetic fields

This is the documentation of GalMag, a Python package for computing galactic magnetic fields based on mean field dynamo theory.

For a detailed description of the physics, please refer to the code paper.

GalMag basic tutorial

Welcome to the GalMag tutorial and basic usage guide.

Quick start

The code uses an object named B_field to store all the magnetic field information.

Associated with this, there is a coordinate grid which has to be set up at the time of the B_field initialization.

import galmag
from galmag.B_field import B_field
import numpy as np
import matplotlib.pyplot as plt

# Maximum and minimum values of the coordinates for each direction
box_limits = [[-15, 15],[-15, 15],[-15, 15]] # kpc
box_resolution = [100,100,100]

B = B_field(box_limits, box_resolution)

The grid - which is an uniform cartesian grid in the default case - can be accessed through the attribute B_field.grid.

This attribute contains a Grid object, which itself contains numpy arrays in the coordinate attributes x, y, z, r_cylindrical, r_spherical, theta and phi.

The coordinate systems place the galactic mid-plane at \(z=0\), \(\theta = \pi/2\).

# Example access to the coordinate grid
print('r_cylindrical type:',type(B.grid.r_spherical))
print('r_cylindrical-grid shape:', B.grid.r_cylindrical.shape)

r_cylindrical type: <class 'numpy.ndarray'>
r_cylindrical-grid shape: (100, 100, 100)

Once the B_field object is initialized, one can add a disk field component. One can do this by specifying 1. The strength of \(B_\phi\) (the azimutal component of the magnetic field) at the solar radius (on the midplade), 2. The positions of the field reversals (the code will try to find solutions with at least these reversals, but there can be more of them), 3. The minimum number of modes to be used for this.

Itens 1 and 3 are actually optional: if absent, a strength of \(-3\,\mu{\rm G}\) is assumed for the \(\phi\) component at the solar radius. There are, of course, many other parameters involved in this calculation which for now are kept with their default values.

The following code snippet exemplifies the use of B_field.add_disk_field method.

reversals_positions = [4.7,12.25] # Requires 2 reversals at 4.1 kpc and 12.25
B_phi_solar_radius = -3 # muG
number_of_modes = 3

B.add_disk_field(reversals=reversals_positions,
                 B_phi_solar_radius=B_phi_solar_radius,
                 number_of_modes=number_of_modes)

One can also include a halo field, using the B_field.add_halo_field method.

Here we keep all the parameters with their standard values and choose to use only adjust the value of \(B_{\phi, {\rm halo}}\) at a reference position which roughly corresponds to the position of the Sun in the Galaxy.

Bphi_sun = 0.1 # muG at the Sun's position

B.add_halo_field(halo_ref_Bphi=Bphi_sun)

Now the galactic magnetic field has been computed and is stored in the B_field object. Specific field components can be accessed in the same way as grid components were accessed. The following snippet exemplifies this showing the \(x\)-depence of \(B_\phi\).

x = B.grid.x[50:,50,50]
Bphi = B.phi[50:,50,50]
y = B.grid.x[50,50,50]
z = B.grid.x[50,50,50]

import galmag.analysis.visualization as visu
visu.std_setup() # Sets matplotlib's rc parameters (fonts, colors, sizes, etc)

plt.figure(figsize=(12,3))
plt.title(r'$y={0:.2f}\,\rm kpc$, $z={1:.2f}\,\rm kpc$'.format(y,z))
plt.plot(x, Bphi)
plt.xlabel(r'$x\,\,[{\rm kpc}]$')
plt.ylabel(r'$B_\phi\,\,[\mu{\rm G}]$')
plt.grid()

plt.figure(figsize=(8,7))
plt.title(r'$z={0:.2f} \,\rm kpc$'.format(z))
visu.plot_x_y_uniform(B.disk, iz=50, field_lines=False);

Obs. Note that the slicing used to create the previous plot exploits the fact that the indices in the cartesian follow the coordinates (i.e. the first index corresponds to \(x\), the second to \(y\) and the third to \(z\)).

Separate field information can be accessed at any point through the attributes B_field.disk and B_field.halo. For example, it is possible to plot the \(x\)-dependence of the halo field in the following way.

x = B.grid.x[50:,50,50]
Bphi_halo = B.halo.phi[50:,50,50]
Bphi_disk = B.disk.phi[50:,50,50]
y = B.grid.x[50,50,50]
z = B.grid.x[50,50,50]

plt.figure(figsize=(12,5))

plt.subplot(2,1,1)
plt.title('$y={0:.2f}$, $z={1:.2f}$'.format(y,z), fontsize=14)
plt.plot(x, Bphi_disk, )
plt.ylabel(r'$B_{\phi,\rm disk}\,\,[\mu{\rm G}]$', fontsize=14)

plt.subplot(2,1,2)
plt.plot(x, Bphi_halo, )
plt.xlabel(r'$x\,\,[{\rm kpc}]$', fontsize=14)
plt.ylabel(r'$B_{\phi,\rm halo}\,\,[\mu{\rm G}]$', fontsize=14);

All the used parameters can be accessed at any point through the parameters attribute, which contain a dictionary.

B.parameters

{'disk_modes_normalization': array([ 1.26205786,  0.20083939, -1.54587971]),
 'disk_height': 0.5,
 'disk_radius': 17.0,
 'disk_turbulent_induction': 0.386392513,
 'disk_dynamo_number': -20.4924192,
 'disk_shear_function': <function galmag.disk_profiles.Clemens_Milky_Way_shear_rate(R, R_d=1.0, Rsun=8.5, normalize=True)>,
 'disk_rotation_function': <function galmag.disk_profiles.Clemens_Milky_Way_rotation_curve(R, R_d=1.0, Rsun=8.5, normalize=True)>,
 'disk_height_function': <function galmag.disk_profiles.exponential_scale_height(R, h_d=1.0, R_HI=5, R_d=1.0, Rsun=8.5)>,
 'disk_regularization_radius': None,
 'disk_ref_r_cylindrical': 8.5,
 'disk_field_decay': True,
 'disk_newman_boundary_condition_envelope': False,
 'halo_symmetric_field': True,
 'halo_n_free_decay_modes': 4,
 'halo_dynamo_type': 'alpha2-omega',
 'halo_rotation_function': <function galmag.halo_profiles.simple_V(rho, theta, phi, r_h=1.0, Vh=220, fraction=0.2, normalize=True, fraction_z=None, legacy=False)>,
 'halo_alpha_function': <function galmag.halo_profiles.simple_alpha(rho, theta, phi, alpha0=1.0)>,
 'halo_Galerkin_ngrid': 501,
 'halo_growing_mode_only': False,
 'halo_compute_only_one_quadrant': True,
 'halo_turbulent_induction': 4.331,
 'halo_rotation_induction': 203.65472,
 'halo_radius': 15.0,
 'halo_ref_radius': 8.5,
 'halo_ref_z': 0.02,
 'halo_ref_Bphi': 0.1,
 'halo_rotation_characteristic_radius': 3.0,
 'halo_rotation_characteristic_height': 1000,
 'halo_manually_specified_coefficients': None,
 'halo_do_not_normalize': False,
 'halo_field_growth_rate': (0.002791486037391211+7.133179987229059j),
 'halo_field_coefficients': array([ 0.14206043,  0.85703689,  0.09889403, -0.41846486])}

All the parameter names can be used as keyword arguments for the methods add_disk_field and add_halo_field previously shown.

Also, it is possible (and preferable) to access separately the parameters associated with a specific component – i.e. halo or disc.

B.disk.parameters

{'disk_modes_normalization': array([ 1.26205786,  0.20083939, -1.54587971]),
 'disk_height': 0.5,
 'disk_radius': 17.0,
 'disk_turbulent_induction': 0.386392513,
 'disk_dynamo_number': -20.4924192,
 'disk_shear_function': <function galmag.disk_profiles.Clemens_Milky_Way_shear_rate(R, R_d=1.0, Rsun=8.5, normalize=True)>,
 'disk_rotation_function': <function galmag.disk_profiles.Clemens_Milky_Way_rotation_curve(R, R_d=1.0, Rsun=8.5, normalize=True)>,
 'disk_height_function': <function galmag.disk_profiles.exponential_scale_height(R, h_d=1.0, R_HI=5, R_d=1.0, Rsun=8.5)>,
 'disk_regularization_radius': None,
 'disk_ref_r_cylindrical': 8.5,
 'disk_field_decay': True,
 'disk_newman_boundary_condition_envelope': False}

(This is actually safer than just using ``B_field.parameters``)

Grid object

All the calculations are done on top of a coordinate grid, which is stored as a `Grid <http://galmag.readthedocs.io/en/latest/galmag.html#galmag.Grid.Grid>`__ object. This class can be directly accessed from the module Grid, and instantiated using a similar syntax as the one used for the B_field.

# A cartesian grid can be constructed as follows
cartesian_grid = galmag.Grid(box=[[-15, 15],[-15, 15],[-15, 15]],
                      resolution = [12,12,12],
                      grid_type = 'cartesian' # Optional
                      )
# For cylindrical grid, the limits are specified assuming
# the order: r (cylindrical), phi, z
cylindrical_grid = galmag.Grid(box=[[0.25, 15],[-np.pi,np.pi],[-15,15]],
                      resolution = [9,12,9],
                      grid_type = 'cylindrical')

# For spherical grid, the limits are specified assuming
# the order: r (spherical), theta, phi (azimuth)
spherical_grid = galmag.Grid(box=[[0, 15],[0, np.pi], [-np.pi,np.pi],],
                      resolution = [12,10,10],
                      grid_type = 'spherical')

The geometry of the different grids can be easily understood from the following plot (note the automatic conversion to cartesian coordinates through the attributes x, y, z).

plt.figure(figsize=(12,5))
plt.subplot(1,2,1)
plt.scatter(cartesian_grid.x, cartesian_grid.y, color='y', label='cartesian', alpha=0.5)
plt.scatter(cylindrical_grid.x, cylindrical_grid.y, label='cylindrical', alpha=0.5)
plt.scatter(spherical_grid.x, spherical_grid.y, color='g', label='spherical', alpha=0.5)
plt.xlabel('x'); plt.ylabel('y')
plt.legend()

plt.subplot(1,2,2)
plt.scatter(cartesian_grid.x, cartesian_grid.z, color='y', label='cartesian', alpha=0.5)
plt.scatter(cylindrical_grid.x, cylindrical_grid.z, label='cylindrical', alpha=0.5)
plt.scatter(spherical_grid.x, spherical_grid.z, label='spherical', color='g', alpha=0.5)
plt.xlabel('x'); plt.ylabel('z')
plt.legend();

Disc field

The disc field in GalMag is constructed by expanding the dynamo solution (associated with a given choice of parameters) in a series of modes and assigning an adequate choice of coefficients \(C_n\!\)’s to the \(n\) dominant modes.

Direct choice of normalizations

The choice of the \(C_n\)’ can be done explicitly setting up the parameter disk_modes_normalization, which should contain an array of coefficients. Each mode is pre-normalized so that its absolute value at a reference cylindrical radius is unit, i.e. \(| \mathbf{B}_{\rm mode}(R_{\rm ref},0,0)|=1\,\). \(\;R_{\rm ref}\) is set by the parameter disk_ref_r_cylindrical, whose default value is \(8.5\,{\rm kpc}\), roughly the Sun’s radius.

To exemplify this, we will generate a new B_field object based on an uniform cylindrical grid and add to it a disc magnetic field corresponding to the second mode only.

plt.figure(figsize=(12,5))
box_limits = [[1e-8, 15],# r [kpc]
              [0,0], # phi
              [0,0]] # z [kpc]
box_resolution = [100,1,1] # No phi or z variation
Bcyl = B_field(box_limits, box_resolution, grid_type='cylindrical')

Bcyl.add_disk_field(disk_modes_normalization=[0,2]) # Skips first mode, normalization 2 for the second mode

y = np.sqrt(Bcyl.r_cylindrical[:,0,0]**2 + Bcyl.phi[:,0,0]**2 +Bcyl.z[:,0,0]**2)
x = Bcyl.grid.r_cylindrical[:,0,0]
plt.plot(x,y)
plt.ylabel(r'$| \mathbf{B}_{\rm mode}|$')
plt.xlabel(r'$R\,[{\rm kpc}]$')
plt.grid();

For a more interesting application, let us examine how \(B_\phi\) components of the different modes behave.

plt.figure(figsize=(12,5))
box_limits = [[1e-8, 15],# r [kpc]
              [0,0], # phi
              [0,0]] # z [kpc]
box_resolution = [100,1,1] # No phi or z variation
Bcyl = B_field(box_limits, box_resolution, grid_type='cylindrical')

nmodes = 6

for i in range(nmodes):
    mode_norm = np.zeros((i+1))
    mode_norm[i]=1

    # Overwrites the disk field with another mode
    Bcyl.add_disk_field(disk_modes_normalization=mode_norm)

    x = Bcyl.grid.r_cylindrical[:,0,0]
    y = Bcyl.phi[:,0,0]
    plt.plot(x,y/abs(y).max(), label='mode {}'.format(i+1))
    plt.ylabel(r'$B_\phi/\max(B_\phi)$')
    plt.xlabel(r'$R\,[{\rm kpc}]$')
plt.legend(loc='upper left')
plt.grid();

Reversals

As it was discussed in the Quick Start section, the disk field can be constructed by specifying the position of the magnetic field reversals (positions \(R_{\rm rev}\) along the mid-plane where \(B_\phi\) changes sign), the strength of the field at a reference point (generally \(R_{\rm Sun}\approx 8.5\,{\rm kpc}\)) and the number of modes to be used.

This is done using a least squares method to find the choice of \(C_n\)’s which offers best fit to these constraints. It is worth noting that such solution is not necessarily unique and, for some special cases, may not be exact.

Thus, while this is a good tool for the quick generation of examples, the user should use it carefully.

plt.figure(figsize=(12,5))
reversals_positions = [3,10]
B_phi_ref = 3 # muG
number_of_modes = 4

Bcyl.add_disk_field(reversals=reversals_positions,
                    B_phi_ref=B_phi_ref,
                    number_of_modes=number_of_modes)

y = Bcyl.phi[:,0,0]
x = Bcyl.grid.r_cylindrical[:,0,0]
plt.title(r'Reversals at $3$ and $10\,\rm kpc$ using 4 modes.')
plt.plot(x,x*0, 'r:', alpha=0.75)
plt.plot(x,y)
plt.ylabel(r'$B_\phi$')
plt.xlabel(r'$R\,[{\rm kpc}]$')
plt.grid();

GalMag’s choice of \(C_n\)’s can found in the parameters dictionary, under the name disk_modes_normalization.

Bcyl.parameters

{'disk_modes_normalization': array([-0.91287645, -0.82464857,  0.83137772,  0.43987423]),
 'disk_height': 0.5,
 'disk_radius': 17.0,
 'disk_turbulent_induction': 0.386392513,
 'disk_dynamo_number': -20.4924192,
 'disk_shear_function': <function galmag.disk_profiles.Clemens_Milky_Way_shear_rate(R, R_d=1.0, Rsun=8.5, normalize=True)>,
 'disk_rotation_function': <function galmag.disk_profiles.Clemens_Milky_Way_rotation_curve(R, R_d=1.0, Rsun=8.5, normalize=True)>,
 'disk_height_function': <function galmag.disk_profiles.exponential_scale_height(R, h_d=1.0, R_HI=5, R_d=1.0, Rsun=8.5)>,
 'disk_regularization_radius': None,
 'disk_ref_r_cylindrical': 8.5,
 'disk_field_decay': True,
 'disk_newman_boundary_condition_envelope': False}

Disc parameters

The following parameters control the disc magnetic field:

1. disk_dynamo_number - Default: -20.49

   Dynamo number, \(D= R_\omega R_\alpha\), associated with the disc component.
2. disk_field_decay - Default: True

   Sets behaviour for \(|z|>h(R)\). If True, the field will decay with \(|z|^{-3}\), otherwise, the field will be assumed to be constant.
3. disk_height - Default: 0.5 \([\rm kpc ]\),

   Disc height at the solar radius, \(h_\odot\).
4. disk_height_function - Default: galmag.disk_profiles.exponential_scale_height

   Disc height as a function of the cylindrical radius, \(h(R)\). The default function corresponds to an exponential disc. A function of a constant height disc can be found in galmag.disk_profiles.constant_scale_height
5. disk_modes_normalization.

   An n-array containing the the coefficients of the first n modes (see Direct choice of normalizations section).
6. disk_radius - Default: 17 \([\rm kpc]\)

   Radius of the dynamo active region of the disc.
7. disk_rotation_function - Default: galmag.disk_profiles.Clemens_Milky_Way_rotation_curve

   Rotation curve of the disc. By default, assumes the rotation curve obtained by Clemens (1985). One alternative, mostly for testing is galmag.disk_profiles.solid_body_rotation_curve
8. disk_shear_function - Default: galmag.disk_profiles.Clemens_Milky_Way_shear_rate

   Shear rate as a function of the cylindrical radius of the disc. By default, assumes a shear rate compatible with the rotation curve obtained by Clemens (1985). Also available is (again, mostly for testing) is a constant shear profile: galmag.disk_profiles.constant_shear_rate .
9. disk_turbulent_induction - Default: 0.386,

   Dimensionless intensity of helical turbulence, \(R_\alpha\).
10. disk_ref_r_cylindrical: 8.5 \([\rm kpc]\)

    Reference cylindrical radius corresponding, \(s_0\). This is used both for the normalization of the magnetic field and for the normalization of the height profile.

Halo field

To specify the halo field, one chooses the field strength at a reference position, whether the field is symmetric or anti-symmetric and the number of free-decay modes used to construct the solution. All these (and other parameters listed later) can be used as arguments for the B_field.add_halo_field method.

Generally, more than one solution is found for a given set of parameters. The B_field.add_halo_field method will include the one with the largest growth rate.

It is possible to access the growth/decay rate of the solution found through the attribute B_field_component.growth_rate.

Similarly, the coefficients used can be accessed using B_field_component.coefficients.

# Using the previously defined B object
B.add_halo_field(
    # Number of free decay modes to be used in the Galerkin expansion
    halo_n_free_decay_modes = 4,
    # cylindrical radius of the reference point
    halo_ref_radius = 9,
    # z coordinate of the reference point
    halo_ref_z = 0.02,
    # B_phi value at the reference point
    halo_ref_Bphi = 4.1, # [muG]
    # Chooses between symmetric and anti-symmetric solution
    halo_symmetric_field = True # Symmetric
    )

plt.figure(figsize=(13.1,11))
plt.subplot(2,2,1)
plt.title('Halo only')
visu.plot_x_z_uniform(B.halo, skipx=4,skipz=4, iy=49)
plt.subplot(2,2,2)
plt.title('Halo + disk')
visu.plot_x_z_uniform(B, skipx=4,skipz=4, iy=49)
plt.subplot(2,2,3)
plt.title('Halo only')
visu.plot_x_y_uniform(B.halo, skipy=5,skipx=5, iz=49, field_lines=False)
plt.subplot(2,2,4)
plt.title('Halo + disk')
visu.plot_x_y_uniform(B, skipy=5,skipx=5, iz=49, field_lines=False)
plt.tight_layout()

print('Growth rate: ', B.halo.growth_rate)
print('Coefficients: ', B.halo.coefficients)
B.halo.parameters

Growth rate:  (0.002791486037391211+7.133179987229059j)
Coefficients:  [ 0.14206043  0.85703689  0.09889403 -0.41846486]

{'halo_symmetric_field': True,
 'halo_n_free_decay_modes': 4,
 'halo_dynamo_type': 'alpha2-omega',
 'halo_rotation_function': <function galmag.halo_profiles.simple_V(rho, theta, phi, r_h=1.0, Vh=220, fraction=0.2, normalize=True, fraction_z=None, legacy=False)>,
 'halo_alpha_function': <function galmag.halo_profiles.simple_alpha(rho, theta, phi, alpha0=1.0)>,
 'halo_Galerkin_ngrid': 501,
 'halo_growing_mode_only': False,
 'halo_compute_only_one_quadrant': True,
 'halo_turbulent_induction': 4.331,
 'halo_rotation_induction': 203.65472,
 'halo_radius': 15.0,
 'halo_ref_radius': 9,
 'halo_ref_z': 0.02,
 'halo_ref_Bphi': 4.1,
 'halo_rotation_characteristic_radius': 3.0,
 'halo_rotation_characteristic_height': 1000,
 'halo_manually_specified_coefficients': None,
 'halo_do_not_normalize': False,
 'halo_field_growth_rate': (0.002791486037391211+7.133179987229059j),
 'halo_field_coefficients': array([ 0.14206043,  0.85703689,  0.09889403, -0.41846486])}

Galerkin expansion coefficients

It is likely that one may want to study directly which choices of parameters lead to which growth rates. This can be done using the module galmag.galerkin, which contains the function Galerkin_expansion_coefficients.

from galmag.galerkin import Galerkin_expansion_coefficients

The function requires a dictionary containing the following parameters (described individually in the comments):

p = { # Square root of the number of grid points used in the calculations
      'halo_Galerkin_ngrid': 500,
      # Function specifying the alpha profile of the halo
      'halo_alpha_function': galmag.halo_profiles.simple_alpha,
      # R_\alpha
      'halo_turbulent_induction': 7.7,
      # R_\omega
     'halo_rotation_induction': 200.0,
      # Number of free decay modes to be used for the expansion
     'halo_n_free_decay_modes': 4,
      # Whether the field is symmetric or antisymmetric
      'halo_symmetric_field': True,
      # Defines the dynamo type
      'halo_dynamo_type': 'alpha-omega',
      # Function specifying the rotation curve of the halo
      'halo_rotation_function': galmag.halo_profiles.simple_V,
      # Parameters used in this function
      'halo_rotation_characteristic_height': 1e3,
      'halo_rotation_characteristic_radius': 3.0,
      # Halo radius
      'halo_radius': 20.,
     }

The function returns the growth rates and the coefficients, \(a_i\) of the free decay modes. Optionally, the intermediate matrix \(W_{ij} = \int \mathbf{B_j} \cdot \hat{W} \,\mathbf{B_i} \;{\rm d}^3\mathbf{r}\;\) is

vals, vecs, Wij = Galerkin_expansion_coefficients(p, return_matrix=True)

print('The resultant matrix is:\n', Wij)

for i, (val, vec) in enumerate(zip(vals, vecs)):
    print('\nSolution ',i)
    print(' Growth rate =', val)
    for i, c in enumerate(vec):
        print('    a{0} = {1}'.format(i+1, c))

The resultant matrix is:
 [[-2.01907286e+01  1.25802759e+01  1.98763272e+01  1.49637059e+01]
 [ 1.38178883e+02 -2.01907286e+01 -1.12133344e+01  0.00000000e+00]
 [-7.44585112e+00  1.11181174e-03 -4.88311936e+01 -2.27594189e+01]
 [-7.28554698e+01  0.00000000e+00 -1.54579219e+02 -4.88311936e+01]]

Solution  0
 Growth rate = (-82.5472249922034+9.72793002473695j)
    a1 = (-0.3304728627932634+0.04590311536614698j)
    a2 = (-0.3304728627932634-0.04590311536614698j)
    a3 = (0.19556281135575382+0.039281681425336674j)
    a4 = (0.19556281135575382-0.039281681425336674j)

Solution  1
 Growth rate = (-82.5472249922034-9.72793002473695j)
    a1 = (0.7728996311795923+0j)
    a2 = (0.7728996311795923-0j)
    a3 = (0.7782745573919754+0j)
    a4 = (0.7782745573919754-0j)

Solution  2
 Growth rate = (13.525302792151383+9.587979337319839j)
    a1 = (0.2257082569733093-0.10486375376800403j)
    a2 = (0.2257082569733093+0.10486375376800403j)
    a3 = (0.06976706192021101-0.18140737243306812j)
    a4 = (0.06976706192021101+0.18140737243306812j)

Solution  3
 Growth rate = (13.525302792151383-9.587979337319839j)
    a1 = (0.39769866232974604-0.26683684329612717j)
    a2 = (0.39769866232974604+0.26683684329612717j)
    a3 = (-0.3315122408562829+0.45477951371599235j)
    a4 = (-0.3315122408562829-0.45477951371599235j)

Halo parameters

1. halo_Galerkin_ngrid

   Number of grid points used for the Galerkin expansion calculations.
2. halo_alpha_function - Default: galmag.halo_profiles.simple_alpha

   alpha-effect profile. By default, it is assumed simply: \(\alpha \propto \cos(\theta)\,\).
3. halo_dynamo_type - Default: 'alpha2-omega'

   Chooses the dominant type of dynamo. Available options: 'alpha-omega' and 'alpha2-omega'.
4. halo_growing_mode_only - Default: False.

   If true, returns a zero field if all modes are decaying for a particular choice of parameters.
5. halo_n_free_decay_modes - Default: 4

   Number of free decay modes to be used for the Galerkin expansion.
6. halo_radius - Default: 15 [\({\rm kpc}\)]

   Radius of the dynamo active region of the halo.
7. halo_ref_Bphi - Default: -0.5 [\({\mu\rm G}\)]

   Strength of the halo magnetic field at the halo reference point.
8. halo_ref_radius - Default: 8.5 [\({\rm kpc}\)]

   Cylindrical radius of the halo reference point.
9. halo_ref_z - Default: 0.02 [\({\rm kpc}\)]

   z-coordinate of the halo reference point.
10. halo_rotation_function - Default: galmag.halo_profiles.simple_V

    Rotation curve of the halo. By default a simple rotation curve with no z dependence, exponentially decaying with the cylindrical radius, is assumed.
11. halo_rotation_induction - Default: 203.65

    Dimensionless induction by differential rotation, \(R_\omega\), for the halo.
12. halo_symmetric_field - Default: True

    If True, the field is assumed to be symmetric over the midplane. If False, it is assumed to be anti-symmetric.
13. halo_turbulent_induction: 4.331

    Dimensionless intensity of helical turbulence, \(R_\alpha\), for the halo.

Plotting tools

GalMag comes with a small selection of plotting tools, particularly to facilitate the preparation of 2D slices for diagnostic. These can be found in the galmag.analysis.visualization module.

help(galmag.analysis.visualization)

Help on module galmag.analysis.visualization in galmag.analysis:

NAME
    galmag.analysis.visualization

FUNCTIONS
    plot_r_z_uniform(B, skipr=3, skipz=5, quiver=True, contour=True, quiver_color='0.25', cmap='viridis', field_lines=True, vmin=None, vmax=None, levels=None, **kwargs)
        Plots a r-z slice of the field. Assumes B is created using a cylindrical
        grid - for a more sophisticated/flexible plotting script which does not
        rely on the grid structure check the plot_slice.

        The plot consists of:
          1) a coloured contourplot of :math:`B_\phi`
          2) quivers showing the x-z projection of the field

        Parameters
        ----------
        B : B_field
            a B_field or B_field_component object
        quiver : bool
            If True, shows quivers. Default True
        contour : bool
            If True, shows contours. Default: True
        skipx/skipz : int
            Tweaks the display of quivers (Default: skipz=5, skipr=3)

    plot_slice()
        Plots slice of arbitrary orientation

        Note
        ----
        Not implemented yet

    plot_x_y_uniform(B, skipx=5, skipy=5, iz=0, field_lines=True, quiver=True, vmin=None, vmax=None, contour=True, levels=None, quiver_color='0.25', cmap='viridis', **kwargs)
        Plots a x-y slice of the field. Assumes B is created using a cartesian
        grid - for a more sophisticated/flexible plotting script which does not
        rely on the grid structure check the plot_slice.

        The plot consists of:
          1) a coloured contourplot of :math:`|B|^2`
          2) Field lines of the :math:`B_x` and :math:`B_y` field
          3) Quivers showing the :math:`B_x` and :math:`B_y` field

        Parameters
        ----------
        B : B_field
            a B_field or B_field_component object
        field_lines : bool
            If True, shows field lines. Default: True
        quiver : bool
            If True, shows quivers. Default True
        contour : bool
            If True, shows contours. Default: True
        skipx/skipy : int
            Tweaks the display of quivers (Default: skipx=5, skipy=5)

    plot_x_z_uniform(B, skipx=1, skipz=5, iy=0, quiver=True, contour=True, quiver_color='0.25', cmap='viridis', vmin=None, vmax=None, no_colorbar=False, **kwargs)
        Plots a x-z slice of the field. Assumes B is created using a cartesian
        grid - for a more sophisticated/flexible plotting script which does not
        rely on the grid structure check the plot_slice.

        The plot consists of:
          1) a coloured contourplot of :math:`B_\phi`
          2) quivers showing the x-z projection of the field

        Parameters
        ----------
        B : B_field
            a B_field or B_field_component object
        quiver : bool
            If True, shows quivers. Default True
        contour : bool
            If True, shows contours. Default: True
        skipx/skipz : int
            Tweaks the display of quivers (Default: skipz=5, skipx=1)

    plot_y_z_uniform(B, skipy=5, skipz=5, ix=0, quiver=True, contour=True, quiver_color='0.25', cmap='viridis', vmin=None, vmax=None, **kwargs)
        Plots a y-z slice of the field. Assumes B is created using a cartesian
        grid - for a more sophisticated/flexible plotting script which does not
        rely on the grid structure check the plot_slice.

        The plot consists of:
          1) a coloured contourplot of :math:`B_\phi`
          2) Quivers showing the y-z projection of the field

        Parameters
        ----------
        B : B_field
            a B_field or B_field_component object
        quiver : bool
            If True, shows quivers. Default True
        contour : bool
            If True, shows contours. Default: True
        skipy/skipz : int
            Tweaks the display of quivers (Default: skipz=5, skipy=5)

    std_setup()
        Adjusts matplotlib default settings

FILE
    /home/lrodrigues/GalMag/galmag/analysis/visualization.py

Advanced

Adding extra major components: Gaussian random field example

If one explores the B_field object a little further, one finds there is a method called B_field.set_field_component which takes a component name and a B_field_component object as inputs.

We exemplify, below, how to add an extra personalised Gaussian random field component to our B_field object.

We start by generating the noisy field. To keep things divergence-free we first construct a random vector potential, \(\mathbf A\), and then take the curl of it.

from galmag.util import derive

mu = 0
sigma = 0.0778 # magic number to generate Brms~0.5
print(sigma)
# Defines a random vector potential
A_rnd = {} # Dictionary of vector components
for i, c in enumerate(('x','y','z')):
    A_rnd[c] = np.random.normal(mu, sigma, B.resolution.prod())
    A_rnd[c] = A_rnd[c].reshape(B.resolution)

# Prepares the derivatives to compute the curl
dBi_dj ={}
for i, c in enumerate(['x','y','z']):
    for j, d in enumerate(['x','y','z']):
        dj = (B.box[j,-1]-B.box[j,0])/float(B.resolution[j])
        dBi_dj[c,d] = derive(A_rnd[c],dj, axis=j)

# Computes the curl of A_rnd
tmpBrnd = {}
tmpBrnd['x'] = dBi_dj['z','y'] -  dBi_dj['y','z']
tmpBrnd['y'] = dBi_dj['x','z'] -  dBi_dj['z','x']
tmpBrnd['z'] = dBi_dj['y','x'] -  dBi_dj['x','y']

del dBi_dj; del A_rnd

for c in tmpBrnd:
    print('B{0}: mean = {1}  std = {2}'.format(c,
                                                tmpBrnd[c].mean(),
                                                tmpBrnd[c].std()))

0.0778
Bx: mean = -0.00017563838972479574  std = 0.288666276254484
By: mean = 0.0001620638666346189  std = 0.2886417232520821
Bz: mean = 5.918382704113669e-05  std = 0.2882613118574375

Note that the format is compatible with the grid of the previously generated field object.

Now, we can use these data create a B_field_component object. Note that we included the parameters used for generating the data for future reference.

from galmag.B_field import B_field_component
Brnd_component = B_field_component(B.grid,
                                     x=tmpBrnd['x'],
                                     y=tmpBrnd['y'],
                                     z=tmpBrnd['z'],
                                     parameters={'mu': mu,
                                                 'sigma':sigma})
del tmpBrnd

One can now access the x,y,z vector components through the Brnd_component’s attributes.

One can also access the same information in a cylindrical or spherical coordinate system.

print('A sample of Brnd_x')
print(Brnd_component.x[1:3,1:3,1:3])
print('\nA sample of Brnd_theta')
print(Brnd_component.theta[1:3,1:3,1:3])

A sample of Brnd_x
[[[-0.11884636 -0.11828462]
  [ 0.32860367 -0.43862623]]

 [[-0.15642111 -0.34380286]
  [-0.0559142   0.37455865]]]

A sample of Brnd_theta
[[[ 0.00549858  0.0563936 ]
  [ 0.48793848 -0.00555552]]

 [[-0.65227692  0.31251592]
  [-0.06968609  0.14013689]]]

Finally, one can add this component to the previously defined B_field object using

B.set_field_component('random',Brnd_component)

print('Brms', (B.random.x**2+B.random.y**2+B.random.z**2).mean()**0.5)

Brms 0.4997368387247418

Now the new component can be accessed easily using B.random.

plt.figure(figsize=(13.1,11))
plt.subplot(2,2,1)
plt.title('random only')
visu.plot_x_z_uniform(B.random, skipx=4,skipz=4, iy=49)

plt.subplot(2,2,2)
plt.title('random + halo + disk')
visu.plot_x_z_uniform(B, skipx=4,skipz=4, iy=49)

plt.subplot(2,2,3)
plt.title('random only')
visu.plot_x_y_uniform(B.random, skipy=4,skipx=4, iz=49, field_lines=False)

plt.subplot(2,2,4)
plt.title('random + halo + disk')
visu.plot_x_y_uniform(B, skipy=4,skipx=4, iz=49, field_lines=False)
plt.tight_layout()

Generators

To produce the built-in field and halo field components, a generator class is used. The generator is actually stored in the B_component_field object itself. E.g.

B.halo.generator?

Type:        B_generator_halo
String form: <galmag.B_generators.B_generator_halo.B_generator_halo object at 0x7fd8dbdb6090>
File:        ~/GalMag/galmag/B_generators/B_generator_halo.py
Docstring:
Generator for the halo field

Parameters
----------
box : 3x2-array_like
     Box limits
resolution : 3-array_like
     containing the resolution along each axis.
grid_type : str, optional
    Choice between 'cartesian', 'spherical' and 'cylindrical' *uniform*
    coordinate grids. Default: 'cartesian'
dtype : numpy.dtype, optional
    Data type used. Default: np.dtype(np.float)

Most of the physics is coded in the generator classes \(-\) while the B_field and B_field_component classes deal only with coordinate transformations and convenience.

The reader is invited to inspect these modules

galmag.B_generators.B_generator_halo.B_generator_halo and galmag.B_generators.B_generator_disk.B_generator_disk for further details on the computation of the magnetic fields.

One can use the generators directly to produce the magnetic field components (and then possibly include them using the set_field_component method). Also, for the sake of consistency, long term inclusion of new components (e.g. a random/turbulent field component) should be implemented using similar generators.

galmag

galmag package

Subpackages

galmag.B_generators package

Submodules

galmag.B_generators.B_generator module

class galmag.B_generators.B_generator.B_generator(grid=None, box=None, resolution=None, grid_type='cartesian', default_parameters={}, dtype=<class 'float'>)

Bases: object

Base class for B-field generators

Note

This class does not work on its own.

get_B_field()

galmag.B_generators.B_generator_disk module

class galmag.B_generators.B_generator_disk.B_generator_disk(grid=None, box=None, resolution=None, grid_type='cartesian', default_parameters={}, dtype=<class 'float'>)

Bases: galmag.B_generators.B_generator.B_generator

Generator for the disk field

Parameters:

• box (3x2-array_like) – Box limits
• resolution (3-array_like) – containing the resolution along each axis.
• grid_type (str, optional) – Choice between ‘cartesian’, ‘spherical’ and ‘cylindrical’ uniform coordinate grids. Default: ‘cartesian’
• dtype (numpy.dtype, optional) – Data type used. Default: np.dtype(np.float)

find_B_field(B_phi_ref=-3.0, reversals=None, number_of_modes=0, **kwargs)

Constructs B_field objects for the disk field based on constraints.

Parameters:

• B_phi_ref (float) – Magnetic field intensity at the reference radius. Default: -3
• reversals (list_like) – a list containing the r-positions of field reversals over the midplane (units consitent with the grid).
• number_of_modes (int, optional) – Minimum of modes to be used. NB: modes_count = max(number_of_modes, len(reversals)+1)

Note

Other disc parameters should be specified as keyword arguments.

Returns:The computed disc component. Return type:B_field_component

get_B_field(**kwargs)

Computes a B_field object containing the specified disk field.

Parameters:disk_modes_normalization (array_like) – array of coefficients for the disc modes

Note

Other disc parameters should be specified as keyword arguments.

Returns:result_field_obj – The computed disc field component. Return type:galmag.B_field.B_field_component

galmag.B_generators.B_generator_halo module

class galmag.B_generators.B_generator_halo.B_generator_halo(grid=None, box=None, resolution=None, grid_type='cartesian', default_parameters={}, dtype=<class 'float'>)

Bases: galmag.B_generators.B_generator.B_generator

Generator for the halo field

Parameters:

• box (3x2-array_like) – Box limits
• resolution (3-array_like) – containing the resolution along each axis.
• grid_type (str, optional) – Choice between ‘cartesian’, ‘spherical’ and ‘cylindrical’ uniform coordinate grids. Default: ‘cartesian’
• dtype (numpy.dtype, optional) – Data type used. Default: np.dtype(np.float)

get_B_field(**kwargs)

Constructs a B_field component object containing a solution of the dynamo equation for the halo field.

Module contents

galmag.analysis package

Submodules

galmag.analysis.visualization module

galmag.analysis.visualization.plot_r_z_uniform(B, skipr=3, skipz=5, quiver=True, contour=True, quiver_color='0.25', cmap='viridis', field_lines=True, vmin=None, vmax=None, levels=None, **kwargs)

Plots a r-z slice of the field. Assumes B is created using a cylindrical grid - for a more sophisticated/flexible plotting script which does not rely on the grid structure check the plot_slice.

The plot consists of:

1. a coloured contourplot of \(B_\phi\)
2. quivers showing the x-z projection of the field

Parameters:

• B (B_field) – a B_field or B_field_component object
• quiver (bool) – If True, shows quivers. Default True
• contour (bool) – If True, shows contours. Default: True
• skipx/skipz (int) – Tweaks the display of quivers (Default: skipz=5, skipr=3)

galmag.analysis.visualization.plot_slice()

Plots slice of arbitrary orientation

Note

Not implemented yet

galmag.analysis.visualization.plot_x_y_uniform(B, skipx=5, skipy=5, iz=0, field_lines=True, quiver=True, vmin=None, vmax=None, contour=True, levels=None, quiver_color='0.25', cmap='viridis', **kwargs)

Plots a x-y slice of the field. Assumes B is created using a cartesian grid - for a more sophisticated/flexible plotting script which does not rely on the grid structure check the plot_slice.

The plot consists of:

1. a coloured contourplot of \(|B|^2\)
2. Field lines of the \(B_x\) and \(B_y\) field
3. Quivers showing the \(B_x\) and \(B_y\) field

Parameters:

• B (B_field) – a B_field or B_field_component object
• field_lines (bool) – If True, shows field lines. Default: True
• quiver (bool) – If True, shows quivers. Default True
• contour (bool) – If True, shows contours. Default: True
• skipx/skipy (int) – Tweaks the display of quivers (Default: skipx=5, skipy=5)

galmag.analysis.visualization.plot_x_z_uniform(B, skipx=1, skipz=5, iy=0, quiver=True, contour=True, quiver_color='0.25', cmap='viridis', vmin=None, vmax=None, no_colorbar=False, **kwargs)

Plots a x-z slice of the field. Assumes B is created using a cartesian grid - for a more sophisticated/flexible plotting script which does not rely on the grid structure check the plot_slice.

The plot consists of:

1. a coloured contourplot of \(B_\phi\)
2. quivers showing the x-z projection of the field

Parameters:

• B (B_field) – a B_field or B_field_component object
• quiver (bool) – If True, shows quivers. Default True
• contour (bool) – If True, shows contours. Default: True
• skipx/skipz (int) – Tweaks the display of quivers (Default: skipz=5, skipx=1)

galmag.analysis.visualization.plot_y_z_uniform(B, skipy=5, skipz=5, ix=0, quiver=True, contour=True, quiver_color='0.25', cmap='viridis', vmin=None, vmax=None, **kwargs)

Plots a y-z slice of the field. Assumes B is created using a cartesian grid - for a more sophisticated/flexible plotting script which does not rely on the grid structure check the plot_slice.

The plot consists of:

1. a coloured contourplot of \(B_\phi\)
2. Quivers showing the y-z projection of the field

Parameters:

• B (B_field) – a B_field or B_field_component object
• quiver (bool) – If True, shows quivers. Default True
• contour (bool) – If True, shows contours. Default: True
• skipy/skipz (int) – Tweaks the display of quivers (Default: skipz=5, skipy=5)

galmag.analysis.visualization.std_setup()

Adjusts matplotlib default settings

Module contents

Submodules

galmag.B_field module

Contains the definition of the B_field and B_field_component classes.

class galmag.B_field.B_field(box, resolution, grid_type='cartesian', dtype=<class 'float'>, **kwargs)

Bases: object

Galactic magnetic field

B_field objects store data for the whole galactic magnetic field.

Individual galactic magnetic field components (e.g. disk or halo) are stored as attributes containing B_field_component objects.

The properties x, y, z, r_spherical, r_cylindrical, theta and phi return the magnetic field along each of these axis. These are actually the sum of all the galactic components available.

A disk or halo component can be added using the add_disk_field and add_halo_field methods.

A custom new magnetic field component can be set using the set_field_component method. Alternatively, the custom component can be added during initialization as an extra keyword (i.e. B_field.set_field_component(name, component) is equivalent to B_field(…, name=component).

Parameters:

• box (3x2-array_like) – Box limits
• resolution (3-array_like) – containing the resolution along each axis.
• grid_type (str, optional) – Choice between ‘cartesian’, ‘spherical’ and ‘cylindrical’ uniform coordinate grids. Default: ‘cartesian’
• dtype (numpy.dtype, optional) – Data type used. Default: np.dtype(np.float)

add_disk_field(name='disk', **kwargs)

Includes a disc magnetic field component.

Parameters:

• name (str, optional) – Name of the disc component. Default value: ‘disk’
• **kwargs – See galmag.B_generator_disk for the list of disc parameters

add_halo_field(name='halo', **kwargs)

Parameters:

• name – (Default value = ‘halo’)
• **kwargs –

phi

Azimuthal coordinate component, \(\phi\)

r_cylindrical

Cylindrical radial coordinate component, \(s\)

r_spherical

Spherical radial coordinate component, \(r\)

reset_cache()

Erases cached total field values. Major components as ‘disk’ or ‘halo’ are preserved and coordinate component values will be generated on next time they are called.

set_field_component(name, component)

Parameters:

• name –
• component –

theta

Polar coordinate component, \(\theta\)

x

Horizontal coordinate component, \(x\)

y

Horizontal coordinate component, \(y\)

z

Vertical coordinate component, \(z\)

class galmag.B_field.B_field_component(grid, x=None, y=None, z=None, r_spherical=None, r_cylindrical=None, theta=None, phi=None, copy=True, dtype=dtype('float64'), generator=None, parameters={})

Bases: object

A single galactic magnetic field component.

B_field_component objects store data for a single galactic magnetic field component (e.g. the galactic halo magnetic field).

The properties x, y, z, r_spherical, r_cylindrical, theta and phi return the magnetic field along each of these axis (with any required coordinate transformation being performed on the fly).

Parameters:

• grid (grid_object) – See :class: galmag.Grid
• y, z, r_spherical, r_cylindrical, theta, phi (x,) – Values of the magnetic field in coordinates compatible with the grid parameter should be supplied at initialization. Unused coordinates should be set to None. E.g. r_cylindrical=[[..],[..],[..]], z=[[..],[..],[..]], phi=[[..],[..],[..]]
• copy (bool) – whether a fresh copy of the field data should be made upon initialization. Default: True

phi

Azimuthal coordinate component, \(\phi\)

r_cylindrical

Cylindrical radial coordinate component, \(s\)

r_spherical

Spherical radial coordinate component, \(r\)

set_field_data(name, data, copy=True)

Includes fresh data in a particular component

Parameters:

• name (str) – Name of the coordinate component
• data (array_like) – Data
• copy (bool, optional) – Copy (usual) or reference the data? Default: True.

theta

Polar coordinate component, \(\theta\)

x

Horizontal coordinate component, \(x\)

y

Horizontal coordinate component, \(y\)

z

Vertical coordinate component, \(z\)

galmag.Grid module

Contains the definition of the Grid class.

class galmag.Grid.Grid(box, resolution, grid_type='cartesian')

Bases: object

Defines a 3D grid object for a given choice of box dimensions and resolution.

Calling the attributes does the conversion between different coordinate systems automatically.

Parameters:

• box (3x2-array_like) – Box limits
• resolution (3-array_like) – containing the resolution along each axis.
• grid_type (str, optional) – Choice between ‘cartesian’, ‘spherical’ and ‘cylindrical’ uniform coordinate grids. Default: ‘cartesian’

coordinates

A dictionary contaning all the coordinates

cos_phi

\(\cos(\phi)\)

cos_theta

\(\cos(\theta)\)

get_prototype(dtype=None)

phi

Azimuthal coordinate, \(\phi\)

r_cylindrical

Cylindrical radial coordinate, \(s\)

r_spherical

Spherical radial coordinate, \(r\)

sin_phi

\(\sin(\phi)\)

sin_theta

\(\sin(\theta)\)

theta

Polar coordinate, \(\theta\)

x

Horizontal coordinate, \(x\)

y

Horizontal coordinate, \(y\)

z

Vertical coordinate, \(z\)

galmag.Observables module

galmag.disk_profiles module

Contains the definitions of the disk rotation curve, radial shear, alpha profile and disk scale height.

galmag.disk_profiles.Clemens_Milky_Way_rotation_curve(R, R_d=1.0, Rsun=8.5, normalize=True)

Rotation curve of the Milky Way obtained by Clemens (1985)

Parameters:

• R (float or array) – radial coordinate
• Rsun (float or array) – sun’s radius in kpc. Default: 8.5 kpc
• R_d (float) – unit of R in kpc [e.g. R_d=(disk radius in kpc) for r=0..1 within the disk]. Default: 1.0
• normalize (bool) – If True , with result normalized to unit at solar radius, if False, the result will be in km/s for R and Rsun in kpc

Returns:

rotation curve

Return type:

same as R

galmag.disk_profiles.Clemens_Milky_Way_shear_rate(R, R_d=1.0, Rsun=8.5, normalize=True)

Shear rate of the Milky Way based on the rotation curve obtained by Clemens (1985)

Parameters:

• R (float or array) – radial coordinate
• Rsun (float or array) – sun’s radius in kpc. Default: 8.5 kpc
• R_d (float) – unit of R in kpc [e.g. R_d=(disk radius in kpc) for r=0..1 within the disk]. Default: 1.0
• normalize (bool) – If True , with result normalized to unit at solar radius, if False, the result will be in km/s for R and Rsun in kpc

Returns:

shear rate profile, with

Return type:

same as R

galmag.disk_profiles.Omega(rotation_curve, R, Rsun=8.5, R_d=1.0, normalize=True, **kwargs)

Simple wrapper to avoid making mistakes when using dimensionless (aka normalized) quantities.

galmag.disk_profiles.constant_scale_height(R, h_d=1.0, R_d=1.0, Rsun=8.5)

Constant scale height for testing.

galmag.disk_profiles.constant_shear_rate(R, R_d=1.0, Rsun=8.5, S0=25, normalize=True)

Constant shear for testing.

\[V(R) = cte\]

Parameters:R (float or array) – cylindrical radius

galmag.disk_profiles.exponential_scale_height(R, h_d=1.0, R_HI=5, R_d=1.0, Rsun=8.5)

Exponential disk scale-heigh profile profile

\[h(R)=\exp\left(\frac{R-R_\odot}{R_{\rm sh}}\right)\]

Parameters:

• R (float or array) – radial coordinate
• h_d (float) – normalization of the scaleheight. Default: 1.0
• R_d (float) – unit of R in kpc [e.g. R_d=(disk radius in kpc) for r=0..1 within the disk]. Default: 1.0
• R_HI (float) – Parameter \(R_{\rm sh}\), characterizes how “flared” the disc is.
• Rsun (float) – Sun’s radius in kpc. Default: 8.5 kpc

Returns:

scale height normalized to h_d at the solar radius

Return type:

same as R

galmag.disk_profiles.regularize(r, Om, S, r_reg, Om_reg, k=4)

Avoids unphysically large values of Omega and Shear near the origin applying an exponential cutoff on Omega to prevent it.

\[\Omega(r) = \exp\left[-(r_\xi/r)^k\right]\left[ \tilde\Omega(r)-\Omega_\xi\right] + \Omega_\xi\,,\]

and

\[S(r) = e^{-(r_\xi/r)^k}\left\{k\left(\frac{r_\xi}{r}\right)^k\left[\tilde\Omega(r) -\Omega_\xi \right] +\tilde S \right\}\]

Parameters:

• R (array) – radial coordinate
• Om (float) – Angular velocity profile
• r_reg (float) – regularization radius
• Om_reg (float) – Value of Omega to be used for math:r lesssim r_{reg
• k (float) – How sharp is the cutoff. Default: 4

galmag.disk_profiles.simple_rotation_curve(R, R_d=1.0, Rsun=8.5, V0=220, normalize=True, fraction=0.029411764705882353)

Simple flat rotation curve

\[V = V_0 \left[1-\exp(-R/(f R_\odot)\right]\]

with fraction \(= f\), Rsun \(= R_\odot\), V0 \(= V_0\)

Parameters:

• R (float or array) – radial coordinate
• Rsun (float or array) – sun’s radius in kpc. Default: 8.5 kpc
• R_d (float) – unit of R in kpc [e.g. R_d=(disk radius in kpc) for r=0..1 within the disk]. Default: 1.0
• V0 (float) – Circular velocity at infinity (i.e. at the flat part). Default: 220 (km/s)
• fraction (float) – Fraction of the solar radius at which the profile decays exponentially. Default: 0.03 (i.e. 250 pc for Rsun=8.5).
• normalize (bool) – If True , with result normalized to unit at solar radius, if False, the result will be in units of V0.

Returns:

rotation curve

Return type:

same as R

galmag.disk_profiles.simple_shear_rate(R, R_d=1.0, Rsun=8.5, V0=220, normalize=True, fraction=0.029411764705882353)

A simple shear rate profile, compatible with the simple_rotation_curve() .

Parameters:

• R (float or array) – radial coordinate
• Rsun (float or array) – sun’s radius in kpc. Default: 8.5 kpc
• R_d (float) – unit of R in kpc [e.g. R_d=(disk radius in kpc) for r=0..1 within the disk]. Default: 1.0
• V0 (float) – Circular velocity at infinity (i.e. at the flat part). Default: 220 (km/s)
• fraction (float) – Fraction of the solar radius at which the profile decays exponentially. Default: 0.03 (i.e. 250 pc for Rsun=8.5).
• normalize (bool) – If True , with result normalized to unit at solar radius, if False, the result will be in km/s for R and Rsun in kpc

Returns:

shear rate

Return type:

same as R

galmag.disk_profiles.solid_body_rotation_curve(R, R_d=1.0, Rsun=8.5, V0=220, normalize=True)

Solid body rotation curve for testing.

\[V(R) = R\]

Parameters:R (float or array) – cylindrical radius

galmag.electron_profiles module

Distributions of electron densities and cosmic ray densities

galmag.electron_profiles.constant_ncr(rho, theta, phi, ncr0=1.0)

Dummy for constant cosmic ray electron density

galmag.electron_profiles.constant_ne(rho, theta, phi, ne0=1.0)

Dummy for constant cosmic ray electron density

galmag.electron_profiles.simple_ne(rho, theta, phi, ne0=1.0, Rs=2.6, R_d=1.0, Rsun=8.5, **kwargs)

Simple electron density profile

Comprises a exponential flared disc.

\[n = n_0 \exp\left[-\frac{r \sin(\theta)-s_\odot}{s_s}\right] \exp\left[\frac{-r \cos(\theta)}{h(s)}\right]\]

where \(h(s)\) is given by galmag.disk_profiles.exponential_scale_height()

Parameters:

• rho (array) – Spherical radial coordinate, \(r\)
• theta (array) – Polar coordinate, \(\theta\)
• phi (array) – Azimuthal coordinate, \(\phi\)
• ne0 (float) – Number density at Rsun. Default: 1.0
• R_d (float) – unit of rho in kpc [e.g. R_d=(disk radius in kpc) for r=0..1 within the disk]. Default: 1.0
• Rs (float) – Scale radius, \(s_s\), of the disc in kpc. Default: 2.6 kpc
• Rsun (float) – Reference radius, \(s_\odot\), in kpc. Default: 8.5 kpc

galmag.galerkin module

Contains functions which deal with the computation the coefficients associated with a Galerkin expansion of the solution of the mean field dynamo equation.

galmag.galerkin.Galerkin_expansion_coefficients(parameters, return_matrix=False, dtype=<class 'numpy.float64'>)

Calculates the Galerkin expansion coefficients.

First computes the transformation M defined by:

\[M_{ij} = W_{ij}, \text{ for } i \neq j\]

\[M_{ij} = \gamma_j, \text{ for } i=j\]

where:

\[W_{ij} = \int B_j \cdot \hat{W} B_i\]

Then, solves the eigenvalue/eigenvector problem using numpy.linalg.eig() , computing thus all the coefficients \(a_i\) and the growth rates (eigenvalues) \(\Gamma_i\).

Parameters:

• return_matrix (bool, optional) – If True, the matrix \(W_{ij}\) will be returned as well.
• p (dict) –

  A dictionary of parameters dictionary of parameters containing the parameters:

  • halo_Galerkin_ngrid -> Number of grid points used in the

    calculation of the Galerkin expansion

  • halo_rotation_function -> a function specifying the halo rotation

    curve

  • halo_alpha_function -> a function specifying the dependence of

    the alpha effect

  • halo_turbulent_induction -> \(R_{\alpha}\)
  • halo_rotation_induction -> \(R_{\omega}\)

  • halo_n_free_decay_modes -> number of free decay modes to be used

    in the expansion

  • halo_symmetric_field -> Symmetric or anti-symmetric field

    solution

  • halo_rotation_characteristic_radius -> turn-over radius of the flat

    rotation curve

  • halo_rotation_characteristic_height -> characteristic z used in some

    rotation curve prescriptions

Returns:

• val (array) – n-array containing growth rates (the eigenvalues of \(Mij\)).
• vec (array) – \(a_i\)’s: nx3 array containing the Galerkin coefficients associated.
• Wij (array) – The matrix \(W_{ij}\). Only present if return_matrix is True.

galmag.galerkin.perturbation_operator(r, theta, phi, Br, Bt, Bp, Vr, Vt, Vp, alpha, Ra, Ro, dynamo_type='alpha-omega')

Applies the perturbation operator associated with an dynamo to a magnetic field in uniform spherical coordinates.

Parameters:

• r/B/alpha/V (array) – position vector (not radius!), magnetic field, alpha and rotation curve, respectively, expressed as 3xNxNxN arrays containing the \(r\), \(\theta\) and \(\phi\) components in [0,…], [1,…] and [2,…], respectively.
• p (dict) – A dictionary of parameters containing ‘Ralpha_h’.

Returns:

A list containing NxNxN arrays corresponding to the 3 components of \(\hat W(\mathbf{B})\)

Return type:

list

galmag.halo_free_decay_modes module

Contains functions to compute the free decay modes of the magnetic of the halo of a galaxy.

They should be accessed through the function get_mode().

galmag.halo_free_decay_modes.get_B_a_1(r, theta, phi, C=0.346, k=3.141592653589793)

Computes the first (pure poloidal) antisymmetric free decay mode. Purely poloidal.

Parameters:

• r (array_like) – NxNxN array containing azimuthal coordinates.
• theta (array_like) – NxNxN array containing azimuthal coordinates.
• phi (array_like) – NxNxN array containing azimuthal coordinates.
• C (float, optional) – Mode normalization.
• k (float, optional) – \(k\)

Returns:

List containing three numpy arrays corresponding to: \(B_r\), \(B_\theta\), \(B_\phi\)

Return type:

list

galmag.halo_free_decay_modes.get_B_a_2(r, theta, phi, C=0.25, k=5.763)

Computes the second antisymmetric free decay mode - one of a degenerate pair, with eigenvalue \(\gamma_2=-(5.763)^2\). Purely poloidal.

Parameters:

• r (array_like) – NxNxN array containing azimuthal coordinates.
• theta (array_like) – NxNxN array containing azimuthal coordinates.
• phi (array_like) – NxNxN array containing azimuthal coordinates.
• C (float, optional) – Mode normalization.
• k (float, optional) – \(k\)

Returns:

List containing three numpy arrays corresponding to: \(B_r\), \(B_\theta\), \(B_\phi\)

Return type:

list

galmag.halo_free_decay_modes.get_B_a_3(r, theta, phi, C=3.445, k=5.763)

Computes the third antisymmetric free decay mode - one of a degenerate pair, with eigenvalue \(\gamma_2=-(5.763)^2\). Purely toroidal.

Parameters:

• r (array_like) – NxNxN array containing azimuthal coordinates.
• theta (array_like) – NxNxN array containing azimuthal coordinates.
• phi (array_like) – NxNxN array containing azimuthal coordinates.
• C (float, optional) – Mode normalization.
• k (float, optional) – \(k\)

Returns:

List containing three numpy arrays corresponding to: \(B_r\), \(B_\theta\), \(B_\phi\)

Return type:

list

galmag.halo_free_decay_modes.get_B_a_4(r, theta, phi, C=0.244, k=6.283185307179586)

Computes the forth antisymmetric free decay mode. Purely poloidal.

Parameters:

• r (array_like) – NxNxN array containing azimuthal coordinates.
• theta (array_like) – NxNxN array containing azimuthal coordinates.
• phi (array_like) – NxNxN array containing azimuthal coordinates.
• C (float, optional) – Mode normalization.
• k (float, optional) – \(k\)

Returns:

List containing three numpy arrays corresponding to: \(B_r\), \(B_\theta\), \(B_\phi\)

Return type:

list

galmag.halo_free_decay_modes.get_B_s_1(r, theta, phi, C=0.653646562698, k=4.493409457909)

Computes the first (poloidal) symmetric free decay mode. Purely poloidal.

Parameters:

• r (array_like) – NxNxN array containing azimuthal coordinates.
• theta (array_like) – NxNxN array containing azimuthal coordinates.
• phi (array_like) – NxNxN array containing azimuthal coordinates.
• C (float, optional) – Mode normalization.
• k (float, optional) – \(k\)

Returns:

List containing three numpy arrays corresponding to: \(B_r\), \(B_\theta\), \(B_\phi\)

Return type:

list

galmag.halo_free_decay_modes.get_B_s_2(r, theta, phi, C=1.32984358196, k=4.493409457909)

Computes the second symmetric free decay mode (one of a degenerate pair, with eigenvalue gamma_2=-(5.763)^2 . Purely toroidal.

Parameters:

• r (array_like) – NxNxN array containing azimuthal coordinates.
• theta (array_like) – NxNxN array containing azimuthal coordinates.
• phi (array_like) – NxNxN array containing azimuthal coordinates.
• C (float, optional) – Mode normalization.
• k (float, optional) – \(k\)

Returns:

List containing three numpy arrays corresponding to: \(B_r\), \(B_\theta\), \(B_\phi\)

Return type:

list

galmag.halo_free_decay_modes.get_B_s_3(r, theta, phi, C=0.0169610298034, k=6.987932000501)

Computes the third symmetric free decay mode. Purely poloidal

Parameters:

• r (array_like) – NxNxN array containing azimuthal coordinates.
• theta (array_like) – NxNxN array containing azimuthal coordinates.
• phi (array_like) – NxNxN array containing azimuthal coordinates.
• C (float, optional) – Mode normalization.
• k (float, optional) – \(k\)

Returns:

List containing three numpy arrays corresponding to: \(B_r\), \(B_\theta\), \(B_\phi\)

Return type:

list

galmag.halo_free_decay_modes.get_B_s_4(r, theta, phi, C=0.539789362061, k=6.987932000501)

Computes the fourth symmetric free decay mode. Purely toroidal.

Parameters:

• r (array_like) – NxNxN array containing azimuthal coordinates.
• theta (array_like) – NxNxN array containing azimuthal coordinates.
• phi (array_like) – NxNxN array containing azimuthal coordinates.
• C (float, optional) – Mode normalization.
• k (float, optional) – \(k\)

Returns:

List containing three numpy arrays corresponding to: \(B_r\), \(B_\theta\), \(B_\phi\)

Return type:

list

galmag.halo_free_decay_modes.get_mode(r, theta, phi, n_mode, symmetric)

Computes the n_mode’th free decay mode.

(This is actually a wrapper invoking other functions)

Parameters:

• r (array_like) – NxNxN array containing azimuthal coordinates.
• theta (array_like) – NxNxN array containing azimuthal coordinates.
• phi (array_like) – NxNxN array containing azimuthal coordinates.
• nmode (int) – index of the mode to be computed
• symmetric – If True, the mode is assumed to be symmetric over the midplane (i.e. quadrupolar). If False the mode is assumed to be antisymmetric (dipolar). If None (or anything else), than the mode is drawn from a mixture of symmetric and antisymmetric modes.

Note

At the moment, only the first 4 symmetric and first 4 antisymmetric modes are available. Later versions will support arbitrary choice of mode.

Returns:List containing three numpy arrays corresponding to: \(B_r\), \(B_\theta\), \(B_\phi\) Return type:list

class galmag.halo_free_decay_modes.xi_lookup_table(filepath='.xilookup.npy', regenerate=False, **kwargs)

Bases: object

Stores a look-up table of the roots of the equation

\[J_{n-1/2}(\xi_{nl}) J_{n+1/2}(\xi_{nl}) = 0\]

which can be accessed through the method get_xi().

These are related to the decay rates through: \(\gamma_{nl} = -(\xi_{nl})^2\)

which can be access through the method get_gamma().

Parameters:

• filepath (str) – Path to the file where the lookup table will be stored. Default: ‘.xilookup.npy’
• regenerate (bool) – If True, the original file on hard disk will be overwritten.

generate_xi_lookup_table(max_n=4, max_l=4, number_of_guesses=150, max_guess=25, save=True)

Produces a (max_n,max_l)-array containing containing the roots of

\[J_{n-1/2}(\xi_{nl}) J_{n+1/2}(\xi_{nl}) = 0\]

These are related to the decay rates through:

\[\gamma_{nl} = -(\xi_{nl})^2\]

The root are found using the mpmath.findroot function. The search for roots supplying findroot with number_of_guesses initial guesses uniformly distributed in the interval [3,max_guess].

Parameters:

• max_l (max_n,) – Maximum values for n and l in the table. Default: 4
• number_of_guesses (int) – Number of uniformly distributed guesses. Default 150
• max_guess (float) – Maximum guess value for root.
• save – Saves the lookup table to file

Returns:

A (max_n,max_l)-array containing containing the roots of

Return type:

numpy.ndarray

get_gamma(n, l, **kwargs)

Returns the free-decay-mode decay rate, \(\gamma_{nl}\)

Parameters:l (n,) – indices of the free decay mode Returns:\(\gamma_{nl}\) Return type:float

get_xi(n, l, regenerate=False, **kwargs)

Computes or reads the root of \(J_{n-1/2}(\xi_{nl}) J_{n+1/2}(\xi_{nl}) = 0\)

Parameters:

• l (n,) – indices of the free decay mode
• regenerate – (Default value = False)
• **kwargs –

Returns:

\(\xi_{nl}\)

Return type:

float

galmag.halo_profiles module

GalMag

Contains the definitions of the halo rotation curve and alpha profile.

galmag.halo_profiles.simple_V(rho, theta, phi, r_h=1.0, Vh=220, fraction=0.2, normalize=True, fraction_z=None, legacy=False)

Simple form of the rotation curve to be used for the halo

\[V(r,\theta,\phi) \propto [1-\exp(-r \sin(\theta) / s_v) ]\]

Note

This simple form has no z dependence

Parameters:

• rho (array) – Spherical radial coordinate, \(r\)
• theta (array) – Polar coordinate, \(\theta\)
• phi (array) – Azimuthal coordinate, \(\phi\)
• fraction – fraction of the halo radius corresponding to the turnover of the rotation curve.
• r_h – halo radius in the same units as rho. Default: 1.0
• Vh – Value of the rotation curve at rho=r_h. Default: 220 km/s
• normalize (bool, optional) – if True, the rotation curve will be normalized to one at rho=r_h

Returns:

List containing three numpy arrays corresponding to: \(V_r\), \(V_\theta\), \(V_\phi\)

Return type:

list

galmag.halo_profiles.simple_V_exp(rho, theta, phi, r_h=1.0, Vh=220, fraction=0.2, fraction_z=0.7333333333333333, normalize=True, legacy=False)

Variation on simple_V which decays exponentially with z

\[V(r,\theta,\phi) \propto (1-\exp(-r \sin(\theta) / s_v)) \exp(-r \cos(\theta)/z_v)\]

Parameters:

• rho (array) – Spherical radial coordinate, \(r\)
• theta (array) – Polar coordinate, \(\theta\)
• phi (array) – Azimuthal coordinate, \(\phi\)
• fraction – fraction of the halo radius corresponding to the turnover of the rotation curve.
• fraction_z – fraction of the halo radius corresponding to the characteristic vertical decay length of the rotation
• r_h – halo radius in the same units as rho. Default: 1.0
• Vh – Value of the rotation curve at rho=r_h. Default: 220 km/s
• normalize (bool, optional) – if True, the rotation curve will be normalized to one at rho=r_h

Returns:

List containing three numpy arrays corresponding to: \(V_r\), \(V_\theta\), \(V_\phi\)

Return type:

list

galmag.halo_profiles.simple_V_legacy(rho, theta, phi, r_h=1.0, Vh=220, fraction=0.5, fraction_z=None, normalize=True)

Rotation curve employed in version 0.1 and in the MMath Final Report of James Hollins. Same as simple_V but with a slight change in the way it is normalized.

galmag.halo_profiles.simple_V_linear(rho, theta, phi, r_h=1.0, Vh=220, fraction=0.2, fraction_z=0.7333333333333333, normalize=True, legacy=False)

Variation on simple_V which decays linearly with z, reaching 0 at z=(halo radius), and V_h at z=0

\[V(r,\theta,\phi) \propto [1-\exp(-r \sin(\theta) / s_v)] (1-z/z_v)\]

Parameters:

• rho (array) – Spherical radial coordinate, \(r\)
• theta (array) – Polar coordinate, \(\theta\)
• phi (array) – Azimuthal coordinate, \(\phi\)
• fraction – fraction of the halo radius corresponding to the turnover of the rotation curve. (s_v = fraction*r_h)
• fraction_z – fraction of the halo radius controling the “lag” of the rotation curve. (z_v = fraction_z*r_h)
• r_h – halo radius in the same units as rho. Default: 1.0
• Vh – Value of the rotation curve at rho=r_h. Default: 220 km/s
• normalize (bool, optional) – if True, the rotation curve will be normalized to one at rho=r_h
• Returns – List containing three numpy arrays corresponding to: \(V_r\), \(V_\theta\), \(V_\phi\)

galmag.halo_profiles.simple_alpha(rho, theta, phi, alpha0=1.0)

Simple profile for alpha

\[\alpha(\mathbf{r}) = \alpha_0\cos(\theta)\]

Parameters:

• rho (array) – Spherical radial coordinate, \(r\)
• theta (array) – Polar coordinate, \(\theta\)
• phi (array) – Azimuthal coordinate, \(\phi\)
• alpha0 (float, optional) – Normalization. Default: 1.0

galmag.util module

Auxiliary functions.

galmag.util.curl_spherical(rr, tt, pp, Br, Bt, Bp, order=4)

Computes the curl of a vector in spherical coordinates.

Parameters:

• rr/tt/pp (array_like) – NxNxN arrays containing the \(r\), \(\theta\) and \(\phi\) coordinates
• Br/Bt/Bp (array_like) – NxNxN arrays \(r\), \(\theta\) and \(\phi\) components of the vector in the same coordinate grid.

Returns:

three NxNxN arrays containing the components of the curl.

Return type:

list

galmag.util.derive(V, dx, axis=0, order=2)

Computes the numerical derivative of a function specified over a 3 dimensional uniform grid. Uses second order finite differences.

Obs: extremities will use forward or backwards finite differences.

Parameters:

• V (array_like) – NxNxN array
• dx (float) – grid spacing
• axis (int) – specifies over which axis the derivative should be performed. Default: 0.
• order (int) – order of the finite difference method. Default: 2

Returns:

The derivative, dV/dx

Return type:

same as V

galmag.util.get_max_jobs()

galmag.util.simpson(f, r)

Integrates over the last axis

galmag.version module

GalMag version

Stores version number

Module contents

GalMag - A Python tool for computing realistic galactic magnetic fields

Generates realistic galactic magnetic field based on the dynamo equation

Example

Initializes a B_field object on a uniform cartesian 100x100x100 grid:

from galmag import B_field

box_limits = [[-15, 15],[-15, 15],[-15, 15]] # kpc
box_resolution = [100,100,100]
B = B_field(box_limits, box_resolution)

Adds a disc magnetic field with reversals at 4.7 kpc and 12.25 kpc and a field strength of -3 microgauss close to the Sun:

B.add_disk_field(reversals=[4.7,12.25] ,
                 B_phi_solar_radius=-3, # muG
                 number_of_modes=number_of_modes)

Adds a halo magnetic field with intensity 0.1 microgauss close to the Sun:

B.add_halo_field(halo_ref_Bphi=Bphi_sun)

Each coordinate component of the produced magnetic can then be accessed through the relevant coordinate attributes, e.g.:

B.r_spherical
B.theta

For a quick introduction to usage, please check the tutorial jupyter notebook.

• Index
• Module Index
• Search Page