Merge pull request #172 from dstansby/fline-construct

Revert to faster FieldLine specification

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

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

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

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
@@ -116,25 +117,35 @@ 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.
-
-    Attributes
-    ----------
-    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):
+        """
+        Field line `~astropy.coordinates.SkyCoord`.
+        """
+        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 +187,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 +210,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()

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)

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)