pfsspy is a python package for carrying out Potential Field Source Surface modelling. For more information on the actually PFSS calculation see this document.

Note

pfsspy is a very new package, so elements of the API are liable to change with the first few releases. If you find any bugs or have any suggestions for improvement, please raise an issue here: https://github.com/dstansby/pfsspy/issues

Installing

pfsspy can be installed from PyPi using

pip install pfsspy

Improving performance

pfsspy automatically detects an installation of numba, which compiles some of the numerical code to speed up pfss calculations. To enable this simply install numba and use pfsspy as normal.

Citing

If you use pfsspy in work that results in publication, please cite the archived code at both

• https://zenodo.org/record/2566462
• https://zenodo.org/record/1472183

Citation details can be found at the lower right hand of each web page.

Code reference

For the main user-facing code and a changelog see

pfsspy Package

Functions

load_output(file)	Load a saved output file.
pfss(input)	Compute PFSS model.

Classes

FieldLine(*args[, copy])	A single magnetic field line.
Grid(ns, nphi, nr, rss)	Grid on which the solution is calculated.
Input(br, nr, rss[, dtime])	Input to PFSS modelling.
Output(alr, als, alp, grid[, dtime])	Output of PFSS modelling.

Changelog

0.2.0

• pfsspy.Input and pfsspy.Output now take the optiona keyword argument dtime, which stores the datetime on which the magnetic field measurements were made. This is then propagated to the obstime attribute of computed field lines, allowing them to be transformed in to coordinate systems other than Carrignton frames.
• pfsspy.FieldLine no longer overrrides the SkyCoord __init__; this should not matter to users, as FieldLine objects are constructed internally by calling pfsspy.Output.trace()

0.1.5

• Output.plot_source_surface now accepts keyword arguments that are given to Matplotlib to control the plotting of the source surface.

0.1.4

• Added more explanatory comments to the examples
• Corrected the dipole solution calculation
• Added pfsspy.coords.sph2cart() to transform from spherical to cartesian coordinates.

0.1.3

• pfsspy.Output.plot_pil() now accepts keyword arguments that are given to Matplotlib to control the style of the contour.
• pfsspy.FieldLine.expansion_factor is now cached, and is only calculated once if accessed multiple times.

for usage examples see

pfsspy examples

Dipole source solution

Dipole source solution

A simple example showing how to use pfsspy to compute the solution to a dipole source field.

First, import required modules

import astropy.constants as const
import matplotlib.pyplot as plt
import matplotlib.patches as mpatch
import numpy as np
import pfsspy
import pfsspy.coords as coords

To start with we need to construct an input for the PFSS model. To do this, first set up a regular 2D grid in (phi, s), where s = cos(theta) and (phi, theta) are the standard spherical coordinate system angular coordinates. In this case the resolution is (360 x 180).

nphi = 360
ns = 180

phi = np.linspace(0, 2 * np.pi, nphi)
s = np.linspace(-1, 1, ns)
s, phi = np.meshgrid(s, phi)

Now we can take the grid and calculate the boundary condition magnetic field.

def dipole_Br(r, s):
    return 2 * s / r**3


br = dipole_Br(1, s).T

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 = 50
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(br, nrho, rss)

Using the Input object, plot the input field

fig, ax = plt.subplots()
mesh = input.plot_input(ax)
fig.colorbar(mesh)
ax.set_title('Input dipole field')

Now calculate the PFSS solution.

output = pfsspy.pfss(input)

Using the Output object we can plot the source surface field, and the polarity inversion line.

fig, ax = plt.subplots()
mesh = output.plot_source_surface(ax)
fig.colorbar(mesh)
output.plot_pil(ax)
ax.set_title('Source surface magnetic field')

Finally, using the 3D magnetic field solution we can trace some field lines. In this case 32 points equally spaced in theta are chosen and traced from the source surface outwards.

fig, ax = plt.subplots()
ax.set_aspect('equal')

# Take 32 start points spaced equally in theta
r = 1.01
phi = np.pi / 2
for theta in np.linspace(0, np.pi, 33):
    x0 = coords.sph2cart(r, theta, phi)
    field_line = output.trace(np.array(x0))
    color = {0: 'black', -1: 'tab:blue', 1: 'tab:red'}.get(field_line.polarity)
    ax.plot(field_line.y / const.R_sun,
            field_line.z / const.R_sun, color=color)

# Add inner and outer boundary circles
ax.add_patch(mpatch.Circle((0, 0), 1, color='k', fill=False))
ax.add_patch(mpatch.Circle((0, 0), input.grid.rss, color='k', linestyle='--',
                           fill=False))
ax.set_title('PFSS solution for a dipole source field')
plt.show()

GONG PFSS extrapolation

GONG PFSS extrapolation

Calculating PFSS solution for a GONG synoptic magnetic field map.

First, import required modules

import os
import astropy.constants as const
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import sunpy.map

import pfsspy
import pfsspy.coords as coords

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

if not os.path.exists('190310t0014gong.fits') and not os.path.exists('190310t0014gong.fits.gz'):
    import urllib.request
    urllib.request.urlretrieve(
        'https://gong2.nso.edu/oQR/zqs/201903/mrzqs190310/mrzqs190310t0014c2215_333.fits.gz',
        '190310t0014gong.fits.gz')

if not os.path.exists('190310t0014gong.fits'):
    import gzip
    with gzip.open('190310t0014gong.fits.gz', 'rb') as f:
        with open('190310t0014gong.fits', 'wb') as g:
            g.write(f.read())

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

The mean is subtracted to enforce div(B) = 0 on the solar surface: n.b. it is not obvious this is the correct way to do this, so use the following lines at your own risk!

[[br, header]] = sunpy.io.fits.read('190310t0014gong.fits')
br = br - np.mean(br)

GONG maps have their LH edge at -180deg in Carrington Longitude, so roll to get it at 0deg. This way the input magnetic field is in a Carrington frame of reference, which matters later when lining the field lines up with the AIA image.

br = np.roll(br, header['CRVAL1'] + 180, axis=1)

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 = 60
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(br, nrho, rss)

Using the Input object, plot the input field

fig, ax = plt.subplots()
mesh = input.plot_input(ax)
fig.colorbar(mesh)
ax.set_title('Input field')

Now calculate the PFSS solution, and plot the polarity inversion line.

output = pfsspy.pfss(input)
output.plot_pil(ax)

Using the Output object we can plot the source surface field, and the polarity inversion line.

fig, ax = plt.subplots()
mesh = output.plot_source_surface(ax)
fig.colorbar(mesh)
output.plot_pil(ax)
ax.set_title('Source surface magnetic field')

Finally, using the 3D magnetic field solution we can trace some field lines. In this case 256 points equally gridded in theta and phi are chosen and traced from the source surface outwards.

fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

# Loop through 16 values in theta and 16 values in phi
r = 1.01
for theta in np.linspace(0, np.pi, 17):
    for phi in np.linspace(0, 2 * np.pi, 17):
        x0 = np.array(coords.sph2cart(r, theta, phi))
        field_line = output.trace(x0)
        color = {0: 'black', -1: 'tab:blue', 1: 'tab:red'}.get(field_line.polarity)
        ax.plot(field_line.x / const.R_sun,
                field_line.y / const.R_sun,
                field_line.z / const.R_sun,
                color=color, linewidth=1)

ax.set_title('PFSS solution')
plt.show()

# sphinx_gallery_thumbnail_number = 3

Overplotting field lines on AIA maps

Overplotting field lines on AIA maps

This example shows how to take a PFSS solution, trace some field lines, and overplot the traced field lines on an AIA 193 map.

First, we import the required modules

from datetime import datetime
import os

import astropy.units as u
import matplotlib.pyplot as plt
import numpy as np
import sunpy.map
import sunpy.io.fits

import pfsspy
import pfsspy.coords as coords

Load a GONG magnetic field map. The map date is 10/03/2019

if not os.path.exists('190310t0014gong.fits') and not os.path.exists('190310t0014gong.fits.gz'):
    import urllib.request
    urllib.request.urlretrieve(
        'https://gong2.nso.edu/oQR/zqs/201903/mrzqs190310/mrzqs190310t0014c2215_333.fits.gz',
        '190310t0014gong.fits.gz')

if not os.path.exists('190310t0014gong.fits'):
    import gzip
    with gzip.open('190310t0014gong.fits.gz', 'rb') as f:
        with open('190310t0014gong.fits', 'wb') as g:
            g.write(f.read())

Load the corresponding AIA 193 map

if not os.path.exists('AIA20190310.fits'):
    import urllib.request
    urllib.request.urlretrieve(
        'http://jsoc2.stanford.edu/data/aia/synoptic/2019/03/10/H0000/AIA20190310_0000_0193.fits',
        'AIA20190310.fits')

aia = sunpy.map.Map('AIA20190310.fits')
dtime = aia.date

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

The mean is subtracted to enforce div(B) = 0 on the solar surface: n.b. it is not obvious this is the correct way to do this, so use the following lines at your own risk!

[[br, header]] = sunpy.io.fits.read('190310t0014gong.fits')
br = br - np.mean(br)

GONG maps have their LH edge at -180deg in Carrington Longitude, so roll to get it at 0deg. This way the input magnetic field is in a Carrington frame of reference, which matters later when lining the field lines up with the AIA image.

br = np.roll(br, header['CRVAL1'] + 180, axis=1)

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 grid points in rho, and the source surface radius.

nrho = 60
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(br, nrho, rss, dtime=dtime)

Using the Input object, plot the input photospheric magnetic field

fig, ax = plt.subplots()
mesh = input.plot_input(ax)
fig.colorbar(mesh)
ax.set_title('Input field')

We can also plot the AIA map to give an idea of the global picture. There is a nice active region in the top right of the AIA plot, that can also be seen in the top left of the photospheric field plot above.

ax = plt.subplot(1, 1, 1, projection=aia)
aia.plot(ax)

Now we construct a 10 x 10 grid of footpoitns to trace some magnetic field lines from.

The figure shows a zoom in of the magnetic field map, with the footpoints overplotted. The footpoints are centered around the active region metnioned above.

s, phi = np.meshgrid(np.linspace(0.1, 0.2, 5),
                     np.deg2rad(np.linspace(55, 65, 5)))

fig, ax = plt.subplots()
mesh = input.plot_input(ax)
fig.colorbar(mesh)
ax.scatter(np.rad2deg(phi), s, color='k', s=1)

ax.set_xlim(50, 70)
ax.set_ylim(0, 0.35)
ax.set_title('Field line footpoints')

Compute the PFSS solution from the GONG magnetic field input

output = pfsspy.pfss(input)

Trace field lines from the footpoints defined above. pfsspy.coords is used to convert the s, phi cooridnates into the cartesian coordinates that are needed by the tracer.

flines = []
for s, phi in zip(s.ravel(), phi.ravel()):
    x0 = np.array(pfsspy.coords.strum2cart(0.01, s, phi))
    flines.append(output.trace(x0, atol=1e-8))

Plot the input GONG magnetic field map, along with the traced mangetic field lines.

fig, ax = plt.subplots()
mesh = input.plot_input(ax)
for fline in flines:
    fline.representation_type = 'spherical'
    ax.plot(fline.lon / u.deg, np.sin(fline.lat), color='black', linewidth=1)

ax.set_xlim(55, 65)
ax.set_ylim(0.1, 0.25)
ax.set_title('Photospheric field and traced field lines')

Plot the AIA map, along with the traced magnetic field lines. Inside the loop the field lines are converted to the AIA observer coordinate frame, and then plotted on top of the map.

fig = plt.figure()
ax = plt.subplot(1, 1, 1, projection=aia)
transform = ax.get_transform('world')
aia.plot(ax)
for fline in flines:
    fline = fline.transform_to(aia.coordinate_frame)
    Tx = fline.Tx.to(u.deg)
    Ty = fline.Ty.to(u.deg)
    ax.plot(Tx, Ty, transform=transform,
            alpha=0.8, linewidth=1, color='black')

ax.set_xlim(500, 900)
ax.set_ylim(400, 800)
plt.show()

# sphinx_gallery_thumbnail_number = 5

and for the helper modules (behind the scense!) see

Helper modules

pfsspy.plot Module

Functions

contour(phi, costheta, field, levels[, ax])	Parameters:
radial_cut(phi, costheta, field[, ax])

pfsspy.coords Module

Helper functions for coordinate transformations used in the PFSS domain.

The PFSS solution is calculated on a “strumfric” grid defined by

• \(\rho = \log (r)\)
• \(s = \cos (\theta )\)
• \(\phi\)

where \(r, \theta, \phi\) are spherical cooridnates that have ranges

• \(1 < r < r_{ss}\)
• \(0 < \theta < \pi\)
• \(0 < \phi < 2\pi\)

The transformation between cartesian coordinates used by the tracer and the above coordinates is given by

• \(x = r\sin (\theta) \cos (\phi)\)
• \(y = r\sin (\theta) \sin (\phi)\)
• \(z = r \cos (\theta)\)

Functions

cart2strum(x, y, z)	Convert cartesian coordinates to strumfric coordinates.
sph2cart(r, theta, phi)	Convert spherical coordinates to cartesian coordinates.
strum2cart(rho, s, phi)	Convert strumfric coordinates to cartesian coordinates.

Indices and tables

• Index
• Module Index
• Search Page