<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 Temperature Response Functions
In this notebook, we'll use the code already on my PR branch to first calculate the wavelength response functions. Then, we'll try to use CHIANTI and ChiantiPy to calculate the temperature response functions for a few ions. 
We need to come up with a way to easily constrain how many wavelengths (or more precisely which ions) we need for each channel. We can certainly just search through every ion but this takes time!
Step1:  Wavelength Response
Step2:  Temperature Response
Set a temperature and density range. This is the range of temperatures and densities over which the contribution function for each ion will be calculated. Note that this is not a grid of temperatures and densities, but rather a list of $T$ and $n$ pairs.
According to Boerner et al. (2012), we should use a constant pressure of $10^{15}$ cm$^{-3}$ K.
Step3:  The main question is
Step4:  Boerner et al. (2012) use the coronal abundances of Feldman and Widing (1993) and the ionization balances of Dere et al. (2009). We also want to make sure we are calculating all of the emissivities in units of photons rather than ergs as this makes it easier when multiplying by the instrument response function.
Step5:  Now iterate through the ion list, calculating the emission and subsequently the contribution function at each stage and then interpolating the wavelength response function to the appropriate wavelengths.
Not that for the line emission, we are using the expression,
$$
G_X(\lambda,T) = \frac{1}{4\pi}\epsilon_X(\lambda,T)\mathrm{Ab}(X)\frac{N(X^{+m}}{N(X)}\frac{1}{n_e}
$$
which has units of photons cm$^{3}$ s$^{-1}$ sr$^{-1}$. The factor of $1/4\pi$ is included in the expression return by the emiss() method on the ChiantiPy ion object.
<ASSISTANT_TASK:>
Python Code:
import json
import numpy as np
import h5py
import seaborn as sns
from scipy.interpolate import splev,splrep
import matplotlib.pyplot as plt
import astropy.units as u
from sunpy.instr import aia
import ChiantiPy.core as ch
import ChiantiPy.tools.data as ch_data
%matplotlib inline
Explanation: AIA Temperature Response Functions
In this notebook, we'll use the code already on my PR branch to first calculate the wavelength response functions. Then, we'll try to use CHIANTI and ChiantiPy to calculate the temperature response functions for a few ions. 
We need to come up with a way to easily constrain how many wavelengths (or more precisely which ions) we need for each channel. We can certainly just search through every ion but this takes time!
End of explanation
response = aia.Response(ssw_path='/Users/willbarnes/Documents/Rice/Research/ssw/',
                        #channel_list=[131,171,193,211,304]
                       )
response.calculate_wavelength_response(include_crosstalk=True)
response.peek_wavelength_response()
data = np.loadtxt('../aia_sample_data/aia_wresponse_raw.dat')
channels = sorted(list(response.wavelength_response.keys()))
ssw_results = {}
for i in range(len(channels)):
    ssw_results[channels[i]] = {'wavelength':data[:,0],
                                'response':data[:,i+1]}
fig,axes = plt.subplots(3,3,figsize=(12,12))
for c,ax in zip(channels,axes.flatten()):
    #ssw
    ax.plot(ssw_results[c]['wavelength'],ssw_results[c]['response'],
            #color=response.channel_colors[c],
            label='ssw')
    #sunpy
    ax.plot(response.wavelength_response[c]['wavelength'],response.wavelength_response[c]['response'],
            #color=response.channel_colors[c],
            #marker='.',ms=6,markevery=5,
            label='SunPy',linestyle=':',alpha=0.95,lw=2)
    if c!=335 and c!=304:
        ax.set_xlim([c-20,c+20])
    #if c==335:
    #    ax.set_xlim([120,140])
    #if c==304:
    #    ax.set_xlim([80,100])
    ax.set_title('{} $\mathrm{{\mathring{{A}}}}$'.format(c),fontsize=20)
    ax.set_xlabel(r'$\lambda$ ({0:latex})'.format(response.wavelength_response[c]['wavelength'].unit),fontsize=20)
    ax.set_ylabel(r'$R_i(\lambda)$ ({0:latex})'.format(response.wavelength_response[c]['response'].unit),fontsize=20)
# contamination plots
#304
#ssw
ax = axes.flatten()[-2]
ax.plot(ssw_results[304]['wavelength'],ssw_results[304]['response'],
        #color=response.channel_colors[c],
        label='ssw')
#sunpy
ax.plot(response.wavelength_response[304]['wavelength'],response.wavelength_response[304]['response'],
        label='SunPy',linestyle='',alpha=0.95,lw=2,
        #color=response.channel_colors[c],
        marker='.',ms=8,markevery=2,
       )
ax.set_xlim([80,100])
ax.set_title('{} $\mathrm{{\mathring{{A}}}}$ contamination from 94'.format(304),fontsize=14)
ax.set_xlabel(r'$\lambda$ ({0:latex})'.format(response.wavelength_response[c]['wavelength'].unit),fontsize=20)
ax.set_ylabel(r'$R_i(\lambda)$ ({0:latex})'.format(response.wavelength_response[c]['response'].unit),fontsize=20)
#335
ax = axes.flatten()[-1]
ax.plot(ssw_results[335]['wavelength'],ssw_results[335]['response'],
        #color=response.channel_colors[c],
        label='ssw')
#sunpy
ax.plot(response.wavelength_response[335]['wavelength'],response.wavelength_response[335]['response'],
        label='SunPy',linestyle='',alpha=0.95,lw=2,
        #color=response.channel_colors[c],
        marker='.',ms=8,markevery=2,
       )
ax.set_xlim([120,140])
ax.set_title('{} $\mathrm{{\mathring{{A}}}}$ contamination from 131'.format(335),fontsize=14)
ax.set_xlabel(r'$\lambda$ ({0:latex})'.format(response.wavelength_response[c]['wavelength'].unit),fontsize=20)
ax.set_ylabel(r'$R_i(\lambda)$ ({0:latex})'.format(response.wavelength_response[c]['response'].unit),fontsize=20)
axes[0,0].legend(loc='best')
plt.tight_layout()
fig,axes = plt.subplots(3,3,figsize=(12,12),sharey=True,sharex=True)
for c,ax in zip(channels,axes.flatten()):
    #ssw
    ax2 = ax.twinx()
    ssw_interp = ssw_results[c]['response']*response.wavelength_response[c]['response'].unit
    delta_response = np.fabs(response.wavelength_response[c]['response'] - ssw_interp)/(ssw_interp)
    ax.plot(response.wavelength_response[c]['wavelength'],delta_response,
            #color=response.channel_colors[c]
           )
    ax2.plot(response.wavelength_response[c]['wavelength'],response.wavelength_response[c]['response'],
            color='k',linestyle='--')
    ax.set_title('{} $\mathrm{{\mathring{{A}}}}$'.format(c),fontsize=20)
    ax.set_xlabel(r'$\lambda$ ({0:latex})'.format(response.wavelength_response[c]['wavelength'].unit),fontsize=20)
    ax.set_ylabel(r'$\frac{|\mathrm{SSW}-\mathrm{SunPy}|}{\mathrm{SSW}}$',fontsize=20)
    ax2.set_ylabel(r'$R_i(\lambda)$ ({0:latex})'.format(response.wavelength_response[c]['response'].unit))
    ax.set_ylim([-1.1,1.1])
plt.tight_layout()
info_table = aia.response.aia_instr_properties_to_table([94,131,171,193,211,335],
                        ['/Users/willbarnes/Documents/Rice/Research/ssw/sdo/aia/response/aia_V6_all_fullinst.genx'])
Explanation: Wavelength Response
End of explanation
temperature = np.logspace(5,8,50)*u.K
pressure = 1e15*u.K*u.cm**(-3)
density = pressure/temperature
Explanation: Temperature Response
Set a temperature and density range. This is the range of temperatures and densities over which the contribution function for each ion will be calculated. Note that this is not a grid of temperatures and densities, but rather a list of $T$ and $n$ pairs.
According to Boerner et al. (2012), we should use a constant pressure of $10^{15}$ cm$^{-3}$ K.
End of explanation
ion_list = (['fe_{}'.format(i) for i in np.arange(6,26)] 
            + ['ca_{}'.format(i) for i in np.arange(10,20)])
Explanation: The main question is: how exactly do we calculate the temperature response functions? According to Boerner et al. (2012), the response function $K$ for channel $i$ is given by,
$$
K_i(T) = \int_0^{\infty}\mathrm{d}\lambda\,G(\lambda,T)R_i(\lambda)
$$
where $K$ has units DN cm$^{-5}$ s$^{-1}$ pix$^{-1}$. So what is the right expression for $G(\lambda,n,T)$, the contribution function? 
Let's divide the contribution function into a line emission part and a continuum part such that,
$$
G(\lambda,T) = G_{continuum}(\lambda,T) + \sum_XG_X(\lambda,T)
$$
where $X$ denotes an ion in the CHIANTI database. In this way, the response for each channel becomes,
$$
\begin{align}
K_i(T) &= \int_0^{\infty}\mathrm{d}\lambda\,\big(G_{continuum}(\lambda,T) + \sum_XG_X(\lambda,T)\big)R_i(\lambda) \
       &= \int_0^{\infty}\mathrm{d}\lambda\,G_{continuum}(\lambda,T)R_i(\lambda) + \sum_X\int_0^{\infty}\mathrm{d}\lambda\,G_X(\lambda,T)R_i(\lambda)
\end{align}
$$
So there is a contribution from the continuum plus a contribution from each ion at every wavelength where there is a spectral line. In particular, the continuum includes three different contributions: free-free losses, free-bound losses, and two-photon losses,
$$
G_{continuum}(\lambda,T) = G_{ff}(\lambda,T) + G_{fb}(\lambda,T) + G_{tp}(\lambda,T)
$$
The expressions for these losses can be found in Landi et al. (1999). These are calculated and summed over each ion as well. Thus, the full expression for the temperature response is given by,
$$
K_i(T) = \int_0^{\infty}\mathrm{d}\lambda\,(G_{ff}(\lambda,T) + G_{fb}(\lambda,T)R_i(\lambda) + G_{tp}(\lambda,T)) + \sum_X\int_0^{\infty}\mathrm{d}\lambda\,G_X(\lambda,T)R_i(\lambda)
$$
Make a list of ions. Exactly what ions should be included is not really clear. Certainly, at least all ions of Fe. Really though, just look at all the ions in CHIANTI. Though, the user should have a chance to select which ions are included so that they can easily calculate the response functions, given only a few lines/ions.
End of explanation
ch_data.Defaults['flux'] = 'photon'
ch_data.Defaults['abundfile'] = 'sun_coronal_1992_feldman'
ch_data.Defaults['ioneqfile'] = 'chianti'
Explanation: Boerner et al. (2012) use the coronal abundances of Feldman and Widing (1993) and the ionization balances of Dere et al. (2009). We also want to make sure we are calculating all of the emissivities in units of photons rather than ergs as this makes it easier when multiplying by the instrument response function.
End of explanation
temperature_responses = {k:np.zeros(len(temperature)) for k in response.wavelength_response}
for ion in ch_data.MasterList:
    #if ion.split('_')[0] != 'fe':
    #    continue
    print('{}: Calculating contribution function for {}'.format(ch_data.MasterList.index(ion),ion))
    #declare ion object
    tmp = ch.ion(ion,temperature=temperature.value,eDensity=density.value,
                 abundance='sun_coronal_1992_feldman')
    #calculate emissivity
    tmp.emiss()
    em = tmp.Emiss['emiss'][np.argsort(tmp.Emiss['wvl']),:]
    wvl = np.sort(tmp.Emiss['wvl'])
    #calculate contribution function.
    gofnt = tmp.Abundance*em*tmp.IoneqOne/tmp.EDensity
    #iterate over channels
    for channel in response.wavelength_response:
        #print('Adding to channel {}'.format(channel))
        #interpolate response function to transitions
        rsp = splev(wvl,splrep(response.wavelength_response[channel]['wavelength'].value,
                               response.wavelength_response[channel]['response'].value))
        rsp = np.where(rsp<0,0,rsp)*response._channel_info[channel]['plate_scale'].value
        #weighted sum over wavelength
        #add to temperature response
        temperature_responses[channel] += np.dot(rsp,gofnt)
# ssw responses
precalculated_responses_data = np.loadtxt('../aia_sample_data/aia_tresponse_raw.dat')
precalc_channels = [94,131,171,193,211,304,335]
precalculated_responses = {c: precalculated_responses_data[:,i+1] for i,c in enumerate(precalc_channels)}
precalculated_responses['temperature'] = precalculated_responses_data[:,0]
# ssw responses with chiantifix and evenorm fix
precalculated_responses_data = np.loadtxt('../aia_sample_data/aia_tresponse_fix.dat')
precalculated_responses_fix = {c: precalculated_responses_data[:,i+1] for i,c in enumerate(precalc_channels)}
precalculated_responses_fix['temperature'] = precalculated_responses_data[:,0]
channel_colors = {c: sns.color_palette('Set2',7)[i] for i,c in enumerate(response.wavelength_response)}
fig,axes = plt.subplots(4,2,figsize=(15,30),sharex=True)
for channel,ax in zip(sorted(list(temperature_responses.keys())),axes.flatten()):
    ax.plot(temperature,temperature_responses[channel]/(0.83*(1./4./np.pi)),
             label=r'lines',
             color=channel_colors[channel])
    ax.plot(temperature,(temperature_responses[channel]/(0.83*(1./4./np.pi))
                          + continuum_contributions[channel]),
             linestyle=':',
             label='lines + continuum',
             color=channel_colors[channel])
    ax.set_xscale('log')
    ax.set_yscale('log')
    ax.set_ylim([1e-30,2e-24])
    ax.set_xlim([1e5,1e8])
    ax.set_title(r'{} $\AA$'.format(channel))
    
for i,pc in enumerate(precalc_channels):
    axes.flatten()[i].plot(10**precalculated_responses['temperature'],
                            precalculated_responses[pc],
                            linestyle='--',
                            label=r'SSW',
                            color=channel_colors[pc])
    axes.flatten()[i].plot(10**precalculated_responses_fix['temperature'],
                        precalculated_responses_fix[pc],
                        linestyle='-.',
                        label=r'SSW with chiantifix',
                        color=channel_colors[pc])
axes[0,1].legend(loc='best')
continuum_contributions = {k:np.zeros(len(temperature)) for k in response.wavelength_response}
wvl = response.wavelength_response[94]['wavelength'].value
for ion in ch_data.MasterList:
    #if ion.split('_')[0] != 'fe':
    #    continue
    print('{}: Calculating contribution function for {}'.format(ch_data.MasterList.index(ion),ion))
    tmp = ch.ion(ion,temperature=temperature.value,eDensity=density.value,abundance='sun_coronal_1992_feldman')
    #two photon emiss
    tmp.twoPhoton(wvl)
    if 'rate' in tmp.TwoPhoton:
        two_photon = tmp.TwoPhoton['rate']
    else:
        two_photon = tmp.TwoPhoton['emiss']
    #free-free
    tmp_cont = ch.continuum(ion,temperature.value,abundance='sun_coronal_1992_feldman')
    if tmp_cont.Ion > 1:
        tmp_cont.freeFree(wvl)
        if 'rate' in tmp_cont.FreeFree:
            free_free = tmp_cont.FreeFree['rate']
        else:
            free_free = np.zeros((len(temperature),len(wvl)))
    else:
        free_free = np.zeros((len(temperature),len(wvl)))
    #free-bound
    if tmp_cont.Ion > 1:
        tmp_cont.freeBound(wvl)
        if 'rate' in tmp_cont.FreeBound:
            free_bound = tmp_cont.FreeBound['rate']
        else:
            free_bound = np.zeros((len(temperature),len(wvl)))
    else:
        free_bound = np.zeros((len(temperature),len(wvl)))
    #add to channels
    for channel in response.wavelength_response:
        continuum_contributions[channel] += np.dot((two_photon + free_free + free_bound),
                                                   (response.wavelength_response[channel]['response'].value
                                                    *response._channel_info[channel]['plate_scale'].value))
plt.figure(figsize=(8,8))
for channel in continuum_contributions:
    plt.plot(temperature,
             continuum_contributions[channel],
             label=channel,color=channel_colors[channel])
plt.xscale('log')
plt.yscale('log')
plt.xlim([1e5,1e8])
plt.ylim([1e-32,1e-27])
plt.legend(loc='best')
Explanation: Now iterate through the ion list, calculating the emission and subsequently the contribution function at each stage and then interpolating the wavelength response function to the appropriate wavelengths.
Not that for the line emission, we are using the expression,
$$
G_X(\lambda,T) = \frac{1}{4\pi}\epsilon_X(\lambda,T)\mathrm{Ab}(X)\frac{N(X^{+m}}{N(X)}\frac{1}{n_e}
$$
which has units of photons cm$^{3}$ s$^{-1}$ sr$^{-1}$. The factor of $1/4\pi$ is included in the expression return by the emiss() method on the ChiantiPy ion object.
End of explanation
<END_TASK>