Magnetic field along a field line

How to get the value of the magnetic field along a field line traced through the PFSS solution.

First, import required modules

import astropy.constants as const
import astropy.units as u
from astropy.coordinates import SkyCoord
import matplotlib.pyplot as plt
import numpy as np
import sunpy.map

import pfsspy
from pfsspy import tracing
from pfsspy.sample_data import get_gong_map

Load a GONG magnetic field map. If ‘gong.fits’ is present in the current directory, just use that, otherwise download a sample GONG map.

gong_fname = get_gong_map()

Out:

Files Downloaded:   0%|          | 0/1 [00:00<?, ?file/s]

mrzqs200901t1304c2234_022.fits.gz:   0%|          | 0.00/242k [00:00<?, ?B/s]

mrzqs200901t1304c2234_022.fits.gz:   0%|          | 100/242k [00:00<08:00, 502B/s]

mrzqs200901t1304c2234_022.fits.gz:  83%|########3 | 201k/242k [00:00<00:00, 824kB/s]
Files Downloaded: 100%|##########| 1/1 [00:00<00:00,  1.82file/s]
Files Downloaded: 100%|##########| 1/1 [00:00<00:00,  1.82file/s]


                                                                                    

We can now use SunPy to load the GONG fits file, and extract the magnetic field data.

The mean is subtracted to remove the monopole component

gong_map = sunpy.map.Map(gong_fname)
# Remove the mean
gong_map = sunpy.map.Map(gong_map.data - np.mean(gong_map.data), gong_map.meta)

The PFSS solution is calculated on a regular 3D grid in (phi, s, rho), where rho = ln(r), and r is the standard spherical radial coordinate. We need to define the number of rho grid points, and the source surface radius.

nrho = 35
rss = 2.5

From the boundary condition, number of radial grid points, and source surface, we now construct an Input object that stores this information

input = pfsspy.Input(gong_map, nrho, rss)
output = pfsspy.pfss(input)

Now take a seed point, and trace a magnetic field line through the PFSS solution from this point

tracer = tracing.PythonTracer()
r = 1.2 * const.R_sun
lat = 70 * u.deg
lon = 0 * u.deg

seeds = SkyCoord(lon, lat, r, frame=output.coordinate_frame)
field_lines = tracer.trace(seeds, output)

From this field line we can extract the coordinates, and then use .b_along_fline to get the components of the magnetic field along the field line.

From the plot we can see that the non-radial component of the mangetic field goes to zero at the source surface, as expected.

field_line = field_lines[0]
B = field_line.b_along_fline
r = field_line.coords.radius
fig, ax = plt.subplots()

ax.plot(r.to(const.R_sun), B[:, 0], label=r'$B_{r}$')
ax.plot(r.to(const.R_sun), B[:, 1], label=r'$B_{\theta}$')
ax.plot(r.to(const.R_sun), B[:, 2], label=r'$B_{\phi}$')
ax.legend()
ax.set_xlabel(r'r / $R_{\odot}$')
ax.set_ylabel(f'B ({B.unit})')

plt.show()
plot field line magnetic field

Total running time of the script: ( 0 minutes 5.604 seconds)

Gallery generated by Sphinx-Gallery