'''
Written by Li Feng, Nov. 25, 2016. @NJNU
'''
# ------------------------ python code begins here

# Import necessary modules
import pyfits
import numpy as np
import matplotlib.pyplot as plt
 
#Read data
hdulist = pyfits.open('SDSS_filter_curves.fits')
data1 = hdulist[1].data
data2 = hdulist[2].data
data3 = hdulist[3].data
data4 = hdulist[4].data
data5 = hdulist[5].data
print(data1)
lamda_u = data1.field(0) # wavelength
respt_u = data1.field(1) # quantum efficiency for a point source
lamda_g = data2.field(0) # wavelength
respt_g = data2.field(1) # quantum efficiency for a point source
lamda_r = data3.field(0) # wavelength
respt_r = data3.field(1) # quantum efficiency for a point source
lamda_i = data4.field(0) # wavelength
respt_i = data4.field(1) # quantum efficiency for a point source
lamda_z = data5.field(0) # wavelength
respt_z = data5.field(1) # quantum efficiency for a point source
#print(lamda_u,respt_u)

#create plot

plt.plot(lamda_u, respt_u, 'g', linewidth=4.5)
plt.plot(lamda_g, respt_g, 'b', linewidth=4.5)
plt.plot(lamda_r, respt_r, 'k', linewidth=4.5)
plt.plot(lamda_i, respt_i, 'r', linewidth=4.5)
plt.plot(lamda_z, respt_z, 'c', linewidth=4.5)

# Axis Lables
#plt.ylabel('Flux ($\mathrm{10^{-13} ergs cm^{-2} s^{-1} \AA ^{-1}}$)', fontsize=15)
plt.ylabel('Filter response', fontsize=25)
plt.xlabel('Wavelength ($\mathrm{\AA}$)', fontsize=25)

# Define limits
plt.xlim(2800,11000)
plt.ylim(0,0.55)

plt.show()
# ------------------------ python code ends here