from yt.mods import *

from matplotlib import rc
rc('text', usetex=True)
rc('font', family='serif')
rc('font', size=22)
from matplotlib.font_manager import fontManager, FontProperties
font= FontProperties(size='small')

import pylab as pl
from h5py import File
import os
import string
from matplotlib import pyplot 
from matplotlib import pyplot as plt
from matplotlib import pylab
from matplotlib import colors
import numpy as np
import math
from matplotlib.ticker import NullFormatter
nullfmt = NullFormatter()

import pickle

nion_mass = 1.12e46

filename = 'virialmass_escapefraction'

fig = plt.figure(figsize=[12,12])
plt.subplots_adjust(left=0.07, right=0.95, bottom=0.08, top=0.98, hspace=1e-3)

fesc= [[[],[],[]],[[],[],[]]]

for ii,n in zip([0,1],[13,41]):

 f=open('RD%04i_escape_fraction.txt'%n,'r')
 lines=f.readlines()
 f.close()

 Mass = []
 #Mass1 = []
 #fesc = [[],[],[]]
 #fesc1 = []
 fstar = []
 fgas = []
 for i in xrange(1,len(lines)):
    ftemp = float(lines[i].split()[7])
    if ftemp  < 0.0: continue
    m = float(lines[i].split()[2])
    if m > 10**7.5 and m <= 10**8:
      fesc[ii][0].append(ftemp)
    if m > 10**8.0 and m <= 10**8.5:
      fesc[ii][1].append(ftemp)
    if m > 10**8.5 and m <= 10**9.0:
      fesc[ii][2].append(ftemp)
    #if float(lines[i].split()[7]) <= 0.0: continue
    #Mass1.append(float(lines[i].split()[2]))
    #fesc1.append(float(lines[i].split()[7]))
    #fstar.append(float(lines[i].split()[5])/float(lines[i].split()[3]))    
    #fgas.append(float(lines[i].split()[3])/float(lines[i].split()[2]))

fo=open("escapefraction.dat",'w')
pickle.dump(fesc, fo)
fo.close()


xlabel = [r'$f_{\rm esc}$']
ylabel = [r'$f$']

for iplot in xrange(3):

    ax = fig.add_subplot(3,1,iplot+1)

    # the histogram of the data
    if fesc[0][iplot]!=[]: n, bins, patches = ax.hist(fesc[0][iplot], 50,range=[0,1], normed=1,cumulative=1, histtype='step',color='green',label="z=18.5", alpha=0.75)
    if fesc[1][iplot]!=[]: n, bins, patches = ax.hist(fesc[1][iplot], 50,range=[0,1], normed=1,cumulative=1, histtype='step',color='red',label="z=15", alpha=0.75)   
    #if fesc[2][iplot]!=[]: n, bins, patches = ax.hist(fesc[2][iplot], 50,range=[0,1], normed=1,cumulative=1, histtype='step',color='blue', alpha=0.75)

    if iplot == 0:
       ax.text(0.7,0.2,r'7.5 $<$ log$_{10}$ (M$_{vir}$) $<$ 8.0')
       ax.legend(loc=8)
    elif iplot == 1:
       ax.text(0.7,0.2,r'8.0 $<$ log$_{10}$ (M$_{vir}$) $<$ 8.5')
    elif iplot == 2:
       ax.text(0.7,0.2,r'8.5 $<$ log$_{10}$ (M$_{vir}$) $<$ 9.0')

    # add a 'best fit' line
    #y = mlab.normpdf( bins, mu, sigma)
    #l = plt.plot(bins, y, 'r--', linewidth=1)

    if iplot == 2:pl.xlabel(r'$f_{\rm esc}$')
    pl.ylabel(r'$f$')
    if iplot < 2:ax.xaxis.set_major_formatter(nullfmt)
    pl.ylim(0.01,1.05)
    #plt.title(r'$\mathrm{Histogram\ of\ IQ:}\ \mu=100,\ \sigma=15$')
    #plt.axis([40, 160, 0, 0.03])
    #plt.grid(True)
    #plt.show()


pl.savefig('%s.eps'% (filename),format='eps')
fig.clf()