Field lines tutorial

Generating field lines and plotting them in 3D

Make sure that the pyRich python module is in your PYTHONPATH.

This notebook and the input files for the fortran codes are in the /tutorial/FieldLines folder in the main repository.

import pyRich.field as fl
import matplotlib.pyplot as plt
import numpy as np
%matplotlib notebook

Creating the field lines with the fortran code

In your fortran code, do a make all in the mag/field folder, to create executables for the Build_ programs.

It is a good idea to define an environment variable that contains the path to the executables (in the example below, $RICH points to the repo main folder).

1. Creating the potential field file

The first step is to create a potential field (add refs). The input file for this example is build_potential_dipole.in, which sets up a dipolar field tilted at 45 degrees with respect to the z-axis (which is the rotational axis in later codes).

$RICH/mag/field/build_pot_field < build_potential_dipole.in

will create pf_beta45.h5

2. Creating the field lines file

The second step is to create a set of field lines. The input file for this example is build_field_lines.in. The lines are randomly selected at the surface of the star, with the option of weighting the probability according to the surface field. Other options are available, see (link to the documentation)

$RICH/mag/field/build_field_lines < build_field_lines.in

will create fl_beta45.h5

Loading the field_lines h5 file in a python object

Here we use the field module of pyRich. A field_lines object contains a set of field lines (see pyRich.field documentation).

You can access the field lines informations by accessing objects in the class.

# Read the h5 file containing the set of field lines
lines = fl.load_field_lines('fl_beta45.h5')

# print the number of field lines contained in the field_lines object
print('there are {} field lines'.format(lines.n_lines))

there are 128 field lines

Making an interactive 3D graph of the field lines

The class has a built-in method to create an interactive 3D graph.

The cartesian axes are in the main frame of reference of the fortran code (usually taken to be the rotational frame of reference in the later on codes, with the z-axis being the rotational frame). But if you are using the field tools as stand-alone, the main frame of reference is your choice.

In the example above, the polar axis was tilted by 45\(^\circ\) with respect to the z-axis of the main frame of reference

lines = fl.load_field_lines('fl_beta45.h5')

# Using the build-in function that returns a plot of the field lines in 3D.
fig, ax = lines.plot(star=True)

Making a 2D graph in the plane of the sky

For this nuilt-in function, the z-axis of the main frame of reference is taken to be the rotational axis.

The graph function makes a transformation to the line-of-sight (LOS) frame of reference, for an inclination of the rotational axis \(i\) and a rotational phase \(\phi\).

The rotation is a right-hand rule about the rotational axis.

In the figure, the red line shows the rotational axis, and the green line shows the rotational x-direction.

lines = fl.load_field_lines('fl_beta45.h5')

# The inclination (here 45) is given in degrees
# The rotational phase (here 0.3) is given in the range [0,1]

fig, ax = lines.plot_2D(45, phi_one=0.3)

## the figure is returned, so that you can make some changes
# e.g.:

ax.set_ylim(-6,6)
ax.set_xlim(-6,6)

not animated

(-6.0, 6.0)

Making a 2D animation in the plane of the sky

There is also an option to use plot_2D to create an animation showing the rotation on the plane of the sky.

lines = fl.load_field_lines('fl_beta45.h5')

# The inclination (here 45) is given in degrees

anim = lines.plot_2D(45, animated=True)

# This is to animate on a jupyter notebook.
# If you are running from a script, this is not needed.
# Aslo to run without error, you may want to make sure that ffmeg is installed
# conda install -c conda-forge ffmpeg

#from IPython.display import HTML
#HTML(anim.to_html5_video())

# To save as a mpeg
#from matplotlib.animation import FFMpegWriter
#anim.save('fl_beta45_incl45.mp4', writer=FFMpegWriter(fps=60))

This markdown cell is to render the animation in the documentation and can be ignored if you are using the notebook.

This first line is necessary to copy the video in the doc-src when nbsphinx runs:

fl_beta45_incl45.mp4

Then the html command will render it embedded:

Using the field line data

This is an example on how the line data can be accessed for other science task

Here, let’s make a graph of 1. the magnetic field magnitude along a field line 2. the magnetic field magnetidue along the field line after we have transformed the loop to be in the observer’s frame of reference for i=45 and phase=0.0 3. the magnetic field component along the field loop in the direction of the observer for i=45 and phase=0.0

lines = fl.load_field_lines('fl_beta45.h5')

fig, ax = plt.subplots(1,1)

# Let's select the 42nd line in the ensemble (just because :)
n = 42

# just get the field line of interest
line = lines.field_line[n]

# find the apex of the field line
# the get_r function returns the radial coordiate of the points on the field line.
n = np.where( line.get_r() == np.max(line.get_r()) )

# Thus the apex of the line is at s:
s_apex = line.s[n[0][0]]
# add a dashed vertical line at s=0
ax.axvline(x=0.0, ls='dotted', c='k')

# Let's make a graph of field strength along the loop,
# setting s=0 to be at the loop apex.
# The get_B function returns the magnetitude of the field along the field line.
ax.plot( line.s-s_apex, line.get_B(), label='Field strength' )

ax.set_xlabel('Coordiate along the loop, centered at the loop apex')
ax.set_ylabel('Magnetic field strength (G)')

# Now let's look at the line-of sight component of the magnetic field along the line.
# We can use the rotate function to transform from
# the fortran code main frame of reference (here assuming rotational, see above)
# Let's use here inclination of 45 degree and a rotational phase of zero

LOS_line = line.rotate(45*np.pi/180, 0.0) # this makes a copy

# If we plot the field strength along the loop again, we should get the same curve as before
ax.plot( LOS_line.s-s_apex, LOS_line.get_B(), ls='--', label='Test field strength after transformation' )

# Now let's graph the z component of this rotated loop (where z is now in the direction of the observer)
# The Bx/By/Bz are in the vB array [n,3]
ax.plot( LOS_line.s-s_apex, LOS_line.vB[:,2], label='Longitudinal field for i=90, phase=0.0' )

ax.legend(loc=0)

<matplotlib.legend.Legend at 0x7fd6b4e6c460>

Modifying the object to only have the field loop graphed above

lines = fl.load_field_lines('fl_beta45.h5')

n = 42

# We can use the field_line constructor to make a copy of "lines",
# But with a single field loop of our choice.
one_line = fl.field_lines(1, lines.pot_field, lines.star, [lines.field_line[n]])

anim = one_line.plot_2D(45, animated=True)