more fixes

synthesizAR/instruments/base.py

@@ -105,6 +105,10 @@ def detector(self):
     def observatory(self):
         return self.name
 
+    @property
+    def _expected_unit(self):
+        return None
+
     def get_instrument_name(self, channel):
         return self.name
 
@@ -351,7 +355,7 @@ def integrate_los(self, time, channel, skeleton, coordinates_centers, bins, bin_
             coordinates_centers.Ty.value,
             bins=bins,
             range=((blc.Tx.value, trc.Tx.value), (blc.Ty.value, trc.Ty.value)),
-            weights=kernels.value * visible,
+            weights=kernels.to_value(self._expected_unit) * visible,
         )
         # For some quantities, need to average over all components along a given LOS
         if self.average_over_los:
@@ -363,18 +367,17 @@ def integrate_los(self, time, channel, skeleton, coordinates_centers, bins, bin_
                 weights=visible,
             )
             hist /= np.where(_hist == 0, 1, _hist)
-        new_header = copy.deepcopy(header)
-        new_header['bunit'] = kernels.unit.to_string('fits')
         # NOTE: Purposefully using a nonstandard key to record this time as we do not
         # want this to have the implicit consequence of changing the coordinate frame
         # by changing a more standard time key. However, still want to record this
         # information somewhere in the header.
         # FIXME: Figure out a better way to deal with this.
+        new_header = copy.deepcopy(header)
         new_header['date_sim'] = (self.observer.obstime + time).isot
 
         return Map(hist.T, new_header)
 
-    def get_header(self, channel, coordinates, unit=None):
+    def get_header(self, channel, coordinates):
         """
         Create the FITS header for a given channel and set of loop coordinates
         that define the needed FOV.
@@ -395,7 +398,7 @@ def get_header(self, channel, coordinates, unit=None):
             telescope=self.telescope,
             detector=self.detector,
             wavelength=channel.channel,
-            unit=unit,
+            unit=self._expected_unit,
         )
         return header
 

synthesizAR/instruments/hinode.py

@@ -80,6 +80,10 @@ def detector(self):
     def telescope(self):
         return 'Hinode'
 
+    @property
+    def _expected_unit(self):
+        return u.DN / (u.pix * u.s)
+
     def get_instrument_name(self, channel):
         return self.detector
 

synthesizAR/instruments/physical.py

@@ -10,7 +10,6 @@
 import numpy as np
 from scipy.interpolate import interp1d
 import ndcube
-from sunpy.map import GenericMap
 
 from synthesizAR.instruments import ChannelBase, InstrumentBase
 from synthesizAR.util import los_velocity
@@ -53,6 +52,10 @@ def get_instrument_name(self, channel):
         # This ensures that the temperature bin labels are in the header
         return f'{self.name}_{channel.name}'
 
+    @property
+    def _expected_unit(self):
+        return u.cm**(-5)
+
     @staticmethod
     @return_quantity_as_tuple
     def calculate_intensity_kernel(loop, channel, **kwargs):
@@ -121,71 +124,6 @@ def calculate_intensity(dem, spectra, header, meta=None):
         wcs = astropy.wcs.WCS(header=wave_header)
         return ndcube.NDCube(intensity, wcs=wcs, meta=meta)
 
-    def make_slope_map(self, dem, time_index, temperature_bounds=None, em_threshold=None,
-                       rsquared_tolerance=0.5, full=False):
-        """
-        Calculate emission measure slope :math:`a` in each pixel
-
-        Create map of emission measure slopes by fitting :math:`\mathrm{EM}\sim T^a` for a
-        given temperature range. A slope is masked if a value between the `temperature_bounds`
-        is less than :math:`\mathrm{EM}`. Additionally, the "goodness-of-fit" is evaluated using
-        the correlation coefficient, :math:`r^2=1 - R_1/R_0`, where :math:`R_1` and :math:`R_0`
-        are the residuals from the first and zeroth order polynomial fits, respectively. We mask
-        the slope if :math:`r^2` is less than `rsquared_tolerance`.
-
-        Parameters
-        ----------
-        temperature_bounds : `~astropy.units.Quantity`, optional
-        em_threshold : `~astropy.units.Quantity`, optional
-            Mask slope if any emission measure in the fit interval is below this value
-        rsquared_tolerance : `float`
-            Mask any slopes with a correlation coefficient, :math:`r^2`, below this value
-        full : `bool`
-            If True, return maps of the intercept and :math:`r^2` values as well
-        """
-        if temperature_bounds is None:
-            temperature_bounds = u.Quantity((1e6, 4e6), u.K)
-        if em_threshold is None:
-            em_threshold = u.Quantity(1e25, u.cm**(-5))
-        # Get temperature fit array
-        index_temperature_bounds = np.where(np.logical_and(
-            self.temperature_bin_centers >= temperature_bounds[0],
-            self.temperature_bin_centers <= temperature_bounds[1]))
-        temperature_fit = np.log10(
-            self.temperature_bin_centers[index_temperature_bounds].to(u.K).value)
-        if temperature_fit.size < 3:
-            warnings.warn(f'Fitting to fewer than 3 points in temperature space: {temperature_fit}')
-        # Cut on temperature
-        data = u.Quantity([dem[c.name][time_index].quantity for c in self.channels]).T
-        # Get EM fit array
-        em_fit = np.log10(data.value.reshape((np.prod(data.shape[:2]),) + data.shape[2:]).T)
-        em_fit[np.logical_or(np.isinf(em_fit), np.isnan(em_fit))] = 0.0  # Filter before fitting
-        # Fit to first-order polynomial
-        coefficients, rss, _, _, _ = np.polyfit(temperature_fit, em_fit, 1, full=True,)
-        # Create mask from correlation coefficient EM threshold
-        _, rss_flat, _, _, _ = np.polyfit(temperature_fit, em_fit, 0, full=True,)
-        rss = 0.*rss_flat if rss.size == 0 else rss  # returns empty residual when fit is exact
-        rsquared = 1. - rss/rss_flat
-        rsquared_mask = rsquared.reshape(data.shape[:2]) < rsquared_tolerance
-        em_mask = np.any(data < em_threshold, axis=2)
-        # Rebuild into a map
-        base_meta = self[0].meta.copy()
-        base_meta['temp_a'] = 10.**temperature_fit[0]
-        base_meta['temp_b'] = 10.**temperature_fit[-1]
-        base_meta['bunit'] = ''
-        base_meta['detector'] = 'EM slope'
-        base_meta['comment'] = 'Linear fit to log-transformed LOS EM'
-        plot_settings = self[0].plot_settings.copy()
-        plot_settings['norm'] = None
-        m = GenericMap(
-            coefficients[0, :].reshape(data.shape[:2]),
-            base_meta,
-            mask=np.stack((em_mask, rsquared_mask), axis=2).any(axis=2),
-            plot_settings=plot_settings
-        )
-
-        return (m, coefficients, rsquared) if full else m
-
 
 class InstrumentQuantityBase(InstrumentBase):
 
@@ -206,6 +144,10 @@ def calculate_intensity_kernel(loop, *args, **kwargs):
             raise ValueError('Must pass in observer to compute LOS velocity.')
         return los_velocity(loop.velocity_xyz, observer)
 
+    @property
+    def _expected_unit(self):
+        return u.km / u.s
+
 
 class InstrumentTemperature(InstrumentQuantityBase):
     name = 'temperature'
@@ -214,3 +156,7 @@ class InstrumentTemperature(InstrumentQuantityBase):
     @return_quantity_as_tuple
     def calculate_intensity_kernel(loop, *args, **kwargs):
         return loop.electron_temperature
+
+    @property
+    def _expected_unit(self):
+        return u.K

synthesizAR/instruments/sdo.py

@@ -75,6 +75,10 @@ def detector(self):
     def telescope(self):
         return 'SDO/AIA'
 
+    @property
+    def _expected_unit(self):
+        return u.DN / (u.pix * u.s)
+
     def get_instrument_name(self, channel):
         return f'{self.detector}_{channel.telescope_number}'
 

synthesizAR/models/geometry.py

@@ -16,7 +16,7 @@ def semi_circular_loop(length: u.cm=None,
                        observer=None,
                        obstime=None,
                        n_points=1000,
-                       offset: u.cm = 0*u.cm,
+                       offset: u.cm = None,
                        gamma: u.deg = 0*u.deg,
                        inclination: u.deg = 0*u.deg,
                        ellipticity=0):
@@ -38,8 +38,8 @@ def semi_circular_loop(length: u.cm=None,
     n_points : `int`, optional
         Number of points in the coordinate. Only used if `s` is not specified.
     offset : `~astropy.units.Quantity`
-        Offset in the direction perpendicular to loop, convenient for simulating
-        arcade geometries.
+        Offset in the direction perpendicular and parallel to loop, convenient for simulating
+        arcade geometries. This should always be of length 2.
     gamma : `~astropy.units.Quantity`
         Angle between the loop axis and the HCC x-axis. This defines the orientation
         of the loop in the HCC x-y plane. `gamma=0` corresponds to a loop who's
@@ -82,10 +82,12 @@ def semi_circular_loop(length: u.cm=None,
         observer=observer,
         obstime=observer.obstime if obstime is None else obstime,
     )
-    return SkyCoord(x=-(offset + z * np.sin(inclination)) * np.sin(gamma) + x * np.cos(gamma),
-                    y=(offset + z * np.sin(inclination)) * np.cos(gamma) + x * np.sin(gamma),
-                    z=z * np.cos(inclination) + const.R_sun,
-                    frame=hcc_frame)
+    if offset is None:
+        offset = [0, 0] * u.cm
+    x_hcc = -(offset[1] + z * np.sin(inclination)) * np.sin(gamma) + (x + offset[0]) * np.cos(gamma)
+    y_hcc = (offset[1] + z * np.sin(inclination)) * np.cos(gamma) + (x + offset[0]) * np.sin(gamma)
+    z_hcc = z * np.cos(inclination) + const.R_sun
+    return SkyCoord(x=x_hcc, y=y_hcc, z=z_hcc, frame=hcc_frame)
 
 
 @u.quantity_input
@@ -114,7 +116,7 @@ def semi_circular_bundle(length: u.cm, radius: u.cm, num_strands, **kwargs):
     # and nominal length
     lengths = length + np.pi * r * np.cos(theta)
     # The resulting Y Cartesian coordinate is the offset from the HCC origin
-    offset = r * np.sin(theta)
+    offset = u.Quantity([np.zeros(r.shape)*r.unit, r*np.sin(theta)]).T
     return [semi_circular_loop(length=l, offset=o, **kwargs) for l, o in zip(lengths, offset)]