In [5]:
import os
import matplotlib.pyplot as plt
import seaborn as sns
import numpy as np
import astropy.units as u
import ChiantiPy.tools.data as ch_data
import hissw
from sunpy.instr.aia.response import Response
from sunpy.instr.aia.response_utils import EmissTableInterface,aia_instr_properties_to_table,make_emiss_table
%matplotlib inline
In [6]:
ssw_script="""
; wavelength response
wresp = aia_get_response(/area,/dn,/timedepend_date,/evenorm)
wavelength = wresp.wave
wresponse_94 = wresp.a94.ea
wresponse_131 = wresp.a131.ea
wresponse_171 = wresp.a171.ea
wresponse_193 = wresp.a193.ea
wresponse_211 = wresp.a211.ea
wresponse_304 = wresp.a304.ea
wresponse_335 = wresp.a335.ea
; temperature response
tresp = aia_get_response(/temp,/dn,/timedepend_date,/evenorm)
temperature = 10.^tresp.logte
tresponse_94 = tresp.a94.tresp
tresponse_131 = tresp.a131.tresp
tresponse_171 = tresp.a171.tresp
tresponse_193 = tresp.a193.tresp
tresponse_211 = tresp.a211.tresp
tresponse_304 = tresp.a304.tresp
tresponse_335 = tresp.a335.tresp
"""
In [7]:
save_vars = ['temperature',
'wavelength',
'tresponse_94',
'tresponse_131',
'tresponse_171',
'tresponse_193',
'tresponse_211',
'tresponse_304',
'tresponse_335',
'wresponse_94',
'wresponse_131',
'wresponse_171',
'wresponse_193',
'wresponse_211',
'wresponse_304',
'wresponse_335']
In [8]:
ssw_runner = hissw.ScriptMaker(ssw_pkg_list=['sdo/aia'],ssw_path_list=['aia'])
In [9]:
ssw_response_functions = ssw_runner.run([(ssw_script,{})],save_vars=save_vars,cleanup=True,verbose=True)
In [10]:
info_table = aia_instr_properties_to_table([94,131,171,193,211,335],
[os.path.join('/usr/local/ssw/','sdo/aia/response/aia_V6_all_fullinst.genx')])
In [11]:
info_table
Out[11]:
In [12]:
aia_response = Response()
When the Response object is instantiated, if no instrument files are supplied, it checks the the default SunPy downloads directory for the EUV and FUV instrument files. If they are not there, they are downloaded from the Goddard SSW mirror. This way, a local installation of SSW is not needed. However, to calculate the temperature response functions, you'll need CHIANTI and ChiantiPy installed.
In [13]:
wavelength_response = aia_response.calculate_wavelength_response()
In [14]:
fig,axes = plt.subplots(3,3,figsize=(12,12))
for channel,ax in zip(sorted(wavelength_response.keys()),axes.flatten()):
# sunpy
ax.plot(wavelength_response[channel].wavelength,wavelength_response[channel].response,
linestyle=':',lw=3,color=sns.color_palette('deep')[0],label='SunPy')
# SSW
ax.plot(ssw_response_functions['wavelength'],ssw_response_functions['wresponse_{}'.format(channel)],
linestyle='-',color=sns.color_palette('deep')[2],label='SSW')
ax.set_xlim([channel-20,channel+20])
ax.set_xlabel(r'$\lambda$ [{:latex}]'.format(wavelength_response[channel].wavelength.unit))
ax.set_ylabel(r'$R_{{{}}}$ [{:latex}]'.format(channel,wavelength_response[channel].response.unit))
axes[0,0].legend(loc='best')
plt.tight_layout()
Note that this shows the SSW AIA wavelength response functions with the time-dependent corrections so I guess we need to account for this somehow...
Now, calculate the temperature response functions. The first time you do this, it will build an emissivity table which may take around 30 minutes, but this only needs to be done once. Alternatively, you can also build your own emissivity table using the built-in function. This is necessary if, for example, you want to change the abundances or ionization equilibrium file, or calculate over a different temperature/density range.
In [15]:
temperature_response = aia_response.calculate_temperature_response()
In [16]:
fig,axes = plt.subplots(3,3,figsize=(12,12),sharex=True,sharey=True)
for channel,ax in zip(sorted(wavelength_response.keys()),axes.flatten()):
# sunpy
ax.plot(temperature_response[channel].temperature,temperature_response[channel].response,
linestyle=':',lw=3,color=sns.color_palette('deep')[0],label='SunPy')
# SSW
ax.plot(ssw_response_functions['temperature'],ssw_response_functions['tresponse_{}'.format(channel)],
linestyle='-',color=sns.color_palette('deep')[2],label='SSW')
ax.set_xlabel(r'$T$ [{:latex}]'.format(temperature_response[channel].temperature.unit))
ax.set_ylabel(r'$K_{{{}}}$ [{:latex}]'.format(channel,temperature_response[channel].response.unit))
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ylim([1e-30,2e-24])
ax.set_xlim([1e5,1e8])
axes[0,0].legend(loc='best')
plt.tight_layout()
Note that there is an issue with the 304 angstrom channel (even when the time-dependent correction to the wavelength response is not included). Perhaps there is some correction made in the 304 angstrom case, or ions that contribute in this wavlength range, that I've not accounted for here.
Oddly, the lack of a time-dependent correction in the SunPy case actually gives better agreement between the two in the 304 angstrom case.
You can also choose the number of ions that you'd like to include in the temperature response function calculation. So for example if wanted to see how only including Fe ions affected the response functions,
In [17]:
fe_ions = [ion for ion in ch_data.MasterList if ion.split('_')[0]=='fe']
Or no Fe at all.
In [18]:
no_fe_ions = [ion for ion in ch_data.MasterList if ion.split('_')[0]!='fe']
In [19]:
temperature_response_fe_only = aia_response.calculate_temperature_response(ion_list=fe_ions)
In [20]:
temperature_response_no_fe = aia_response.calculate_temperature_response(ion_list=no_fe_ions)
In [21]:
fig,axes = plt.subplots(3,3,figsize=(12,12),sharex=True,sharey=True)
for channel,ax in zip(sorted(wavelength_response.keys()),axes.flatten()):
# all_ions
ax.plot(temperature_response[channel].temperature,temperature_response[channel].response,
linestyle='-',color=sns.color_palette('deep')[0],label='All ions')
# Fe only
ax.plot(temperature_response_fe_only[channel].temperature,
temperature_response_fe_only[channel].response,
linestyle='-',color=sns.color_palette('deep')[2],label='Fe only')
# No Fe
ax.plot(temperature_response_no_fe[channel].temperature,
temperature_response_no_fe[channel].response,
linestyle='-',color=sns.color_palette('deep')[1],label='No Fe')
ax.set_xlabel(r'$T$ [{:latex}]'.format(temperature_response[channel].temperature.unit))
ax.set_ylabel(r'$K_{{{}}}$ [{:latex}]'.format(channel,temperature_response[channel].response.unit))
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ylim([1e-30,2e-24])
ax.set_xlim([1e5,1e8])
axes[0,0].legend(loc='best')
plt.tight_layout()
If you want to change the temperature/density/wavelength over which the emissivities are calculated or change things like the abundance, you can calculate your own emissivity table and then pass that file name in when calculating the temperature response functions.
In this short example, we'll calculate a new table using photospheric abundances, with constant density, and using only handful of ions with very high ionization states. To show the progress bar in the notebook, we need to pass in the notebook=True kwarg. Default is false and nothing will be displayed in notebook. For this to actually work though, you seem to need to have ipywidgets installed.
In [22]:
temperature = np.logspace(5.5,8,100)*u.K
density = 1e9*np.ones(temperature.shape)*u.cm**(-3)
abundance_file = 'sun_photospheric_1998_grevesse'
high_ionization_state_ions = [ion for ion in ch_data.MasterList if ion[-1] != 'd' and int(ion.split('_')[1]) >= 20]
emiss_table_filename = '/home/wtb2/Desktop/test_emiss_table.h5'
In [23]:
make_emiss_table(emiss_table_filename,ion_list=high_ionization_state_ions,temperature=temperature,density=density,
notebook=True, abundance_file=abundance_file)
In [24]:
custom_temperature_response = aia_response.calculate_temperature_response(ion_list=high_ionization_state_ions,
emiss_table_file='/home/wtb2/Desktop/test_emiss_table.h5')
In [25]:
fig,axes = plt.subplots(3,3,figsize=(12,12),sharex=True,sharey=True)
for channel,ax in zip(sorted(wavelength_response.keys()),axes.flatten()):
# all_ions
ax.plot(temperature_response[channel].temperature,temperature_response[channel].response,
linestyle='-',color=sns.color_palette('deep')[0],label='Original')
# Fe only
ax.plot(custom_temperature_response[channel].temperature,
custom_temperature_response[channel].response,
linestyle='-',color=sns.color_palette('deep')[2],label='Custom')
ax.set_xlabel(r'$T$ [{:latex}]'.format(temperature_response[channel].temperature.unit))
ax.set_ylabel(r'$K_{{{}}}$ [{:latex}]'.format(channel,temperature_response[channel].response.unit))
ax.set_xscale('log')
ax.set_yscale('log')
ax.set_ylim([1e-30,2e-24])
ax.set_xlim([1e5,1e8])
axes[0,0].legend(loc='best')
plt.tight_layout()
In [26]:
emiss_table = EmissTableInterface('/home/wtb2/sunpy/data/aia_emiss_table.h5')
In [27]:
emiss_table.temperature
Out[27]:
In [28]:
emiss_table.density
Out[28]:
In [29]:
for i in range(emiss_table['fe_12'].contribution_function.shape[0]):
plt.plot(emiss_table.temperature,emiss_table['fe_12'].contribution_function[:,i])
plt.xlabel(r'$T$ [{:latex}]'.format(emiss_table.temperature.unit))
plt.ylabel(r'$G(T)$ [{:latex}], Fe XII'.format(emiss_table['fe_12'].contribution_function.unit))
plt.ylim([1e-30,1e-18])
plt.yscale('log')
plt.xscale('log')
In [ ]: