In [104]:
import sys
import os
import numpy as np
import scipy.io
from astropy.table import Table
import matplotlib.pyplot as plt
import seaborn as sns
sns.set_context('notebook')
%matplotlib inline
Next, open up the instrument file. We'll use the most recent version of the files, V6. The _fullinst file should contain most of the information that we need. These .genx files have been read into IDL and then resaved as normal IDL .sav files.
In [122]:
v6_all_fullinst_genx = scipy.io.readsav('/Users/willbarnes/Documents/Projects/gsoc/aia_response/ssw_aia_response_data/aia_V6_all_fullinst')
Because of the way the files were saved, all of the information is inside the data keyword.
In [123]:
v6_all_fullinst_genx = v6_all_fullinst_genx['data']
The data structure returned here is a numpy.recarray. Here is a brief tutorial on how to use these data structures.
In [124]:
type(v6_all_fullinst_genx)
Out[124]:
What are the members and how do we access them?
In [125]:
v6_all_fullinst_genx.dtype.names
Out[125]:
Let's look at the 94 Å data. It turns out that each member is a numpy array of length one and that the first (and only) entry is again a numpy.recarray. So let's see what this data structure looks like.
NOTE: A<CHANNEL_NUM> just contains wavelength, effective area, and the platescale while A<CHANNEL_NUM>_FULL contains many more pieces of data, such as information about the primary and secondary mirrors and the CCD. Furthermore, A<CHANNEL_NUM>_FILE is just a single filename, maybe where this data lived originally though it is not clear where this file actually lives or if it is accessible at all.
In [126]:
v6_all_fullinst_genx['a94_full'][0].dtype.names
Out[126]:
In [127]:
v6_all_fullinst_genx['a94'][0]['wave']
Out[127]:
Try plotting the effective area as a function of $\lambda$ for each channel.
In [156]:
fig = plt.figure(figsize=(8,8))
ax = fig.gca()
cp = sns.color_palette('hls',int(len(v6_all_fullinst_genx['channels'][0])))
for channel,i in zip(v6_all_fullinst_genx['channels'][0],range(len(v6_all_fullinst_genx['channels'][0]))):
#skip the thick channels
if b'thick' in channel or b'THICK' in channel:
continue
print('Plotting channel %s'%(channel.decode('utf8')))
wavelength,a_eff = v6_all_fullinst_genx[channel.decode('utf8')][0]['wave'][0],v6_all_fullinst_genx[str(channel.decode('utf8'))][0]['ea'][0]
ax.plot(wavelength,a_eff/np.max(a_eff),label=channel.decode('utf8')[1::]+' angstrom',color=cp[i])
ax.set_xlabel(r'Wavelength (angstroms)')
ax.set_ylabel(r'Effective Area (normalized)')
ax.set_xlim([80,350])
ax.legend(loc='best')
Out[156]:
Similarly, we should be able to to calculate $A_{eff}(\lambda)$ according to the equation in section 2 of Boerner et al. (2012), $$ A_eff = A_{geo}R_pR_sT_eT_fD(t)Q $$ where
geoareaprimary, secondaryent_filter, fp_filtercontamccdThen, we are able to calculate the wavelength response function, $$ R_i(\lambda) = A_{eff,i}(\lambda,t)G(\lambda) $$ where $G(\lambda)=(12398/\lambda/3.65)g$ is the gain of the CCD-camera system in DN per photon. $g$ is the camera gain in DN per electron.
$g$ does not seem to be included in any of the data files. From Table 2 of Boerner et al. (2012),
In [129]:
chan = np.array([93.9, 131.2, 171.1, 195.1, 211.3, 303.8, 335.4])
G = np.array([2.128, 1.523, 1.168, 1.024, 0.946, 0.658, 0.596])
g = G/12398.0*3.65*chan
boerner_table_2 = Table([chan,G,g],names=('Channel', '$G$', '$g$'))
boerner_table_2['Channel'].unit = 'angstrom'
boerner_table_2['$G$'].unit = 'DN/photon'
boerner_table_2['$g$'].unit = 'DN/electron'
boerner_table_2
Out[129]:
So $g$ is roughly constant and is approximately $g\approx0.0588$
In [130]:
1.0/v6_all_fullinst_genx['a94_full'][0]['elecperdn'][0]
Out[130]:
So it seems that this parameter is included in the data file as $1/g$? However, it is not clear why this number is off by about $4\times10^{-3}$. For now let's just make a gain function that sets $g=0.0588$.
In [131]:
def ccd_gain(wavelength):
g = 0.0588
return 12398.0/wavelength/3.65*g
Now lets make a function that will calculate the effective area and return both the calculate version and that read from the file.
In [146]:
def wvl_response(data_struct):
response = {}
for channel in data_struct['channels'][0]:
if b'thick' in channel or b'THICK' in channel:
continue
full_channel = data_struct[channel.decode('utf8')+'_FULL'][0]
wave = full_channel['wave'][0]
effective_area_file = full_channel['effarea'][0]
a_geo,rp,rs,te,tf,d,q = full_channel['geoarea'][0],full_channel['primary'][0],full_channel['secondary'][0],full_channel['ent_filter'][0],full_channel['fp_filter'][0],full_channel['contam'][0],full_channel['ccd'][0]
G = ccd_gain(wave)
response_calc = a_geo*rp*rs*te*tf*d*q*G
response_file = effective_area_file*G
response[channel.decode('utf8')] = {'wave':wave,'file':response_file,'calc':response_calc}
return response
Now do the calculations and plot the functions on top of each other.
In [147]:
response = wvl_response(v6_all_fullinst_genx)
In [155]:
fig = plt.figure(figsize=(8,8))
ax = fig.gca()
cp = sns.color_palette('hls',len(response))
for key,i in zip(response,range(len(response))):
ax.plot(response[key]['wave'],response[key]['calc']/np.max(response[key]['calc']),color=cp[i],label=key)
ax.plot(response[key]['wave'],response[key]['file']/np.max(response[key]['file']),color=cp[i],linestyle='dashed')
ax.set_xlim([80,350])
ax.legend(loc='best')
ax.set_ylabel(r'$R_i(\lambda)$')
ax.set_xlabel(r'$\lambda$')
Out[155]:
Calculating the temperature response functions is now just a matter of pulling the relevant contribution functions from CHIANTI via ChiantiPy.
In [ ]: