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(x, y, z, output)	A single magnetic field line.
Grid(ns, nphi, nr, rss)	Grid on which the solution is calculated.
Input(br, nr, rss)	Input to PFSS modelling.
Output(alr, als, alp, grid)	Output of PFSS modelling.

Changelog

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('gong.fits') and not os.path.exists('gong.fits.gz'):
    import urllib.request
    urllib.request.urlretrieve(
        'https://gong2.nso.edu/oQR/zqs/201901/mrzqs190108/mrzqs190108t1114c2212_050.fits.gz',
        'gong.fits.gz')

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

We can now use SunPy to load the .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!

map = sunpy.map.Map('gong.fits')
br = map.data
br = br - np.mean(br)

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')
ax.set_aspect('equal')

# 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

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