# This program reads in data from Henley and Shelton (2012) and plots OVII line intensity versus Galactic Longitude. It also plots line intensity on an Aitoff projection in Galactic coordinates.
import pyfits
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

hdulist = pyfits.open('table2.fit') 	# Open the fits file
tbdata = hdulist[1].data		# Load the data from the table
cols = hdulist[1].columns		# Load column names

# Load individual columns of data
glon = tbdata.field('GLON')		# Galactic longitude
glat = tbdata.field('GLAT')		# Galactic latitude
o7, eo7 = tbdata.field('I7'), tbdata.field('e7')	# OVII intensity, error
o8, eo8 = tbdata.field('I8'), tbdata.field('eI8')	# OVIII intensity, error

# Filter out bad/null data
good = np.where((o7 >= 0.0))
lon_rad, lat_rad = glon[good], glat[good]

# Wrap data to conform to Aitoff projection
wrap1 = np.where(lon_rad<=180.0)
wrap2 = np.where(lon_rad>180.0)

lon_rad1, lat_rad1 = lon_rad[wrap1]*np.pi/180.0, lat_rad[wrap1]*np.pi/180.0
lon_rad2, lat_rad2 = (lon_rad[wrap2]-360.0)*np.pi/180.0, lat_rad[wrap2]*np.pi/180.0

plt.ion()
# Plotting figures
plt.figure()
plt.plot(glon[good],o7[good],'go',label='OVII')
plt.plot(glon[good],o8[good],'bo',label='OVIII')
plt.xlabel('Galactic Longitude (deg)')
plt.ylabel('Line intensity (L.U.)')
plt.legend()

fig = plt.figure()
ax = plt.subplot(111,projection="aitoff")
c1 = plt.scatter(lon_rad1,lat_rad1,c=o7[good][wrap1],s=36,edgecolors='none',cmap=plt.cm.rainbow)
c2 = plt.scatter(lon_rad2,lat_rad2,c=o7[good][wrap2],s=36,edgecolors='none',cmap=plt.cm.rainbow)
plt.colorbar(c1,orientation='horizontal',aspect=40,pad=0.05)
plt.title('OVII Line Intensities',y=1.1)
plt.grid(True)
plt.show()