From 4335fd442d240bcf20b8642bbba12697fd3974fe Mon Sep 17 00:00:00 2001 From: Will Barnes Date: Thu, 3 Oct 2024 12:09:58 -0400 Subject: [PATCH] more fixes --- synthesizAR/instruments/base.py | 13 +++-- synthesizAR/instruments/hinode.py | 4 ++ synthesizAR/instruments/physical.py | 78 +++++------------------------ synthesizAR/instruments/sdo.py | 4 ++ synthesizAR/models/geometry.py | 18 ++++--- 5 files changed, 38 insertions(+), 79 deletions(-) diff --git a/synthesizAR/instruments/base.py b/synthesizAR/instruments/base.py index e826f55..65d6522 100644 --- a/synthesizAR/instruments/base.py +++ b/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 diff --git a/synthesizAR/instruments/hinode.py b/synthesizAR/instruments/hinode.py index e6310e7..92f2690 100644 --- a/synthesizAR/instruments/hinode.py +++ b/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 diff --git a/synthesizAR/instruments/physical.py b/synthesizAR/instruments/physical.py index 55a92a1..daec196 100644 --- a/synthesizAR/instruments/physical.py +++ b/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 diff --git a/synthesizAR/instruments/sdo.py b/synthesizAR/instruments/sdo.py index 8f8438c..871f2bc 100644 --- a/synthesizAR/instruments/sdo.py +++ b/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}' diff --git a/synthesizAR/models/geometry.py b/synthesizAR/models/geometry.py index 2aa533c..d44a86e 100644 --- a/synthesizAR/models/geometry.py +++ b/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)]