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

Code reference

For the main user-facing code 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.

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

Set up a 1degree by 1degree grid in theta and phi

nphi = 360
ntheta = 180

phi = np.linspace(0, 2 * np.pi, nphi)
theta = np.linspace(-np.pi / 2, np.pi / 2, ntheta)
theta, phi = np.meshgrid(theta, phi)

Define the number of radial grid points and the source surface radius

nr = 50
rss = 2.5

Compute radial component ofa dipole field

def dipole_Br(r, theta):
    return 2 * np.sin(theta) / r**3


br = dipole_Br(1, theta).T

Create PFSS input object

input = pfsspy.Input(br, nr, rss)

Plot input magnetic field

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

Calculate PFSS solution

output = pfsspy.pfss(input)

Plot output field

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

Trace some field lines

br, btheta, bphi = output.bg

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

# Take 32 start points spaced equally in theta
r = 1.01
for theta in np.linspace(0, np.pi, 33):
    x0 = np.array([0, r * np.sin(theta), r * np.cos(theta)])
    field_line = output.trace(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 pfsspy
import sunpy.map

If a gong magnetic field map isn’t present, download one

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

Use SunPy to read the .fits file with the data

map = sunpy.map.Map('gong.fits')
nr = 60
rss = 2.5

Extract the data, and remove the mean to enforce div(B) = 0 on the solar surface

br = map.data
br = br - np.mean(br)

Create PFSS input object

input = pfsspy.Input(br, nr, rss)

Plot input magnetic field

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

Calculate PFSS solution

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

Plot output field

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

Trace some field lines

br, btheta, bphi = output.bg

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([r * np.cos(phi),
                       r * np.sin(theta) * np.sin(phi),
                       r * np.cos(theta) * np.sin(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)

# Add inner and outer boundary circles
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])
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.
strum2cart(rho, s, phi)	Convert strumfric coordinates to cartesian coordinates.

Indices and tables

• Index
• Module Index
• Search Page