from yt.mods import *

from matplotlib import rc
rc('text', usetex=True)
rc('font', family='serif')

import pylab as pl
#from scipy import *
#import scipy.signal
from h5py import File
#import pyfits
import os
import string
from matplotlib import pyplot 
from matplotlib import pyplot as plt
from matplotlib import pylab
from matplotlib import colors
#from pyreadcol import *
import numpy as np
from matplotlib.ticker import NullFormatter
import pickle

n=93

#f=open('../../Ionization/FOF/RD%04i_escape_fraction.txt'%n,'r')
f=open('../../Ionization/FOF/RedshiftOutput%04i-halos.dat'%n,'r')
lines=f.readlines()
f.close()

f=open('Luminosity.out','r')
lines1=f.readlines()
f.close()

Mass = []
Lum = []
Mag = []
ii = 0
for line in lines1:
    line1=lines[int(line.split()[0])+1]
    if float(line1.split()[0])==-1.0: continue
    Mass.append(float(line1.split()[4]))
    Lum.append(float(line.split()[1]))
    Mag.append(float(line.split()[3]))
    

#for i in xrange(1,len(lines)):
# if i-1 != int(lines1[ii].split()[0]): continue
# Mass.append(float(lines[i].split()[2]))
# Lum.append(float(lines1[ii].split()[1]))
# Mag.append(float(lines1[ii].split()[3])) 
# ii=ii+1
# if ii >= len(lines1): break

Mass=np.array(Mass)
Lum=np.array(Lum)
Mag=np.array(Mag)
x_max = -1 # int(Mass.max()+1)
x_min = -22

y_max = 5
y_min = -1

x_bins = 63    #60
y_bins = 60  #200

filename = 'Magnitude_MLratio'
fields = ["HaloNumber"]

fields_label = ['Count']
fields_min = [0]

#prof2d = BinnedProfile2D(sphere, d_bins, "Density", d_min, d_max, True,beta_bins, "beta", beta_min, beta_max, True)
data,xvals,yvals = pylab.histogram2d(Mag,np.log10(Mass/Lum),bins=[x_bins,y_bins],range=[[x_min,x_max],[y_min,y_max]])

f=open('/Users/haoxu/Desktop/Void_region/no_BG/RD0028/SED/Mass_Magnitude_Luminosity.aj','r')
lines=f.readlines()
f.close()
Mag1=[]
MLratio1=[]
for line in lines:
    Mag1.append(float(line.split()[2]))
    MLratio1.append(float(line.split()[1])/float(line.split()[3]))

for field, cbarlabel,field_min in zip(fields,fields_label,fields_min):
    ylabel = r'log$_{10}$ (mass-to-light ratio [$M_\odot/L_\odot]$)'
    xlabel = r'$M_V$'

    dy = yvals[1]-yvals[0]
    dx = xvals[1]-xvals[0]
    xvals1 = xvals[0:-1]+dx/2.0
    yvals1 = yvals[0:-1]+dy/2.0

    #data /= dlogx*dlogy

    intmin = 0
    intmax = data[0:,0:].max()

    fig = pl.figure(figsize=[8,6])
    ax = fig.add_subplot(111)

    mynorm = colors.Normalize([intmin, intmax])
    #im = ax.contourf(xvals1, yvals1, data.T, np.linspace(intmin,intmax,3), cmap=pl.cm.binary, norm=mynorm)

    #if xvals.size == data.T.shape[0]:
    #    xvals = np.append(xvals, 2.*10.**np.log10(xvals[-1]) - 10.**np.log10(xvals[-2]))
    #if yvals.size == data.T.shape[0]:
    #    yvals = np.append(yvals, 2.*10.**np.log10(yvals[-1]) - 10.**np.log10(yvals[-2]))

    fontsize=24
    extent = [xvals[0], xvals[-1], yvals[0], yvals[-1]]

    # put in ,norm=colors.LogNorm() if you want colorbar to log
    im = ax.imshow(data.T, cmap=pl.cm.binary,extent=extent, interpolation='nearest',\
               origin='lower',aspect='auto')

    ax.scatter(Mag1,np.log10(MLratio1),alpha=0.1)

    #X, Y = np.meshgrid(xvals,yvals)
    #im = ax.imshow(data,cmap=pl.cm.binary, norm=mynorm)
    #ax.set_xscale('log')
    #ax.set_yscale('log')
    cbar = fig.colorbar(im, fraction=0.05, format = '%i', ticks=np.linspace(intmin,intmax, (intmax-intmin)+1), pad=0.01, shrink=0.9)
    cbar.set_label('%s' % cbarlabel)

    pl.xlim(xvals.max(), xvals.min())
    #pl.ylim(yvals.min(), yvals.max())
    pl.ylim(yvals.min(), yvals.max())

    pl.xlabel('%s' % xlabel)
    pl.ylabel('%s'% ylabel)
    pl.savefig('RD%04i_%s.eps'% (n,filename),format='eps',dpi=160)
    fig.clf()