diff --git a/.gitignore b/.gitignore index b90ffee..990ef0c 100644 --- a/.gitignore +++ b/.gitignore @@ -1,5 +1,6 @@ # Sphinx gallery docs/generated +docs/sg_execution_times.rst # Coverage reports coverage.xml diff --git a/examples/loop-bundle-rtv.py b/examples/loop-bundle-rtv.py index 6127a7d..23422c4 100644 --- a/examples/loop-bundle-rtv.py +++ b/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 diff --git a/setup.cfg b/setup.cfg index 0de2932..8b86fe4 100644 --- a/setup.cfg +++ b/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 diff --git a/synthesizAR/interfaces/hydrad/hydrad.py b/synthesizAR/interfaces/hydrad/hydrad.py index e03ebe3..1171b09 100644 --- a/synthesizAR/interfaces/hydrad/hydrad.py +++ b/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): diff --git a/synthesizAR/interfaces/martens.py b/synthesizAR/interfaces/martens.py index 0e9aa5d..871d3a1 100644 --- a/synthesizAR/interfaces/martens.py +++ b/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, diff --git a/synthesizAR/models/scaling_laws.py b/synthesizAR/models/scaling_laws.py index 3948e3c..05cfd5d 100644 --- a/synthesizAR/models/scaling_laws.py +++ b/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,9 +93,10 @@ 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 @@ -101,14 +104,12 @@ def __init__(self, s: u.cm, heating_constant, alpha=0, beta=0, gamma=0.5, 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