Time-average Synthesized EIS Intensities

In this notebook, pull in intensity cubes ($N_x\times N_y\times N_{\lambda}$) for a given simulation time interval, $t=7500-12500$ s, and average it. Then save these intensity cubes to be analyzed later.



In [1]:

    
import os

from scipy.interpolate import interp1d,splev,splrep
import numpy as np
import astropy.units as u
import astropy.constants as const

from synthesizAR.instruments import InstrumentHinodeEIS
from synthesizAR.util import EISCube









    



/home/wtb2/anaconda3/envs/systematic-ar-study/lib/python3.5/site-packages/IPython/html.py:14: ShimWarning: The `IPython.html` package has been deprecated. You should import from `notebook` instead. `IPython.html.widgets` has moved to `ipywidgets`.
  "`IPython.html.widgets` has moved to `ipywidgets`.", ShimWarning)
WARNING: AstropyDeprecationWarning: astropy.utils.compat.odict.OrderedDict is now deprecated - import OrderedDict from the collections module instead [astropy.utils.compat.odict]
/home/wtb2/anaconda3/envs/systematic-ar-study/lib/python3.5/site-packages/matplotlib/__init__.py:1401: UserWarning:  This call to matplotlib.use() has no effect
because the backend has already been chosen;
matplotlib.use() must be called *before* pylab, matplotlib.pyplot,
or matplotlib.backends is imported for the first time.

  warnings.warn(_use_error_msg)



In [2]:

    
base1 = '/data/datadrive1/ar_forward_modeling/systematic_ar_study/noaa1109_tn{}'
base2 = '/data/datadrive2/ar_viz/systematic_ar_study/noaa1109_tn{}/'



In [3]:

    
eis = InstrumentHinodeEIS([7.5e3,1.25e4]*u.s)



In [4]:

    
frequencies = [250,750,'750-ion',2500,5000]



In [5]:

    
time_averaged_intensity_cubes = {'{}'.format(freq):{chan['name']:None for chan in eis.channels} for freq in frequencies}

Iterate over each time step and add intensity contribution from each timestep.



In [6]:

    
for freq in frequencies:
    print('t_N = {}'.format(freq))
    if type(freq) == int:
        base = base1
    else:
        base = base2
    for channel in eis.channels:
        print('channel = {}'.format(channel['name']))
        for i,t in enumerate(eis.observing_time):
            path2file = os.path.join(base.format(freq),
                                     'Hinode_EIS',channel['name'],
                                     'map_t{:06d}.h5'.format(i+750))
            if time_averaged_intensity_cubes['{}'.format(freq)][channel['name']] is None:
                time_averaged_intensity_cubes['{}'.format(freq)][channel['name']] = EISCube(path2file)
            else:
                time_averaged_intensity_cubes['{}'.format(freq)][channel['name']] += EISCube(path2file)









    



t_N = 250






    



WARNING: UnitsWarning: 'cm ct g(1/2) / (Angstrom erg(1/2) pix s2)' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]






    



channel = 170-210
channel = 250-290
t_N = 750
channel = 170-210
channel = 250-290
t_N = 750-ion
channel = 170-210
channel = 250-290
t_N = 2500
channel = 170-210






    



WARNING: UnitsWarning: 'cm cm ct g(1/2) / (Angstrom cm erg(1/2) pix s2)' contains multiple slashes, which is discouraged by the FITS standard [astropy.units.format.generic]






    



channel = 250-290
t_N = 5000
channel = 170-210
channel = 250-290

Convert to appropriate physical units and time-average.



In [7]:

    
for freq in frequencies:
    for channel in eis.channels:
        tmp = time_averaged_intensity_cubes['{}'.format(freq)][channel['name']]
        conversion = splev(tmp.wavelength.value,
                           splrep(channel['response']['x'].value,channel['response']['y'].value))
        conversion = conversion*channel['response']['y'].unit
        conversion /= const.h.cgs*const.c.cgs/tmp.wavelength.to(u.cm)/u.photon
        time_averaged_intensity_cubes['{}'.format(freq)][channel['name']] = tmp*(1./conversion)*(1./eis.observing_time.shape[0])

Save the results.



In [10]:

    
for freq in time_averaged_intensity_cubes.keys():
    for channel in time_averaged_intensity_cubes[freq].keys():
        time_averaged_intensity_cubes[freq][channel].save('../data/eis_intensity_{}_tn{}_t7500-12500.h5'.format(channel,freq))



In [ ]: