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


Populating the interactive namespace from numpy and matplotlib

In [53]:
import importlib as imp
bf = imp.reload(bf)

everything in units of h! (h= 0.7)


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 )


0.433059160264 0.441616986353 0.125323853382

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()


/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/fromnumeric.py:2889: RuntimeWarning: Mean of empty slice.
  out=out, **kwargs)
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:80: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:135: RuntimeWarning: Degrees of freedom <= 0 for slice
  keepdims=keepdims)
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:105: RuntimeWarning: invalid value encountered in true_divide
  arrmean, rcount, out=arrmean, casting='unsafe', subok=False)
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:127: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)

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)


[-0.21300785  0.06711576  0.03710117  0.07841519  0.07478181 -0.08878566
 -0.03220624]
[ 0.02037728  0.00549812  0.00224463  0.00191571  0.00260292  0.00447742
  0.0105819 ]
[   10.            19.30697729    37.2759372     71.9685673    138.94954944
   268.26957953   517.94746792  1000.        ]
(1110, 32)

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]:
812

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')


/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/fromnumeric.py:2889: RuntimeWarning: Mean of empty slice.
  out=out, **kwargs)
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:80: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:135: RuntimeWarning: Degrees of freedom <= 0 for slice
  keepdims=keepdims)
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:105: RuntimeWarning: invalid value encountered in true_divide
  arrmean, rcount, out=arrmean, casting='unsafe', subok=False)
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:127: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)
Out[627]:
<matplotlib.legend.Legend at 0x10c087208>

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()


/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/fromnumeric.py:2889: RuntimeWarning: Mean of empty slice.
  out=out, **kwargs)
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:80: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:135: RuntimeWarning: Degrees of freedom <= 0 for slice
  keepdims=keepdims)
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:105: RuntimeWarning: invalid value encountered in true_divide
  arrmean, rcount, out=arrmean, casting='unsafe', subok=False)
/Library/Frameworks/Python.framework/Versions/3.6/lib/python3.6/site-packages/numpy/core/_methods.py:127: RuntimeWarning: invalid value encountered in double_scalars
  ret = ret.dtype.type(ret / rcount)

In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]:


In [ ]: