notebook.community

Edit and run 





In [56]:



    


import numpy as np
import matplotlib.pyplot as plt
import os
import sys
import readsubfHDF5
import emcee



    











In [57]:



    


# Read catalog
basedir = '../../../../../../../media/Janus/Illustris-3'
snapid = 135
catalog = readsubfHDF5.subfind_catalog(basedir,snapid)



    











In [58]:



    


# Set BH Masses and Accretion Rate
bm = catalog.GroupBHMass
#bm *= 1.e7
bmdot=catalog.GroupBHMdot
#bmdot *= 1.e
print len(bm)
print len(bmdot)



    


















    






131727
131727

















In [59]:



    


# Plot m vs mdot with metallicity as color

bm = catalog.GroupBHMass
#bm *= 1.e7
bmdot=catalog.GroupBHMdot
#bmdot *= 1.e7

fig = plt.figure(figsize = (10,10))

sc = plt.scatter (np.log10(bm), np.log10(bmdot), c = catalog.GroupStarMetallicity, s = 8, cmap = 'jet')
plt.colorbar(sc)
plt.show()



    


















    
























In [60]:



    


# Plot m vs mdot colored by sfr

bm = catalog.GroupBHMass
#bm *= 1.e7
bmdot = catalog.GroupBHMdot
#bmdot *= 1.e7

sfr = catalog.GroupSFR
fig = plt.figure(figsize = (10,10))
#sc = plt.scatter (np.log10(bm), np.log10(bmdot), s = 8, alpha = 0.5)

sc = plt.scatter (np.log10(bm), np.log10(bmdot), c = np.log10(sfr), cmap = 'jet', s = 8, alpha = 0.5)
plt.colorbar(sc)
plt.show()



    


















    
























In [61]:



    


bm = catalog.GroupBHMass
#bm *= 1.e7
bmdot = catalog.GroupBHMdot
#bmdot *= 1.e7

print 'Number of BHs:', len(bm)

plt.hist(bm, bins=30, log = True)
plt.xlabel('log BH mass [$10^{10} M_{\odot}$]')
plt.show()



    


















    






Number of BHs: 131727










    
























In [62]:



    


bm = np.asarray(catalog.GroupBHMass)
#bm *= 1.e7
bmdot = np.asarray(catalog.GroupBHMdot)
#bmdot *= 1.e7

print np.sum(bm==0), np.sum(bm!=0), len(bm)
print np.sum(bmdot==0), np.sum(bmdot!=0), len(bmdot)

cond = (bm != 0) * (bmdot != 0)
bm = bm[cond]
bmdot = bmdot[cond]

plt.hist(bm, log = True)
plt.show()

plt.hist(bmdot, log = True)



    


















    






107215 24512 131727
107215 24512 131727










    

















    
Out[62]:








(array([  2.44960000e+04,   1.00000000e+01,   3.00000000e+00,
         0.00000000e+00,   1.00000000e+00,   0.00000000e+00,
         0.00000000e+00,   0.00000000e+00,   1.00000000e+00,
         1.00000000e+00]),
 array([  9.47833023e-24,   1.99749351e-02,   3.99498701e-02,
         5.99248052e-02,   7.98997402e-02,   9.98746753e-02,
         1.19849610e-01,   1.39824545e-01,   1.59799480e-01,
         1.79774415e-01,   1.99749351e-01]),
 <a list of 10 Patch objects>)










    
























In [63]:



    


fig = plt.figure(figsize = (8,4), dpi = 400)

plt.loglog(bm, bmdot, '.', alpha = 0.1)
plt.ylabel('$\.{M}_{BH}$', size = 14)
plt.xlabel('$M_{BH}$', size = 14)
plt.title('Illustris-2 Black hole population', size = 18)



    


















    
Out[63]:








<matplotlib.text.Text at 0x7f552e61a050>










    
























In [64]:



    


bm = np.asarray(catalog.GroupBHMass)
bmdot = np.asarray(catalog.GroupBHMdot)

np.savetxt('bm.txt', bm)
np.savetxt('bmdot.txt', bmdot)



    











In [65]:



    


from matplotlib.colors import LogNorm

bm = np.asarray(catalog.GroupBHMass)
bmdot = np.asarray(catalog.GroupBHMdot)

cond = (bm != 0) * (bmdot != 0)
bm = bm[cond]
bmdot = bmdot[cond]

'''
now transform m and mdot from simulation units to real units
(Msol and Msol/s)
'''

bm_phys = bm * 4.3e10
bmdot_phys = bmdot * 9.72e-7

from sklearn import mixture
import matplotlib.pyplot
import matplotlib.mlab
import numpy as np
clf = mixture.GMM(n_components=2, covariance_type='full')
clf.fit(np.log10(bmdot_phys))
m1, m2 = clf.means_
print m1, m2
w1, w2 = clf.weights_
c1, c2 = clf.covars_
histdist = matplotlib.pyplot.hist(np.log10(bmdot_phys), 100, normed=True)
plotgauss1 = lambda x: plot(x,w1*matplotlib.mlab.normpdf(x,m1,np.sqrt(c1))[0], linewidth=2)
plotgauss2 = lambda x: plot(x,w2*matplotlib.mlab.normpdf(x,m2,np.sqrt(c2))[0], linewidth=2)
plotgauss1(histdist[1])
plotgauss2(histdist[1])

plt.xlabel('Accretion Rate [$log_{10}(M_{\odot}/s)$]', size = 14)
plt.ylabel('fraction of sample', size = 14)
plt.title('$\.{M}_{BH}$ distribution', size = 18)

plt.show()

fig = plt.figure(figsize = (8,4), dpi = 400)

plt.loglog(bm_phys, bmdot_phys, '.', alpha = 0.1)
plt.ylabel('$\.{M}_{BH}$[$log_{10}(M_{\odot}/s)$]', size = 14)
plt.xlabel('$M_{BH}$[$log_{10}(M_{\odot})$]', size = 14)
plt.title('Illustris-2 Black hole population', size = 18)
plt.show()

fig = plt.figure(figsize = (12,8), dpi = 400)

bmdot_phys = bmdot_phys[ np.log10(bmdot_phys) > m2]
bm_phys = bm_phys[ np.log10(bmdot_phys) > m2 ]

p, residual, rank, singular_values, rcond = np.polyfit( np.log10(bm_phys), np.log10(bmdot_phys), 1, full = True )
plt.plot(np.linspace(5., 11., 100), np.polyval(p, np.linspace(5., 11., 100)), c = 'g', linestyle = '--')

relation_params = p

H, xedges, yedges = np.histogram2d(np.log10(bm_phys), np.log10(bmdot_phys), 30)
plt.imshow(H.T, origin = 'lower', extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]], cmap = 'Oranges', norm=LogNorm(), aspect = 0.667, interpolation = 'None')
plt.colorbar(shrink = 0.8)
plt.scatter(np.log10(bm_phys[::50]), np.log10(bmdot_phys[::50]), marker = '.', c = 'grey')

plt.text(6.5,-7.5, '$log_{10}(\.{M})$ = ' + str(np.round(p[0], 3)) + ' $log_{10}(M)$ ' + str(np.sign(p[1]))[0] + str(np.abs(np.round(p[1], 3))), size = 16, color = 'k' )

plt.ylabel('Accretion Rate [$log_{10}(M_{\odot}/s)$]', size = 16)
plt.xlabel('BH Mass [$log_{10}(M_{\odot})$]', size = 16)
plt.title('Illustris-2 Black hole population',  size = 16)
plt.show()



    


















    






[-12.53454338] [-15.68895702]










    

















    

















    
























In [66]:



    


residuals = (np.log10(bmdot_phys)) - (np.polyval(p, np.log10(bm_phys)))
s = np.sqrt( np.sum(residuals**2.) / (len(residuals) - 2) )
se = s / np.sqrt( np.sum( (np.log10(bm_phys) - np.mean(np.log10(bm_phys)))**2. ) )
print se



    


















    






0.0067397834713

Lx proportional to M





In [20]:



    


import scipy.stats

LHS_case1 = np.log10(623.04*bm_phys + 1.656e-15) - np.log10(bm_phys)
RHS_case1 = np.log10(bmdot_phys)

slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(RHS_case1, LHS_case1)
print slope, intercept, r_value**2.

fig = plt.figure(figsize = (12,8), dpi = 400)

#H, xedges, yedges = np.histogram2d(RHS_case1, LHS_case1, (40, 20))
#plt.imshow(H.T, origin = 'lower', extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]], cmap = 'Blues', norm=LogNorm(), interpolation = 'None', zorder = 0)
plt.plot(np.linspace(-7., 5., 100), np.polyval([slope, intercept], np.linspace(-7., 5., 100)), c = 'r', linestyle = '--', zorder = 1)
#plt.scatter(RHS_case1, LHS_case1)

#plt.scatter(RHS_case1, LHS_case1, marker = '.', s = 1, zorder = 0)
plt.ylabel('$log_{10}(L_X) - \,log_{10}(M)$', size = 16)
plt.xlabel('$log_{10}(\.{M})$', size = 16)
#plt.colorbar(shrink=.6)
plt.text(-6.5, 2.9, '$log_{10}(L_X) - log_{10}(M)$ = ' + str(np.round(slope, 3)) + ' $(log_{10}(\.{M}))$ + ' + str(np.abs(np.round(intercept, 3))), size = 16, color = 'k' )
plt.title('$L_X \propto M$', size = 18)
plt.show()



    


















    






-1.57353789987e-18 2.79451592984 4.72554241666e-11










    
























In [21]:



    


plt.ylabel('$log_{10}(L_X) - \,log_{10}(M)$', size = 16)
plt.xlabel('$log_{10}(\.{M})$', size = 16)
#plt.colorbar(shrink=.6)
#plt.text(-6.5, .005, '$log_{10}(L_X) - log_{10}(M)$ = ' + str(np.round(slope, 3)) + ' $(log_{10}(\.{M}))$ + ' + str(np.abs(np.round(intercept, 3))), size = 16, color = 'k' )
plt.title('$L_X \propto M$', size = 18)
plt.scatter(RHS_case1, LHS_case1)



    


















    
Out[21]:








<matplotlib.collections.PathCollection at 0x7f5541049490>










    
























In [22]:



    


import scipy.stats

LHS_case1 = np.log10(623.04*bm_phys + 1.656e-15)
RHS_case1 = np.log10(bmdot_phys) + np.log10(bm_phys)

slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(RHS_case1, LHS_case1)
print slope, intercept, r_value**2.

fig = plt.figure(figsize = (12,8), dpi = 400)

H, xedges, yedges = np.histogram2d(RHS_case1, LHS_case1, (40, 20))
plt.imshow(H.T, origin = 'lower', extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]], cmap = 'Blues', norm=LogNorm(), interpolation = 'None', zorder = 0)
plt.plot(np.linspace(-7., 5., 100), np.polyval([slope, intercept], np.linspace(-7., 5., 100)), c = 'r', linestyle = '--', zorder = 1)

#plt.scatter(RHS_case1, LHS_case1, marker = '.', s = 1, zorder = 0)
plt.ylabel('$log_{10}(L_X) = log_{10}(623.04 \, M + 1.656 \\times 10^{-15})$', size = 16)
plt.xlabel('$log_{10}(\.{M}) + \,log_{10}(M)$', size = 16)
plt.colorbar(shrink=.6)
plt.text(-6.5, 14., '$log_{10}(L_X)$ = ' + str(np.round(slope, 3)) + ' ($log_{10}(\.{M}) + log_{10}(M))$ + ' + str(np.abs(np.round(intercept, 3))), size = 16, color = 'k' )
plt.title('$L_X \propto M$', size = 18)



    


















    






0.410738676998 11.5640719191 0.767644746866










    
Out[22]:








<matplotlib.text.Text at 0x7f5540f3bcd0>

Lx proportional to Mdot





In [15]:



    


import scipy.stats

LHS_case2 = np.log10(4.64e19 * bmdot_phys) - np.log10(bm_phys) # y axis
RHS_case2 = np.log10(bmdot_phys) # x axis

slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(RHS_case2, LHS_case2)
print slope, intercept, r_value**2.

fig = plt.figure(figsize = (12,8), dpi = 400)

H, xedges, yedges = np.histogram2d(RHS_case2, LHS_case2, (40, 20))
plt.imshow(H.T, origin = 'lower', extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]], cmap = 'Blues', norm=LogNorm(), interpolation = 'None', zorder = 0)
plt.plot(np.linspace(-14., -7., 100), np.polyval([slope, intercept], np.linspace(-14., -7., 100)), c = 'r', linestyle = '--', zorder = 1)

#plt.scatter(RHS_case1, LHS_case1, marker = '.', s = 1, zorder = 0)
plt.ylabel('$log_{10}(L_X) - log_{10}(M)$', size = 16)
plt.xlabel('$log_{10}(\.{M})$', size = 16)
plt.colorbar(shrink=.6)
plt.text(-13., 4., '$log_{10}(L_X) - log_{10}(M)$ = ' + str(np.round(slope, 3)) + ' $log_{10}(\.{M})$ + ' + str(np.abs(np.round(intercept, 3))), size = 16, color = 'k' )
plt.title('$L_X \propto \.{M}$', size = 18)



    


















    






0.520536469162 7.50099453526 0.457039931768










    
Out[15]:








<matplotlib.text.Text at 0x3ea9710>










    
























In [16]:



    


import scipy.stats

LHS_case2 = np.log10(4.64e19 * bmdot_phys)
RHS_case2 = np.log10(bmdot_phys) + np.log10(bm_phys)

slope, intercept, r_value, p_value, std_err = scipy.stats.linregress(RHS_case2, LHS_case2)
print slope, intercept, r_value**2.

fig = plt.figure(figsize = (12,8), dpi = 400)

H, xedges, yedges = np.histogram2d(RHS_case2, LHS_case2, (40, 20))
plt.imshow(H.T, origin = 'lower', extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]], cmap = 'Blues', norm=LogNorm(), interpolation = 'None', zorder = 0)
plt.plot(np.linspace(-7., 5., 100), np.polyval([slope, intercept], np.linspace(-7., 5., 100)), c = 'r', linestyle = '--', zorder = 1)

#plt.scatter(RHS_case1, LHS_case1, marker = '.', s = 1, zorder = 0)
plt.ylabel('$log_{10}(L_X) = log_{10}(4.64 \\times 10^{19} \, \.{M})$', size = 16)
plt.xlabel('$log_{10}(\.{M}) + \,log_{10}(M)$', size = 16)
plt.colorbar(shrink=.6)
plt.text(-6.5, 12., '$log_{10}(L_X)$ = ' + str(np.round(slope, 3)) + ' ($log_{10}(\.{M}) + log_{10}(M))$ + ' + str(np.abs(np.round(intercept, 3))), size = 16, color = 'k' )
plt.title('$L_X \propto \.{M}$', size = 18)



    


















    






0.589261323002 10.8969619913 0.871790637515










    
Out[16]:








<matplotlib.text.Text at 0x7f5541919ad0>

finding q and K using convergence

General idea: For Lx ~ M Define function f(M,q,k) = log10(a + b/M) - 0.694 q log10(M) + 16.3769 - k = 0 a = 623.04 b = 1.656e-15 find zeros For Lx ~ \.{M} Define function f(\.{M},q,k) = log10(a + b/\.{M}) - 0.694 q log10(\.{M}(M)} + 16.3769 - k = 0 a = 9.03e18 b = 1.656e-15 \.{M}(M) is given by the orange plot solve system of equations: log10(a1 + b/M) - d*q*log10(M) + e*q - log10(a2*Mdot + b) + 1/d * log10(Mdot) - e/d * q * log10(Mdot) = 0 




In [77]:



    


import scipy.optimize as spopt

def f(q, m, mdot):
    a1 = 623.04
    b = 1.656e-15
    d = relation_params[0]
    a2 = 9.03e18
    e = -relation_params[1]
    return np.log10(a1 + b/m) - d*q*np.log(m) + e*q - np.log10(a2*mdot + b) + (1. - e*q)/d * np.log10(mdot)

sols = []
#sol = spopt.root(f, [1.], args = (bm_phys[0], bmdot_phys[0]))
#print sol.x[0]

for i in range(len(bmdot_phys)):
    sol = spopt.root(f, [1.], args = (bm_phys[i], bmdot_phys[i]))
    sols.append(sol.x[0])
    #print sol.x[0]    

plt.figure(figsize=(8,6), dpi = 400)
#print sols
plt.hist(sols, bins = 50)
plt.axvline(np.mean(sols), c = 'r', linestyle = '--')
print str(np.round(np.mean(sols), 4))
plt.xlabel('best-fit $q$ value', size = 14)
plt.ylabel('# of BHs', size = 14)
plt.title('best-fit $q$ values', size = 18)

sol_perc = np.percentile(sols, [14,50,86])
print sol_perc[1] - sol_perc[0], sol_perc[1], sol_perc[2] - sol_perc[1]

plt.show()

hist, bin_edges = np.histogram(sols, bins = 50)
bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])
print np.max(hist), bin_centers[np.argmax(hist)]



    


















    






0.0693
0.00384899395363 0.068533819202 0.00496553036837










    

















    






3911 0.0684625466785

Test M and Mdot





In [68]:



    


alpha = 3.2e4 #propto M
beta = 4.65e19 #propto Mdot

plt.hist(np.log10((beta*bm_phys)/(alpha*bmdot_phys)), bins = 20)



    


















    
Out[68]:








(array([  4.00000000e+00,   4.00000000e+00,   3.10000000e+01,
         1.38000000e+02,   5.34000000e+02,   1.49900000e+03,
         3.43300000e+03,   4.94500000e+03,   4.72500000e+03,
         3.50900000e+03,   2.12600000e+03,   1.14600000e+03,
         5.73000000e+02,   3.37000000e+02,   1.87000000e+02,
         6.80000000e+01,   1.00000000e+01,   2.00000000e+00,
         5.00000000e+00,   1.00000000e+00]),
 array([ 30.67309172,  31.04875633,  31.42442094,  31.80008555,
        32.17575016,  32.55141477,  32.92707938,  33.30274399,
        33.6784086 ,  34.05407321,  34.42973782,  34.80540243,
        35.18106704,  35.55673165,  35.93239626,  36.30806087,
        36.68372548,  37.05939009,  37.4350547 ,  37.81071931,  38.18638392]),
 <a list of 20 Patch objects>)

New Analysis





In [96]:



    


alpha = 3.2e3 #propto Mdot
beta = 4.65e19 #propto M
epsilon = 0.866
eta = 4.705

m_mdot = [bm_phys, bmdot_phys]

def fm(m_mdot, q, k):
    m = m_mdot[0]
    mdot = m_mdot[1]
    return epsilon * np.log10(m*beta) + eta - np.log10(m) - q*np.log10(mdot) - k

def fmdot(m_mdot, q, k):
    m = m_mdot[0]
    mdot = m_mdot[1]
    return epsilon * np.log10(mdot*alpha) + eta - np.log10(m) - q*np.log10(mdot) - k

New Analysis - M





In [103]:



    


opt, cov = spopt.curve_fit(fm, m_mdot, np.zeros(len(bm_phys)))
print opt, np.sqrt(np.diagonal(cov))

x_pts = np.log10(bmdot_phys)
y_pts = np.log10(beta*bm_phys) - np.log10(bm_phys)

plt.scatter(x_pts[:75], y_pts[:75], marker = 'x', color = 'b')



    


















    






[ -0.06424811  20.10683412] [ 0.00049833  0.00624176]










    
Out[103]:








<matplotlib.collections.PathCollection at 0x104fe890>

New Analysis - Mdot





In [114]:



    


opt, cov = spopt.curve_fit(fmdot, m_mdot, np.zeros(len(bm_phys)))
print opt, np.sqrt(np.diagonal(cov))

x_pts = np.log10(bmdot_phys)
y_pts = np.log10(alpha*bmdot_phys) - np.log10(bm_phys)

plt.scatter(x_pts[:75], y_pts[:75], marker = 'x', color = 'b')
plt.plot(np.linspace(np.min(x_pts), np.max(x_pts), 50), opt[0] * np.linspace(np.min(x_pts), np.max(x_pts), 50) + opt[1], linestyle = '--', c = 'g')



    


















    






[ 0.38653647 -4.42506356] [ 0.00371889  0.04658031]










    
Out[114]:








[<matplotlib.lines.Line2D at 0x111d8b10>]

New Analysis - Summary





In [ ]:

Content source: illustris-bh-wisc/illustris-bh

Similar notebooks:

notebook.community | gallery | about