|
| 1 | +""" |
| 2 | +Simulating Time-dependent Emission from Impulsively Heated Loops with EBTEL |
| 3 | +=========================================================================== |
| 4 | +
|
| 5 | +This example demonstrates how to model the resulting AIA emission from an |
| 6 | +arcade of loops heated impulsively and modeled using the `ebtelplusplus` code. |
| 7 | +""" |
| 8 | +import astropy.time |
| 9 | +import astropy.units as u |
| 10 | +import matplotlib.pyplot as plt |
| 11 | +import numpy as np |
| 12 | +import sunpy.map |
| 13 | + |
| 14 | +from astropy.coordinates import SkyCoord |
| 15 | +from astropy.visualization import AsinhStretch, ImageNormalize, quantity_support |
| 16 | +from sunpy.coordinates import get_horizons_coord |
| 17 | + |
| 18 | +import synthesizAR |
| 19 | + |
| 20 | +from synthesizAR.instruments import InstrumentSDOAIA |
| 21 | +from synthesizAR.interfaces.ebtel import EbtelInterface |
| 22 | +from synthesizAR.interfaces.ebtel.heating_models import PowerLawNanoflareTrain |
| 23 | +from synthesizAR.models import semi_circular_arcade |
| 24 | + |
| 25 | +# sphinx_gallery_thumbnail_number = -1 |
| 26 | + |
| 27 | +########################################################################### |
| 28 | +# First, set up the coordinates for the arcade. The structure we will model |
| 29 | +# is an arcade of longer overlying loops with an arcade of successively |
| 30 | +# shorter loops underneath. |
| 31 | + |
| 32 | +obstime = astropy.time.Time('2022-11-14T22:00:00') |
| 33 | +pos = SkyCoord(lon=15*u.deg, |
| 34 | + lat=25*u.deg, |
| 35 | + radius=1*u.AU, |
| 36 | + obstime=obstime, |
| 37 | + frame='heliographic_stonyhurst') |
| 38 | +arcade_coords = [] |
| 39 | +delta_s = 0.3 * u.Mm |
| 40 | +for l in np.arange(25,150,25)*u.Mm: |
| 41 | + n_points = int(np.ceil((l/delta_s).decompose())) |
| 42 | + arcade_coords += semi_circular_arcade(l, 5*u.deg, 50, pos, n_points=n_points) |
| 43 | + |
| 44 | +########################################################################### |
| 45 | +# Next, build a `~synthesizAR.Skeleton` from the coordinates of the strands |
| 46 | +# in our arcade. |
| 47 | +strands = [synthesizAR.Strand(f'strand{i}', c) for i, c in enumerate(arcade_coords)] |
| 48 | +arcade = synthesizAR.Skeleton(strands) |
| 49 | + |
| 50 | +########################################################################### |
| 51 | +# We can visualize what this structure would look like as observed from |
| 52 | +# the Solar Dynamics Observatory. |
| 53 | +sdo_observer = get_horizons_coord('SDO', time=obstime) |
| 54 | +arcade.peek(observer=sdo_observer, |
| 55 | + axes_limits=[(175, 300)*u.arcsec, (300, 450)*u.arcsec]) |
| 56 | + |
| 57 | +########################################################################### |
| 58 | +# Next, we will model the hydrodynamic response to an impulsive heating event |
| 59 | +# on each strand using the `ebtelplusplus` code. We will simulate a total of |
| 60 | +# 3 h of simulation time where each loop is heated by a single event with an |
| 61 | +# energy chosen from a powerlaw distribution. |
| 62 | +event_model = PowerLawNanoflareTrain( |
| 63 | + [0,200]*u.s, 200*u.s, 0*u.s, [1e-3,1e-1]*u.Unit('erg cm-3 s-1'), -1.5 |
| 64 | +) |
| 65 | +ebtel = EbtelInterface(3*u.h, event_builder=event_model) |
| 66 | + |
| 67 | +########################################################################### |
| 68 | +# To attach the results of our loop simulation to each strand, we pass the |
| 69 | +# interface to the geometric model of our arcade we built above. |
| 70 | +arcade.load_loop_simulations(ebtel) |
| 71 | + |
| 72 | +########################################################################### |
| 73 | +# We can then visualize the temperature and density evolution of each strand |
| 74 | +# as a function of time. Note that because EBTEL is a spatially-averaged model, |
| 75 | +# it is assumed that eadch point along the strand has the same temperature and |
| 76 | +# density. |
| 77 | +with quantity_support(): |
| 78 | + fig = plt.figure(figsize=(10,5)) |
| 79 | + ax1 = fig.add_subplot(121) |
| 80 | + ax2 = fig.add_subplot(122) |
| 81 | + for s in arcade.strands: |
| 82 | + ax1.plot(s.time, s.electron_temperature[:,0], color='k', alpha=0.25) |
| 83 | + ax2.plot(s.time, s.density[:,0], color='k', alpha=0.25) |
| 84 | + |
| 85 | +########################################################################### |
| 86 | +# The last step is to use the temperature and density along each strand to |
| 87 | +# compute the emission as observed by the AIA instrument. We'll model the |
| 88 | +# emission from 500 s to 6000 s at a cadence of 50 s for the 193 Å channel. |
| 89 | +aia = InstrumentSDOAIA(np.arange(500,6e3,50)*u.s, |
| 90 | + sdo_observer, |
| 91 | + pad_fov=(20, 20)*u.arcsec) |
| 92 | +maps = aia.observe(arcade, channels=aia.channels[3:4]) |
| 93 | + |
| 94 | +########################################################################### |
| 95 | +# We can easily visualize this time-dependent emission using a |
| 96 | +# `~sunpy.map.MapSequence`. |
| 97 | +mseq = sunpy.map.MapSequence(maps['193'], sequence=True) |
| 98 | +fig = plt.figure() |
| 99 | +ax = fig.add_subplot(projection=mseq[0]) |
| 100 | +ani = mseq.plot(axes=ax, norm=ImageNormalize(vmin=0, vmax=5, stretch=AsinhStretch())) |
| 101 | + |
| 102 | +plt.show() |
0 commit comments