<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 AIA Response Functions
Step2:  SSW Response Functions
Calculate the wavelength and temperature response functions using SSW.
Step3:  Detailed Channel Information
Note first that we can easily read detailed information about each of the AIA channels. Here, I've used instrument files from local SSW installation.
Step4:  Wavelength Response Functions
Next, we'll actually instantiate the response object and calculate the wavelength response functions. We'll go ahead and calculate the response for all 7 EUV channels.
Step5:  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.
Step6:  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...
Temperature Response Functions
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.
Step7:  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,
Step8:  Or no Fe at all.
Step9:  Calculating a Custom Emission Table
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.
Step10:  Detailed Information about Emissivity Tables
The EmissTableInterface object provides a convenient way to inspect the different parts of the atomic data used to calculate the response functions, including the temperature, density, and line and continuum contributions.
<ASSISTANT_TASK:>
Python Code:
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
Explanation: AIA Response Functions: Comparing SunPy and SSW Results
End of explanation
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
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']
ssw_runner = hissw.ScriptMaker(ssw_pkg_list=['sdo/aia'],ssw_path_list=['aia'])
ssw_response_functions = ssw_runner.run([(ssw_script,{})],save_vars=save_vars,cleanup=True,verbose=True)
Explanation: SSW Response Functions
Calculate the wavelength and temperature response functions using SSW.
End of explanation
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')])
info_table
Explanation: Detailed Channel Information
Note first that we can easily read detailed information about each of the AIA channels. Here, I've used instrument files from local SSW installation.
End of explanation
aia_response = Response()
Explanation: Wavelength Response Functions
Next, we'll actually instantiate the response object and calculate the wavelength response functions. We'll go ahead and calculate the response for all 7 EUV channels.
End of explanation
wavelength_response = aia_response.calculate_wavelength_response()
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()
Explanation: 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.
End of explanation
temperature_response = aia_response.calculate_temperature_response()
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()
Explanation: 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...
Temperature Response Functions
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.
End of explanation
fe_ions = [ion for ion in ch_data.MasterList if ion.split('_')[0]=='fe']
Explanation: 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,
End of explanation
no_fe_ions = [ion for ion in ch_data.MasterList if ion.split('_')[0]!='fe']
temperature_response_fe_only = aia_response.calculate_temperature_response(ion_list=fe_ions)
temperature_response_no_fe = aia_response.calculate_temperature_response(ion_list=no_fe_ions)
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()
Explanation: Or no Fe at all.
End of explanation
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'
make_emiss_table(emiss_table_filename,ion_list=high_ionization_state_ions,temperature=temperature,density=density,
                 notebook=True, abundance_file=abundance_file)
custom_temperature_response = aia_response.calculate_temperature_response(ion_list=high_ionization_state_ions,
                                                    emiss_table_file='/home/wtb2/Desktop/test_emiss_table.h5')
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()
Explanation: Calculating a Custom Emission Table
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.
End of explanation
emiss_table = EmissTableInterface('/home/wtb2/sunpy/data/aia_emiss_table.h5')
emiss_table.temperature
emiss_table.density
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')
Explanation: Detailed Information about Emissivity Tables
The EmissTableInterface object provides a convenient way to inspect the different parts of the atomic data used to calculate the response functions, including the temperature, density, and line and continuum contributions.
End of explanation
<END_TASK>