%pylab inline %load_ext autoreload %autoreload 2 from make2DHistogram import * from generateInfo import * # Global paths and constants data_path = "../LG_Kinematics/data/" G_GRAV = 4.54E-48 #units of Mpc^3 Msun^-1 s^-2 KM_TO_MPC = 3.2E-20 # HUBBLE = 0.70 # Original March 13, 2013 # Planck 2013 http://arxiv.org/abs/1303.5076 20 Mar 2013 # H0=67.3+/-1.2 km/s/Mpc HUBBLE = 0.673 E_UNITS = 1.0E-36 L_UNITS = 1.0 halo_finder="BDM" narrow_data = False # The reference is van der Marel et al. ApJ,753:8 V_M31 = array([66.1,-76.3,45.1]) #km/s V_M31_sigma = array([26.7, 19.0, 26.5]) R_M31 = array([-378.9,612.7,-283.1]) # kpc R_M31_sigma = abs(R_M31/20.0) #This is an approximated value, not the actual observational uncertainty print R_M31_sigma, sqrt(sum(R_M31*R_M31)) R_M31 = R_M31 / 1000.0 # in Mpc R_M31_sigma = R_M31_sigma /1000.0 M31_Mass = 1.6E12 M31_Mass_sigma = 0.5E12 MW_Mass = 1.6E12 MW_Mass_sigma = 0.5E12 Mass_Total = 3.14E12 # in units of Msun Mass_Total_sigma = 0.58E12 #generate the velocities, positions and masses n_points = 1000000 V_M31_MonteCarlo = zeros((n_points,3)) R_M31_MonteCarlo = zeros((n_points,3)) M31_Mass_MonteCarlo = (random.normal(M31_Mass,M31_Mass_sigma, n_points)) MW_Mass_MonteCarlo = (random.normal(MW_Mass,MW_Mass_sigma, n_points)) negative_M31 = where(M31_Mass_MonteCarlo < 0.0) negative_MW = where(MW_Mass_MonteCarlo < 0.0) Total_Mass_MonteCarlo = M31_Mass_MonteCarlo + MW_Mass_MonteCarlo Reduced_Mass_MonteCarlo = (M31_Mass_MonteCarlo * MW_Mass_MonteCarlo)/Total_Mass_MonteCarlo for i in range(3): V_M31_MonteCarlo[:,i] = random.normal(V_M31[i],V_M31_sigma[i],n_points) R_M31_MonteCarlo[:,i] = random.normal(R_M31[i],R_M31_sigma[i],n_points) #don't use the items with negative values of the mass positive = where((M31_Mass_MonteCarlo > 0.0) & (MW_Mass_MonteCarlo > 0.0))[0] V_M31_MonteCarlo = V_M31_MonteCarlo[positive,:] R_M31_MonteCarlo = R_M31_MonteCarlo[positive,:] M31_Mass_MonteCarlo = M31_Mass_MonteCarlo[positive] MW_Mass_MonteCarlo = MW_Mass_MonteCarlo[positive] Total_Mass_MonteCarlo = M31_Mass_MonteCarlo + MW_Mass_MonteCarlo Reduced_Mass_MonteCarlo = (M31_Mass_MonteCarlo * MW_Mass_MonteCarlo)/Total_Mass_MonteCarlo # Compute the angular momentum, energy and lambda # array of angular momentum (vector) and the norm J = cross(R_M31_MonteCarlo, V_M31_MonteCarlo) J_norm = sqrt(J[:,0]**2 + J[:,1]**2 + J[:,2]**2) # Mechanical energy R_M31_MonteCarlo_norm = sqrt(R_M31_MonteCarlo[:,0]**2 + R_M31_MonteCarlo[:,1]**2 + R_M31_MonteCarlo[:,2]**2) E_kin = 0.5*(V_M31_MonteCarlo[:,0]**2 + V_M31_MonteCarlo[:,1]**2 + V_M31_MonteCarlo[:,2]**2) E_kin = E_kin * KM_TO_MPC * KM_TO_MPC E_pot = -G_GRAV * Total_Mass_MonteCarlo / R_M31_MonteCarlo_norm E_kin = E_kin/E_UNITS E_pot = E_pot/E_UNITS E = E_kin + E_pot # Lambda lambda_obs = (Reduced_Mass_MonteCarlo**1.5) * (J_norm * L_UNITS * KM_TO_MPC)* sqrt(abs(E * E_UNITS)) / (G_GRAV * (Total_Mass_MonteCarlo)**(5.0/2.0)) # Total mass histogram n, bins, patches = hist(Total_Mass_MonteCarlo, 20, normed=0, color = 'silver') # M31 mass histogram n, bins, patches = hist(M31_Mass_MonteCarlo, 20, normed=0, color = 'silver') #surface density plots in the E-J space. This defines the 1-sigma surface H, xedges,yedges = histogram2d(J_norm,E, bins=(30,38)) X,Y=meshgrid(yedges[0:-1],xedges[0:-1]) th=4000 # 4000, corresponds to 1-sigma, 200 to 3-sigma print 'This is the fraction of points included inside the contour:',sum(H[where(H>th)])/sum(H) levels=array([th]) C=contour(X,Y,H, levels, colors='k', linewidths=2) ax = axes() ax.set_xlabel("Energy") ax.set_ylabel("Angular Momentum") # Generate the pair data for Bolshoi FOF_file=data_path+"Bolshoi_catalog_"+halo_finder+".dat" ID_file=data_path+"IsoPairs_Bolshoi_"+halo_finder+"_ID.dat" output_file=data_path+"pair_physical_vweb_Bolshoi"+halo_finder+".dat" generate_info(FOF_file, ID_file, output_file=output_file, hubble=HUBBLE) # Generate the pair data for CLUES 1 FOF_file=data_path+"IC2710_halos_catalog.dat" ID_file=data_path+"index_LG_IC2710" output_file=data_path+"pair_physical_vweb_IC2710.dat" generate_info(FOF_file, ID_file, output_file=output_file, hubble=0.73) # generate_info(FOF_file, ID_file, output_file=output_file, hubble=HUBBLE) # Generate the pair data for CLUES 2 FOF_file=data_path+"IC16953_halos_catalog.dat" ID_file=data_path+"index_LG_IC16953" output_file=data_path+"pair_physical_vweb_IC16953.dat" generate_info(FOF_file, ID_file, output_file=output_file, hubble=0.73) # generate_info(FOF_file, ID_file, output_file=output_file, hubble=HUBBLE) # Generate the pair data for CLUES 3 FOF_file=data_path+"IC10909_halos_catalog.dat" ID_file=data_path+"index_LG_IC10909" output_file=data_path+"pair_physical_vweb_IC10909.dat" generate_info(FOF_file, ID_file, output_file=output_file, hubble=0.73) # generate_info(FOF_file, ID_file, output_file=output_file, hubble=HUBBLE) #load the data from Bolshoi and make selections on distances and mass bolshoi_file=data_path+"pair_physical_vweb_Bolshoi"+halo_finder+".dat" physical_prop=np.loadtxt(bolshoi_file) #define the main quantities we will use later to prepare the plots radial_vel = physical_prop[:,0] tangential_vel = physical_prop[:,1] distance = physical_prop[:,8] total_mass = physical_prop[:,9]/ 1E12 energy = physical_prop[:,5] + physical_prop[:,6] angular_momentum = physical_prop[:,7] mu = physical_prop[:,10] lambda_orbit = mu**1.5 * (angular_momentum * L_UNITS * KM_TO_MPC)* sqrt(abs(energy * E_UNITS)) / (G_GRAV * (total_mass *1E12 )**(5.0/2.0)) if(narrow_data): index_good_mass = where((total_mass>1.0) & (total_mass<4.0) & (distance>0.70) & (distance<0.80)) index_good_mass = index_good_mass[0] physical_prop = physical_prop[index_good_mass,:] radial_vel = physical_prop[:,0] tangential_vel = physical_prop[:,1] distance = physical_prop[:,8] total_mass = physical_prop[:,9]/1E12 energy = physical_prop[:,5] + physical_prop[:,6] angular_momentum = physical_prop[:,7] mu = physical_prop[:,10] lambda_orbit = mu**1.5 * (angular_momentum * L_UNITS * KM_TO_MPC)* sqrt(abs(energy * E_UNITS)) / (G_GRAV * (total_mass *1E12 )**(5.0/2.0)) #Figure with the separation distribution fig = figure(1, figsize=(9.5,6.5)) ax = axes() ax.set_xlabel("$\mathrm{Pair\ Separation\ [Mpc]}$",fontsize=25) ax.set_ylabel("$\mathrm{N_{pairs}}$",fontsize=25) ticklabels_x = ax.get_xticklabels() ticklabels_y = ax.get_yticklabels() for label_x in ticklabels_x: label_x.set_fontsize(18) label_x.set_family('serif') for label_y in ticklabels_y: label_y.set_fontsize(18) label_y.set_family('serif') n, bins, patches = ax.hist(distance, 20, normed=0, color = 'silver') ax.set_xlim([0.4,1.0]) filename='separation_'+halo_finder # if(narrow_data): # filename=filename+'_narrow' # savefig(filename + '.eps',format = 'eps', transparent=True) # savefig(filename + '.pdf',format = 'pdf', transparent=True) # savefig(filename + '.png',format = 'png', transparent=True) #Figure with the total mass distribution print amin(total_mass), amax(total_mass) fig = figure(1, figsize=(9.5,6.5)) ax = axes() ax.set_xlabel("$\mathrm{Total\ Mass\ [10^{12}\ M_{\odot}]}$",fontsize=25) ax.set_ylabel("$\mathrm{N_{pairs}}$",fontsize=25) ticklabels_x = ax.get_xticklabels() ticklabels_y = ax.get_yticklabels() for label_x in ticklabels_x: label_x.set_fontsize(18) label_x.set_family('serif') for label_y in ticklabels_y: label_y.set_fontsize(18) label_y.set_family('serif') ax.set_xlim([1.0,12.0]) n, bins, patches = ax.hist(total_mass, 20, normed=0, color = 'silver') filename='total_mass_'+halo_finder if(narrow_data): filename=filename+'_narrow' savefig(filename + '.eps',format = 'eps', transparent=True) savefig(filename + '.pdf',format = 'pdf', transparent=True) savefig(filename + '.png',format = 'png', transparent=True) #observational data mean_point = [-110.0, 0.0] sigma_point = [4.4, 34.3] delta_radial = 4.5 delta_tangential = 17.5 #load the data from the CLUES pairs n_clues = 3 clues1=data_path+"pair_physical_vweb_IC10909.dat" clues2=data_path+"pair_physical_vweb_IC16953.dat" clues3=data_path+"pair_physical_vweb_IC2710.dat" clues=chararray([n_clues], itemsize=1024) clues[0]=clues1;clues[1]=clues2; clues[2]=clues3 clues_data = np.zeros([n_clues,4]) if(narrow_data): n_clues = 1 for i in range(n_clues): tmp_dat = np.loadtxt(clues[i]) clues_data[i,0] = tmp_dat[0]# radial clues_data[i,1] = tmp_dat[1]# tangential clues_data[i,2] = tmp_dat[8]# separation clues_data[i,3] = tmp_dat[9]/1E12 #total mass filename='test_rt_'+halo_finder if(narrow_data): filename=filename+'_narrow' n=MakeSimple2DPlot(radial_vel,tangential_vel,filename=filename, clues=clues_data,clues_label='$\mathrm{Constrained\ Sim.}$', contours=0,xlims=[-240.0,0.0], ylims=[0.0,240.0], nbins=13,nxbins=13,nybins=13,bw=1, xlabel='$\mathrm{Radial\ Velocity\ [km\ s^{-1}]}$', ylabel='$\mathrm{Tangential\ Velocity\ [km\ s^{-1}]}$', mean_obs=mean_point, sigma_obs=sigma_point) #find the highest density peak H2d, r_edges, t_edges = histogram2d(radial_vel, tangential_vel, bins=(12,12), range=[[-240.0,0.0],[0.0,240.0]]) max_radial_id = argmax(H2d, axis=0) max_tangential_id = argmax(H2d, axis=1) max_r_id = argmax(H2d[arange(12), max_tangential_id]) max_t_id = argmax(H2d[max_radial_id, arange(12)]) theory_point = array([r_edges[max_r_id]+10.0, t_edges[max_t_id]+10.0]) print 'total number of pairs', size(radial_vel) print 'location of the highest density peak in Bolshoi', theory_point #selection of the number of galaxies around the observational constraints index_obs = where((radial_vel>(mean_point[0]-sigma_point[0])) & (radial_vel<(mean_point[0]+sigma_point[0])) & (tangential_vel(theory_point[0]-sigma_point[0])) & (radial_vel<(theory_point[0]+sigma_point[0])) & (tangential_vel<(theory_point[1]+sigma_point[1]/2.0)) & (tangential_vel>(theory_point[1]-sigma_point[1]/2.0))) print 'number of points around the highest density peak', size(index_max) #find the number of points with radial velocities consistent with observations tang_to_radial = abs(tangential_vel / radial_vel) index = where(tang_to_radial<0.32) print 'number of points with the right radial to tang', size(index), '. In percentage:', 1.0*size(index)/radial_vel.size #Figure in the E-J plane n_clues=3 clues1=data_path+"pair_physical_vweb_IC10909.dat" clues2=data_path+"pair_physical_vweb_IC16953.dat" clues3=data_path+"pair_physical_vweb_IC2710.dat" clues=chararray([n_clues], itemsize=1024) clues[0]=clues1;clues[1]=clues2; clues[2]=clues3 clues_data = np.zeros([n_clues,2]) if(narrow_data): n_clues = 1 for i in range(n_clues): tmp_dat = np.loadtxt(clues[i]) clues_data[i,0] = tmp_dat[5] + tmp_dat[6] # energy clues_data[i,1] = tmp_dat[7] #angular momentum filename='test_EJ_'+halo_finder if(narrow_data): filename=filename+'_narrow' MakeSimple2DPlot(energy,angular_momentum,filename=filename, clues=clues_data, clues_label='$\mathrm{Constrained\ Sim.}$', xlims=[-25.0,0.0], ylims=[0.0,150.0], nbins=13,nxbins=13,nybins=13,bw=1, xlabel='$\mathrm{e_{tot} [10^{-36}} \mathrm{Mpc^2\ s^{-2}]}$', ylabel='$\mathrm{l_{orb} [Mpc\ } \mathrm{km\ s^{-1}]}$', contours=0, X_field=X, Y_field=Y, Z_field=H, levels=levels) #Figure in the E-J plane # Previous plot with 12 bins rather than 13! n_clues=3 jw_bin_number = 12 clues1=data_path+"pair_physical_vweb_IC10909.dat" clues2=data_path+"pair_physical_vweb_IC16953.dat" clues3=data_path+"pair_physical_vweb_IC2710.dat" clues=chararray([n_clues], itemsize=1024) clues[0]=clues1;clues[1]=clues2; clues[2]=clues3 clues_data = np.zeros([n_clues,2]) if(narrow_data): n_clues = 1 for i in range(n_clues): tmp_dat = np.loadtxt(clues[i]) clues_data[i,0] = tmp_dat[5] + tmp_dat[6] # energy clues_data[i,1] = tmp_dat[7] #angular momentum filename='test_EJ_'+halo_finder if(narrow_data): filename=filename+'_narrow' MakeSimple2DPlot(energy,angular_momentum,filename=filename, clues=clues_data, clues_label='$\mathrm{Constrained\ Sim.}$', xlims=[-25.0,0.0], ylims=[0.0,150.0], nbins=jw_bin_number,nxbins=jw_bin_number,nybins=jw_bin_number,bw=1, xlabel='$\mathrm{e_{tot} [10^{-36}} \mathrm{Mpc^2\ s^{-2}]}$', ylabel='$\mathrm{l_{orb} [Mpc\ } \mathrm{km\ s^{-1}]}$', contours=0, X_field=X, Y_field=Y, Z_field=H, levels=levels) index = where((angular_momentum<40.0) & (energy>-20.0) & (energy<-5.0)) print 'approx. number of points inside contour', size(index), 'in percentage', 1.0 * size(index)/energy.size #Figure with the spin parameter distribution fig = figure(1, figsize=(9.5,6.5)) ax = axes() ax.set_xlabel("$\log_{10}\mathrm{\lambda}$",fontsize=25) ax.set_ylabel("$\mathrm{Normalized\ Histogram}$",fontsize=25) ticklabels_x = ax.get_xticklabels() ticklabels_y = ax.get_yticklabels() for label_x in ticklabels_x: label_x.set_fontsize(18) label_x.set_family('serif') for label_y in ticklabels_y: label_y.set_fontsize(18) label_y.set_family('serif') n_clues=3 clues1=data_path+"pair_physical_vweb_IC10909.dat" clues2=data_path+"pair_physical_vweb_IC16953.dat" clues3=data_path+"pair_physical_vweb_IC2710.dat" clues=chararray([n_clues], itemsize=1024) clues[0]=clues1;clues[1]=clues2; clues[2]=clues3 clues_data = np.zeros([n_clues,4]) if(narrow_data): n_clues = 1 for i in range(n_clues): tmp_dat = np.loadtxt(clues[i]) clues_data[i,0] = tmp_dat[5] + tmp_dat[6] # energy clues_data[i,1] = tmp_dat[7] #angular momentum clues_data[i,2] = tmp_dat[9] # total mass clues_data[i,3] = tmp_dat[10] # reduced_mass e_CLUES = clues_data[i,0] l_CLUES = clues_data[i,1] total_mass_CLUES = clues_data[i,2] mu_CLUES = clues_data[i,3] lambda_CLUES = mu_CLUES**1.5 * (l_CLUES * L_UNITS * KM_TO_MPC) * \ sqrt(abs(e_CLUES * E_UNITS)) / (G_GRAV * (total_mass_CLUES )**(5.0/2.0)) lambda_CLUES = log10(lambda_CLUES) y_value=1.9 ax.axvline(x=lambda_CLUES,ymin=0.8,ymax=0.95,c="black",linewidth=2) if(i==0): ax.scatter(lambda_CLUES, y_value, color='white', edgecolors='black', s=200, label='$\mathrm{Constrained\ Sim.}$') else: ax.scatter(lambda_CLUES, y_value, color='white', edgecolors='black', s=200) hist_A, bin_edges_A = histogram(log10(lambda_obs), bins=arange(-3.0,0.5,0.1), density=True) hist_B, bin_edges_B = histogram(log10(lambda_orbit), bins=arange(-3.0,0.5,0.1), density=True) plot(bin_edges_A[:-1],hist_A,color='black',ms=1,linewidth=3.0, label="$\mathrm{Obs.}$") plot(bin_edges_B[:-1],hist_B,'r--',color='black',ms=2,linewidth=3.0, label="$\mathrm{Bolshoi}$") ax.legend(loc=0, scatterpoints=1, prop={'size':20}) ax.set_xlim([-3.0,0.5]) ax.set_ylim([0.0,2.2]) filename='test_lambda_'+halo_finder if(narrow_data): filename=filename+'_narrow' savefig(filename + '.eps',format = 'eps', transparent=True) savefig(filename + '.pdf',format = 'pdf', transparent=True) savefig(filename + '.png',format = 'png', transparent=True) mean_lambda_obs = mean(log10(lambda_obs)) variance_obs = var(log10(lambda_obs)) mean_lambda_sim = mean(log10(lambda_orbit)) variance_sim = var(log10(lambda_orbit)) print 'mean lambda, observations', mean_lambda_obs, variance_obs print 'mean lambda, bolshoi', mean_lambda_sim, variance_sim print 'lambda CLUES', lambda_CLUES # calculate the number of points in the simulation that are consistent with mean value and sigma index_sim = where((log10(lambda_orbit)<(mean_lambda_obs+variance_obs)) & (log10(lambda_orbit) > (mean_lambda_obs-variance_obs))) print 'there are ', size(index_sim), 'pairs consistent with the observations' # mean lambda, observations -1.72013855372 0.0797602386996 # mean lambda, bolshoi -1.470309254 0.137777720569 # lambda CLUES -1.72803661534 # there are 257 pairs consistent with the observations