Merge pull request #89 from wtbarnes/hydrad-interface-updates

Clean up HYDRAD model interface

.gitignore

@@ -1,5 +1,6 @@
 # Sphinx gallery
 docs/generated
+docs/sg_execution_times.rst
 
 # Coverage reports
 coverage.xml

examples/loop-bundle-rtv.py

@@ -11,7 +11,7 @@
 import astropy.units as u
 from astropy.visualization import quantity_support
 from astropy.coordinates import SkyCoord
-from sunpy.map import extract_along_coord
+from sunpy.map import pixelate_coord_path, sample_at_coords
 
 import synthesizAR
 from synthesizAR.instruments import InstrumentSDOAIA
@@ -59,13 +59,15 @@
 
 ###########################################################################
 # Additionally, we can look at intensity profiles  along and across the
-# loop axis using `sunpy.map.extract_along_coord`.
+# loop axis using `sunpy.map.pixelate_coord_path` and `sunpy.map.sample_at_coords`.
 coord_axis = SkyCoord(Tx=[-30, 30]*u.arcsec, Ty=0*u.arcsec,
                       frame=m_171.coordinate_frame)
+coord_axis = pixelate_coord_path(m_171, coord_axis)
+profile_axis = sample_at_coords(m_171, coord_axis)
 coord_xs = SkyCoord(Tx=0*u.arcsec, Ty=[-10, 10]*u.arcsec,
                     frame=m_171.coordinate_frame)
-profile_axis, coord_axis = extract_along_coord(m_171, coord_axis)
-profile_xs, coord_xs = extract_along_coord(m_171, coord_xs)
+coord_xs = pixelate_coord_path(m_171, coord_xs)
+profile_xs = sample_at_coords(m_171, coord_xs)
 
 ###########################################################################
 # Note that the intensity is highest at the footpoints because we are
@@ -80,9 +82,9 @@
 with quantity_support():
     plt.figure(figsize=(11, 5))
     plt.subplot(121)
-    plt.plot(coord_axis.Tx, profile_axis, color='C0')
+    plt.plot(coord_axis.separation(coord_axis[0]).to('arcsec'), profile_axis, color='C0')
     plt.subplot(122)
-    plt.plot(coord_xs.Ty, profile_xs, color='C1')
+    plt.plot(coord_xs.separation(coord_xs[0]).to('arcsec'), profile_xs, color='C1')
 
 ###########################################################################
 # Finally, we can also compute the AIA 171 intensity as viewed from the

setup.cfg

@@ -29,10 +29,10 @@ aia =
 xrt =
     xrtpy @ git+https://github.com/HinodeXRT/xrtpy.git@main
 atomic =
-    plasmapy[theory]
+    plasmapy
     fiasco @ git+https://github.com/wtbarnes/fiasco.git@main
 hydrad =
-    pydrad @ git+https://github.com/rice-solar-physics/pydrad.git@master
+    pydrad @ git+https://github.com/rice-solar-physics/pydrad.git@main
 parallel =
     dask[complete]
     distributed

synthesizAR/interfaces/hydrad/hydrad.py

@@ -2,12 +2,11 @@
 Model interface for the HYDrodynamics and RADiation (HYDRAD) code
 """
 import os
+import pathlib
 
 import numpy as np
 from scipy.interpolate import splrep, splev
 import astropy.units as u
-import astropy.constants as const
-import sunpy.sun.constants as sun_const
 
 from pydrad.configure import Configure
 from pydrad.parse import Strand
@@ -16,23 +15,28 @@
 __all__ = ['HYDRADInterface']
 
 
-class HYDRADInterface(object):
+class HYDRADInterface:
     """
     Interface to the HYDrodynamics and RADiation (HYDRAD) code
 
     Configure, interpolate, and load results for HYDRAD simulations
-    for each loop in the magnetic skeleton.
+    for each loop in the magnetic skeleton. Note that if you want to
+    just use a preexisting HYDRAD simulation, the simulations should
+    just be placed in the `output_dir` directory with the name of
+    each subdirectory corresponding to the appropriate loop. If you do
+    this, all other input parameters besides `interpolate_to_norm` can
+    be omitted.
 
     Parameters
     ----------
-    base_config: `dict`
-        Dictionary of configuration options for the HYDRAD model.
-    hydrad_dir: `str`
-        Path to a "clean copy" of HYDRAD.
-    output_dir: `str`
+    output_dir: pathlike
         Root directory to place all of the HYDRAD results in. Subdirectories
         will be named according to the name of each loop.
-    heating_model: object
+    base_config: `dict`, optional
+        Dictionary of configuration options for the HYDRAD model.
+    hydrad_dir: pathlike, optional
+        Path to a "clean copy" of HYDRAD.
+    heating_model: object, optional
         Instance of a heating model class that describes when and
         where the heating should occur along the loop
     use_gravity: `bool`, optional
@@ -42,7 +46,7 @@ class HYDRADInterface(object):
     use_initial_conditions: `bool`, optional
         If true, use only the hydrostatic initial conditions as the model
     maximum_chromosphere_ratio: `float`, optional
-        Maximum allowed ratio between the loop length and the total chromsphere depth.
+        Maximum allowed ratio between the loop length and the total chromosphere depth.
         If specified, `general.footpoint_height` will be set to this ratio times the
         loop length if ``2 * general.footpoint_height / length`` is greater than this
         ratio.
@@ -56,18 +60,18 @@ class HYDRADInterface(object):
 
     @u.quantity_input
     def __init__(self,
-                 base_config,
-                 hydrad_dir,
                  output_dir,
-                 heating_model,
+                 base_config=None,
+                 hydrad_dir=None,
+                 heating_model=None,
                  use_gravity=True,
                  use_magnetic_field=True,
                  use_initial_conditions=False,
                  maximum_chromosphere_ratio=None,
                  interpolate_to_norm=False):
+        self.output_dir = pathlib.Path(output_dir)
         self.base_config = base_config
-        self.hydrad_dir = hydrad_dir
-        self.output_dir = output_dir
+        self.hydrad_dir = hydrad_dir if hydrad_dir is None else pathlib.Path(hydrad_dir)
         self.heating_model = heating_model
         self.use_gravity = use_gravity
         self.use_magnetic_field = use_magnetic_field
@@ -77,16 +81,11 @@ def __init__(self,
 
     def configure_input(self, loop):
         config = self.base_config.copy()
-        # NOTE: Truncate precision in loop length here as passing a loop length
-        # with too many significant figures to HYDRAD can cause issues with the
-        # initical conditions calculation. Here, we truncate the loop length such
-        # that it has 1 significant figure when expressed in Mm
-        length = float(f'{loop.length.to("Mm").value:.1f}') * u.Mm
-        config['general']['loop_length'] = length
-        config['initial_conditions']['heating_location'] = length / 2.
+        config['general']['loop_length'] = loop.length
+        config['initial_conditions']['heating_location'] = loop.length / 2.
         if self.maximum_chromosphere_ratio:
-            config['general']['footpoint_height'] = min(config['general']['footpoint_height'],
-                                                        self.maximum_chromosphere_ratio * length / 2)
+            config['general']['footpoint_height'] = min(
+                config['general']['footpoint_height'], self.maximum_chromosphere_ratio * loop.length / 2)
         if self.use_gravity:
             config['general']['poly_fit_gravity'] = self.configure_gravity_fit(loop)
         if self.use_magnetic_field:
@@ -98,9 +97,10 @@ def configure_input(self, loop):
                            verbose=False)
 
     def load_results(self, loop):
+        strand = Strand(self.output_dir / loop.name)
         return self._load_results_from_strand(
             loop,
-            Strand(os.path.join(self.output_dir, loop.name)),
+            strand,
             use_initial_conditions=self.use_initial_conditions,
             interpolate_to_norm=self.interpolate_to_norm,
         )
@@ -110,38 +110,37 @@ def _load_results_from_strand(loop,
                                   strand,
                                   use_initial_conditions=False,
                                   interpolate_to_norm=False):
-        loop_coord_center = loop.field_aligned_coordinate_center.to(u.cm).value
+        loop_coord_center = loop.field_aligned_coordinate_center.to_value('cm')
         if interpolate_to_norm:
-            loop_coord_center /= loop.length.to(u.cm).value
+            loop_coord_center = loop.field_aligned_coordinate_center_norm.value
         if use_initial_conditions:
             time = strand.initial_conditions.time.reshape((1,))
             strand = [strand.initial_conditions]
         else:
             time = strand.time
         shape = time.shape + loop_coord_center.shape
-        electron_temperature = np.zeros(shape)
-        ion_temperature = np.zeros(shape)
-        density = np.zeros(shape)
-        velocity = np.zeros(shape)
+        quantities = {
+            'electron_temperature': np.zeros(shape) * u.K,
+            'ion_temperature': np.zeros(shape) * u.K,
+            'density': np.zeros(shape) * u.cm**(-3),
+            'velocity': np.zeros(shape) * u.cm/u.s,
+        }
         for i, p in enumerate(strand):
-            coord = p.coordinate.to(u.cm).value
+            coord = p.coordinate.to('cm').value
             if interpolate_to_norm:
-                coord /= strand.loop_length.to(u.cm).value
-            tsk = splrep(coord, p.electron_temperature.to(u.K).value,)
-            electron_temperature[i, :] = splev(loop_coord_center, tsk, ext=0)
-            tsk = splrep(coord, p.ion_temperature.to(u.K).value,)
-            ion_temperature[i, :] = splev(loop_coord_center, tsk, ext=0)
-            tsk = splrep(coord, p.electron_density.to(u.cm**(-3)).value)
-            density[i, :] = splev(loop_coord_center, tsk, ext=0)
-            tsk = splrep(coord, p.velocity.to(u.cm/u.s).value,)
-            velocity[i, :] = splev(loop_coord_center, tsk, ext=0)
+                coord /= strand.loop_length.to('cm').value
+            for k in quantities:
+                q = getattr(p, 'electron_density' if k == 'density' else k)
+                tsk = splrep(coord, q.to_value(quantities[k].unit))
+                q_interp = splev(loop_coord_center, tsk, ext=0)
+                quantities[k][i, :] = u.Quantity(q_interp, quantities[k].unit)
 
         return (
             time,
-            electron_temperature*u.K,
-            ion_temperature*u.K,
-            density*u.cm**(-3),
-            velocity*u.cm/u.s,
+            quantities['electron_temperature'],
+            quantities['ion_temperature'],
+            quantities['density'],
+            quantities['velocity'],
         )
 
     def configure_gravity_fit(self, loop):

synthesizAR/interfaces/martens.py

@@ -43,7 +43,7 @@ def load_results(self, loop):
         time = u.Quantity([0, ], 's')
         s_half = np.linspace(0, 1, 1000)*loop.length/2
         H = self.get_heating_constant(loop)
-        msl = MartensScalingLaws(s_half, H, **self.model_parameters)
+        msl = MartensScalingLaws(s_half, loop.length/2, H, **self.model_parameters)
         # Make sure there are no temperatures below specified cutoff
         msl_temperature = np.where(msl.temperature < self.temperature_cutoff,
                                    self.temperature_cutoff,

synthesizAR/models/scaling_laws.py

@@ -70,6 +70,8 @@ class MartensScalingLaws(object):
     ----------
     s : `~astropy.units.Quantity`
         Field-aligned loop coordinate for half of symmetric, semi-circular loop
+    loop_length : `~astropy.units.Quantity`
+        Loop half-length
     heating_constant : `astropy.units.Quantity`
         Constant of proportionality that relates the actual heating rate to the
         scaling with temperature and pressure. The actual units will depend on
@@ -91,24 +93,23 @@ class MartensScalingLaws(object):
     """
 
     @u.quantity_input
-    def __init__(self, s: u.cm, heating_constant, alpha=0, beta=0, gamma=0.5,
+    def __init__(self, s: u.cm, loop_length: u.cm, heating_constant, alpha=0, beta=0, gamma=0.5,
                  chi=10**(-18.8) * u.erg * u.cm**3 / u.s * u.K**(0.5)):
         self.s = s
+        self.loop_length = loop_length
         self.heating_constant = heating_constant
         self.alpha = alpha
         self.beta = beta
         self.gamma = gamma
         self.chi = chi
         self.chi_0 = self.chi/(4.*(const.k_B**2))
 
-    @property
-    @u.quantity_input
-    def loop_length(self,) -> u.cm:
-        return np.diff(self.s).sum()
-
     @property
     def x(self,):
-        return (self.s/self.loop_length).decompose()
+        x = (self.s/self.loop_length).decompose()
+        if (x > 1).any():
+            raise ValueError()
+        return x
 
     @property
     @u.quantity_input