From 0e6e054905a0c5d98790cd76fc2c77ea90939619 Mon Sep 17 00:00:00 2001 From: David Stansby Date: Tue, 28 Apr 2020 16:47:22 +0100 Subject: [PATCH 1/4] Revert to faster FieldLine specification --- doc/source/changes.rst | 6 ------ pfsspy/fieldline.py | 31 +++++++++++++++++++-------- pfsspy/tests/test_fieldline.py | 39 ++++++---------------------------- pfsspy/tracing.py | 16 ++------------ 4 files changed, 30 insertions(+), 62 deletions(-) diff --git a/doc/source/changes.rst b/doc/source/changes.rst index 57ec31e8..2bd2bbfc 100644 --- a/doc/source/changes.rst +++ b/doc/source/changes.rst @@ -70,12 +70,6 @@ pfsspy now comes with built in `sunpy` map sources for GONG and ADAPT synoptic maps, which automatically fix some non-compliant FITS header values. To use these, just import ``pfsspy`` and load the .FITS files as normal with sunpy. -Using `~pfsspy.fieldline.FieldLine` objects -~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -The input to `~pfsspy.fieldline.FieldLine` objects is now a set of -`~astropy.coordinates.SkyCoord` along with the `~pfsspy.Output` through which -the field lines were traced. - Tracing seeds ~~~~~~~~~~~~~ `pfsspy.tracing.Tracer` no longer has a ``transform_seeds`` helper method, which diff --git a/pfsspy/fieldline.py b/pfsspy/fieldline.py index 2f12af0c..ab25a025 100644 --- a/pfsspy/fieldline.py +++ b/pfsspy/fieldline.py @@ -1,6 +1,7 @@ import astropy.coordinates as coord import astropy.constants as const import sunpy.coordinates.frames as frames +import astropy.units as u import numpy as np import scipy.interpolate @@ -126,15 +127,27 @@ class FieldLine: coords : astropy.coordinates.SkyCoord Field line coordinates. """ - def __init__(self, coords, output): - self.coords = coords + def __init__(self, x, y, z, output): + self._x = np.array(x) + self._y = np.array(y) + self._z = np.array(z) + self._r = np.sqrt(self._x**2 + self._y**2 + self._z**2) self._output = output # Set _is_open - r = coords.radius atol = 0.1 - self._is_open = np.abs(r[0] - r[-1]) > const.R_sun * atol + self._is_open = np.abs(self._r[0] - self._r[-1]) > atol # Set _polarity - self._polarity = -np.sign(r[0] - r[-1]) * self._is_open + self._polarity = -np.sign(self._r[0] - self._r[-1]) * self._is_open + + @property + def coords(self): + r, lat, lon = coord.cartesian_to_spherical( + self._x, self._y, self._z) + r *= const.R_sun + lon += self._output._lon0 + 180 * u.deg + coords = coord.SkyCoord( + lon, lat, r, frame=self._output.coordinate_frame) + return coords @property def is_open(self): @@ -176,9 +189,9 @@ def solar_footpoint(self): this case. """ if self.polarity == 1 or not self.is_open: - return coord.SkyCoord(self.coords[0]) + return self.coords[0] else: - return coord.SkyCoord(self.coords[-1]) + return self.coords[-1] @property def source_surface_footpoint(self): @@ -199,9 +212,9 @@ def source_surface_footpoint(self): this case. """ if self.polarity == 1 or not self.is_open: - return coord.SkyCoord(self.coords[-1]) + return self.coords[-1] else: - return coord.SkyCoord(self.coords[0]) + return self.coords[0] @property @functools.lru_cache() diff --git a/pfsspy/tests/test_fieldline.py b/pfsspy/tests/test_fieldline.py index 769aafe9..63e2211e 100644 --- a/pfsspy/tests/test_fieldline.py +++ b/pfsspy/tests/test_fieldline.py @@ -1,25 +1,19 @@ from pfsspy.fieldline import FieldLine, FieldLines, OpenFieldLines, ClosedFieldLines import astropy.units as u -import astropy.constants as const from astropy.time import Time -from astropy.coordinates import SkyCoord from sunpy.coordinates import frames as sunframes import pytest -rsun = const.R_sun - -@pytest.mark.parametrize('r, open, pol', +@pytest.mark.parametrize('x, open, pol', [[[1, 2.5], True, 1], [[2.5, 1], True, -1], [[1, 1], False, 0], ]) -def test_open(r, open, pol): - coord = SkyCoord(0 * u.deg, 0 * u.deg, r * rsun, - frame='heliographic_carrington') - fline = FieldLine(coord, None) +def test_open(x, open, pol): + fline = FieldLine(x, [0, 0], [0, 0], None) assert (fline.is_open == open) assert (fline.polarity == pol) @@ -30,32 +24,11 @@ def test_open(r, open, pol): assert len(flines.closed_field_lines) == int(not open) -@pytest.mark.parametrize('r, cls', +@pytest.mark.parametrize('x, cls', [[[1, 2.5], ClosedFieldLines], [[1, 1], OpenFieldLines], ]) -def test_flines_errors(r, cls): - coord = SkyCoord(0 * u.deg, 0 * u.deg, r * rsun, - frame='heliographic_carrington') - fline = FieldLine(coord, None) +def test_flines_errors(x, cls): + fline = FieldLine(x, [0, 0], [0, 0], None) with pytest.raises(ValueError): cls([fline]) - - -def test_transform(): - # Check that field lines can be transformed into different coordinates - obstime = Time('1992-12-21') - stonyhurst = sunframes.HeliographicStonyhurst( - 12 * u.deg, 0 * u.deg, obstime=obstime) - coord = SkyCoord([0 * u.deg, 0 * u.deg], - [0 * u.deg, 0 * u.deg], - [1 * rsun, 2.5 * rsun], - frame='heliographic_carrington', - observer='earth', - obstime=obstime) - fline = FieldLine(coord, None) - # Check field line transform - fline.coords.transform_to(stonyhurst) - - # Check footpoint transform - fline.solar_footpoint.transform_to(stonyhurst) diff --git a/pfsspy/tracing.py b/pfsspy/tracing.py index f96679cc..9bce47dd 100644 --- a/pfsspy/tracing.py +++ b/pfsspy/tracing.py @@ -73,15 +73,6 @@ def coords_to_xyz(seeds, output): z = z.to_value(const.R_sun) return x, y, z - @staticmethod - def xyz_to_coords(x, y, z, output): - r, lat, lon = astrocoords.cartesian_to_spherical(x, y, z) - r *= const.R_sun - lon += output._lon0 + 180 * u.deg - coords = astrocoords.SkyCoord( - lon, lat, r, frame=output.coordinate_frame) - return coords - class FortranTracer(Tracer): r""" @@ -168,8 +159,7 @@ def trace(self, seeds, output): f'(currently set to {self.max_steps}) and try again.') xs = [np.stack(pfsspy.coords.strum2cart(x[:, 2], x[:, 1], x[:, 0]), axis=-1) for x in xs] - xs = [self.xyz_to_coords(x[:, 0], x[:, 1], x[:, 2], output) for x in xs] - flines = [fieldline.FieldLine(coords, output) for coords in xs] + flines = [fieldline.FieldLine(x[:, 0], x[:, 1], x[:, 2], output) for x in xs] return fieldline.FieldLines(flines) @@ -200,9 +190,7 @@ def trace(self, seeds, output): xback = output._integrate_one_way(-1, seed, self.rtol, self.atol) xback = np.flip(xback, axis=1) xout = np.row_stack((xback.T, xforw.T)) - # Hacky way to roate back by 180deg - coords = self.xyz_to_coords(xout[:, 0], xout[:, 1], xout[:, 2], output) - fline = fieldline.FieldLine(coords, output) + fline = fieldline.FieldLine(xout[:, 0], xout[:, 1], xout[:, 2], output) flines.append(fline) return fieldline.FieldLines(flines) From 782289d1da4bcaf564f898b9549fa4ce1841ed60 Mon Sep 17 00:00:00 2001 From: David Stansby Date: Tue, 28 Apr 2020 17:03:33 +0100 Subject: [PATCH 2/4] Fix examples --- examples/using_pfsspy/plot_dipole.py | 7 ++++--- examples/using_pfsspy/plot_gong.py | 9 +++++---- 2 files changed, 9 insertions(+), 7 deletions(-) diff --git a/examples/using_pfsspy/plot_dipole.py b/examples/using_pfsspy/plot_dipole.py index 63e3ab40..63179295 100644 --- a/examples/using_pfsspy/plot_dipole.py +++ b/examples/using_pfsspy/plot_dipole.py @@ -103,10 +103,11 @@ def dipole_Br(r, s): field_lines = tracer.trace(seeds, output) for field_line in field_lines: - field_line.coords.representation_type = 'cartesian' + coords = field_line.coords + coords.representation_type = 'cartesian' color = {0: 'black', -1: 'tab:blue', 1: 'tab:red'}.get(field_line.polarity) - ax.plot(field_line.coords.y / const.R_sun, - field_line.coords.z / const.R_sun, color=color) + ax.plot(coords.y / const.R_sun, + coords.z / const.R_sun, color=color) # Add inner and outer boundary circles ax.add_patch(mpatch.Circle((0, 0), 1, color='k', fill=False)) diff --git a/examples/using_pfsspy/plot_gong.py b/examples/using_pfsspy/plot_gong.py index 908fd59d..89e0dd38 100644 --- a/examples/using_pfsspy/plot_gong.py +++ b/examples/using_pfsspy/plot_gong.py @@ -108,10 +108,11 @@ def set_axes_lims(ax): for field_line in field_lines: color = {0: 'black', -1: 'tab:blue', 1: 'tab:red'}.get(field_line.polarity) - field_line.coords.representation_type = 'cartesian' - ax.plot(field_line.coords.x / const.R_sun, - field_line.coords.y / const.R_sun, - field_line.coords.z / const.R_sun, + coords = field_line.coords + coords.representation_type = 'cartesian' + ax.plot(coords.x / const.R_sun, + coords.y / const.R_sun, + coords.z / const.R_sun, color=color, linewidth=1) ax.set_title('PFSS solution') From a484e8c9fd7135f3c892b70b0622809698b30a26 Mon Sep 17 00:00:00 2001 From: David Stansby Date: Wed, 29 Apr 2020 16:36:10 +0100 Subject: [PATCH 3/4] Update docstring --- pfsspy/fieldline.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pfsspy/fieldline.py b/pfsspy/fieldline.py index ab25a025..2b6957bf 100644 --- a/pfsspy/fieldline.py +++ b/pfsspy/fieldline.py @@ -117,8 +117,8 @@ class FieldLine: Parameters ---------- - coords : astropy.coordinates.SkyCoord - Field line coordinates. + x, y, z : + Field line coordinates in cartesian coordinates. output : Output The PFSS output through which this field line was traced. From e8d2d6e3ba6aac03a1938097993176136ee81ffe Mon Sep 17 00:00:00 2001 From: David Stansby Date: Wed, 29 Apr 2020 17:01:46 +0100 Subject: [PATCH 4/4] Fix doc --- pfsspy/fieldline.py | 8 +++----- 1 file changed, 3 insertions(+), 5 deletions(-) diff --git a/pfsspy/fieldline.py b/pfsspy/fieldline.py index 2b6957bf..a047dbc7 100644 --- a/pfsspy/fieldline.py +++ b/pfsspy/fieldline.py @@ -121,11 +121,6 @@ class FieldLine: Field line coordinates in cartesian coordinates. output : Output The PFSS output through which this field line was traced. - - Attributes - ---------- - coords : astropy.coordinates.SkyCoord - Field line coordinates. """ def __init__(self, x, y, z, output): self._x = np.array(x) @@ -141,6 +136,9 @@ def __init__(self, x, y, z, output): @property def coords(self): + """ + Field line `~astropy.coordinates.SkyCoord`. + """ r, lat, lon = coord.cartesian_to_spherical( self._x, self._y, self._z) r *= const.R_sun