Repository: wtbarnes/synthesizAR

Issue:
Add a way to compute a velocity DEM (VDEM)
Add an instrument class to compute a VDEM or a velocity differential emission measure. This can be a generalization of the DEM calculation, but binning in 4 dimensions instead of 3.

Discussion:

Patch:
diff --git a/.rtd-environment.yml b/.rtd-environment.yml
index e58b607..da83b8a 100644
--- a/.rtd-environment.yml
+++ b/.rtd-environment.yml
@@ -5,3 +5,4 @@ dependencies:
   - python=3.12
   - pip
   - graphviz!=2.42.*,!=2.43.*
+  - sqlite<3.49.0
diff --git a/docs/references.bib b/docs/references.bib
index d560c59..c530b1f 100644
--- a/docs/references.bib
+++ b/docs/references.bib
@@ -49,3 +49,19 @@
   doi = {10.1088/0004-637X/764/2/193},
   urldate = {2014-07-23}
 }
+
+@article{mikic_importance_2013,
+  title = {The {{Importance}} of {{Geometric Effects}} in {{Coronal Loop Models}}},
+  author = {Miki{\'c}, Zoran and Lionello, Roberto and Mok, Yung and Linker, Jon A. and Winebarger, Amy R.},
+  year = {2013},
+  month = aug,
+  journal = {The Astrophysical Journal},
+  volume = {773},
+  number = {2},
+  pages = {94},
+  issn = {0004-637X},
+  doi = {10.1088/0004-637X/773/2/94},
+  url = {http://iopscience.iop.org.ezproxy.rice.edu/0004-637X/773/2/94},
+  urldate = {2013-09-16},
+  langid = {english}
+}
diff --git a/examples/ebtel-nanoflares.py b/examples/ebtel-nanoflares.py
index 7443291..2b039d4 100644
--- a/examples/ebtel-nanoflares.py
+++ b/examples/ebtel-nanoflares.py
@@ -94,7 +94,7 @@ maps = aia.observe(arcade, channels=aia.channels[3:4])
 ###########################################################################
 # We can easily visualize this time-dependent emission using a
 # `~sunpy.map.MapSequence`.
-mseq = sunpy.map.MapSequence(maps['193'], sequence=True)
+mseq = sunpy.map.Map(maps['193'], sequence=True)
 fig = plt.figure()
 ax = fig.add_subplot(projection=mseq[0])
 ani = mseq.plot(axes=ax, norm=ImageNormalize(vmin=0, vmax=5, stretch=AsinhStretch()))
diff --git a/examples/multi-instrument.py b/examples/multi-instrument.py
index 6a275cd..b6044dc 100644
--- a/examples/multi-instrument.py
+++ b/examples/multi-instrument.py
@@ -87,7 +87,7 @@ with quantity_support():
 # We'll select a field of view by specifying the center of the field of view
 # as well as the width and height.
 center = SkyCoord(Tx=0*u.arcsec, Ty=-550*u.arcsec, frame=Helioprojective(observer=sdo, obstime=sdo.obstime))
-aia = InstrumentSDOAIA([0, 1]*u.s, sdo, fov_center=center, fov_width=(250, 250)*u.arcsec)
+aia = InstrumentSDOAIA([0]*u.s, sdo, fov_center=center, fov_width=(250, 250)*u.arcsec)
 aia_images = aia.observe(skeleton)
 for k in aia_images:
     aia_images[k][0].peek()
@@ -95,7 +95,7 @@ for k in aia_images:
 ###############################################################################
 # We can carry out this same procedure for *Hinode* XRT for the same field of view.
 # We'll look just at the Be-thin and Al-poly channels.
-xrt = InstrumentHinodeXRT([0, 1]*u.s, hinode, ['Be-thin', 'Al-poly'],
+xrt = InstrumentHinodeXRT([0]*u.s, hinode, ['Be-thin', 'Al-poly'],
                           fov_center=center, fov_width=(250, 250)*u.arcsec)
 xrt_images = xrt.observe(skeleton)
 for k in xrt_images:
@@ -111,13 +111,13 @@ for k in xrt_images:
 class InstrumentSTEREOEUVI(InstrumentSDOAIA):
     name = 'STEREO_EUVI'
 
-    @property
-    def resolution(self) -> u.Unit('arcsec / pixel'):
-        return u.Quantity([1.58777404, 1.58777404], 'arcsec / pixel')
-
-    @property
-    def cadence(self) -> u.s:
-        return 1 * u.h
+    def __init__(self, *args, **kwargs):
+        super().__init__(
+            *args,
+            resolution=u.Quantity([1.58777404, 1.58777404], 'arcsec / pixel'),
+            cadence=1*u.h,
+            **kwargs,
+        )
 
     @property
     def observatory(self):
@@ -138,6 +138,6 @@ class InstrumentSTEREOEUVI(InstrumentSDOAIA):
 # We can then use our custom instrument class in the exact same way as our
 # predefined classes to model the emission from EUVI. Note that we'll only do
 # this for the 171 Å channel.
-euvi = InstrumentSTEREOEUVI([0, 1]*u.s, stereo_a, fov_center=center, fov_width=(250, 250)*u.arcsec)
+euvi = InstrumentSTEREOEUVI([0]*u.s, stereo_a, fov_center=center, fov_width=(250, 250)*u.arcsec)
 euvi_images = euvi.observe(skeleton, channels=euvi.channels[2:3])
 euvi_images['171'][0].peek()
diff --git a/synthesizAR/instruments/base.py b/synthesizAR/instruments/base.py
index fa209de..a71af8d 100644
--- a/synthesizAR/instruments/base.py
+++ b/synthesizAR/instruments/base.py
@@ -37,11 +37,15 @@ class InstrumentBase:
     Parameters
     ----------
     observing_time : `~astropy.units.Quantity`
-        Tuple of start and end observing times
-    observer_coordinate : `~astropy.coordinates.SkyCoord`
+        If cadence is also provided and this has a length of 2, this is interpreted as
+        the start and end time of the observation period and an observing time is
+        constructed based on this interval and the cadence. Otherwise, this is interpreted
+        as the times at which the observations should be forward modeled.
+    observer : `~astropy.coordinates.SkyCoord`
         Coordinate of the observing instrument
-    cadence : `~astropy.units.Quantity`
     resolution : `~astropy.units.Quantity`
+    cadence : `~astropy.units.Quantity`, optional
+        If specified, this is used to construct the observing time
     pad_fov : `~astropy.units.Quantity`, optional
         Two-dimensional array specifying the padding to apply to the field of view of the synthetic
         image in both directions. If None, no padding is applied and the field of view is defined
@@ -55,13 +59,17 @@ class InstrumentBase:
     def __init__(self,
                  observing_time: u.s,
                  observer,
+                 resolution: u.Unit('arcsec/pix'),
+                 cadence: u.s = None,
                  pad_fov: u.arcsec = None,
                  fov_center=None,
                  fov_width: u.arcsec = None,
                  average_over_los=False):
         self.observer = observer
+        self.cadence = cadence
         self.observing_time = observing_time
-        self.pad_fov = (0, 0) * u.arcsec if pad_fov is None else pad_fov
+        self.resolution = resolution
+        self.pad_fov = pad_fov
         self.fov_center = fov_center
         self.fov_width = fov_width
         self.average_over_los = average_over_los
@@ -72,19 +80,27 @@ class InstrumentBase:
 
     @observing_time.setter
     def observing_time(self, value):
-        if self.cadence is None or len(value) > 2:
-            self._observing_time = value
-        else:
+        if self.cadence is not None and len(value) == 2:
             self._observing_time = np.arange(*value.to_value('s'),
                                              self.cadence.to_value('s')) * u.s
+        else:
+            self._observing_time = value
 
     @property
     def cadence(self):
-        return None
+        return self._cadence
+
+    @cadence.setter
+    def cadence(self, value):
+        self._cadence = value
 
     @property
     def resolution(self) -> u.arcsec/u.pix:
-        return (1, 1) * u.arcsec / u.pix
+        return self._resolution
+
+    @resolution.setter
+    def resolution(self, value):
+        self._resolution = value
 
     @property
     def observer(self):
@@ -94,6 +110,16 @@ class InstrumentBase:
     def observer(self, value):
         self._observer = value
 
+    @property
+    def pad_fov(self) -> u.arcsec:
+        return self._pad_fov
+
+    @pad_fov.setter
+    def pad_fov(self, value):
+        if value is None:
+            value = [0, 0] * u.arcsec
+        self._pad_fov = value
+
     @property
     def telescope(self):
         return self.name
@@ -108,7 +134,11 @@ class InstrumentBase:
 
     @property
     def _expected_unit(self):
-        return None
+        raise NotImplementedError
+
+    @property
+    def channels(self):
+        raise NotImplementedError
 
     def get_instrument_name(self, channel):
         return self.name
diff --git a/synthesizAR/instruments/hinode.py b/synthesizAR/instruments/hinode.py
index 6014f3b..3d5331d 100644
--- a/synthesizAR/instruments/hinode.py
+++ b/synthesizAR/instruments/hinode.py
@@ -2,21 +2,16 @@
 Class for Hinode/EIS instrument. Holds information about spectral, temporal, and spatial resolution
 and other instrument-specific information.
 """
-import astropy.constants as const
 import astropy.units as u
-import json
 import numpy as np
-import pkg_resources
 import xrtpy
 
 from dataclasses import dataclass
-from scipy.ndimage import gaussian_filter
 
 from synthesizAR.instruments import ChannelBase, InstrumentBase
-from synthesizAR.util import SpatialPair
 from synthesizAR.util.decorators import return_quantity_as_tuple
 
-__all__ = ['InstrumentHinodeEIS', 'InstrumentHinodeXRT']
+__all__ = ['InstrumentHinodeXRT']
 
 
 @dataclass
@@ -35,12 +30,21 @@ class InstrumentHinodeXRT(InstrumentBase):
     name = 'Hinode_XRT'
 
     def __init__(self, observing_time, observer, filters, **kwargs):
-        super().__init__(observing_time, observer, **kwargs)
-        self.channels = self._setup_channels(filters)
+        resolution = kwargs.pop('resolution', [2.05719995499, 2.05719995499] * u.arcsec/u.pixel)
+        cadence = kwargs.pop('cadence', 20 * u.s)
+        self._filters = filters
+        super().__init__(
+            observing_time=observing_time,
+            observer=observer,
+            resolution=resolution,
+            cadence=cadence,
+            **kwargs
+        )
 
-    def _setup_channels(self, filters):
+    @property
+    def channels(self):
         channels = []
-        for f in filters:
+        for f in self._filters:
             trf = xrtpy.response.TemperatureResponseFundamental(f, self.observer.obstime)
             # Assign filter wheel
             if f in xrtpy.response.effective_area.index_mapping_to_fw1_name:
@@ -60,14 +64,6 @@ class InstrumentHinodeXRT(InstrumentBase):
             channels.append(c)
         return channels
 
-    @property
-    def cadence(self) -> u.s:
-        return 20 * u.s
-
-    @property
-    def resolution(self) -> u.arcsec / u.pix:
-        return [2.05719995499, 2.05719995499] * u.arcsec/u.pixel
-
     @property
     def observatory(self):
         return 'Hinode'
@@ -100,150 +96,3 @@ class InstrumentHinodeXRT(InstrumentBase):
                         channel.temperature,
                         channel.response)
         return K_T * loop.density**2
-
-
-class InstrumentHinodeEIS(InstrumentBase):
-    """
-    Class for Extreme-ultraviolet Imaging Spectrometer (EIS) instrument on the Hinode spacecraft.
-    Converts emissivity calculations for each loop into detector units based on the spectral,
-    spatial, and temporal resolution along with the instrument response functions.
-    """
-
-    def __init__(self, observing_time, observer_coordinate, window=None, apply_psf=True):
-        self.name = 'Hinode_EIS'
-        self.cadence = 10.0*u.s
-        self.resolution = SpatialPair(x=1.0*u.arcsec/u.pixel, y=2.0*u.arcsec/u.pixel, z=None)
-        self.fits_template['telescop'] = 'Hinode'
-        self.fits_template['instrume'] = 'EIS'
-        self.fits_template['detector'] = 'EIS'
-        self.fits_template['waveunit'] = 'angstrom'
-        self.apply_psf = apply_psf
-        self.window = 0.5*u.angstrom if window is None else window
-        super().__init__(observing_time, observer_coordinate=observer_coordinate)
-        self._setup_channels()
-
-    def _setup_channels(self):
-        """
-        Read instrument properties from files. This is a temporary solution and requires that the
-        detector files all be collected into the same directory and be formatted in a specific way.
-
-        .. warning:: This method will be modified once EIS response functions become
-                    available in a different format.
-        """
-        hinode_fn = pkg_resources.resource_filename('synthesizAR',
-                                                    'instruments/data/hinode_eis.json')
-        with open(hinode_fn) as f:
-            eis_info = json.load(f)
-
-        self.channels = []
-        for key in eis_info:
-            if key != 'name' and key != 'description':
-                self.channels.append({
-                    'wavelength': eis_info[key]['wavelength']*u.Unit(eis_info[key]['wavelength_units']),
-                    'name': key,
-                    'response': {
-                        'x': eis_info[key]['response_x']*u.Unit(eis_info[key]['response_x_units']),
-                        'y': eis_info[key]['response_y']*u.Unit(eis_info[key]['response_y_units'])},
-                    'spectral_resolution': eis_info[key]['spectral_resolution']*u.Unit(eis_info[key]['spectral_resolution_units']),
-                    'gaussian_width': {'x': (3.*u.arcsec)/self.resolution.x,
-                                       'y': (3.*u.arcsec)/self.resolution.y},
-                    'instrument_width': eis_info[key]['instrument_width']*u.Unit(eis_info[key]['instrument_width_units']),
-                    'wavelength_range': [eis_info[key]['response_x'][0],
-                                         eis_info[key]['response_x'][-1]]*u.Unit(eis_info[key]['response_x_units'])
-                })
-
-        self.channels = sorted(self.channels, key=lambda x: x['wavelength'])
-
-    def make_fits_header(self, field, channel):
-        """
-        Extend base method to include extra wavelength dimension.
-        """
-        header = super().make_fits_header(field, channel)
-        header['wavelnth'] = channel['wavelength'].value
-        header['naxis3'] = len(channel['response']['x'])
-        header['ctype3'] = 'wavelength'
-        header['cunit3'] = 'angstrom'
-        header['cdelt3'] = np.fabs(np.diff(channel['response']['x']).value[0])
-        return header
-
-    def build_detector_file(self, file_template, chunks, field):
-        """
-        Build HDF5 files to store detector counts
-        """
-        additional_fields = [f'{line.value}' for line in field.strands[0].resolved_wavelengths]
-        super().build_detector_file(file_template, chunks, additional_fields=additional_fields)
-
-    def flatten(self, loop, interp_s, hf, start_index):
-        """
-        Flatten loop emission to HDF5 file for given number of wavelengths
-        """
-        for wavelength in loop.resolved_wavelengths:
-            emiss, ion_name = loop.get_emission(wavelength, return_ion_name=True)
-            dset = hf[f'{str(wavelength.value)}']
-            dset.attrs['ion_name'] = ion_name
-            self.interpolate_and_store(emiss, loop, interp_s, dset, start_index)
-
-    @staticmethod
-    def compute_emission():
-        """
-        Compute intensity for transitions observed by EIS
-        """
-        pass
-
-    def flatten_parallel(self, loops, interpolated_loop_coordinates, save_path, emission_model):
-        """
-        Build task graph for computing EIS spectra
-        """
-        pass
-
-    def detect(self, hf, channel, i_time, header, temperature, los_velocity):
-        """
-        Calculate response of Hinode/EIS detector for given loop object.
-        """
-        import plasmapy
-
-        # trim the instrument response to the appropriate wavelengths
-        trimmed_indices = []
-        for w in channel['model_wavelengths']:
-            indices = np.where(np.logical_and(channel['response']['x'] >= w-self.window,
-                                              channel['response']['x'] <= w+self.window))
-            trimmed_indices += indices[0].tolist()
-        trimmed_indices = list(sorted(set(trimmed_indices+[0, len(channel['response']['x'])-1])))
-        response_x = channel['response']['x'][trimmed_indices]
-        response_y = channel['response']['y'][trimmed_indices]
-
-        # compute the response
-        counts = np.zeros(temperature.shape+response_x.shape)
-        for wavelength in channel['model_wavelengths']:
-            # thermal width + instrument width
-            ion_name = hf[f'{str(wavelength.value)}'].attrs['ion_name']
-            ion_mass = plasmapy.atomic.ion_mass(ion_name.split(' ')[0].capitalize()).cgs
-            thermal_velocity = 2.*const.k_B.cgs*temperature/ion_mass
-            thermal_velocity = np.expand_dims(thermal_velocity, axis=2)*thermal_velocity.unit
-            line_width = ((wavelength**2)/(2.*const.c.cgs**2)*thermal_velocity
-                          + (channel['instrument_width']/(2.*np.sqrt(2.*np.log(2.))))**2)
-            # doppler shift due to LOS velocity
-            doppler_shift = wavelength*los_velocity/const.c.cgs
-            doppler_shift = np.expand_dims(doppler_shift, axis=2)*doppler_shift.unit
-            # combine emissivity with instrument response function
-            dset = hf[f'{str(wavelength.value)}']
-            hist, edges = np.histogramdd(self.total_coordinates.value,
-                                         bins=[self.bins.x, self.bins.y, self.bins.z],
-                                         range=[self.bin_range.x, self.bin_range.y, self.bin_range.z],
-                                         weights=np.array(dset[i_time, :]))
-            emiss = np.dot(hist, np.diff(edges[2])).T
-            emiss = np.expand_dims(emiss, axis=2) * u.Unit(dset.attrs.get('unit')) *self.total_coordinates.unit
-            intensity = emiss*response_y/np.sqrt(2.*np.pi*line_width)
-            intensity *= np.exp(-((response_x - wavelength - doppler_shift)**2)/(2.*line_width))
-            if not hasattr(counts, 'unit'):
-                counts = counts*intensity.unit
-            counts += intensity
-
-        header['bunit'] = counts.unit.to_string()
-        if self.apply_psf:
-            counts = (gaussian_filter(counts.value, (channel['gaussian_width']['y'].value,
-                                                     channel['gaussian_width']['x'].value, 0))
-                      * counts.unit)
-
-        from synthesizAR.analysis import EISCube
-        return EISCube(data=counts, header=header, wavelength=response_x)
diff --git a/synthesizAR/instruments/physical.py b/synthesizAR/instruments/physical.py
index daec196..be4cd99 100644
--- a/synthesizAR/instruments/physical.py
+++ b/synthesizAR/instruments/physical.py
@@ -2,22 +2,22 @@
 Instruments for calculating physical quantities, rather than
 observed counts, projected along a LOS
 """
-from dataclasses import dataclass
-import warnings
-
 import astropy.units as u
 import astropy.wcs
+import ndcube
 import numpy as np
+
+from dataclasses import dataclass
 from scipy.interpolate import interp1d
-import ndcube
 
 from synthesizAR.instruments import ChannelBase, InstrumentBase
+from synthesizAR.instruments.util import add_wave_keys_to_header, extend_celestial_wcs
 from synthesizAR.util import los_velocity
 from synthesizAR.util.decorators import return_quantity_as_tuple
-from synthesizAR.instruments.util import add_wave_keys_to_header, extend_celestial_wcs
 
 __all__ = [
     'InstrumentDEM',
+    'InstrumentVDEM',
     'InstrumentLOSVelocity',
     'InstrumentTemperature'
 ]
@@ -25,12 +25,19 @@ __all__ = [
 
 @dataclass
 class ChannelDEM(ChannelBase):
-    bin_edges: u.Quantity = None
+    # NOTE: These have to have a default value because of the
+    # keyword arguments on the base class.
+    bin_edges: u.Quantity[u.K] = None
 
     def __post_init__(self):
-        self.log_bin_edges = np.log10(self.bin_edges.to('K').value)
+        self.log_bin_edges = np.log10(self.bin_edges.to_value('K'))
         self.name = f'{self.log_bin_edges[0]:.2f}-{self.log_bin_edges[1]:.2f}'
 
+    @property
+    @u.quantity_input
+    def bin_center(self) -> u.K:
+        return self.bin_edges.sum()/2
+
 
 class InstrumentDEM(InstrumentBase):
     name = 'DEM'
@@ -38,11 +45,14 @@ class InstrumentDEM(InstrumentBase):
     @u.quantity_input
     def __init__(self, *args, temperature_bin_edges: u.K, **kwargs):
         self.temperature_bin_edges = temperature_bin_edges
-        n_bins = temperature_bin_edges.shape[0]-1
-        bin_edges = [temperature_bin_edges[[i, i+1]] for i in range(n_bins)]
-        self.channels = [ChannelDEM(bin_edges=be) for be in bin_edges]
         super().__init__(*args, **kwargs)
 
+    @property
+    def channels(self):
+        n_bins = self.temperature_bin_edges.shape[0]-1
+        bin_edges = [self.temperature_bin_edges[[i, i+1]] for i in range(n_bins)]
+        return [ChannelDEM(bin_edges=be) for be in bin_edges]
+
     @property
     @u.quantity_input
     def temperature_bin_centers(self) -> u.K:
@@ -66,9 +76,9 @@ class InstrumentDEM(InstrumentBase):
         kernel = n**2 * bin_mask
         return kernel
 
-    def dem_maps_to_cube(self, dem, time_index):
+    def maps_to_cube(self, dem, time_index):
         """
-        Convert a list of DEM maps to a DEM NDCube
+        Transform a set of DEM maps at a single time step to a `~ndcube.NDCube`
         """
         # NOTE: this is the format that .observe returns
         return type(self).dem_maps_list_to_cube(
@@ -78,14 +88,11 @@ class InstrumentDEM(InstrumentBase):
 
     @staticmethod
     def dem_maps_list_to_cube(dem_maps, temperature_bin_centers):
-        # Construct WCS
-        compound_wcs = extend_celestial_wcs(dem_maps[0].wcs,
-                                            temperature_bin_centers,
-                                            'temperature',
-                                            'phys.temperature')
-        # Stack arrays
+        compound_wcs = extend_celestial_wcs(
+            dem_maps[0].wcs,
+            (temperature_bin_centers, 'temperature', 'phys.temperature'),
+        )
         dem_array = u.Quantity([d.quantity for d in dem_maps])
-
         return ndcube.NDCube(dem_array, wcs=compound_wcs, meta=dem_maps[0].meta)
 
     @staticmethod
@@ -107,14 +114,14 @@ class InstrumentDEM(InstrumentBase):
         temperature_bin_centers = dem.axis_world_coords(0)[0]
         wavelength_spectra = spectra.axis_world_coords(1)[0]
         temperature_spectra = spectra.axis_world_coords(0)[0].to(temperature_bin_centers.unit)
-        # Interpolate spectral cube to DEM temperatures
-        f_interp = interp1d(temperature_spectra.value, spectra.data,
-                            axis=0, bounds_error=False, fill_value=0.0)
+        f_interp = interp1d(temperature_spectra.value,
+                            spectra.data,
+                            axis=0,
+                            bounds_error=False,
+                            fill_value=0.0)
         spectra_interp = f_interp(temperature_bin_centers.value)
         spectra_interp = spectra_interp * spectra.unit
-        # Take dot product between DEM and spectra
         intensity = np.tensordot(spectra_interp, u.Quantity(dem.data, dem.unit), axes=([0], [0]))
-        # Construct cube
         wave_header = add_wave_keys_to_header(wavelength_spectra, header)
         wave_header['BUNIT'] = intensity.unit.to_string()
         wave_header['NAXIS'] = len(intensity.shape)
@@ -125,13 +132,132 @@ class InstrumentDEM(InstrumentBase):
         return ndcube.NDCube(intensity, wcs=wcs, meta=meta)
 
 
+@dataclass
+class ChannelVDEM(ChannelBase):
+    # NOTE: These have to have a default value because of the
+    # keyword arguments on the base class.
+    temperature_bin_edges: u.Quantity[u.K] = None
+    velocity_bin_edges: u.Quantity[u.km/u.s] = None
+
+    def __post_init__(self):
+        self.name = self._make_channel_label(self.temperature_bin_edges,
+                                             self.velocity_bin_edges)
+
+    @staticmethod
+    def _make_channel_label(temperature_bin_edges, velocity_bin_edges):
+        # This is a static method so it can be used externally
+        log_tbin_edges = np.log10(temperature_bin_edges.to_value('K'))
+        vbin_edges = velocity_bin_edges.to_value('km/s')
+        tname = f'{log_tbin_edges[0]:.2f}-{log_tbin_edges[1]:.2f}'
+        vname = f'{vbin_edges[0]:.2f}-{vbin_edges[1]:.2f}'
+        return f'logT:{tname}_v:{vname}'
+
+    @property
+    @u.quantity_input
+    def temperature_bin_center(self) -> u.K:
+        return self.temperature_bin_edges.sum()/2
+
+    @property
+    @u.quantity_input
+    def velocity_bin_center(self) -> u.Unit('km/s'):
+        return self.velocity_bin_edges.sum()/2
+
+
+class InstrumentVDEM(InstrumentBase):
+    name = 'VDEM'
+
+    @u.quantity_input
+    def __init__(self,
+                 *args,
+                 temperature_bin_edges: u.K,
+                 velocity_bin_edges: u.Unit('km/s'),
+                 **kwargs):
+        self.temperature_bin_edges = temperature_bin_edges
+        self.velocity_bin_edges = velocity_bin_edges
+        super().__init__(*args, **kwargs)
+
+    @property
+    def channels(self):
+        channels = []
+        for i in range(self.temperature_bin_centers.shape[0]):
+            for j in range(self.velocity_bin_centers.shape[0]):
+                channels.append(
+                    ChannelVDEM(temperature_bin_edges=self.temperature_bin_edges[[i,i+1]],
+                                velocity_bin_edges=self.velocity_bin_edges[[j,j+1]])
+                )
+        return channels
+
+    @property
+    @u.quantity_input
+    def temperature_bin_centers(self) -> u.K:
+        return (self.temperature_bin_edges[1:] + self.temperature_bin_edges[:-1])/2
+
+    @property
+    @u.quantity_input
+    def velocity_bin_centers(self) -> u.Unit('km/s'):
+        return (self.velocity_bin_edges[1:] + self.velocity_bin_edges[:-1])/2
+
+    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)
+
+    def observe(self, *args, **kwargs):
+        kwargs['observer'] = self.observer
+        return super().observe(*args, **kwargs)
+
+    @staticmethod
+    @return_quantity_as_tuple
+    def calculate_intensity_kernel(loop, channel, **kwargs):
+        observer = kwargs.get('observer')
+        T = loop.electron_temperature
+        n = loop.density
+        v_los = los_velocity(loop.velocity_xyz, observer)
+        in_temperature_bin = np.logical_and(T>=channel.temperature_bin_edges[0],
+                                            T<channel.temperature_bin_edges[1])
+        in_velocity_bin = np.logical_and(v_los>=channel.velocity_bin_edges[0],
+                                         v_los<channel.velocity_bin_edges[1])
+        bin_mask = np.where(np.logical_and(in_temperature_bin, in_velocity_bin), 1, 0)
+        kernel = n**2 * bin_mask
+        return kernel
+
+    def maps_to_cube(self, vdem, time_index):
+        """
+        Transform a set of DEM maps at a single time step to a `~ndcube.NDCube`
+        """
+        arrays = []
+        for j in range(self.velocity_bin_centers.shape[0]):
+            _arrays = []
+            for i in range(self.temperature_bin_centers.shape[0]):
+                key = self.channels[0]._make_channel_label(
+                    self.temperature_bin_edges[[i,i+1]],
+                    self.velocity_bin_edges[[j,j+1]]
+                )
+                _map = vdem[key][time_index]
+                _arrays.append(_map.data)
+            arrays.append(_arrays)
+        arrays = np.array(arrays)
+        compound_wcs = extend_celestial_wcs(
+            _map.wcs,
+            (self.temperature_bin_centers, 'temperature', 'phys.temperature'),
+            (self.velocity_bin_centers, 'velocity', 'phys.velocity'),
+        )
+        return ndcube.NDCube(arrays, wcs=compound_wcs, meta=_map.meta, unit=_map.unit)
+
+
 class InstrumentQuantityBase(InstrumentBase):
 
     @u.quantity_input
     def __init__(self, *args, **kwargs):
-        self.channels = [ChannelBase(name=self.name)]
         super().__init__(*args, average_over_los=True, **kwargs)
 
+    @property
+    def channels(self):
+        return [ChannelBase(name=self.name)]
+
 
 class InstrumentLOSVelocity(InstrumentQuantityBase):
     name = 'los_velocity'
@@ -140,10 +266,12 @@ class InstrumentLOSVelocity(InstrumentQuantityBase):
     @return_quantity_as_tuple
     def calculate_intensity_kernel(loop, *args, **kwargs):
         observer = kwargs.get('observer')
-        if observer is None:
-            raise ValueError('Must pass in observer to compute LOS velocity.')
         return los_velocity(loop.velocity_xyz, observer)
 
+    def observe(self, *args, **kwargs):
+        kwargs['observer'] = self.observer
+        return super().observe(*args, **kwargs)
+
     @property
     def _expected_unit(self):
         return u.km / u.s
diff --git a/synthesizAR/instruments/sdo.py b/synthesizAR/instruments/sdo.py
index e408485..7d28406 100644
--- a/synthesizAR/instruments/sdo.py
+++ b/synthesizAR/instruments/sdo.py
@@ -40,28 +40,15 @@ class InstrumentSDOAIA(InstrumentBase):
     name = 'SDO_AIA'
 
     def __init__(self, observing_time, observer, **kwargs):
-        self.channels = [
-            Channel(94*u.angstrom),
-            Channel(131*u.angstrom),
-            Channel(171*u.angstrom),
-            Channel(193*u.angstrom),
-            Channel(211*u.angstrom),
-            Channel(335*u.angstrom),
-        ]
-        # Add the Gaussian width for the PSF convolution
-        psf_params = filter_mesh_parameters(use_preflightcore=True)
-        for c in self.channels:
-            psf_width = psf_params[c.channel]['width']
-            c.psf_width = u.Quantity([psf_width, psf_width])
-        super().__init__(observing_time, observer, **kwargs)
-
-    @property
-    def cadence(self) -> u.s:
-        return 12.0 * u.s
-
-    @property
-    def resolution(self) -> u.arcsec / u.pix:
-        return [0.600698, 0.600698] * u.arcsec/u.pixel
+        resolution = kwargs.pop('resolution', [0.600698, 0.600698] * u.arcsec/u.pixel)
+        cadence = kwargs.pop('cadence', 12.0 * u.s)
+        super().__init__(
+            observing_time=observing_time,
+            observer=observer,
+            resolution=resolution,
+            cadence=cadence,
+            **kwargs,
+        )
 
     @property
     def observatory(self):
@@ -75,6 +62,23 @@ class InstrumentSDOAIA(InstrumentBase):
     def telescope(self):
         return 'SDO/AIA'
 
+    @property
+    def channels(self):
+        channels = [
+            Channel(94*u.angstrom),
+            Channel(131*u.angstrom),
+            Channel(171*u.angstrom),
+            Channel(193*u.angstrom),
+            Channel(211*u.angstrom),
+            Channel(335*u.angstrom),
+        ]
+        # Add the Gaussian width for the PSF convolution
+        psf_params = filter_mesh_parameters(use_preflightcore=True)
+        for c in channels:
+            psf_width = psf_params[c.channel]['width']
+            c.psf_width = u.Quantity([psf_width, psf_width])
+        return channels
+
     @property
     def _expected_unit(self):
         return u.DN / (u.pix * u.s)
diff --git a/synthesizAR/instruments/util.py b/synthesizAR/instruments/util.py
index ee7fcd6..fdda6fc 100644
--- a/synthesizAR/instruments/util.py
+++ b/synthesizAR/instruments/util.py
@@ -42,19 +42,25 @@ def add_wave_keys_to_header(wavelength_array: u.angstrom, header):
     return header_copy
 
 
-def extend_celestial_wcs(celestial_wcs, array, name, physical_type):
+def extend_celestial_wcs(celestial_wcs, *extra_coords, **kwargs):
     """
-    Add an additional 3rd axis corresponding to a `~astropy.units.Quantity`
-    to a 2D celestial WCS
+    Add additional axes as extra coords to an existing WCS
+
+    Parameters
+    ----------
+    celestial_wcs: `~astropy.wcs.WCS`
+    extra_coords: `list`
+        Each member of the list can be an extra coord or a tuple.
+        If the latter, it is assumed this is a three-tuple specifying
+        the array, name, and physical type of the `QuantityTableCoordinate`
     """
-    temp_table = QuantityTableCoordinate(array,
-                                         names=name,
-                                         physical_types=physical_type)
-    mapping = list(range(celestial_wcs.pixel_n_dim))
-    mapping.extend(
-        [celestial_wcs.pixel_n_dim] * temp_table.wcs.pixel_n_dim
-    )
-    return CompoundLowLevelWCS(celestial_wcs, temp_table.wcs, mapping=mapping)
+    wcses = []
+    for ec in extra_coords:
+        if isinstance(ec, tuple):
+            array, name, physical_type = ec
+            ec = QuantityTableCoordinate(array, names=name, physical_types=physical_type)
+        wcses.append(ec.wcs)
+    return CompoundLowLevelWCS(celestial_wcs, *wcses, **kwargs)
 
 
 def read_cube_from_dataset(filename, axis_name, physical_type):
diff --git a/synthesizAR/models/magnetic.py b/synthesizAR/models/magnetic.py
new file mode 100644
index 0000000..bd64abe
--- /dev/null
+++ b/synthesizAR/models/magnetic.py
@@ -0,0 +1,37 @@
+"""
+Analytic magnetic field profiles
+"""
+import astropy.units as u
+import numpy as np
+
+__all__ = ['mikic_magnetic_field_profile']
+
+
+@u.quantity_input
+def mikic_magnetic_field_profile(s: u.cm,
+                                 length: u.Mm,
+                                 ell: u.Mm=14*u.Mm,
+                                 B_0: u.G=1*u.G,
+                                 B_1: u.G=10*u.G):
+    r"""
+    Magnetic field strength along a strand as prescribed by Eq. 4 of :cite:t:`mikic_importance_2013`
+
+    The analytic magnetic field profile along the strand as a function of the field-aligned coordinate
+    :math:`s` is given by,
+
+    .. math::
+
+        B(s) = B_0 + B_1(\exp{-s/\ell} + \exp{-(L-s)/\ell})
+
+    Parameters
+    ----------
+    s: `~astropy.units.Quantity`
+        Field-aligned coordinate
+    length: `~astropy.units.Quantity`
+        Strand length
+    ell: `~astropy.units.Quantity`, optional
+        Scale height of the magnetic field
+    B_0: `~astropy.units.Quantity`, optional
+    B_1: `~astropy.units.Quantity`, optional
+    """
+    return B_0 + B_1*(np.exp(-s/ell) + np.exp(-(length-s)/ell))
diff --git a/synthesizAR/strand.py b/synthesizAR/strand.py
index 874684e..54e15b7 100644
--- a/synthesizAR/strand.py
+++ b/synthesizAR/strand.py
@@ -177,7 +177,7 @@ Simulation Type: {self.simulation_type}'''
     @property
     def n_s(self):
         """
-        Numb of points along the field-aligned coordinate
+        Number of points along the field-aligned coordinate
         """
         return self.field_aligned_coordinate_center.shape[0]
 
diff --git a/synthesizAR/util/util.py b/synthesizAR/util/util.py
index 81f56c7..a7824c9 100644
--- a/synthesizAR/util/util.py
+++ b/synthesizAR/util/util.py
@@ -9,13 +9,11 @@ import sunpy.coordinates
 import warnings
 
 from astropy.coordinates import SkyCoord
-from collections import namedtuple
 from scipy.interpolate import RegularGridInterpolator
 
 import synthesizAR
 
 __all__ = [
-    'SpatialPair',
     'los_velocity',
     'coord_in_fov',
     'find_minimum_fov',
@@ -27,9 +25,6 @@ __all__ = [
 ]
 
 
-SpatialPair = namedtuple('SpatialPair', 'x y z')
-
-
 @u.quantity_input
 def los_velocity(v_xyz: u.cm/u.s, observer):
     """


Test patch:
diff --git a/synthesizAR/conftest.py b/synthesizAR/conftest.py
index 9d0853c..1b83d01 100644
--- a/synthesizAR/conftest.py
+++ b/synthesizAR/conftest.py
@@ -13,9 +13,13 @@ from synthesizAR.models import semi_circular_arcade, semi_circular_loop
 
 
 @pytest.fixture
-def bare_skeleton():
-    observer = get_earth(time='2020-01-01T00:00:00')
-    arcade = semi_circular_arcade(100*u.Mm, 20*u.deg, 10, observer)
+def earth_observer():
+    return get_earth(time='2020-01-01T00:00:00')
+
+
+@pytest.fixture
+def bare_skeleton(earth_observer):
+    arcade = semi_circular_arcade(100*u.Mm, 20*u.deg, 10, earth_observer)
     loops = [synthesizAR.Strand(f'{i}', c) for i, c in enumerate(arcade)]
     return synthesizAR.Skeleton(loops)
 
@@ -27,6 +31,18 @@ def skeleton_with_model(bare_skeleton):
     return bare_skeleton
 
 
+@pytest.fixture
+def skeleton_with_time_dependent_model(bare_skeleton):
+    from ebtelplusplus.models import HeatingModel, TriangularHeatingEvent
+
+    from synthesizAR.interfaces.ebtel import EbtelInterface
+    heating_model = HeatingModel(partition=1)
+    heating_model.events = [TriangularHeatingEvent(0*u.s, 200*u.s, 5e-3*u.Unit('erg cm-3 s-1'))]
+    interface = EbtelInterface(5e3*u.s, heating_model=heating_model)
+    bare_skeleton.load_loop_simulations(interface)
+    return bare_skeleton
+
+
 @pytest.fixture
 def semi_circle_strand(bare_skeleton):
     coords = semi_circular_loop(length=100*u.Mm)
diff --git a/synthesizAR/instruments/tests/test_physical.py b/synthesizAR/instruments/tests/test_physical.py
new file mode 100644
index 0000000..4b029cb
--- /dev/null
+++ b/synthesizAR/instruments/tests/test_physical.py
@@ -0,0 +1,71 @@
+"""
+Tests for instrument classes that compute physical quantities
+"""
+import astropy.units as u
+import ndcube
+import numpy as np
+import sunpy.map
+
+from synthesizAR.instruments import (
+    InstrumentDEM,
+    InstrumentLOSVelocity,
+    InstrumentTemperature,
+    InstrumentVDEM,
+)
+
+# smoke tests to ensure physical instrument classes work
+# NOTE: This does not test correctness of the output
+
+
+def test_dem_instrument(skeleton_with_model, earth_observer):
+    temperature_bin_edges = np.geomspace(1e5, 1e7, num=21,endpoint=True)*u.K
+    dem = InstrumentDEM([0,]*u.s,
+                        earth_observer,
+                        [1,1]*u.arcsec/u.pixel,
+                        temperature_bin_edges=temperature_bin_edges)
+    dem_maps = dem.observe(skeleton_with_model)
+    assert len(dem_maps) == temperature_bin_edges.shape[0]-1
+    for k, v in dem_maps.items():
+        assert len(v) == 1
+        assert isinstance(v[0], sunpy.map.GenericMap)
+        assert v[0].unit.is_equivalent('cm-5')
+    dem_cube = dem.maps_to_cube(dem_maps, 0)
+    assert isinstance(dem_cube, ndcube.NDCube)
+
+
+def test_vdem_instrument(skeleton_with_time_dependent_model, earth_observer):
+    temperature_bin_edges=10**np.arange(5.6,6.8,0.2)*u.K
+    velocity_bin_edges=np.arange(-20,30,10)*u.km/u.s
+    vdem = InstrumentVDEM([50]*u.s,
+                          earth_observer,
+                          [1,1]*u.arcsec/u.pixel,
+                          temperature_bin_edges=temperature_bin_edges,
+                          velocity_bin_edges=velocity_bin_edges)
+    vdem_maps = vdem.observe(skeleton_with_time_dependent_model)
+    assert len(vdem_maps) == len(vdem.channels)
+    for k, v in vdem_maps.items():
+        assert len(v) == 1
+        assert isinstance(v[0], sunpy.map.GenericMap)
+        assert v[0].unit.is_equivalent('cm-5')
+    vdem_cube = vdem.maps_to_cube(vdem_maps, 0)
+    assert isinstance(vdem_cube, ndcube.NDCube)
+
+
+def test_velocity_instrument(skeleton_with_model, earth_observer):
+    los_vel = InstrumentLOSVelocity([0,]*u.s, earth_observer, [1,1]*u.arcsec/u.pixel)
+    vel_maps = los_vel.observe(skeleton_with_model)
+    assert len(vel_maps) == 1
+    for k, v in vel_maps.items():
+        assert len(v) == 1
+        assert isinstance(v[0], sunpy.map.GenericMap)
+        assert v[0].unit.is_equivalent('km s-1')
+
+
+def test_temperature_instrument(skeleton_with_model, earth_observer):
+    temperature = InstrumentTemperature([0,]*u.s, earth_observer, [1,1]*u.arcsec/u.pixel)
+    vel_maps = temperature.observe(skeleton_with_model)
+    assert len(vel_maps) == 1
+    for k, v in vel_maps.items():
+        assert len(v) == 1
+        assert isinstance(v[0], sunpy.map.GenericMap)
+        assert v[0].unit.is_equivalent('K')
diff --git a/synthesizAR/instruments/tests/test_sdo.py b/synthesizAR/instruments/tests/test_sdo.py
new file mode 100644
index 0000000..801e8a2
--- /dev/null
+++ b/synthesizAR/instruments/tests/test_sdo.py
@@ -0,0 +1,41 @@
+"""
+Tests for SDO instrument classes
+"""
+import astropy.units as u
+import numpy as np
+import pytest
+import sunpy.map
+
+from synthesizAR.instruments.sdo import InstrumentSDOAIA
+
+
+@pytest.mark.remote_data
+def test_aia_instrument(skeleton_with_model, earth_observer):
+    aia = InstrumentSDOAIA([0,]*u.s, earth_observer)
+    aia_maps = aia.observe(skeleton_with_model)
+    for k, v in aia_maps.items():
+        assert len(v) == len(aia.observing_time)
+        for m in v:
+            assert isinstance(m, sunpy.map.GenericMap)
+            assert m.unit.is_equivalent('DN pix-1 s-1')
+
+@pytest.mark.remote_data
+@pytest.mark.parametrize(('input', 'result'), [
+    ([0,]*u.s, [0,]*u.s),
+    ([0,1]*u.s, [0,]*u.s),
+    ([0,25]*u.s, [0,12,24]*u.s),
+    (np.arange(2000,3000,60)*u.s, np.arange(2000,3000,60)*u.s),
+])
+def test_aia_observing_times(earth_observer, input, result):
+    aia = InstrumentSDOAIA(input, earth_observer)
+    assert u.allclose(aia.observing_time, result)
+
+
+@pytest.mark.remote_data
+def test_aia_resolution(earth_observer):
+    aia = InstrumentSDOAIA([0,]*u.s, earth_observer)
+    assert u.allclose(aia.resolution, [0.600698, 0.600698]*u.arcsec/u.pix)
+    aia = InstrumentSDOAIA([0,]*u.s,
+                           earth_observer,
+                           resolution=[1.2,1.2] * u.arcsec / u.pix)
+    assert u.allclose(aia.resolution, [1.2,1.2] * u.arcsec / u.pix)