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
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. |
for usage examples see
pfsspy examples¶
Note
Click here to download the full example code
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()
Total running time of the script: ( 0 minutes 15.111 seconds)
Note
Click here to download the full example code
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
Total running time of the script: ( 0 minutes 30.047 seconds)
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. |