# 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
import matplotlib.pyplot as plt
import sunpy.map
from astropy.coordinates import SkyCoord

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


Load a GONG magnetic field map

gong_fname = get_gong_map()
gong_map = sunpy.map.Map(gong_fname)


Out:

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

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

pfsspy.mrzqs200901t1304c2234_022.fits.gz:   0%|          | 100/242k [00:00<07:05, 568B/s]



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


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

pfss_in = pfsspy.Input(gong_map, nrho, rss)
pfss_out = pfsspy.pfss(pfss_in)


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

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

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


Out:

INFO: Missing metadata for solar radius: assuming the standard radius of the photosphere. [sunpy.map.mapbase]


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
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()


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

Gallery generated by Sphinx-Gallery