New Feature: field line integrals

- Tracers objects now return with an "integral" field.
- By default the integral is field line length, but like in regular
  mapfl, additional arguments can make it do scalar field integrals:
  - integrate_along_fl (logical that turns it on, default False)
  - scalar_input_file (3D RTP H5 file path with scalar field to be integrated).
  - weight_integral_by_area (logical to weight integral by 1/b, default False)
  - max_along_fl (logical to take the max of the scalar field along the
    field line, defualt False)
- Update mapfl to reflect current GitHub version and add this capability.
- Update mapflpy fortran for compatitbility.
- Update Tracers namedtuple definition.
- Compatibility updates in other python files.
- WARNING: Currently the scalar_input_file is the only way to read a
  scalar (i.e. from a file on disk). Eventually it will be possible to
  pass this by memory.
- An example for the docs pages will be added soon (TBD).

Author: cooperdowns

environment.yml

@@ -30,6 +30,7 @@ dependencies:
   - ruff
   - coverage
   - nox
+  - pre-commit
 
   # -------- Docs --------
   - sphinx

mapflpy/globals.py

@@ -30,7 +30,7 @@
 # It contains the geometry of the traces, their start and end positions,
 # and whether they were traced to a boundary.
 #------------------------------------------------------------------------------
-Traces = namedtuple('Traces', ['geometry', 'start_pos', 'end_pos', 'traced_to_boundary'])
+Traces = namedtuple('Traces', ['geometry', 'start_pos', 'end_pos', 'traced_to_boundary', 'integral'])
 Traces.__doc__ = (
     """Named tuple for storing magnetic fieldline trace information.
     
@@ -182,6 +182,8 @@ class Polarity(IntEnum):
 DEFAULT_PARAMS = MappingProxyType({
     'integrate_along_fl_': False,
     'scalar_input_file_': '',
+    'weight_integral_by_area_': False,
+    'max_along_fl_': False,
     'verbose_': False,
     'cubic_': False,
     'trace_fwd_': False,

mapflpy/tracer.py

@@ -54,9 +54,6 @@
 _BS1 = np.zeros(3, order='F').astype(np.float64)
 """Sentinel array for :meth:`mapfl.trace` **bs1** parameter."""
 
-_S = np.zeros(1, order='F').astype(np.float64)
-"""Sentinel array for :meth:`mapfl.trace` **s** parameter."""
-
 
 def state_modifier(method):
     """
@@ -289,19 +286,20 @@ def trace(lps, buffer_size):
             traces = np.full((buffer_size, *lps.shape), np.nan, order='F').astype(np.float64)
             s1 = np.zeros(lps.shape, np.float64, order='F')
             mask = np.full((1, lps.shape[1]), False, order='F')
+            integral = np.full((1, lps.shape[1]), 0.0, order='F')
             for i in range(lps.shape[1]):
                 trace_args = dict(
                     s0=lps[:, i],
                     s1=s1[:, i],
                     bs0=_BS0,
                     bs1=_BS1,
-                    s=_S,
+                    s=integral[:,i],
                     traced_to_r_boundary=mask[:, i],
                     svec=traces[:, :, i],
                     svec_n=buffer_size
                 )
                 mapfl.trace(**trace_args)
-            return Traces(traces, lps, s1, mask[0, :])
+            return Traces(traces, lps, s1, mask[0, :], integral[0,:])
 
         while True:
             method, *args = pipe.recv()
@@ -866,27 +864,29 @@ def _trace(self,
             Structured container with traced fieldline geometry, starting points, ending points,
             and whether the fieldlines reached the radial boundary of the provided mesh.
         """
-        # `traces`, `s1`, and `mask` are initialized to Fortran-contiguous arrays
-        # and are used to store the traced fieldlines, their end points, and boundary masks.
+        # `traces`, `s1`, `mask`, and `integral` are initialized to Fortran-contiguous arrays
+        # and are used to store the traced fieldlines, their end points, boundary masks, and
+        # the scalar field integral along the field line (length by default).
         # `mapflpy_fortran` alters these arrays in-place during tracing.
         # Each array must have a dimensionality > 1 so that a numpy view is passed.
         traces = np.full((buff, *lps.shape), np.nan, order='F').astype(np.float64)
         s1 = np.zeros(lps.shape, np.float64, order='F')
         mask = np.full((1, lps.shape[1]), False, order='F')
+        integral = np.full((1, lps.shape[1]), 0.0, order='F')
 
-        # Since the parameters _bs0, _bs1, and _s are discarded after the trace,
-        # we can use the globals `_BS0`, `_BS1`, and `_S` as placeholders to avoid having
+        # Since the parameters _bs0 and _bs1 are discarded after the trace,
+        # we can use the globals `_BS0` and `_BS1` as placeholders to avoid having
         # to repeatedly create and destruct arrays.
         for i in range(lps.shape[1]):
             self._mapfl.trace(s0=lps[:, i],
                               s1=s1[:, i],
                               bs0=_BS0,
                               bs1=_BS1,
-                              s=_S,
+                              s=integral[:,i],
                               traced_to_r_boundary=mask[:, i],
                               svec=traces[:, :, i],
                               svec_n=buff)
-        return Traces(traces, lps, s1, mask[0, :])
+        return Traces(traces, lps, s1, mask[0, :], integral[0,:])
 
 
 class TracerMP(_Tracer):

mapflpy/utils.py

@@ -220,7 +220,8 @@ def combine_fwd_bwd_traces(fwd_traces: Traces,
                                              fwd_traces.geometry]),
                              bwd_traces.end_pos,
                              fwd_traces.end_pos,
-                             fwd_traces.traced_to_boundary & bwd_traces.traced_to_boundary)
+                             fwd_traces.traced_to_boundary & bwd_traces.traced_to_boundary,
+                             fwd_traces.integral + bwd_traces.integral)
     return combined_traces
 
 
@@ -423,7 +424,8 @@ def trim_fieldline_nan_buffer(traces: Traces | np.ndarray) -> List[np.ndarray]:
 
 def combine_and_pad_fieldlines(arrs: list | tuple,
                                to_trace: bool = False,
-                               traced_to_boundary: Optional[NDArray[np.bool_]] = None
+                               traced_to_boundary: Optional[NDArray[np.bool_]] = None,
+                               integral: Optional[NDArray[np.bool_]] = None
                                ) -> np.ndarray | Traces:
     """
     NaN-pad a list of variable-length 3D point arrays into a single dense array.
@@ -499,7 +501,7 @@ def combine_and_pad_fieldlines(arrs: list | tuple,
             geometry[:a.shape[1], :, i] = a.T   # (3, ni) -> (ni, 3)
             start_pos[:, i] = a[:,0]
             end_pos[:, i] = a[:, -1]
-        return Traces(geometry, start_pos, end_pos, traced_to_boundary)
+        return Traces(geometry, start_pos, end_pos, traced_to_boundary, integral)