import numpy as np
from scipy.stats import rv_continuous
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit


num_events = 10000

"""
right ascension (alpha) -- uniform
goes from 0 to 2*pi
"""
'''

y_range_ra = 2*np.pi
grad_ra = 100

xbins = np.linspace(0,y_range_ra,grad_ra*(y_range_ra))

#random distribution
x=np.array(np.random.uniform(0,2*np.pi,size=num_events))

#actual distribution
y_act= (1.0/(2*np.pi))

#FIGURE 1 - LINEAR AXES
fig1 = plt.figure(1)
ax1 = fig1.add_subplot(111)
plt.hist(x,normed=True,bins=xbins,label="probability")
plt.plot([0,2*np.pi],[y_act,y_act],label="expected function")
fig1.suptitle("Expected Right Ascension Distribution for BBH")
plt.ylabel("p")
plt.xlabel("RA (radians)")

# curve fitting
n, bins = np.histogram(x,normed=True,bins=xbins)
bin_centers = xbins[:-1] + 0.5 * (xbins[1:] - xbins[:-1])

def ra_guess(x,a):
    return a*((x-x)+1)

popt, pcov = curve_fit(ra_guess, bin_centers, n, p0 = [.1])

text = "constant, c="+str(popt[0])+" vs. c="+str(1.0/(2*np.pi))

plt.figure(1)
plt.plot(bin_centers,ra_guess(bin_centers, *popt),label="fit function")
ax1.set_title(text,fontsize=9)
plt.legend(loc='upper left')


plt.show()

'''
"""
declination (delta) -- uniform in cos(delta)
"""
'''

y_range_dec = np.pi
grad_dec= 100

xbins = np.linspace(0,y_range_dec,grad_dec*(y_range_dec))

#random distribution
x=np.array(np.random.uniform(-1,1,size=num_events))
y_rand= np.arccos(x)

#actual distribution
y_act= .5*np.sin(xbins)

#FIGURE 2 - LINEAR AXES
fig2 = plt.figure(2)
ax2 = fig2.add_subplot(111)
plt.hist(y_rand,normed=True,bins=xbins,label="probability")
plt.plot(xbins,y_act,label="expected function")
fig2.suptitle("Expected Declination Angle for BBH")
plt.ylabel("p")
plt.xlabel("delta (radians)")

# curve fitting
n, bins = np.histogram(y_rand,normed=True,bins=xbins)
bin_centers = xbins[:-1] + 0.5 * (xbins[1:] - xbins[:-1])

def dec_guess(x,a,b):
    return a*np.sin(b*x)

popt, pcov_= curve_fit(dec_guess, bin_centers, n, p0 = [.3,.9])

text = "a*sin(bx), a="+str(popt[0])+",b="+str(popt[1])+"vs. a=.5,b=1"

plt.figure(2)
plt.plot(bin_centers,dec_guess(bin_centers, *popt),label="fit function")
ax2.set_title(text,fontsize=9)
plt.legend(loc='upper left')

plt.show()

'''
"""
inclination angle (iota) -- uniform in cos(i)
"""
'''

y_range_iota = np.pi
grad_iota= 100

xbins = np.linspace(0,y_range_iota,grad_iota*(y_range_iota))

#random distribution
x=np.array(np.random.uniform(-1,1,size=num_events))
y_rand= np.arccos(x)

#actual distribution
y_act= .5*np.sin(xbins)

#FIGURE 3 - LINEAR AXES
fig3 = plt.figure(3)
ax3 = fig3.add_subplot(111)
plt.hist(y_rand,normed=True,bins=xbins,label="probability")
plt.plot(xbins,y_act,label="expected function")
fig3.suptitle("Expected Inclintation Angle for BBH")
plt.ylabel("p")
plt.xlabel("iota (radians)")

# curve fitting
n, bins = np.histogram(y_rand,normed=True,bins=xbins)
bin_centers = xbins[:-1] + 0.5 * (xbins[1:] - xbins[:-1])

def iota_guess(x,a,b):
    return a*np.sin(b*x)

popt, pcov_= curve_fit(iota_guess, bin_centers, n, p0 = [.3,.9])

text = "a*sin(bx), a="+str(popt[0])+",b="+str(popt[1])+"vs. a=.5,b=1"

plt.figure(3)
plt.plot(bin_centers,iota_guess(bin_centers, *popt),label="fit function")
ax3.set_title(text,fontsize=9)
plt.legend(loc='upper left')

plt.show()
'''

"""
phase (psi) -- uniform  (ex does the system look like an ellipse? How tilted?)
"""
'''
y_range_phase = 2*np.pi
grad_phase= 100

xbins = np.linspace(0,y_range_phase,grad_phase*(y_range_phase))

#random distribution
x=np.array(np.random.uniform(0,2*np.pi,size=num_events))

#actual distribution
y_act= (1.0/(2*np.pi))

#FIGURE 4 - LINEAR AXES
fig4 = plt.figure(4)
ax4 = fig4.add_subplot(111)
plt.hist(x,normed=True,bins=xbins,label="probability")
plt.plot([0,2*np.pi],[y_act,y_act],label="expected function")
fig4.suptitle("Expected Phase Distribution for BBH")
plt.ylabel("p")
plt.xlabel("psi (radians)")

# curve fitting
n, bins = np.histogram(x,normed=True,bins=xbins)
bin_centers = xbins[:-1] + 0.5 * (xbins[1:] - xbins[:-1])

def phase_guess(x,a):
    return a*((x-x)+1)

popt, pcov = curve_fit(phase_guess, bin_centers, n, p0 = [.1])

text = "constant, c="+str(popt[0])+" vs. c="+str(1.0/(2*np.pi))

plt.figure(4)
plt.plot(bin_centers,phase_guess(bin_centers, *popt),label="fit function")
ax4.set_title(text,fontsize=9)
plt.legend(loc='upper left')


plt.show()
'''


"""
coalesence phase -- uniform
"""
'''
y_range_phase_coal = 2*np.pi
grad_phase_coal= 100

xbins = np.linspace(0,y_range_phase_coal,grad_phase_coal*(y_range_phase_coal))

#random distribution
x=np.array(np.random.uniform(0,2*np.pi,size=num_events))

#actual distribution
y_act= (1.0/(2*np.pi))

#FIGURE 5 - LINEAR AXES
fig5 = plt.figure(5)
ax5 = fig5.add_subplot(111)
plt.hist(x,normed=True,bins=xbins,label="probability")
plt.plot([0,2*np.pi],[y_act,y_act],label="expected function")
fig5.suptitle("Expected Coalesence Phase Distribution for BBH")
plt.ylabel("p")
plt.xlabel("beta (radians)")

# curve fitting
n, bins = np.histogram(x,normed=True,bins=xbins)
bin_centers = xbins[:-1] + 0.5 * (xbins[1:] - xbins[:-1])

def phase_coal_guess(x,a):
    return a*((x-x)+1)

popt, pcov = curve_fit(phase_coal_guess, bin_centers, n, p0 = [.1])

text = "constant, c="+str(popt[0])+" vs. c="+str(1.0/(2*np.pi))

plt.figure(5)
plt.plot(bin_centers,phase_coal_guess(bin_centers, *popt),label="fit function")
ax5.set_title(text,fontsize=9)
plt.legend(loc='upper left')


plt.show()


'''


"""
spin magnitudes -- gaussian, mean = .7, stddev = .1
"""
'''

spin_mean=.7
spin_std_dev=.1
y_range_spin = 1
grad_spin=300

xbins = np.linspace(0,y_range_spin,grad_spin*(y_range_spin))

class Spin_mag_pdf(rv_continuous):
    def _pdf(self, x):
        return (1.0/(2*np.pi*spin_std_dev**2)**.5)*np.exp(-(x-spin_mean)**2/(2*spin_std_dev**2))

#random distribution
spin_mag_pdf = Spin_mag_pdf(name="Spin mag distribution", a=0,b=1)
x_1=np.array(np.random.rand(num_events))
x_2=np.array(np.random.rand(num_events))
y_rand_1=spin_mag_pdf.ppf(x_1)
y_rand_2=spin_mag_pdf.ppf(x_2)

#actual distribution
y_act= (1.0/(2*np.pi*spin_std_dev**2)**.5)*np.exp(-(xbins-spin_mean)**2/(2*spin_std_dev**2))

#FIGURE 6 - LINEAR AXES
fig6 = plt.figure(6)
ax6 = fig6.add_subplot(111)
plt.hist(y_rand_1,normed=True,bins=xbins,label="probability")
plt.plot(xbins,y_act,label="expected function")
fig6.suptitle("Expected Spin Magnitude Distribution for BH_1")
plt.ylabel("p")
plt.xlabel("spin")

#FIGURE 7 - LINEAR AXES
fig7 = plt.figure(7)
ax7 = fig7.add_subplot(111)
plt.hist(y_rand_2,normed=True,bins=xbins,label="probability")
plt.plot(xbins,y_act,label="expected function")
fig7.suptitle("Expected Spin Magnitude Distribution for BH_2")
plt.ylabel("p")
plt.xlabel("spin")

# curve fitting
n_1, bins_1 = np.histogram(y_rand_1,normed=True,bins=xbins)
bin_centers = xbins[:-1] + 0.5 * (xbins[1:] - xbins[:-1])
n_2, bins_2 = np.histogram(y_rand_2,normed=True,bins=xbins)

def spin_mag_guess(x, mean, std_dev,a):
    return a*np.exp(-(x-mean)**2/(2.*std_dev**2))

popt_1, pcov_1 = curve_fit(spin_mag_guess, bin_centers, n_1, p0 = [.5, .01,1])
popt_2, pcov_2 = curve_fit(spin_mag_guess, bin_centers, n_2, p0 = [.5, .01,1])

text_1 = "gaussian, mean="+str(popt_1[0])+",std_dev="+str(popt_1[1])+" vs. mean=.7,std_dev=.1"
text_2 = "gaussian, mean="+str(popt_2[0])+",std_dev="+str(popt_2[1])+" vs. mean=.7,std_dev=.1"

plt.figure(6)
plt.plot(bin_centers,spin_mag_guess(bin_centers, *popt_1),label="fit function")
ax6.set_title(text_1,fontsize=9)
plt.legend(loc='upper left')
plt.figure(7)
plt.plot(bin_centers,spin_mag_guess(bin_centers, *popt_2),label="fit function")
ax7.set_title(text_2,fontsize=9)
plt.legend(loc='upper left')

plt.show()
'''

"""
spin azimuthal angles -- uniform
"""

'''
y_range_spin_azi = 2*np.pi
grad_spin_azi= 100

xbins = np.linspace(0,y_range_spin_azi,grad_spin_azi*(y_range_spin_azi))

#random distribution
x_1=np.array(np.random.uniform(0,2*np.pi,size=num_events))
x_2=np.array(np.random.uniform(0,2*np.pi,size=num_events))
#actual distribution
y_act= (1.0/(2*np.pi))

#FIGURE 8 - LINEAR AXES
fig8 = plt.figure(8)
ax8 = fig8.add_subplot(111)
plt.hist(x_1,normed=True,bins=xbins,label="probability")
plt.plot([0,2*np.pi],[y_act,y_act],label="expected function")
fig8.suptitle("Expected Spin Azimuthal Angle for BH_1")
plt.ylabel("p")
plt.xlabel("gamma (radians)")

#FIGURE 9 - LINEAR AXES
fig9 = plt.figure(9)
ax9 = fig9.add_subplot(111)
plt.hist(x_2,normed=True,bins=xbins,label="probability")
plt.plot([0,2*np.pi],[y_act,y_act],label="expected function")
fig9.suptitle("Expected Spin Azimuthal Angle for BH_2")
plt.ylabel("p")
plt.xlabel("gamma (radians)")


# curve fitting
n_1, bins = np.histogram(x_1,normed=True,bins=xbins)
n_2, bins = np.histogram(x_2,normed=True,bins=xbins)
bin_centers = xbins[:-1] + 0.5 * (xbins[1:] - xbins[:-1])

def spin_azi_guess(x,a):
    return a*((x-x)+1)

popt_1, pcov_1 = curve_fit(spin_azi_guess, bin_centers, n_1, p0 = [.1])
popt_2, pcov_2 = curve_fit(spin_azi_guess, bin_centers, n_2, p0 = [.1])

text_1 = "constant, c="+str(popt_1[0])+" vs. c="+str(1.0/(2*np.pi))
text_2 = "constant, c="+str(popt_2[0])+" vs. c="+str(1.0/(2*np.pi))

plt.figure(8)
plt.plot(bin_centers,spin_azi_guess(bin_centers, *popt_1),label="fit function")
ax8.set_title(text_1,fontsize=9)
plt.legend(loc='upper left')
plt.figure(9)
plt.plot(bin_centers,spin_azi_guess(bin_centers, *popt_2),label="fit function")
ax9.set_title(text_2,fontsize=9)
plt.legend(loc='upper left')


plt.show()
'''
"""
spin polar angles -- uniform in cos(theta)
"""
'''
y_range_spin_pol = np.pi
grad_spin_pol= 100

xbins = np.linspace(0,y_range_spin_pol,grad_spin_pol*(y_range_spin_pol))

#random distribution
x_1=np.array(np.random.uniform(-1,1,size=num_events))
x_2=np.array(np.random.uniform(-1,1,size=num_events))
y_rand_1= np.arccos(x_1)
y_rand_2=np.arccos(x_2)

#actual distribution
y_act= .5*np.sin(xbins)

#FIGURE 10 - LINEAR AXES
fig10= plt.figure(10)
ax10 = fig10.add_subplot(111)
plt.hist(y_rand_1,normed=True,bins=xbins,label="probability")
plt.plot(xbins,y_act,label="expected function")
fig10.suptitle("Expected Spin Polar Angle for BH_1")
plt.ylabel("p")
plt.xlabel("kappa (radians)")

#FIGURE 11 - LINEAR AXES
fig11 = plt.figure(11)
ax11 = fig11.add_subplot(111)
plt.hist(y_rand_2,normed=True,bins=xbins,label="probability")
plt.plot(xbins,y_act,label="expected function")
fig11.suptitle("Expected Spin Polar Angle for BH_2")
plt.ylabel("p")
plt.xlabel("kappa (radians)")


# curve fitting
n_1, bins = np.histogram(y_rand_1,normed=True,bins=xbins)
n_2, bins = np.histogram(y_rand_2,normed=True,bins=xbins)
bin_centers = xbins[:-1] + 0.5 * (xbins[1:] - xbins[:-1])

def spin_pol_guess(x,a,b):
    return a*np.sin(b*x)

popt_1, pcov_1 = curve_fit(spin_pol_guess, bin_centers, n_1, p0 = [.3,.9])
popt_2, pcov_2 = curve_fit(spin_pol_guess, bin_centers, n_2, p0 = [.3,.9])

text_1 = "a*sin(bx), a="+str(popt_1[0])+",b="+str(popt_1[1])+"vs. a=.5,b=1"
text_2 = "a*sin(bx), a="+str(popt_2[0])+",b="+str(popt_2[1])+"vs. a=.5,b=1"

plt.figure(10)
plt.plot(bin_centers,spin_pol_guess(bin_centers, *popt_1),label="fit function")
ax10.set_title(text_1,fontsize=9)
plt.legend(loc='upper left')
plt.figure(11)
plt.plot(bin_centers,spin_pol_guess(bin_centers, *popt_2),label="fit function")
ax11.set_title(text_2,fontsize=9)
plt.legend(loc='upper left')

plt.show()
'''
"""
coalesence time -- uniform
70% +/ 5% chance of being on
assume you're on / off in hour chunks
"""
num_events=10000
t_obs = 1
grad_t = 200 #minutes; should be 3600...????
percent_on = [.7,.6,.8]
error = [.05,.1,.05]
num_detectors=len(percent_on)

xbins = np.linspace(0,t_obs,grad_t*(t_obs))
event_time = np.array(np.random.uniform(0,1,size=num_events))

detected={}
for i in range(num_detectors):
    event = np.array(np.random.uniform(0,1,size=num_events))
    y_detected=event<np.random.normal(percent_on[i], error[i], size=num_events)
    detected[i+1]=y_detected*1
    detected[i+1].astype(float)

num_detectors_per_event=np.zeros((1,num_events))
for i in range(num_detectors):
    num_detectors_per_event+=detected[i+1]
num_detectors_per_event=num_detectors_per_event[0]

xbins = range(0,num_detectors+1)

y = []
for i in range(num_detectors+1):
    y.append(list(num_detectors_per_event).count(i))

fig12 = plt.figure(12)
ax12 = fig12.add_subplot(111)
plt.bar(xbins,y)
fig12.suptitle("Number of Detectors Detecting Each Event")
plt.ylabel("# events")
plt.xlabel("# detectors")
#text = "num events="+str(num_events)+"; num events detected="+str(num_detected)
#ax12.set_title(text,fontsize=9)
plt.legend(loc='upper left')

num_bars = 20.

x=np.linspace(0.03,.97,20)
fig13 = plt.figure(13)


for i in range(num_detectors):
    plt.subplot(num_detectors, 1, i+1)
    plt.tight_layout()
    title="Detector "+str(i+1)+"; Time on: mean="+str(percent_on[i])+",error="+str(error[i])
    plt.title(title)

    y = np.mean(detected[i+1].reshape(-1, int(num_events/num_bars)), axis=1)
    plt.bar(x,y,width=.04,color="g",label="prob. detected")
    plt.bar(x,1-y,width=.04,color="r",bottom=y,label="prob. undetected")
    plt.legend(loc='upper left')








#detected_both = (detected_1*1)*(detected_2*1)
'''
detected_events = detected_both*event
detected_events = list(filter(lambda a: a != 0, detected_events))
not_detected_events = (1-detected_both)*event
not_detected_events = list(filter(lambda a: a != 0, not_detected_events))
'''
#num_detected = sum(detected_both)

'''
fig12 = plt.figure(12)
ax12 = fig12.add_subplot(111)
plt.hist(detected_events, bins = xbins,label="detected events",color="g",)
plt.hist(not_detected_events,bins=xbins,label="undetected events",color="r",)
fig12.suptitle("Expected Time Distribution for BBH")
plt.ylabel("counts")
plt.xlabel("t (in total runtime)")
text = "num events="+str(num_events)+"; num events detected="+str(num_detected)
ax12.set_title(text,fontsize=9)
plt.legend(loc='upper left')
'''

plt.show()







"""
masses -- Salpeter mass distribution
"""
'''

min_bh_mass = .5 #default: .5, realistically ~2.5
y_range_mass = 100
grad_mass=3
xbins_mass = np.linspace(1,y_range_mass,grad_mass*(y_range_mass))

y_range_eta = .25
grad_eta=100
xbins_eta = np.linspace(0,y_range_eta,grad_eta*(y_range_eta))

class Mass_pdf(rv_continuous):
    def _pdf(self, x):
        return 1.35*((2*min_bh_mass)**(1.35))*(x**(-2.35))
#random distribution
mass_pdf = Mass_pdf(name="Mass distribution", a=2*min_bh_mass)
x=np.array(np.random.rand(num_events))
y_rand_tot=mass_pdf.ppf(x) #y_rand_tot is the list of random masses

# NOW GENERATE ETA
eta_std_dev=.05
eta_mean=.25

class Eta_pdf(rv_continuous):
    def _pdf(self, x):
        return (2.0/(2*np.pi*eta_std_dev**2)**.5)*np.exp(-(x-eta_mean)**2/(2*eta_std_dev**2))
#random distribution
eta_pdf = Eta_pdf(name="Eta distribution", a=0,b=.25)
x=np.array(np.random.rand(num_events))
y_rand_eta=eta_pdf.ppf(x) #y_rand_eta is the list of random etas

# CALCULATE M_1 AND M_2 AND PLOT

m_1 = .5*y_rand_tot*(1+(1-4*y_rand_eta))
m_2 = .5*y_rand_tot*(1-(1-4*y_rand_eta))

y_act_mass= 1.35*((2*min_bh_mass)**(1.35))*(xbins_mass**(-2.35))

y_range = 100000
xbins = np.linspace(1,y_range,grad_mass*(y_range))

#FIGURE 13 - LINEAR AXES
fig13 = plt.figure(13)
plt.scatter(m_1,m_2,color="b",label="m_1 vs m_2",s=2)
plt.plot(xbins,xbins,color="g",label="m_1=m_2")
fig13.suptitle("m_1 vs m_2 for BH's")
plt.ylabel("m_2 (M_solar)")
plt.xlabel("m_1 (M_solar)")
plt.yscale('log', nonposy='clip')
plt.xlabel("mass (M_solar)")
plt.legend(loc='upper left')

#FIGURE 14 - LINEAR AXES
fig14 = plt.figure(14)
ax14 = fig14.add_subplot(111)
plt.hist(m_2,normed=True,bins=xbins_mass,label="probability")
plt.plot(xbins_mass,y_act_mass,label="expected function")
fig14.suptitle("Expected Mass for BH_2")
plt.ylabel("p")
plt.xlabel("mass (M_solar)")
plt.legend(loc='upper left')

#FIGURE 15 - LINEAR AXES
fig15 = plt.figure(15)
ax15 = fig15.add_subplot(111)
plt.hist(m_1,normed=True,bins=xbins_mass,label="probability")
plt.plot(xbins_mass,y_act_mass,label="expected function")
fig15.suptitle("Expected Mass for BH_1")
plt.ylabel("p")
plt.xlabel("mass (M_solar)")
plt.legend(loc='upper left')


logbins_mass = 10 ** np.linspace(np.log10(1), np.log10(y_range_mass), grad_mass*y_range_mass)

#FIGURE 16 - LOG AXES
fig16 = plt.figure(16)
ax16 = fig16.add_subplot(111)
plt.hist(m_2,normed=True,bins=logbins_mass,label="probability")
plt.plot(xbins_mass,y_act_mass,label="expected function")
fig16.suptitle("Expected Mass for BH_2")
plt.ylabel("p")
plt.gca().set_xscale("log")
plt.yscale('log', nonposy='clip')
plt.xlabel("mass (M_solar)")
plt.legend(loc='upper left')

#FIGURE 17 - LOG AXES
fig17 = plt.figure(17)
ax17 = fig17.add_subplot(111)
plt.hist(m_1,normed=True,bins=logbins_mass,label="probability")
plt.plot(xbins_mass,y_act_mass,label="expected function")
fig17.suptitle("Expected Mass for BH_1")
plt.ylabel("p")
plt.gca().set_xscale("log")
plt.yscale('log', nonposy='clip')
plt.xlabel("mass (M_solar)")
plt.legend(loc='upper left')

# curve fitting
n_2, bins_2 = np.histogram(m_2,normed=True,bins=xbins_mass)
n_1, bins_1 = np.histogram(m_1,normed=True,bins=xbins_mass)
bin_centers = xbins_mass[:-1] + 0.5 * (xbins_mass[1:] - xbins_mass[:-1])

def mass_guess(x, a, b):
    return a*x**b

popt_2, pcov_2 = curve_fit(mass_guess, bin_centers, n_2, p0 = [1, -3])
popt_1, pcov_1 = curve_fit(mass_guess, bin_centers, n_1, p0 = [1, -3])

text_2 = "a*x^b a="+str(popt_2[0])+",b="+str(popt_2[1])+" vs. a=1.35,b=-2.35"
text_1 = "a*x^b a="+str(popt_1[0])+",b="+str(popt_1[1])+" vs. a=1.35,b=-2.35"

plt.figure(14)
plt.plot(bin_centers,mass_guess(bin_centers, *popt_2),label="fit function")
ax14.set_title(text_2,fontsize=9)
plt.legend(loc='upper right')
plt.figure(16)
plt.plot(bin_centers,mass_guess(bin_centers, *popt_2),label="fit function")
ax16.set_title(text_2,fontsize=9)
plt.legend(loc='upper right')
plt.figure(15)
plt.plot(bin_centers,mass_guess(bin_centers, *popt_1),label="fit function")
ax15.set_title(text_1,fontsize=9)
plt.legend(loc='upper right')
plt.figure(17)
plt.plot(bin_centers,mass_guess(bin_centers, *popt_1),label="fit function")
ax17.set_title(text_1,fontsize=9)
plt.legend(loc='upper right')

plt.show()
'''



"""
luminosity distance (d_L) -- quadratic
"""

'''
grad_dist=1000
density = 50000000
horizon_dist = 1

class Distance_pdf(rv_continuous):
    def _pdf(self, x):
        return 3.0*(x**2)*(1.0/horizon_dist)**3

xbins = np.linspace(0,horizon_dist,grad_dist*(horizon_dist))

#random distribution
distance_pdf = Distance_pdf(name="Distance distribution", a=0,b=horizon_dist)
x=np.array(np.random.rand(num_events))
y_rand=distance_pdf.ppf(x)

#actual distribution
y_act = 3.0*(1.0/horizon_dist)**3*(xbins**2)

#FIGURE 18 - LINEAR X AXIS
fig18 = plt.figure(18)
ax18 = fig18.add_subplot(111)
plt.hist(y_rand,normed=True,bins=xbins,label="probability")
plt.plot(xbins,y_act,label="expected function")
fig18.suptitle("Expected Radial Distribution for BH's,normed")
plt.ylabel("p")
plt.xlabel("r (Gpc)")

#FIGURE 19 - LOG ALL AXES
fig19 = plt.figure(19)
ax19 = fig19.add_subplot(111)
logbins = 10 ** np.linspace(np.log10(1), np.log10(horizon_dist), grad_dist*horizon_dist)
plt.hist(y_rand,normed=True,bins=xbins,label="probability")
plt.plot(xbins,y_act,label="expected function")
fig19.suptitle("Expected Radial Distribution for BH's,normed")
plt.ylabel("p")
plt.gca().set_xscale("log")
plt.yscale('log', nonposy='clip')
plt.xlabel("r (Gpc)")

# curve fitting
n, bins = np.histogram(y_rand,normed=True,bins=xbins)
bin_centers = bins[:-1] + 0.5 * (bins[1:] - bins[:-1])
def f(x, a, b):
    return a*x**b
popt, pcov = curve_fit(f, bin_centers, n, p0 = [1,1])

text = "a*x^b a="+str(popt[0])+",b="+str(popt[1])+" vs. a="+str(3.0*(1.0/horizon_dist)**3)+",b=2"

plt.figure(18)
plt.plot(bin_centers,f(bin_centers, *popt),label="fit function")
ax18.set_title(text,fontsize=9)
plt.legend(loc='upper right')
plt.figure(19)
plt.plot(bin_centers,f(bin_centers, *popt),label="fit function")
ax19.set_title(text,fontsize=9)
plt.legend(loc='upper right')

# NOW UNNORMED


#random distribution
normed_value = (4.0/3)*np.pi*density*horizon_dist**3
hist, bins = np.histogram(y_rand, bins=xbins, density=True)
widths = np.diff(bins)
hist *= normed_value

#FIGURE 20 - LINEAR X AXIS
fig20 = plt.figure(20)
ax20 = fig20.add_subplot(111)
plt.bar(bins[:-1], hist, widths,label="probability",color="b")
plt.plot(xbins,y_act,label="expected function")
fig20.suptitle("Expected Radial Distribution for BH's")
plt.ylabel("p")
plt.xlabel("r (Gpc)")

#FIGURE 21 - LOG ALL AXES
fig21 = plt.figure(21)
ax21 = fig21.add_subplot(111)
plt.bar(logbins[1:], hist, widths,label="probability",color="b")
plt.plot(xbins,y_act,label="expected function")
fig21.suptitle("Expected Radial Distribution for BH's")
plt.ylabel("p")
plt.gca().set_xscale("log")
plt.yscale('log', nonposy='clip')
plt.xlabel("r (Gpc)")

# curve fitting
bin_centers = bins[:-1] + 0.5 * (bins[1:] - bins[:-1])
def f(x, a, b):
    return a*x**b
popt, pcov = curve_fit(f, bin_centers, hist, p0 = [50000000,1])

text = "a*x^b a="+str(popt[0])+",b="+str(popt[1])+" vs. a="+str((4.0)*np.pi*density)+",b=2"

plt.figure(20)
plt.plot(bin_centers,f(bin_centers, *popt),label="fit function")
ax20.set_title(text,fontsize=9)
plt.legend(loc='upper right')
plt.figure(21)
plt.plot(bin_centers,f(bin_centers, *popt),label="fit function")
ax21.set_title(text,fontsize=9)
plt.legend(loc='upper right')

plt.show()

'''