import numpy as np
from astropy.io import fits
from astropy.cosmology import FlatLambdaCDM
import astropy.units as u
import math
import matplotlib.pyplot as plt
from pylab import * # This program requires the matplotlib python library for the plots.
from matplotlib import rc
font = {'family': 'Droid Sans',
        'weight': 'normal',
        'size': 20}

cosmo = FlatLambdaCDM(H0=70 * u.km / u.s / u.Mpc, Om0=0.3)

hdul_garik = fits.open('garik7_rcsed.fits')
data_garik = hdul_garik[1].data
ml_ = fits.open('output.fits')
ml_data = ml_[0].data
f = open('list_fits_files_CM.csv')
Msol = 4.8
x_arr = []
y_arr = []
x_arr_lit = []
y_arr_lit = []
for line in f:
    args = line.split(',')
    path = args[2] + '/output_B.fits'
    try:
        hdulist = fits.open(path)
        hdr = hdulist[2].header
        mag_data = hdr['2_MAG'].split('+/-')
        hdulist.close()
        i = 0
        for row in data_garik:
            if str(row['bestObjID']) == args[0]:
                ml = ml_data[10,i]
                Dl = cosmo.luminosity_distance(float(row['z'])).value
                if (hdulist[1].header['FILTER']=='g.MP9401') & ('*' not in str(hdr['2_AR'])):
                    lum = float(mag_data[0].replace("*", "")) - 5*math.log(Dl)/math.log(10) - 25 - row['kcorr_i']
                    L = 10**(0.4*(Msol - lum))
                    print row['MBH'], ',',L,',',ml
                    x_arr.append(L*ml)
                    y_arr.append(row['MBH'])
                if row['literature']:
                    x_arr_lit.append(L*ml)
                    y_arr_lit.append(row['MBH'])   
            i=i+1
    except IOError:
        pass

plot(x_arr,y_arr,'ro')
plot(x_arr_lit,y_arr_lit,'bo')
plt.yscale('log')
plt.xscale('log')
plt.savefig('foo.png', dpi=100)