Add instrument class for computing VDEM (#103)

* minor instr refactoring

* fix tests

* add analytic mikic B-field model

* remove SpatialPair

* add vdem instrument class

* add function for making VDEM cubes

* fix doc build

* more doc fixes

Author: wtbarnes

.rtd-environment.yml

@@ -5,3 +5,4 @@ dependencies:
   - python=3.12
   - pip
   - graphviz!=2.42.*,!=2.43.*
+  - sqlite<3.49.0

docs/references.bib

@@ -49,3 +49,19 @@ @article{reep_diagnosing_2013
   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}
+}

examples/ebtel-nanoflares.py

@@ -94,7 +94,7 @@
 ###########################################################################
 # 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()))

examples/multi-instrument.py

@@ -87,15 +87,15 @@ def get_heating_constant(self, loop):
 # 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()
 
 ###############################################################################
 # 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 @@ def get_heating_constant(self, loop):
 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 @@ def get_instrument_name(self, *args):
 # 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()

synthesizAR/conftest.py

@@ -13,9 +13,13 @@
 
 
 @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)

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 @@ def observing_time(self) -> u.s:
 
     @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 @@ def observer(self):
     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 @@ def observatory(self):
 
     @property
     def _expected_unit(self):
-        return None
+        raise NotImplementedError
+
+    @property
+    def channels(self):
+        raise NotImplementedError
 
     def get_instrument_name(self, channel):
         return self.name