Modeling Intensities from Multiple Instruments

This example shows how to compute the synthetic intensities from three different observatories: AIA, XRT, and EUVI. It also demonstrates how to define a custom instrument class.

import astropy.units as u
from astropy.coordinates import SkyCoord
import numpy as np
import matplotlib.pyplot as plt
from sunpy.coordinates import get_earth, get_horizons_coord, HeliographicStonyhurst, Helioprojective
from astropy.visualization import quantity_support

import synthesizAR
from synthesizAR.models import semi_circular_arcade
from synthesizAR.interfaces import MartensInterface
from synthesizAR.instruments import InstrumentSDOAIA, InstrumentHinodeXRT
# sphinx_gallery_thumbnail_number = -1

First, we’ll set up the geometry for the active region that we are modeling. We’ll use a simple arcade of semi-circular loops all with a length of 150 Mm.

obstime = '2021-10-28T15:00:00'
loc = SkyCoord(HeliographicStonyhurst(lon=0*u.deg,lat=-30*u.deg, radius=1*u.R_sun, obstime=obstime))
arcade = semi_circular_arcade(150*u.Mm, 10*u.deg, 50, loc, gamma=90*u.deg, n_points=5000)
skeleton = synthesizAR.Skeleton([synthesizAR.Loop(f'{i}', c) for i,c in enumerate(arcade)])

We’ll select a few different observer locations for SDO and STERO-A and use sunpy.coordinates.get_horizons_coord to get these locations from the JPL Horizons service. We’ll approximate the location of Hinode as Earth since its location is not available through JPL Horizons.

sdo = get_horizons_coord('SDO', time=obstime)
stereo_a = get_horizons_coord('STEREO-A', time=obstime)
hinode = get_earth(time=obstime)

INFO: Obtained JPL HORIZONS location for Solar Dynamics Observatory (spac [sunpy.coordinates.ephemeris]
INFO: Obtained JPL HORIZONS location for STEREO-A (spacecraft) (-234) [sunpy.coordinates.ephemeris]

We can quickly peek at what the structure of our active region looks like from the viewpoints of SDO and STEREO-A.

skeleton.peek(observer=sdo, axes_limits=[(-500,500)*u.arcsec, (-1000,0)*u.arcsec])
skeleton.peek(observer=stereo_a, axes_limits=[(0,1000)*u.arcsec, (-1000,0)*u.arcsec])

• 
• 

Next, we want to calculate some sort of thermodynamic model for each one of these strands. We’ll use the MartensScalingLaws model with the heating rate chosen from a uniform distribution.

class MartensRandom(MartensInterface):

    def get_heating_constant(self, loop):
        h_a = 1e-5 * u.Unit('erg cm-3 s-1')
        h_b = 100*h_a
        return h_a + np.random.random_sample()*(h_b - h_a)

martens = MartensRandom(None)
skeleton.load_loop_simulations(martens)

We can visualize the electron temperature and density profiles for each loop.

with quantity_support():
    plt.figure(figsize=(11, 5))
    ax1 = plt.subplot(121)
    for l in skeleton.loops:
        plt.plot(l.field_aligned_coordinate_center.to('Mm'), l.electron_temperature[0].to('MK'), color='k')
    plt.subplot(122)
    for l in skeleton.loops:
        plt.plot(l.field_aligned_coordinate_center.to('Mm'), l.density[0], color='k')
    plt.yscale('log')

Let’s compute the emission that would be observed from these loops with this particular model for the temperature and density. First, we’ll compute the emission as observed in all channels of AIA. 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_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'],
                          fov_center=center, fov_width=(250, 250)*u.arcsec)
xrt_images = xrt.observe(skeleton)
for k in xrt_images:
    xrt_images[k][0].peek()

• 
• 

Lastly, we want to compute the emission as observed by EUVI on STEREO-A. Currently, synthesizAR.instruments does not contain an EUVI instrument class. However, we can easily define our own as follows. Note that here, for simplicity, we are using the 171 Å temperature response function for AIA as a proxy for the temperature response of the 171 Å channel on EUVI.

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

    @property
    def observatory(self):
        return 'STEREO A'

    @property
    def telescope(self):
        return 'STEREO'

    @property
    def detector(self):
        return 'EUVI'

    def get_instrument_name(self, *args):
        return 'SECCHI'

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 $AA$ channel.

euvi = InstrumentSTEREOEUVI([0, 1]*u.s, stereo_a, fov_center=center, fov_width=(250, 250)*u.arcsec)
euvi_images = euvi.observe(skeleton, channels=euvi.channels[2:3])
for k in euvi_images:
    euvi_images[k][0].peek()