In [50]:
# Import libaries and illustris library
import numpy as np
import illustris_python as il
import matplotlib.pyplot as plt
%matplotlib inline
import random
%pylab inline
pylab.rcParams['figure.figsize'] = (5, 5)
import data_subroutines.bins as bf
In [53]:
import importlib as imp
bf = imp.reload(bf)
In [28]:
# 0:HaloNumber, 1:SubhaloNumber, 2:StellarMass, 3:GasMass, 4:DarkMatterMass, 5-7:SubhaloPosition x,y,z
# 8-10:CoM Velocity x, CoM Velocity y, CoM Velocity z, 11-13:AngularMomentum x, AngularMomentum y, AngularMomentum z,
# 14:Kappa, 15:Flatness, 16:Ellip, 17:Lambda 18: Index of nearest segment point, 19:absolute distance to filament,
# 20:cos(theta) 21:sin(theta) 22:(theta) 23-25: nearest segement coordinates 26: Index of nearest critical point,
# 27: absolute distanceto crit., 28-30: nearest critical point coordinates, 31: bool in cluster, 32: bool in filament
galaxy_data_file = './data/galaxy_data_full_sim_v3.txt'
galaxy_data = np.loadtxt(galaxy_data_file, comments='#')
# 0-2: CM position, 3: halo mass, 4: index of nearest filament point, 5: filament absolute distance
# 6: index of nearest crit point, 7: crit absolute distance, 8: bool in cluster, 9: bool in filament
halo_data_file = './data/halo_data_full_sim_v2.txt'
halo_data = np.loadtxt(halo_data_file, comments='#')
In [45]:
# mass versus distance to clust or fil
# 0-2: CM position, 3: halo mass, 4: index of nearest filament point, 5: filament absolute distance
# 6: index of nearest crit point, 7: crit absolute distance, 8: bool in cluster, 9: bool in filament
x = bf.log_bin_data(halo_data[halo_data[:,8]>0], 3, 7, 15, 1)
x1 = bf.log_bin_data(halo_data[halo_data[:,9]>0], 3, 5, 15, 1)
plt.xscale('log'); plt.yscale('linear');
plt.xlabel('$Total M /M_*$'); plt.ylabel('Distance [ckpc]'); plt.title('Mass functions');
plt.errorbar(x[:,0], x[:,1], xerr=x[:,2], yerr=x[:,3],
fmt='o', label = 'Distance to critical points')
plt.errorbar(x1[:,0], x1[:,1], xerr=x1[:,2], yerr=x1[:,3],
fmt='o', label = 'Distance to segments')
plt.legend(loc='best')
plt.show()
In [56]:
# recreate mass density in filaments and clusters graph
# halo_data:
# 0-2: CM position, 3: halo mass, 4: index of nearest filament point, 5: filament absolute distance
# 6: index of nearest crit point, 7: crit absolute distance, 8: bool in cluster, 9: bool in filament
V_ill = (75/0.7)**3
crits_mass_array = halo_data[halo_data[:,8]>0]
segs_mass_array = halo_data[halo_data[:,9]>0]
void_mass_array = halo_data[(halo_data[:,8]==0)&(halo_data[:,9]==0)]
total_mass_array = halo_data
#compute mass fraction
total_mass = np.sum(halo_data[:,3]); crits_mass = np.sum(crits_mass_array[:,3]);
segs_mass = np.sum(segs_mass_array[:,3]); voids_mass = np.sum(void_mass_array[:,3]);
print(crits_mass/total_mass, segs_mass/total_mass, voids_mass/total_mass )
In [5]:
import importlib as imp
In [57]:
x = bf.mass_function(crits_mass_array[:,3], 50)
x1 = bf.mass_function(segs_mass_array[:,3], 50)
x2 = bf.mass_function(total_mass_array[:,3], 50)
plt.plot( x[:,1], x[:,2]/(V_ill * 0.021),'.', label = 'Clusters');
plt.plot( x1[:,1], x1[:,2]/(V_ill * 0.33),'.', label = 'Segments');
plt.plot( x2[:,1], x2[:,2]/(V_ill),'.' , label = 'Total');
plt.xscale('log');plt.yscale('log'); plt.legend(loc='lower left');
plt.xlabel('Mass $[M/{M_\odot}]$'); plt.ylabel('dn/dM $[{M_\odot}^{-1} {Mpc}^{-3}]$');
plt.title('Mass function of haloes residing in filaments or clusters \n Abundances in box are rescaled by corresponding volume fractions.')
plt.show()
In [8]:
# 0:HaloNumber, 1:SubhaloNumber, 2:StellarMass, 3:GasMass, 4:DarkMatterMass, 5-7:SubhaloPosition x,y,z
# 8-10:CoM Velocity x, CoM Velocity y, CoM Velocity z, 11-13:AngularMomentum x, AngularMomentum y, AngularMomentum z,
# 14:Kappa, 15:Flatness, 16:Ellip, 17:Lambda 18: Index of nearest segment point, 19:absolute distance to filament,
# 20:cos(theta) 21:sin(theta) 22:(theta) 23-25: nearest segement coordinates 26: Index of nearest critical point,
# 27: absolute distanceto crit., 28-30: nearest critical point coordinates, 31: bool in cluster, 32: bool in filament
In [56]:
galaxy_data_elip = galaxy_data[(galaxy_data[:,14]<0.5) & (galaxy_data[:,32]>0)]
galaxy_data_disk = galaxy_data[(galaxy_data[:,14]>0.5) & (galaxy_data[:,32]>0)]
binned_data_elip = log_bin_data((galaxy_data_elip[galaxy_data_elip[:,2]<10**1.5]), 2, 22, 8,-1)
binned_data_disk = log_bin_data(galaxy_data_disk[galaxy_data_disk[:,2]<10**1.5] , 2, 22, 5,-1)
binned_data_all = log_bin_data((abs(galaxy_data[(galaxy_data[:,32]>0)] )), 2, 22, 10,-1)
x = binned_data_elip[:,0]
y = binned_data_elip[:,1]
x1 = binned_data_all[:,0]
y1 = binned_data_all[:,1]
x2 = binned_data_disk[:,0]
y2 = binned_data_disk[:,1]
plt.xscale('log')
plt.yscale('linear')
plt.xlabel('$Stellar M /M_*$')
plt.ylabel('Mean $\Theta$')
plt.title('Major axis spin alignment vs mass - ellipticals ($\kappa$ < 0.5)')
plt.errorbar(x, y, xerr=binned_data_elip[:,2], yerr=binned_data_elip[:,3], fmt='o', label = 'ETG' )
plt.errorbar(x1, y1, xerr=binned_data_all[:,2], yerr=binned_data_all[:,3] ,fmt='o', label='All')
plt.errorbar(x2, y2, xerr=binned_data_disk[:,2], yerr=binned_data_disk[:,3] ,fmt='o', label='Disk')
plt.axhline(3.141592654/2, color='grey', alpha=0.3, linestyle='--')
plt.legend(loc='best')
plt.show()
x1 = binned_data_all[:,0]
y1 = binned_data_all[:,1]
plt.plot( x1, y1,'.' );
plt.xscale('log')
plt.yscale('linear')
plt.xlabel('$Stellar M /M_*$')
plt.ylabel('Median $\Theta$')
plt.title('Major axis spin alignment vs mass - all')
plt.errorbar(x1, y1, xerr=binned_data_all[:,2], yerr=binned_data_all[:,3] ,fmt='o')
plt.axhline(.5, color='grey', alpha=0.3, linestyle='--')
plt.axhline(3.141592654/2, color='grey', alpha=0.3, linestyle='--')
plt.show()
In [ ]:
In [90]:
ellipticals=galaxy_data[(galaxy_data[:,14]>0.5) & (galaxy_data[:,19]<800),:]
minM=10**1
maxM=10**3
nBins=8
massBins=10**(np.linspace(np.log10(minM),np.log10(maxM),nBins))
mean=np.zeros(nBins-1)
err=np.zeros(nBins-1)
for i in range(nBins-1):
whichGal=np.argwhere( (ellipticals[:,4]>massBins[i]) & (ellipticals[:,4]<massBins[i+1]) )
if (whichGal.size==0):
continue
mean[i]=np.median(ellipticals[whichGal,20])
err[i]=np.std(ellipticals[whichGal,20])/whichGal.size
plt.plot( massBins[:-1], mean,'.' )
plt.errorbar( massBins[:-1], mean,yerr=err )
plt.scatter(ellipticals[:,4],ellipticals[:,20])
plt.xlim(0.9*minM,1.1*maxM)
plt.xscale('log')
plt.yscale('linear')
plt.show()
print(mean)
print(err)
print(massBins)
print(ellipticals.shape)
In [58]:
# 0:HaloNumber, 1:SubhaloNumber, 2:StellarMass, 3:GasMass, 4:DarkMatterMass, 5-7:SubhaloPosition x,y,z
# 8-10:CoM Velocity x, CoM Velocity y, CoM Velocity z, 11-13:AngularMomentum x, AngularMomentum y, AngularMomentum z,
# 14:Kappa, 15:Flatness, 16:Ellip, 17:Lambda 18: Index of nearest segment point, 19:absolute distance to filament,
# 20:cos(theta) 21:sin(theta) 22:(theta) 23-25: nearest segement coordinates 26: Index of nearest critical point,
# 27: absolute distanceto crit., 28-30: nearest critical point coordinates, 31: bool in cluster, 32: bool in filament
def findFrac(gal_data,steps):
lambdaBins = np.linspace(0,0.75,steps)
frac = np.zeros((steps,12))
mass_cut = 10
for i in range (steps):
range_lambda = gal_data[(gal_data[:,17]>i/steps) & (gal_data[:,17]<(i+1)/steps) & (gal_data[:,14]<0.5)]
range_lambda_m1 = range_lambda[range_lambda[:,2]<mass_cut]
fils_range_lambda_m1 = np.argwhere(range_lambda_m1[:,32]==1).size
crits_range_lambda_m1 = np.argwhere(range_lambda_m1[:,31]==1).size
tot_range_lambda_m1 = np.argwhere(range_lambda_m1[:,2]>0).size
range_lambda_m2 = range_lambda[range_lambda[:,2]>mass_cut]
fils_range_lambda_m2 = np.argwhere(range_lambda_m2[:,32]==1).size
crits_range_lambda_m2 = np.argwhere(range_lambda_m2[:,31]==1).size
tot_range_lambda_m2 = np.argwhere(range_lambda_m2[:,2]>0).size
if (tot_range_lambda_m1 != 0):
frac[i,0] = crits_range_lambda_m1/tot_range_lambda_m1
frac[i,1] = fils_range_lambda_m1/tot_range_lambda_m1
frac[i,2] = 1/np.sqrt(crits_range_lambda_m1)
frac[i,3] = 1/np.sqrt(fils_range_lambda_m1)
if (tot_range_lambda_m2 != 0):
frac[i,4] = crits_range_lambda_m2/tot_range_lambda_m2
frac[i,5] = fils_range_lambda_m2/tot_range_lambda_m2
frac[i,6] = 1/np.sqrt(crits_range_lambda_m2)
frac[i,7] = 1/np.sqrt(fils_range_lambda_m2)
return lambdaBins, frac
In [630]:
np.argwhere( (galaxy_data[:,2]<5 ) & (galaxy_data[:,31]==1)).size
Out[630]:
In [65]:
data = findFrac(galaxy_data, 5)
lamda = data[0]
crits_m1 = data[1][:,0]
crits_m2 = data[1][:,4]
fils_m1 = data[1][:,1]
fils_m2 = data[1][:,5]
#plt.errorbar(lamda, crits_m1, yerr=data[1][:,2] ,fmt='o' , label='In clusters with stellar $M/{M_*} < 10$')
#plt.errorbar(lamda, crits_m2, yerr=data[1][:,6] ,fmt='o' , label='In clusters with stellar $M/{M_*}$ > 10')
plt.errorbar(lamda, fils_m1, yerr=data[1][:,3] ,fmt='o' , label='In filaments with stellar $M/{M_*} < 10$')
plt.errorbar(lamda, fils_m2, yerr=data[1][:,7] ,fmt='o' , label='In filaments with stellar $M/{M_*}$ >10')
plt.xscale('linear'); plt.yscale('linear');
plt.ylim([0, 1]); plt.xlim([-.05, .59]);
plt.xlabel(r'$\lambda$'); plt.ylabel('Fraction in cluster/in filament');
plt.title(''); plt.legend(loc='best');
plt.show()
In [ ]:
In [ ]:
# 0:HaloNumber, 1:SubhaloNumber, 2:StellarMass, 3:GasMass, 4:DarkMatterMass, 5-7:SubhaloPosition x,y,z
# 8-10:CoM Velocity x, CoM Velocity y, CoM Velocity z, 11-13:AngularMomentum x, AngularMomentum y, AngularMomentum z,
# 14:Kappa, 15:Flatness, 16:Ellip, 17:Lambda 18: Index of nearest segment point, 19:absolute distance to filament,
# 20:cos(theta) 21:sin(theta) 22:(theta) 23-25: nearest segement coordinates 26: Index of nearest critical point,
# 27: absolute distanceto crit., 28-30: nearest critical point coordinates, 31: bool in cluster, 32: bool in filament
Lamda vs stellar mass
Lamda vs critical distance
In [627]:
data = bin_data(galaxy_data[(galaxy_data[:,14]<0.5) & (galaxy_data[:,32] != 1) & (galaxy_data[:,2] < 5)], 17, 27, 14, 1)
data1 = bin_data(galaxy_data[(galaxy_data[:,14]<0.5) & (galaxy_data[:,32] != 1) & (galaxy_data[:,2] > 5) & (30 > galaxy_data[:,2])], 17, 27, 11, 1)
data2 = bin_data(galaxy_data[(galaxy_data[:,14]<0.5) & (galaxy_data[:,32] != 1) & (galaxy_data[:,2] > 30)], 17, 27, 9, 1)
x = data[:,0]; y = data[:,1]; x1 = data1[:,0]; y1 = data1[:,1]; x2 = data2[:,0]; y2 = data2[:,1]; plt.xscale('linear'); plt.yscale('log'); plt.ylim([10**2, 10**4]);
plt.ylabel('Distance to nearest cluster'); plt.xlabel('$\lambda$'); plt.title('Ellipticals not in filaments ($\kappa$ < 0.5)')
plt.errorbar(x, y, xerr=data[:,2], yerr=data[:,3], fmt='o', c='red', label='Stellar $M/{M_*}$ < 5')
plt.errorbar(x1, y1, xerr=data1[:,2], yerr=data1[:,3], fmt='o', c='green', label='Stellar 30 > $M/{M_*}$ > 5')
plt.errorbar(x2, y2, xerr=data2[:,2], yerr=data2[:,3], fmt='o', c='blue', label='Stellar $M/{M_*}$ > 30')
legend(loc='best')
Out[627]:
In [570]:
# input : array with mass data at mass index, bins: no of bins in log space
# output: binned_mat, 0: median of mass in interval, 1: mass interval step size (dm) 2: dn/dm
data = bin_data(galaxy_data[(galaxy_data[:,14]<0.5) & (galaxy_data[:,32] != 1)], 17, 27, 16, 1)
data1 = log_bin_data(galaxy_data[(galaxy_data[:,14]<0.5) & (galaxy_data[:,32] == 1)], 2, 17, 15, 1)
x = (data[:,0])
y = data[:,1]
#x1 = (data1[:,0])
#y1 = data1[:,1]
plt.plot( x, y,'.', c='red');
#plt.plot( x1, y1,'.', c='blue');
plt.xscale('linear')
plt.yscale('linear')
plt.ylabel('$ y $')
plt.xlabel('$ x $')
plt.title('ellipticals ($\kappa$ < 0.5)')
plt.errorbar(x, y, xerr=data[:,2], yerr=data[:,3], fmt='o' )
#plt.errorbar(x1, y1, xerr=data1[:,2], yerr=data1[:,3], fmt='o' )
#plt.axhline(3.141592654/2, color='grey', alpha=0.3, linestyle='--')
plt.show()
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]:
In [ ]: