diff --git a/maxwellian_kicks_python/kicks.py b/maxwellian_kicks_python/kicks.py
new file mode 100644
index 0000000000000000000000000000000000000000..fcb0f43ca69154380f2d766039acddcf225be7d6
--- /dev/null
+++ b/maxwellian_kicks_python/kicks.py
@@ -0,0 +1,136 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Aug 15 09:20:27 2019
+
+@author: m13239
+"""
+
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Aug 15 09:05:53 2019
+
+@author: m13239
+"""
+
+## KICK VELOCITIES ACCORDING TO A MAXWELLIAN DISTRIBUTION
+
+# We start by importing the necessary
+import math
+
+import matplotlib.pyplot as plt
+
+import numpy as np
+
+from scipy.special import erf
+
+import random
+
+# Hnn
+
+# http://pulsars.info/Thesis.pdf
+
+# https://arxiv.org/pdf/1807.11489.pdf
+
+# http://mathworld.wolfram.com/MaxwellDistribution.html
+
+# Give different formulae
+
+# Especially the first factor is a bit unclear sometimes
+
+​# We define the probability density function and the cumulative probability
+# density
+
+def maxwell_prob_density(vkick, dispersion):
+
+ return np.sqrt(2.0/np.pi)* (vkick**2)/(dispersion**3) * np.exp(-((vkick**2)/(2*(dispersion**2))))
+
+​
+
+def maxwell_cum_density(vkick, dispersion):
+
+ return erf((vkick)/(np.sqrt(2)*dispersion)) - ((vkick * np.exp(-((vkick**2)/(2*(dispersion**2)))))/(dispersion))*(np.sqrt(2.0/math.pi))
+
+# sampling from distribution:
+
+def sample_maxwell(sigma):
+
+ X = random.random()
+
+ Y = random.random()
+
+ W = random.random()
+
+ Z = random.random()
+
+​
+
+ s = sigma*math.sqrt(-2*math.log(1-X))
+
+ p = sigma*math.sqrt(-2*math.log(1-Y))
+
+​
+
+ theta = 2*math.pi*Y
+
+ phi = 2*math.pi*Z
+
+​
+
+ u = s * math.cos(theta)
+
+ v = s * math.sin(theta)
+
+ w = p * math.cos(phi)
+
+​
+
+ kick = math.sqrt((u**2) + (v**2) + (w**2))
+
+ return kick
+
+​
+
+​# We set the velocity dispersion of the maxwellian distribution:
+
+dispersion = 265
+
+​# We set the number of samplings:
+
+samples = 50000
+
+# And then we sample from that distribution according to the number of times set above (and obviously the defined velocity dispersion)
+
+sampled_array = [sample_maxwell(dispersion) for _ in range(samples)]
+
+
+
+v_arr = np.arange(0, 1500, 1) # This arranges values from 0 to 1500 with steps separated by 1 (0,1,2,3,...,1498,1499)
+# We then compute the probability density and the cumulative probability density for the v_arr defined above.
+mw_prob = maxwell_prob_density(v_arr, dispersion)
+
+mw_cum = maxwell_cum_density(v_arr, dispersion)
+
+​# The idea is that now we can compute on top of the sampling of velocities the maxwellian distribution itself.
+
+fig, axes = plt.subplots(nrows=1, ncols=2)
+
+axes[0].plot(v_arr, mw_prob)
+
+axes[0].hist(sampled_array, density=True, bins=100)
+
+axes[0].set_ylabel('Probability density')
+
+axes[0].set_xlabel('v (km s-1)')
+
+​
+
+axes[1].plot(v_arr, mw_cum)
+
+axes[1].set_ylabel('Cumulative density')
+
+axes[1].set_xlabel('v (km s-1)')
+
+​
+
+plt.show()
+
diff --git a/mcluster_script.sh b/mcluster_script.sh
index 985606886bf4ee30bbf8f41e7362f6434b56f043..ecdc06c23c362a2367bf99c628ba1c65b4cbacbe 100755
--- a/mcluster_script.sh
+++ b/mcluster_script.sh
@@ -25,7 +25,8 @@ MCLUSTER_EXECUTABLE='/user/HS103/m13239/Desktop/mcluster/mcluster'
TARGET_DIR=$MCLUSTER_CONFIG_DIR
# Store settings for mcluster in here.
-MCLUSTER_COMMANDS="-C5 -M 5000 -R 0.5 -P 0 -f 1 -m 0.08 -m 150 -B 0 -Z 0.0001 -o"
+MCLUSTER_COMMANDS="-C5 -M 5000 -R 1 -P 0 -f 1 -m 0.08 -m 150 -B 0 -Z 0.001 -o"
+#MCLUSTER_COMMANDS="-C5 -M 5000 -R 0.5 -P 0 -f 1 -m 0.08 -m 150 -B 0 -Z 0.001 -o"
#MCLUSTER_COMMANDS="-C5 -M 50000 -R 2 -P 0 -f 1 -m 0.08 -m 150 -B 0 -Z 0.0001 -o"
#MCLUSTER_COMMANDS="-C5 -M 50000 -R 2 -P 1 -W 3.0 -f 1 -B 0 -o"
diff --git a/results/cluster properties/Cluster properties/CLUSTER_PROPERTIES.m b/results/cluster properties/Cluster properties/CLUSTER_PROPERTIES.m
index d770d6d6081ce16bd50465cf515f65eb7068d090..017dcdfce6119a0c3750f6e71ab986ef16dacc97 100644
--- a/results/cluster properties/Cluster properties/CLUSTER_PROPERTIES.m
+++ b/results/cluster properties/Cluster properties/CLUSTER_PROPERTIES.m
@@ -37,17 +37,37 @@ elseif k==2
end
-%The first thing to compute is the retention fraction for different values
-%of the mass of the cluster and half-mass radius for the kick prescription
-%chosen for different black hole masses.
-M_cl=linspace(10^3,10^6,1000);r_h=linspace(0.2,5,1000);
-[M_cl,r_h]=meshgrid(M_cl,r_h);
-v_esc=fun_v_esc(r_h,M_cl);
+%The first thing to compute is the probability of escape as a function of
+%the BH mass and escape velocity according to the metallicity and the kick
+%prescription chosen. This way we can use it in conjunction with the next
+%plot to know which cluster parameters yield such escape velocity.
+
+%To do that we first compute the minimum ZAMS mass that yields a BH, as now
+%we are only interested in studying those stars that become BHs.
+vec_M_zams_test=linspace(7,40,2000);
+vec_M_bh_test=fun_M_rem(Z,vec_M_zams_test);
+ll=find(abs(vec_M_bh_test-3)<1e-2);
+M_z_min_bh_test=vec_M_zams_test(ll(1));
+
+%Then we can start doing the contour plot:
+vec_M_zams=linspace(M_z_min_bh_test,150,2000);
+vec_v_esc=linspace(5,60,2000);
+[M_zams,v_esc]=meshgrid(vec_M_zams,vec_v_esc);
+M_bh=fun_M_rem(Z,M_zams);
%For the probability we use M_ns=1.4, sig_ns=265km/s. As proposed by Hobbs,
%Lorimer,Lyne and Kramer (2005).
M_ns=1.4;sig_ns=265;
-%We will do this for different masses:
-vec_M_bh=[10 20 30 40];
+prob=prob_escape(Z,M_zams,M_ns,sig_ns,v_esc,k);
+%And then:
+figure(1)
+contour(M_zams,v_esc,prob,'ShowText','on')
+title(['Probability of escape. ' num2str(kick)])
+xlabel('M_{BH} (M_{sun})')
+ylabel('v_{esc} (km/s)')
+
+
+%{
+vec_M_zams=[10 20 30 40];
p=prob_escape(Z,vec_M_bh(1),M_ns,sig_ns,v_esc,k);f_ret=1-p;
figure(1)
subplot(2,2,1)
@@ -85,15 +105,21 @@ ylabel('r_h (pc)')
set(gca, 'YScale', 'log')
legend(['M=' num2str(vec_M_bh(4))])
+%}
+%{
%We then want to plot the escape velocity as a function of the mass of the
%cluster and the radius in parsecs:
+M_cl=linspace(10^3,10^6,1000);r_h=linspace(0.2,5,1000);
+[M_cl,r_h]=meshgrid(M_cl,r_h);
+v_esc=fun_v_esc(r_h,M_cl);
figure(2)
contour(M_cl./(10^5),r_h,v_esc,[10,20,30,40,50,75,100],'ShowText','on')
title('Contour plot for the escape velocity (km/s) as a function of M_5 and r_h')
xlabel('M_5 (M_{cl}/(10^5 M_{sun})')
ylabel('r_h (pc)')
+
%{
%We then plot the probability of escape for different BH masses and escape
%velocities:
@@ -110,7 +136,7 @@ ylabel('Escape velocity (in km/s)')
%of escape and the fallback fraction in one plot by using subplots:
vec_M_zams=linspace(0.07,150,1000);
vec_M_co=zeros(1,length(vec_M_zams));vec_M_rem=zeros(1,length(vec_M_zams));
-vec_f_fb=zeros(1,length(vec_M_zams));vec_M_zams=linspace(10,100,2000);
+vec_f_fb=zeros(1,length(vec_M_zams));vec_M_zams=linspace(10,150,2000);
for ii=1:length(vec_M_zams)
M_zams=vec_M_zams(ii);
vec_M_co(ii)=fun_M_co(Z,M_zams);vec_M_rem(ii)=fun_M_rem(Z,M_zams);
@@ -169,4 +195,5 @@ title(['Fallback fraction. Z=' num2str(Z)])
xlabel('M_{rem} (in M_{sun})')
ylabel('Fallback fraction')
+%}
end
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/GUIDE.m b/results/cluster properties/Cluster properties/GUIDE.m
new file mode 100644
index 0000000000000000000000000000000000000000..cc35783fcff486788cfe07372ba4785a1295dac9
--- /dev/null
+++ b/results/cluster properties/Cluster properties/GUIDE.m
@@ -0,0 +1,112 @@
+%% GUIDE FOR THE USE OF THE CODES
+%{
+In order to compute the different magnitudes of interest there are
+different codes one needs to consider. Below we describe the important
+ones. (There are others, but are very especific, and only used for the
+computation of some plots in some order as to put them in the report).
+%}
+
+%% ESSENTIAL FUNCTIONS
+%{
+The following functions are the essential ones
+%}
+%% prob_escape.m
+%{
+It computes the probability of escape for a given star.
+It requires:
+ fun_M_rem.m
+ fun_fb_rem.m
+And it gives the option to use three different kick prescriptions with the
+change of the variable k, as explained in the function script itself
+%}
+%% prob_esc_ns.m
+%{
+It computes the probability of escape for a neutron star.
+IT does not require any additional function.
+%}
+%% fun_fb.m
+%{
+It computes the fallback fraction after a supernova kick.
+It requires:
+ fun_M_co.m
+%}
+%% fun_M_co.m
+%{
+It computes the mass of the CO core
+It does not require any additional function.
+%}
+%% fun_M_rem.m
+%{
+It computes the remnant mass.
+It requires:
+ fun_M_co.m
+%}
+%% fun_v_esc.m
+%{
+It computes the escape velocity of a cluster.
+It does not require any additional function
+%}
+%% fun_for_contour_fractions_3.m
+%{
+It gives the number fraction of BHs compared to the initial ones (just the
+IMF)
+It requires loading cpdf.mat and it is only usable after a meshgrid
+command. That is, it only works if the matrix for the escape velocity and
+mass of the cluster have the same size.
+It requires:
+ fun_M_rem.m
+ prob_escape.m
+ prob_esc_ns.m
+%}
+%% fun_for_mass_histograms.m
+%{
+It yields a vector of the final and initial masses of BHs in a cluster.
+IMPORTANT: Only works for scalar M_{cl} and v_{esc}, does not work if one
+uses vectors as inputs.
+%}
+%% fun_for_contour_fractions_all_plots
+%{
+This function allows for the plotting of a contour plot for
+the retention fraction of BHs for different stellar cluster masses and
+half-mass radii for the three different kick prescriptions at once.
+%}
+%% ESSENTIAL SCRIPTS
+%{
+And now some functions that use the ones defined above to compute the main
+magnitudes we are interested in.
+%}
+%% CLUSTER_PROPERTIES.m
+%{
+This function does different things. First of all if shows a contour plot
+ for the retention fraction of stars in a cluster ,according to the
+ probability function derived (as compared to the more 'random' approach
+ considered in the SAMPLING.m code where we have defined stars), as a
+ function of the mass of the cluster and the half-mass radii for different
+ black hole masses. It also plots the escape velocity as a function of
+ these same variables.
+ Moreover, it computes the probability of escape as a function of the
+ black hole mass and escape velocity.
+ Finally, in a last plot, it shows the remnant mass (M_rem) and CO core
+ mass (M_co) as a function of M_zams, as well as the probability of
+ escape, and the fallback fraction for different masses according to the
+ kick prescription chosen.
+%}
+%% SAMPLING.m
+%{
+ This function does different things. First of all it samples a
+ cluster with the specifications chosen and then plots a histogram for the
+ different BH masses both before and after supernova kicks. Then it also
+ does 100 different samplings (vec_ii can be changed if one wants a
+ different number of samplings) in order to show how a different sampling
+ affects the results obtained.
+ Then it also plots a contour plot for the chosen kick prescription
+ showing the BH retention fraction for different cluster masses and
+ half-mass radii.
+
+There is also the option to compute the retention fraction and mass
+fraction for the three different kick prescriptions by using:
+ fun_for_contour_fractions_all_plots.m
+And one can also focus on those same results but for lower mass clusters by
+using:
+ fun_for_contour_fractions_low_mass.m
+%}
diff --git a/results/cluster properties/Cluster properties/HYBRID_low_radii.png b/results/cluster properties/Cluster properties/HYBRID_low_radii.png
new file mode 100644
index 0000000000000000000000000000000000000000..9ad833c885ffa9a97357838db1d43b3228bc9bf4
Binary files /dev/null and b/results/cluster properties/Cluster properties/HYBRID_low_radii.png differ
diff --git a/results/cluster properties/Cluster properties/M_rem_M_co_dif_metal.png b/results/cluster properties/Cluster properties/M_rem_M_co_dif_metal.png
new file mode 100644
index 0000000000000000000000000000000000000000..295bfe6e091f2f7e1e4a077eaaea9a64c10874ad
Binary files /dev/null and b/results/cluster properties/Cluster properties/M_rem_M_co_dif_metal.png differ
diff --git a/results/cluster properties/Cluster properties/SAMPLING.m b/results/cluster properties/Cluster properties/SAMPLING.m
index e8b0f994fc53d5a53098cd8845c8b6c127f602b4..6e60bb264515d393efc2be677d8b80bad27b277f 100644
--- a/results/cluster properties/Cluster properties/SAMPLING.m
+++ b/results/cluster properties/Cluster properties/SAMPLING.m
@@ -28,11 +28,11 @@ Description:
%function [f_N_bh_relative,M_cl,r_h]=SAMPLING(Z,M_cl,r_h,k)
function SAMPLING(Z,M_cl,r_h,k)
if k==0
- kick='W/o fallback,reduced sigma';
+ kick='NEUTRINO-DRIVEN NK';
elseif k==1
- kick='W/ fallback,reduced sigma';
+ kick='HYBRID MODEL';
elseif k==2
- kick='W/ fallback,sigma not reduced';
+ kick='STANDARD FALLBACK-CONTROLLED NK';
end
%{
We do the sampling by calculating the cumulative probability distribution
@@ -49,8 +49,9 @@ Below we show what that process is:
%We start by defining the IMF:
alfa1=0.3;alfa2=1.3;alfa3=2.3;
-imf= @(M) (M<0.08).*(M.^(-alfa1)) + (M<0.5 & M>=0.08).*((M.^(-alfa2))*(0.08)) + (M>=0.5).*((M.^(-alfa3))*(0.5));
-%We then define the maximum value for the cumulative probability density
+imf= @(M) ((M<0.08).*(M.^(-alfa1))) + ((M<0.5 & M>=0.08).*((M.^(-alfa2))*(0.08))) + ((M>=0.5).*((M.^(-alfa3))*(0.04)));
+%We then define the maximum value for the cumulative probability
+%distribution
%function so that we can normalise the function itself so it can be treated
%as a probability.
cpdf_max=integral(imf,0.07,150,'RelTol',1e-6,'AbsTol',1e-6);
@@ -59,8 +60,6 @@ cpdf_max=integral(imf,0.07,150,'RelTol',1e-6,'AbsTol',1e-6);
%mass. So that we asign different masses until we reach the maximum value
%we want it to have.
%For example:
-%We initialize the total mass to be 0 before we start "creating" stars:
-M_tot=0;
vec_M_examples=[linspace(0.07,2,15000) linspace(2,150,25000)];
eval_fun=[];
for ii=1:length(vec_M_examples)
@@ -71,7 +70,8 @@ end
%}
load('cpdf')
-%{
+
+
%In order to do the sampling and then compute how many BHs remain in the
%cluster after the supernova kicks we will use fun_for_mass_histograms.m
%For that we first need to compute the escape velocity:
@@ -79,36 +79,39 @@ v_esc=fun_v_esc(r_h,M_cl);
%And then:
[vec_m_bh_initial,vec_m_bh_final]=fun_for_mass_histograms(M_cl,Z,v_esc,k);
figure(1)
-h1=histogram(vec_m_bh_initial,35,'FaceColor','b');
-h1.BinWidth = 1.5;
-title(['Number of BHs for M_{cl}=' num2str(M_cl) ' M_{sun}, r_h=' num2str(r_h) 'pc & Z=' num2str(Z) '. ' num2str(kick)])
+h1=histogram(vec_m_bh_initial,20,'FaceColor','b');
+h1.BinWidth = 3;
+title(['Number of BHs for M_{cl}=' num2str(M_cl) ' M_{sun}, r_h=' num2str(r_h) ' pc & Z=' num2str(Z) '. ' num2str(kick)])
xlabel('M_{BH} (M_{sun})')
ylabel('Number of BHs')
xlim([0 100])
%xticks([0 10 20 30 40 50 60 70 80 90 100])
%set(gca,'YScale','log')
hold on
-h2=histogram(vec_m_bh_final,35,'FaceColor','r');
-h2.BinWidth = 1.5;
+h2=histogram(vec_m_bh_final,20,'FaceColor','r');
+h2.BinWidth = 3;
legend('IMF (before kicks)','after kicks')
%And now we show the mean mass of BHs before and after the SN kicks:
-disp('FOR THE SECTION ON THE MASS HISTOGRAM:')
+disp(['FOR THE SECTION ON THE MASS HISTOGRAM: (' num2str(kick) ')'])
disp(' MEAN MASS OF BHs BEFORE SN KICKS')
mean_m_bh_initial=sum(vec_m_bh_initial)/length(vec_m_bh_initial);
disp(mean_m_bh_initial)
disp(' MEAN MASS OF BHs AFTER SN KICKS')
-mean_m_bh_final=sum(vec_m_bh_initial)/length(vec_m_bh_final);
+mean_m_bh_final=sum(vec_m_bh_final)/length(vec_m_bh_final);
disp(mean_m_bh_final)
+%{
%It is also important to study the stochasticity of the sampling method.
%That is why we compute the number stars, BHs and NSs, as well as some
%important fractions for different samplings in order to study these
%variations.
+%We initialize the total mass to be 0 before we start "creating" stars:
+M_tot=0;
-vec_ii=[1:100];
+vec_ii=[1:10];
mat_m=[];
N_tot=zeros(1,length(vec_ii));
M_tot=zeros(1,length(vec_ii));
@@ -150,6 +153,13 @@ M_stars(ii)=sum(vec_m(find(vec_m<M_z_min_ns_test)));
M_ns(ii)=M_tot(ii)-M_bh(ii)-M_stars(ii);
end
+
+%{
+disp('mean mass of stars')
+mean_mass_stars=sum(vec_m)/length(vec_m);
+disp(mean_mass_stars)
+%}
+
%Now that we have sampled it different times we have to compute the results
%after the supernova kicks according to the prescription chosen:
for ii=1:length(vec_ii)
@@ -211,7 +221,7 @@ vec_bh_ns=fun_M_rem(Z,vec_bh_ns);
M_bh_final(ii)=sum(vec_bh_ns_final(find(vec_bh_ns_final>3)));
M_ns_final(ii)=sum(vec_bh_ns_final(find(vec_bh_ns_final<3)));
M_stars_final(ii)=M_tot_final(ii)-M_bh_final(ii)-M_ns_final(ii);
- f_m_bh(ii)=M_bh_final(ii)/M_bh(ii);f_m_ns(ii)=M_ns_final(ii)/M_ns(ii);f_m_stars=M_stars_final(ii)/M_tot(ii);
+ f_m_bh(ii)=M_bh_final(ii)/M_tot(ii);f_m_ns(ii)=M_ns_final(ii)/M_ns(ii);f_m_stars=M_stars_final(ii)/M_tot(ii);
end
%And now we can plot this:
figure(2)
@@ -221,7 +231,7 @@ xlabel('Different samplings')
ylabel('Number of stars')
hold on
plot(vec_ii,N_tot_final,'-r')
-legend('IMF (before kicks)',['after kicks,' num2str(kick)],'Location','southwest')
+legend('IMF (before kicks)',['after kicks,' num2str(kick)],'Location','southwest(1-feval(@prob_esc_ns,265,mat_v_esc(dd,ff)))')
legend('boxoff')
figure(3)
@@ -270,15 +280,16 @@ var=(f_N_bh-mean_f_N_bh).^2;
sigma_f_N_bh=sqrt(sum(var)/length(vec_ii));
disp('Standard deviation in the number fraction of BHs retained:')
disp(sigma_f_N_bh)
-
+%}
%And now it would be interesting to do a contour plot for different cluster
%masses and half-mass radii. To do that we use the function
%fun_for_contour_fractions_3.m
+%{
%We define the following vectors:
-vec_M_cl=linspace(10^3,10^6,100);
-vec_r_h=linspace(0.1,25,100);
+vec_M_cl=linspace(10^4,10^5,10);
+vec_r_h=linspace(0.1,15,10);
%vec_v_esc=linspace(5,60,100);
%And with this we form the meshgrid:
[M_cl,r_h]=meshgrid(vec_M_cl,vec_r_h);
@@ -293,24 +304,26 @@ figure(4)
contour(M_cl/(10^5),r_h,f_N_bh_relative);
colorbar
%clabel(C,h)
-%title(['BH retention fraction. Z=' num2str(Z) '. ' num2str(kick)])
-title(['BH retention fraction. Z=' num2str(Z) '. W/ fallback , sigma 265'])
+title(['BH retention fraction. Z=' num2str(Z) '. ' num2str(kick)])
xlabel('M_{cl} (10^5 M_{sun})')
ylabel('r_h (in pc)')
set(gca, 'YScale', 'log')
+
%set(gca, 'XScale', 'log')
+
%}
%{
%To obtain the contour plot for all prescriptions at once:
%vec_M_cl=linspace(10^4,10^5,100);
-vec_M_cl=logspace(4,5,100);%100 logarithmically spaced numbers.
-vec_r_h=linspace(0.1,15,100);
+vec_M_cl=logspace(4,6,100);%100 logarithmically spaced numbers.
+vec_r_h=linspace(0.1,1.5,100);
%vec_v_esc=linspace(5,60,100);
%And with this we form the meshgrid:
[M_cl,r_h]=meshgrid(vec_M_cl,vec_r_h);
mat_v_esc=fun_v_esc(r_h,M_cl);
fun_for_contour_fractions_all_plots(M_cl,Z,r_h)
%}
+%{
%To obtain the results for the mean values of low-mass clusters:
vec_M_cl=linspace(10^4,10^5,100);
%vec_M_cl=logspace(4,5,100);%100 logarithmically spaced numbers.
@@ -320,4 +333,5 @@ vec_r_h=linspace(0.1,15,100);
[M_cl,r_h]=meshgrid(vec_M_cl,vec_r_h);
mat_v_esc=fun_v_esc(r_h,M_cl);
fun_for_contour_fractions_low_mass(M_cl,Z,r_h)
+%}
end
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/cpdf.mat b/results/cluster properties/Cluster properties/cpdf.mat
index a23b6c3dd0e23f9c6fc88b526b9c965d47ec6d29..a8d413cebef9d7f1b79bf86d717d201bd5c4090d 100644
Binary files a/results/cluster properties/Cluster properties/cpdf.mat and b/results/cluster properties/Cluster properties/cpdf.mat differ
diff --git a/results/cluster properties/Cluster properties/fun_for_contour_fractions_3.m b/results/cluster properties/Cluster properties/fun_for_contour_fractions_3.m
index c8fb5a15a53eb077e09a06c58b4ab2a79ae3a0c8..6496d482ad2a7ed1c1406d79518d185c795062dd 100644
--- a/results/cluster properties/Cluster properties/fun_for_contour_fractions_3.m
+++ b/results/cluster properties/Cluster properties/fun_for_contour_fractions_3.m
@@ -31,7 +31,7 @@ load('cpdf')
vec_m_initial=[];%We will use this matrix to save the different masses. We initialize this variable so that we store every mass vector in a different row:
vec_M_cl=M_cl(1,:); %We define this vector so that the sampling for every cluster mass is only done once. This way it is faster.
-tic
+
%And then we can sample the IMF:
for s=1:length(vec_M_cl)
%M_cl_it=M_cl(dd,ff);
@@ -61,7 +61,6 @@ M_z_min_ns_test=vec_M_zams_test(ll(1));
%And then we compute the total number of BHs (stars that will become one):
N_bh_initial(s)=sum(vec_m_initial(s,:)>M_z_min_bh_test);
-toc
end
%Now, we take the values we obtained before and use them to determine which
@@ -69,8 +68,7 @@ toc
mat_v_esc=fun_v_esc(r_h,M_cl);
for dd=1:length(r_h(1,:))
for ff=1:length(r_h(:,1))
-vec_m_final=[];%We initialize the vector that will store the retained masses.
-tic
+vec_m_final=[];%We initialize the vector that will store the retained masses.
%In order to not use so many if statements we are going to find which
%positions in the vector vec_m_initial have a mass smaller than the minimum
%neutron star mass. Because this stars are retained and considering them in
@@ -126,7 +124,6 @@ vec_bh_ns=fun_M_rem(Z,vec_bh_ns);
%cluster after considering the supernova kicks
N_bh_final(dd,ff)=sum(vec_bh_ns_final>3);
f_N_bh_relative(dd,ff)=N_bh_final(dd,ff)./N_bh_initial(ff);
- toc
end
end
end
diff --git a/results/cluster properties/Cluster properties/fun_for_contour_fractions_all_plots.m b/results/cluster properties/Cluster properties/fun_for_contour_fractions_all_plots.m
index 52ad520762d8d9ba3d7ab5039eff9d1107e76363..f7a34892633d3f1e200556e10ea52685702ef98b 100644
--- a/results/cluster properties/Cluster properties/fun_for_contour_fractions_all_plots.m
+++ b/results/cluster properties/Cluster properties/fun_for_contour_fractions_all_plots.m
@@ -97,7 +97,9 @@ vec_bh_ns=nonzeros(vec_bh_ns)';%Vector that contains the initial stars that are
%And now we do an if statement for this stars:
ll=length(vec_m_final); %We initialize ll so that we fill the vector after that position.
%We calculate only once the probability for the initial vector of masses:
- f_ret=(1-feval(@prob_escape,Z,vec_bh_ns,1.4,265,mat_v_esc(dd,ff),k)); %As it can be seen by default we choose
+ f_ret=(1-feval(@prob_escape,Z,vec_bh_ns,1.4,265,mat_v_esc(dd,ff),k));
+ f_ret_ns=(1-feval(@prob_esc_ns,265,mat_v_esc(dd,ff)));
+ %As it can be seen by default we choose
%a velocity dispersion of 265km/s and a mean neutron star mass of 1.4
%solar masses.
%It is also important to see that later, in order to really compute the
@@ -121,7 +123,7 @@ vec_bh_ns=fun_M_rem(Z,vec_bh_ns);
end
else
rr=rand;
- if rr<=(1-feval(@prob_esc_ns,265,mat_v_esc(dd,ff)))
+ if rr<=f_ret_ns
vec_m_final(ll)=vec_bh_ns(jj);
vec_bh_ns_final(nn)=vec_bh_ns(jj);
ll=ll+1;nn=nn+1;
@@ -142,27 +144,97 @@ end
end
figure()
subplot(1,2,1)
-contour(M_cl/(10^4),r_h,f_N_bh_relative);
+contour(M_cl/(10^5),r_h,f_N_bh_relative);
+hold on
colorbar
-%clabel(C,h)
title(['BH retention fraction. Z=' num2str(Z) '. ' num2str(kick)])
-%title(['BH retention fraction. Z=' num2str(Z) '. W/ fallback , sigma 265'])
-xlabel('M_{cl} (10^4 M_{sun})')
+xlabel('M_{cl} (10^5 M_{sun})')
ylabel('r_h (in pc)')
-set(gca, 'YScale', 'log')
-%set(gca, 'XScale', 'log')
+%set(gca, 'YScale', 'log')
+
+%And only use this if you want to plot the datapoints for certain clusters: (10^1
+%because we are plotting M/10^5)
hold on
+
+h1(1)=plot(0.67*10^1,0.59,'.k');hold on %2298
+text(0.67*10^1,0.59,'2298');
+h1(2)=plot(0.9*10^1,0.64,'.k');hold on %6218
+text(0.9*10^1,0.64,'6218');
+h1(3)=plot(0.85*10^1,0.42,'.k');hold on %6254
+text(0.85*10^1,0.42,'6254');
+h1(4)=plot(1*10^1,0.23,'.k');hold on %6341
+text(1*10^1,0.23,'6341');
+h1(5)=plot(0.31*10^1,1.08,'.k');hold on %6352
+text(0.31*10^1,1.08,'6352');
+h1(6)=plot(0.26*10^1,0.17,'.k');hold on %6397
+text(0.26*10^1,0.17,'6397');
+h1(7)=plot(0.33*10^1,1.18,'.k');hold on %6496
+text(0.33*10^1,1.18,'6496');
+h1(8)=plot(0.91*10^1,0.42,'.k');hold on %6656
+text(0.91*10^1,0.42,'6656');
+h1(9)=plot(0.55*10^1,0.24,'.k');hold on %6752
+text(0.55*10^1,0.24,'6752');
+h1(10)=plot(0.15*10^1,0.47,'.k');hold on %6838
+text(0.15*10^1,0.47,'6838');
+
+
subplot(1,2,2)
-contour(M_cl/(10^4),r_h,f_M_bh);
+contour(M_cl/(10^5),r_h,f_M_bh);
colorbar
-%clabel(C,h)
title(['BH Mass fraction. Z=' num2str(Z) '. ' num2str(kick)])
-%title(['BH retention fraction. Z=' num2str(Z) '. W/ fallback , sigma 265'])
-xlabel('M_{cl} (10^4 M_{sun})')
+xlabel('M_{cl} (10^5 M_{sun})')
+ylabel('r_h (in pc)')
+%set(gca, 'YScale', 'log')
+
+%Only use this if you want to plot some clusters on top of the cluster
+%fractions
+
+hold on
+
+h1(1)=plot(0.67*10^1,0.59,'.k');hold on %2298
+text(0.67*10^1,0.59,'2298');
+h1(2)=plot(0.9*10^1,0.64,'.k');hold on %6218
+text(0.9*10^1,0.64,'6218');
+h1(3)=plot(0.85*10^1,0.42,'.k');hold on %6254
+text(0.85*10^1,0.42,'6254');
+h1(4)=plot(1*10^1,0.23,'.k');hold on %6341
+text(1*10^1,0.23,'6341');
+h1(5)=plot(0.31*10^1,1.08,'.k');hold on %6352
+text(0.31*10^1,1.08,'6352');
+h1(6)=plot(0.26*10^1,0.17,'.k');hold on %6397
+text(0.26*10^1,0.17,'6397');
+h1(7)=plot(0.33*10^1,1.18,'.k');hold on %6496
+text(0.33*10^1,1.18,'6496');
+h1(8)=plot(0.91*10^1,0.42,'.k');hold on %6656
+text(0.91*10^1,0.42,'6656');
+h1(9)=plot(0.55*10^1,0.24,'.k');hold on %6752
+text(0.55*10^1,0.24,'6752');
+h1(10)=plot(0.15*10^1,0.47,'.k');hold on %6838
+text(0.15*10^1,0.47,'6838');
+
+%{
+%In order to plot the points themselves and not the isolines formed by the
+%contour plot we analyse and plot the matrices row by row:
+for l=1:length(M_cl(1,:))
+figure(k+2)
+subplot(1,2,1)
+scatter3(M_cl(l,:),r_h(l,:),f_N_bh_relative(l,:),'x')
+title(['BH retention fraction. Z=' num2str(Z) '. ' num2str(kick)])
+xlabel('M_{cl} (10^5 M_{sun})')
ylabel('r_h (in pc)')
set(gca, 'YScale', 'log')
-%set(gca, 'XScale', 'log')
-hold off
+colorbar
+hold on
+subplot(1,2,2)
+scatter3(M_cl(l,:),r_h(l,:),f_M_bh(l,:),'x')
+title(['BH Mass fraction. Z=' num2str(Z) '. ' num2str(kick)])
+xlabel('M_{cl} (10^4 M_{sun})')
+ylabel('r_h (in pc)')
+set(gca, 'YScale', 'log')
+colorbar
+hold on
+end
+%}
end
end
diff --git a/results/cluster properties/Cluster properties/fun_for_mass_histograms.m b/results/cluster properties/Cluster properties/fun_for_mass_histograms.m
index 3c4ab6dd8246545b06ce1a35591a263012904135..ed51fbecbcb67b348f1f6b65d164f94b8331ecc1 100644
--- a/results/cluster properties/Cluster properties/fun_for_mass_histograms.m
+++ b/results/cluster properties/Cluster properties/fun_for_mass_histograms.m
@@ -10,7 +10,7 @@ Output:
IMPORTANT: THIS CODE ONLY WORKS FOR A GIVEN MASS AND ESCAPE VELOCITY (OR
HALF-MASS RADIUS BASICALLY)
-For more information on the sampling and how we compute the retained stars see sampling.m
+For more information on the sampling and how we compute the retained stars see SAMPLING.m
%}
function [vec_m_bh_initial,vec_m_bh_final]=fun_for_mass_histograms(M_cl,Z,v_esc,k)
@@ -21,7 +21,7 @@ load('cpdf')
vec_m_initial=[];%We will use this vector to save the different masses. We initialize this variable so that we store every mass.
%And then we can sample the IMF:
- M_tot_initial=0;
+M_tot_initial=0;
jj=1;%We initialize this variable before starting the while loop.
while M_tot_initial<M_cl
r=rand;
diff --git a/results/cluster properties/Cluster properties/fun_for_mass_histograms_all_kicks.m b/results/cluster properties/Cluster properties/fun_for_mass_histograms_all_kicks.m
new file mode 100644
index 0000000000000000000000000000000000000000..c5b6cde4403f05be7e2fbcae6e557260891b67b7
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fun_for_mass_histograms_all_kicks.m
@@ -0,0 +1,151 @@
+%% FUNCTION THAT COMPUTES THE THE MASS HISTOGRAM FOR THE THREE PRESCRIPTIONS AT THE SAME TIME:
+
+function fun_for_mass_histograms_all_kicks(Z,M_cl,r_h)
+v_esc=fun_v_esc(r_h,M_cl);
+%We start by sampling the IMF. And for that we need to load the cpdf
+%(cumulative probability distribution function):
+
+load('cpdf')
+
+vec_m_initial=[];%We will use this vector to save the different masses. We initialize this variable so that we store every mass.
+%And then we can sample the IMF:
+ M_tot_initial=0;
+jj=1;%We initialize this variable before starting the while loop.
+ while M_tot_initial<M_cl
+ r=rand;
+ kk=find(abs(eval_fun-r)<1e-3);
+ m=vec_M_examples(kk(randi(length(kk))));
+ M_tot_initial=M_tot_initial+m;
+ vec_m_initial(jj)=m;% As it can be seen we store in each row the masses of a given sampling (for a given cluster mass)
+ jj=jj+1;
+ end
+N_tot_initial=length(vec_m_initial);%In this vector we save the total number of stars.
+%We now have to define the minimum M_zams for NS and BHs:
+vec_M_zams_test=linspace(7,40,2000);
+vec_M_bh_test=fun_M_rem(Z,vec_M_zams_test);
+ll=find(abs(vec_M_bh_test-3)<1e-2);
+M_z_min_bh_test=vec_M_zams_test(ll(1));
+%And to find the minimum NS mass: The minimum mass of 1.17 solar masses is
+%given by Yudai Suwa, Takashi Yoshida, Masaru Shibata, Hideyuki Umeda & Koh
+%Takahashi in "On the minimum mass of neutron stars" in 2009. The problem
+%is that we cannot get M_rem=1.17 with the function we obtained via the
+%PARSEC paper.That is why we set a minimum value of 1.4:
+ll=find(abs(vec_M_bh_test-1.4)<1e-2);
+M_z_min_ns_test=vec_M_zams_test(ll(1));
+
+
+%Now, we take the values we obtained before and use them to determine which
+%stars are retained and which ones are not:
+vec_m_final=[];%We initialize the vector that will store the retained masses.
+%In order to not use so many if statements we are going to find which
+%positions in the vector vec_m_initial have a mass smaller than the minimum
+%neutron star mass. Because this stars are retained and considering them in
+%an if statement only slows the computation:
+pos_stars=vec_m_initial<M_z_min_ns_test;%Position of such stars in the vector vec_m_initial.
+%(we get 1's in the positions of the stars that verify the requirement, and 0 for the others)
+vec_m_final_i=vec_m_initial.*pos_stars;
+%But we need to eliminate the 0's and only consider those that are
+%different from 0.
+vec_m_final=nonzeros(vec_m_final_i)';
+%And the vector of masses for the NS and BHs before considering the
+%probabilities for them to be kicked out are:
+vec_bh_ns=vec_m_initial-vec_m_final_i;
+vec_bh_ns=nonzeros(vec_bh_ns)';%Vector that contains the initial stars that are (or will be, basically, as these are ZAMS masses) BHs or NSs.
+
+%But we want a vector that contains the initial masses of BHs:
+
+%But we want a vector that contains the initial masses of BHs:
+pos_bh=fun_M_rem(Z,vec_bh_ns)>3;
+vec_m_bh_initial=fun_M_rem(Z,vec_bh_ns).*pos_bh;
+vec_m_bh_initial=nonzeros(vec_m_bh_initial)';
+%And now we do an if statement for this stars:
+ll=length(vec_m_final); %We initialize ll so that we fill the vector after that position.
+%We calculate only once the probability for the initial vector of masses for the different kick prescriptions:
+
+
+
+%And before computing anything else, we make a histogram for these initial
+%masses:
+figure(1)
+subplot(2,2,1)
+h1=histogram(vec_m_bh_initial,35,'FaceColor','k','FaceAlpha',0.6);
+h1.BinWidth = 1.5;
+title(['Number of BHs for M_{cl}=' num2str(M_cl) ' M_{sun}, r_h=' num2str(r_h) 'pc & Z=' num2str(Z) '. '])
+xlabel('M_{BH} (M_{sun})')
+ylabel('Number of BHs')
+xlim([0 100])
+xticks([0 10 20 30 40 50 60 70 80 90 100])
+hold on
+
+
+for k=0:2
+ f_ret(k+1,:)=(1-feval(@prob_escape,Z,vec_bh_ns,1.4,265,v_esc,k));
+end
+%It is also important to see that later, in order to really compute the
+%number of BHs and NSs we need to create a vector that stores the BHs and
+%neutron stars:
+vec_bh_ns_final=[];
+vec_m_final=[];
+
+pos_stars=vec_m_initial<M_z_min_ns_test;%Position of such stars in the vector vec_m_initial.
+%(we get 1's in the positions of the stars that verify the requirement, and 0 for the others)
+vec_m_final_i=vec_m_initial.*pos_stars;
+%But we need to eliminate the 0's and only consider those that are
+%different from 0.
+vec_m_final=nonzeros(vec_m_final_i)';
+%For that we initialize the variable nn
+nn=1;
+%But we need to transform this into the masses that the remnant stars will
+%really have after the explosion(not the ZAMS masses that the BHs or NSs previously had):
+vec_bh_ns=fun_M_rem(Z,vec_bh_ns);
+
+
+%And now we need to study which are retained and which are not for the
+%three different kick prescriptions:
+for k=0:2
+ for jj=1:length(vec_bh_ns)
+ if vec_bh_ns(jj)>3
+ rr=rand;
+ if rr<=f_ret(k+1,jj)
+ vec_m_final(ll)=vec_bh_ns(jj);
+ vec_bh_ns_final(nn)=vec_bh_ns(jj);
+ ll=ll+1;nn=nn+1;
+ else
+ ll=ll;nn=nn;
+ end
+ else
+ rr=rand;
+ if rr<=(1-feval(@prob_esc_ns,265,v_esc))
+ vec_m_final(ll)=vec_bh_ns(jj);
+ vec_bh_ns_final=vec_bh_ns(jj);
+ ll=ll+1;nn=nn+1;
+ else
+ ll=ll;nn=nn;
+ end
+ end
+ end
+ %And now we store the values of the masses of black holes retained
+ %after the kicks:
+ vec_m_bh_final=vec_bh_ns_final(find(vec_bh_ns_final>3));
+ %And now we plot the histogram for the final masses for the different
+ %kick prescriptions on top of the one for the initial masses that we
+ %computed above
+ figure(1)
+ subplot(2,2,k+2)
+ h3=histogram(vec_m_bh_initial,35,'FaceColor','k','FaceAlpha',0.6);
+ h3.BinWidth = 1.5;
+ title(['Number of BHs for M_{cl}=' num2str(M_cl) ' M_{sun}, r_h=' num2str(r_h) 'pc & Z=' num2str(Z) '. '])
+ xlabel('M_{BH} (M_{sun})')
+ ylabel('Number of BHs')
+ xlim([0 100])
+ xticks([0 10 20 30 40 50 60 70 80 90 100])
+ hold on
+ h3=histogram(vec_m_bh_final,35,'FaceAlpha',0.6);
+ h3.BinWidth = 1.5;
+ hold on
+end
+%And then the legend:
+%legend(h1,{'Neutrino-driven NK','W/ fallback & reduced \sigma','Standard fallback-controlled NK','IMF (before kicks)'})
+
+
+
diff --git a/results/cluster properties/Cluster properties/fun_plots_dif_radii.m b/results/cluster properties/Cluster properties/fun_plots_dif_radii.m
new file mode 100644
index 0000000000000000000000000000000000000000..52ed91a81b494bdac5ed488d26fb2732ba7ab864
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fun_plots_dif_radii.m
@@ -0,0 +1,140 @@
+%% FUNCTION THAT COMPUTES THE RETENTION FRACTION AND THE MASS FRACTION FOR DIFFERENT RADII AND A MASS CHOSEN.
+%{
+Inputs:
+ M_cl: Mass of the cluster (in solar masses)
+ Z: Metallicity
+Output:
+ plots the retentrion fraction and mass fraction as functions of the
+ half-mass radius for the cluster mass chosen.
+
+IMPORTANT: THIS CODE ONLY WORKS FOR A SCALAR M_cl and Z.
+For more information on the sampling and how we compute the retained stars see SAMPLING.m
+%}
+
+function f_M_bh=fun_plots_dif_radii(M_cl,Z)
+load('cpdf')
+vec_m_initial=[];%We will use this matrix to save the different masses. We initialize this variable so that we store every mass vector in a different row:
+vec_r_h=linspace(0.1,15,150);
+tic
+%And then we can sample the IMF:
+ M_tot_initial=0;
+ jj=1;%We initialize this variable before starting the while loop.
+ while M_tot_initial<M_cl
+ r=rand;
+ kk=find(abs(eval_fun-r)<1e-3);
+ m=vec_M_examples(kk(randi(length(kk))));
+ M_tot_initial=M_tot_initial+m;
+ vec_m_initial(jj)=m;% As it can be seen we store in each row the masses of a given sampling (for a given cluster mass)
+ jj=jj+1;
+ end
+N_tot_initial=length(vec_m_initial);%In this vector we save the total number of stars.
+%We now have to define the minimum M_zams for NS and BHs:
+vec_M_zams_test=linspace(7,40,2000);
+vec_M_bh_test=fun_M_rem(Z,vec_M_zams_test);
+ll=find(abs(vec_M_bh_test-3)<1e-2);
+M_z_min_bh_test=vec_M_zams_test(ll(1));
+%And to find the minimum NS mass: The minimum mass of 1.17 solar masses is
+%given by Yudai Suwa, Takashi Yoshida, Masaru Shibata, Hideyuki Umeda & Koh
+%Takahashi in "On the minimum mass of neutron stars" in 2009. The problem
+%is that we cannot get M_rem=1.17 with the function we obtained via the
+%PARSEC paper.That is why we set a minimum value of 1.4:
+ll=find(abs(vec_M_bh_test-1.4)<1e-2);
+M_z_min_ns_test=vec_M_zams_test(ll(1));
+
+%And then we compute the total number of BHs (stars that will become one):
+N_bh_initial=sum(vec_m_initial>M_z_min_bh_test);
+toc
+f_N_bh_relative=[];f_M_bh=[];
+for k=0:2
+ if k==0
+ kick='Neutrino-driven NK';
+ elseif k==1
+ kick='Hybrid model';
+ elseif k==2
+ kick='Standard-fallback controlled NK';
+ end
+%Now, we take the values we obtained before and use them to determine which
+%stars are retained and which ones are not:
+
+for dd=1:length(vec_r_h)
+v_esc=fun_v_esc(vec_r_h(dd),M_cl);
+vec_m_final=[];%We initialize the vector that will store the retained masses.
+vec_m_bh_final=[];
+tic
+%In order to not use so many if statements we are going to find which
+%positions in the vector vec_m_initial have a mass smaller than the minimum
+%neutron star mass. Because this stars are retained and considering them in
+%an if statement only slows the computation:
+pos_stars=vec_m_initial<M_z_min_ns_test;%Positioon of such stars in the vector vec_m_initial.
+%(we get 1's in the positions of the stars that verify the requirement, and 0 for the others)
+vec_m_final_i=vec_m_initial.*pos_stars;
+%But we need to eliminate the 0's and only consider those that are
+%different from 0.
+vec_m_final=nonzeros(vec_m_final_i)';
+%And the vector of masses for the NS and BHs before considering the
+%probabilities for them to be kicked out are:
+vec_bh_ns=vec_m_initial-vec_m_final_i;
+vec_bh_ns=nonzeros(vec_bh_ns)';%Vector that contains the initial stars that are (or will be, basically, as these are ZAMS masses) BHs or NSs.
+%And now we do an if statement for this stars:
+ll=length(vec_m_final); %We initialize ll so that we fill the vector after that position.
+%We calculate only once the probability for the initial vector of masses:
+ f_ret=(1-feval(@prob_escape,Z,vec_bh_ns,1.4,265,v_esc,k));
+ f_ret_ns=(1-feval(@prob_esc_ns,265,v_esc));
+ %As it can be seen by default we choose
+ %a velocity dispersion of 265km/s and a mean neutron star mass of 1.4
+ %solar masses.
+%It is also important to see that later, in order to really compute the
+%number of BHs and NSs we need to create a vector that stores the BHs and
+%neutron stars:
+vec_bh_ns_final=[];
+%For that we initialize the variable nn
+nn=1;
+%But we need to transform this into the masses that the remnant stars will
+%really have (not the ZAMS masses that the BHs or NSs previously had):
+vec_bh_ns=fun_M_rem(Z,vec_bh_ns);
+ for jj=1:length(vec_bh_ns)
+ if vec_bh_ns(jj)>3
+ rr=rand;
+ if rr<=f_ret(jj)
+ vec_m_final(ll)=vec_bh_ns(jj);
+ vec_bh_ns_final(nn)=vec_bh_ns(jj);
+ ll=ll+1;nn=nn+1;
+ else
+ ll=ll;nn=nn;
+ end
+ else
+ rr=rand;
+ if rr<=f_ret_ns
+ vec_m_final(ll)=vec_bh_ns(jj);
+ vec_bh_ns_final(nn)=vec_bh_ns(jj);
+ ll=ll+1;nn=nn+1;
+ else
+ ll=ll;nn=nn;
+ end
+ end
+ end
+ %And now we store the values for the fraction of BHs retained in the
+ %cluster after considering the supernova kicks
+ N_bh_final(dd)=sum(vec_bh_ns_final>3);
+ f_N_bh_relative(dd,k+1)=N_bh_final(dd)./N_bh_initial;
+ %And to obtain the mass fractions:
+ vec_m_bh_final=vec_bh_ns_final(find(vec_bh_ns_final>3));
+ f_M_bh(dd,k+1)=sum(vec_m_bh_final)/(M_cl);
+ toc
+end
+ figure(1)
+ h1(k+1)=plot(vec_r_h,f_N_bh_relative(:,k+1));
+ title(['NUMBER FRACTION OF BH RETAINED. M_{cl} ' num2str(M_cl) ' M_{sun} & Z=' num2str(Z)])
+ xlabel('r_{h} (pc)')
+ ylabel('Retention fraction')
+ hold on
+ figure(2)
+ h2(k+1)=plot(vec_r_h,f_M_bh(:,k+1));
+ title(['MASS FRACTION OF BH RETAINED. M_{cl} ' num2str(M_cl) ' M_{sun} & Z=' num2str(Z)])
+ xlabel('r_{h} (pc)')
+ ylabel('Mass fraction')
+ hold on
+end
+legend(h1,{'Neutrino-driven NK','Hybrid model','Standard fallback-controlled NK'});
+legend(h2,{'Neutrino-driven NK','Hybrid model','Standard fallback-controlled NK'});
+end
diff --git a/results/cluster properties/Cluster properties/fun_plots_sig_bh.m b/results/cluster properties/Cluster properties/fun_plots_sig_bh.m
new file mode 100644
index 0000000000000000000000000000000000000000..c3bb7588592a2d95bc24ceedfebdc32253013dd0
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fun_plots_sig_bh.m
@@ -0,0 +1,46 @@
+%% FUNCTION TO OBTAIN THE PLOTS FOR THE VELOCITY DISPERSION FOR DIFFERENT KICKS
+%% FUNCTION THAT YIELDS THE VELOCITY DISPERSION FOR THE DIFFERENT KICK PRESCRIPTIONS
+%{
+Description:
+ This function yields the plots for the velocity dispersion for the different
+ kick prescriptions.
+%}
+%Input:
+% Z: Metallicity
+% M_ns: neutron star mass (in solar masses)
+% sig_ns: neutron star velocity dispersion (in km/s)
+% k: parameter that defines the kick prescription:
+% k=0 ==> neutrino-driven NK
+% k=1 ==> hybrid model
+% K=2 ==> standard kick
+
+%Output:
+% plots for the different velocity dispersions according to
+% the kick
+
+function sig_bh=fun_plots_sig_bh(Z,M_ns,sig_ns)
+vec_M_zams_test=linspace(7,150,50000);
+vec_M_bh_test=fun_M_rem(Z,vec_M_zams_test);
+ll=find(abs(vec_M_bh_test-3)<1e-2);
+M_zams=vec_M_zams_test(ll(1):end);
+M_bh=fun_M_rem(Z,M_zams);
+sig_bh=[];
+ for k=0:2
+ sig_bh(k+1,:)=fun_sigma_bh(Z,M_ns,sig_ns,k);
+ end
+ figure(1)
+ h1(1)=plot(M_bh,sig_bh(1,:));
+ title(['VELOCITY DISPERSION. <M_{NS}>= ' num2str(M_ns) ' M_{sun}, \sigma_{NS}=' num2str(sig_ns) ' km/s & Z=' num2str(Z)])
+ xlabel('M_{BH} (M_{sun})')
+ ylabel('\sigma_{BH} (km/s)')
+ xlim([3 max(M_bh)])
+ hold on
+ h1(2)=plot(M_bh,sig_bh(2,:));
+ hold on
+ h1(3)=plot(M_bh,sig_bh(3,:));
+ legend(h1,{'Neutrino-driven NK','Hybrid model','Standard fallback-controlled NK'})
+end
+
+
+
+
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/fun_sigma_bh.m b/results/cluster properties/Cluster properties/fun_sigma_bh.m
new file mode 100644
index 0000000000000000000000000000000000000000..f4b7bcc89204a02e11187b9609a704a044a953bd
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fun_sigma_bh.m
@@ -0,0 +1,45 @@
+%% FUNCTION THAT YIELDS THE VELOCITY DISPERSION FOR THE DIFFERENT KICK PRESCRIPTIONS
+%{
+Description:
+ This function yields the velocity dispersion for the different
+ kick prescriptions.
+%}
+%Input:
+% Z: Metallicity
+% M_ns: neutron star mass (in solar masses)
+% sig_ns: neutron star velocity dispersion (in km/s)
+% k: parameter that defines the kick prescription:
+% k=0 ==> neutrino-driven NK
+% k=1 ==> hybrid model
+% K=2 ==> standard kick
+
+%Output:
+% sig_bh: the different velocity dispersions according to
+% the kick
+
+function sig_bh=fun_sigma_bh(Z,M_ns,sig_ns,k)
+%In order to compute the BH mass we make use of fun_M_rem.m that gives the
+%mass of the stellar remnant as a function of Z and M.
+%We start by determining which is the minimum ZAMS mass that yields a BH.
+%It is important to do it this way and not just generating the BH masses
+%directly because this way we can also now which is the maximum BH mass
+%that can be formed according to the fitting formulas for the stellar
+%evolution.
+vec_M_zams_test=linspace(7,150,50000);
+vec_M_bh_test=fun_M_rem(Z,vec_M_zams_test);
+ll=find(abs(vec_M_bh_test-3)<1e-2);
+M_zams=vec_M_zams_test(ll(1):end);
+M_bh=fun_M_rem(Z,M_zams);
+%And then we can obtain the BH's velocity dispersion and the fallback
+%depending on the prescription:
+ if k==1
+ f_fb=fun_fb(Z,M_zams);
+ sig_bh=(M_ns./M_bh).*(1-f_fb).*sig_ns;
+ elseif k==0
+ sig_bh=(M_ns./M_bh).*sig_ns;
+ f_fb=0;
+ elseif k==2
+ f_fb=fun_fb(Z,M_zams);
+ sig_bh=(1-f_fb).*sig_ns;
+ end
+end
diff --git a/results/cluster properties/Cluster properties/SAMPLING.m~ b/results/cluster properties/Cluster properties/fun_stochasticy_all_kicks.m
similarity index 50%
rename from results/cluster properties/Cluster properties/SAMPLING.m~
rename to results/cluster properties/Cluster properties/fun_stochasticy_all_kicks.m
index 350e60b3806c56adcc03b4ac8dab30ad7a72f9b4..d444dd80073760f25e6abc4cd6cbb5a0b2c3d1e4 100644
--- a/results/cluster properties/Cluster properties/SAMPLING.m~
+++ b/results/cluster properties/Cluster properties/fun_stochasticy_all_kicks.m
@@ -1,112 +1,23 @@
-%% SCRIPT WITH EVERYTHING REGARDING THE RESULTS OBTAINED VIA SAMPLING
+%% FUNCTION THAT STUDIES THE STOCHASTICITY FOR THE THREE KICK PRESCRIPTIONS WE HAVE WORKED ON
%{
-Description:
- This function does different things. First of all it samples a
- cluster with the specifications chosen and then plots a histogram for the
- different BH masses both before and after supernova kicks. Then it also
- does 100 different samplings (vec_ii can be changed if one wants a
- different number of samplings) in order to show how a different sampling
- affects the results obtained.
- Then it also plots a contour plot for the chosen kick prescription
- showing the BH retention fraction for different cluster masses and
- half-mass radii.
-%}
-%Input:
-% Z: Metallicity
-% M_cl: Mass of the cluster(in solar masses)
-% v_esc: escape velocity (in km/s)
-% k: parameter that defines the kick prescription:
-% k=0 ==> no fallback and reduced velocity dispersion
-% k=1 ==> fallback considered and reduced velocity
-% dispersion
-% K=2 ==> fallback considered and not a reduced
-% velocity dispersion.
-
-%Output:
-% plots explained in the description of the function (see
-% above)
-%function [f_N_bh_relative,M_cl,r_h]=SAMPLING(Z,M_cl,r_h,k)
-function SAMPLING(Z,M_cl,r_h,k)
-if k==0
- kick='W/o fallback,reduced sigma';
-elseif k==1
- kick='W/ fallback,reduced sigma';
-elseif k==2
- kick='W/ fallback,sigma not reduced';
-end
-%{
-We do the sampling by calculating the cumulative probability distribution
-function (cpdf) via the IMF, and then for each star we want to "create" in the
-cluster we first asign a random value r between 0 and 1 and then find
-which mass yields a value of the cpdf closest to that r. In this way we
-generate stars until we reach the mass of the cluster that we are looking
-for.
-As this is a costly process we only computed that once, and then we simply
-load the results that are in the file cpdf.mat saved as eval_fun. (we load
-it by using the command load('cpdf')
-In any case the code is below:
-Below we show what that process is:
-
-%We start by defining the IMF:
-alfa1=0.3;alfa2=1.3;alfa3=2.3;
-imf= @(M) (M<0.08).*(M.^(-alfa1)) + (M<0.5 & M>=0.08).*((M.^(-alfa2))*(0.08)) + (M>=0.5).*((M.^(-alfa3))*(0.5));
-%We then define the maximum value for the cumulative probability density
-%function so that we can normalise the function itself so it can be treated
-%as a probability.
-cpdf_max=integral(imf,0.07,150,'RelTol',1e-6,'AbsTol',1e-6);
-%Once we have this function defined we need to invert it. That it, we need
-%to have it such that we give a random number from 0 to 1 and it returns a
-%mass. So that we asign different masses until we reach the maximum value
-%we want it to have.
-%For example:
-%We initialize the total mass to be 0 before we start "creating" stars:
-M_tot=0;
-vec_M_examples=[linspace(0.07,2,15000) linspace(2,150,25000)];
-eval_fun=[];
-for ii=1:length(vec_M_examples)
- M_example=vec_M_examples(ii);
- fun=@(M) (integral(imf,0.07,M,'RelTol',1e-6,'AbsTol',1e-6)/(cpdf_max));
- eval_fun(ii)=feval(fun,M_example);
-end
+Inputs:
+ M_cl: Mass of the cluster (in solar masses)
+ Z: Metallicity
+ r_h: half-mass radius
+Output:
+ plots the retention fraction for a cluster chosen with the inputs
+ for 100 different samplings (to change this number change vec_ii)
+
+IMPORTANT: THIS CODE ONLY WORKS FOR A SCALAR M_cl and Z.
+For more information on the sampling and how we compute the retained stars see SAMPLING.m
%}
+function fun_stochasticy_all_kicks(Z,M_cl,r_h)
+%We start by loading the cumulative probability distribution function
load('cpdf')
-
-%{
%In order to do the sampling and then compute how many BHs remain in the
%cluster after the supernova kicks we will use fun_for_mass_histograms.m
%For that we first need to compute the escape velocity:
v_esc=fun_v_esc(r_h,M_cl);
-%And then:
-[vec_m_bh_initial,vec_m_bh_final]=fun_for_mass_histograms(M_cl,Z,v_esc,k);
-figure(1)
-h1=histogram(vec_m_bh_initial,35,'FaceColor','b');
-h1.BinWidth = 1.5;
-title(['Number of BHs for M_{cl}=' num2str(M_cl) ' M_{sun}, r_h=' num2str(r_h) 'pc & Z=' num2str(Z) '. ' num2str(kick)])
-xlabel('M_{BH} (M_{sun})')
-ylabel('Number of BHs')
-xlim([0 100])
-%xticks([0 10 20 30 40 50 60 70 80 90 100])
-%set(gca,'YScale','log')
-hold on
-h2=histogram(vec_m_bh_final,35,'FaceColor','r');
-h2.BinWidth = 1.5;
-legend('IMF (before kicks)','after kicks')
-%And now we show the mean mass of BHs before and after the SN kicks:
-disp('FOR THE SECTION ON THE MASS HISTOGRAM:')
-disp(' MEAN MASS OF BHs BEFORE SN KICKS')
-mean_m_bh_initial=sum(vec_m_bh_initial)/length(vec_m_bh_initial);
-disp(mean_m_bh_initial)
-disp(' MEAN MASS OF BHs AFTER SN KICKS')
-mean_m_bh_final=sum(vec_m_bh_initial)/length(vec_m_bh_final);
-disp(mean_m_bh_final)
-
-
-
-%It is also important to study the stochasticity of the sampling method.
-%That is why we compute the number stars, BHs and NSs, as well as some
-%important fractions for different samplings in order to study these
-%variations.
-
vec_ii=[1:100];
mat_m=[];
@@ -149,9 +60,17 @@ M_bh(ii)=sum(feval(@fun_M_rem,Z,vec_m(find(vec_m>M_z_min_bh_test))));
M_stars(ii)=sum(vec_m(find(vec_m<M_z_min_ns_test)));
M_ns(ii)=M_tot(ii)-M_bh(ii)-M_stars(ii);
end
-
%Now that we have sampled it different times we have to compute the results
-%after the supernova kicks according to the prescription chosen:
+%after the supernova kicks according to the three prescriptions we have
+%studied:
+for k=0:2
+ if k==0
+ kick='Neutrino-driven NK.';
+ elseif k==1
+ kick='W/ fallback,reduced sigma NK';
+ elseif k==2
+ kick='Standard falback-controlled NK';
+ end
for ii=1:length(vec_ii)
%In order to not use so many if statements we are going to find which
%positions in the vector vec_m_initial have a mass smaller than the minimum
@@ -211,34 +130,49 @@ vec_bh_ns=fun_M_rem(Z,vec_bh_ns);
M_bh_final(ii)=sum(vec_bh_ns_final(find(vec_bh_ns_final>3)));
M_ns_final(ii)=sum(vec_bh_ns_final(find(vec_bh_ns_final<3)));
M_stars_final(ii)=M_tot_final(ii)-M_bh_final(ii)-M_ns_final(ii);
- f_m_bh(ii)=M_bh_final(ii)/M_bh(ii);f_m_ns(ii)=M_ns_final(ii)/M_ns(ii);f_m_stars=M_stars_final(ii)/M_tot(ii);
+ f_m_bh(ii)=M_bh_final(ii)/M_tot(ii);f_m_ns(ii)=M_ns_final(ii)/M_ns(ii);f_m_stars=M_stars_final(ii)/M_tot(ii);
+ f_m_bh_current(ii)=M_bh_final(ii)/M_tot_final(ii);
end
%And now we can plot this:
-figure(2)
-plot(vec_ii,N_tot,'-k')
+figure(1)
+
+h1(k+1)=plot(vec_ii,N_tot_final);
title(['Number of stars in a cluster of M_{cl}=' num2str(M_cl) ' M_{sun}, r_h=' num2str(r_h) ' pc & Z=' num2str(Z)])
xlabel('Different samplings')
ylabel('Number of stars')
+%legend(['After kicks.' num2str(kick)],'Location','southwest')
+%legend('boxoff')
hold on
-plot(vec_ii,N_tot_final,'-r')
-legend('IMF (before kicks)',['after kicks,' num2str(kick)],'Location','southwest')
-legend('boxoff')
-figure(3)
-subplot(2,1,1)
-plot(vec_ii,N_bh,'-r')
+figure(2)
+subplot(4,1,1)
+h2(k+1)=plot(vec_ii,N_bh_final);
title(['Number of BHs in a cluster of M_{cl}=' num2str(M_cl) ' M_{sun}, r_h=' num2str(r_h) ' pc & Z=' num2str(Z)])
xlabel('Different samplings')
ylabel('Number of BHs')
+%legend(['After kicks.' num2str(kick)],'Location','southwest')
+%legend('boxoff')
hold on
-plot(vec_ii,N_bh_final,'-b')
-legend('IMF (before kicks)',['after kicks,' num2str(kick)],'Location','southwest')
-legend('boxoff')
-subplot(2,1,2)
-plot(vec_ii,f_N_bh,'-k')
+subplot(4,1,2)
+h3(k+1)=plot(vec_ii,f_N_bh);
title('Fraction of BHs retained')
xlabel('Different samplings')
ylabel('Fraction of BHs retained')
+hold on
+subplot(4,1,3)
+h4(k+1)=plot(vec_ii,f_m_bh);
+title('Mass fraction of BHs after NK')
+xlabel('Different samplings')
+ylabel('Mass fraction of BHs')
+hold on
+subplot(4,1,4)
+h5(k+1)=plot(vec_ii,f_m_bh_current);
+title('Fraction of BHs retained')
+xlabel('Different samplings')
+ylabel('Fraction of BHs retained')
+hold on
+
+
%And for these calculations we can study the stochastic nature, that is,
%what is the mean value, standard deviations, etc:
@@ -254,15 +188,16 @@ disp(sigma_N_bh_initial)
%After supernova kicks:
mean_N_bh_final=sum(N_bh_final)/length(vec_ii);
-disp('Mean number of BHs after supernova kicks:')
+disp(['Mean number of BHs after supernova kicks:' ' ( ' num2str(kick) ' )'])
disp(mean_N_bh_final)
var=(N_bh_final-mean_N_bh_final).^2;
sigma_N_bh_final=sqrt(sum(var)/length(vec_ii));
-disp('Standard deviation in the number of BHs after supernova kicks:')
+disp(['Standard deviation in the number of BHs after supernova kicks:' ' ( ' num2str(kick) ' )'] )
disp(sigma_N_bh_final)
%And now we do the same but for the number fraction of black holes:
disp(' VARIATIONS IN THE NUMBER FRACTION OF BHs RETAINED')
+disp([' ' num2str(kick)])
disp('Mean number fraction of BHs retained:')
mean_f_N_bh=sum(f_N_bh)/length(vec_ii);
disp(mean_f_N_bh)
@@ -270,44 +205,16 @@ var=(f_N_bh-mean_f_N_bh).^2;
sigma_f_N_bh=sqrt(sum(var)/length(vec_ii));
disp('Standard deviation in the number fraction of BHs retained:')
disp(sigma_f_N_bh)
+end
+%And we now need to add to the plots above the number of BHs that you find
+%initially (due simply to the IMF)
+figure(1)
+h1(4)=plot(vec_ii,N_tot,'-k');
-
-%And now it would be interesting to do a contour plot for different cluster
-%masses and half-mass radii. To do that we use the function
-%fun_for_contour_fractions_3.m
-
-%We define the following vectors:
-vec_M_cl=linspace(10^3,10^6,100);
-vec_r_h=linspace(0.1,25,100);
-%vec_v_esc=linspace(5,60,100);
-%And with this we form the meshgrid:
-[M_cl,r_h]=meshgrid(vec_M_cl,vec_r_h);
-mat_v_esc=fun_v_esc(r_h,M_cl);
-%So that:
-f_N_bh_relative=fun_for_contour_fractions_3(M_cl,Z,r_h,k);
-%But we want to express it as a function of r_h and nor mat_v_esc. So we
-%basically use the function defined in fun_v_esc.m but rewritten so that it
-%yields the half-mass radius.
-%r_h=(1/40)*((3/(8*pi))^(1/3))*((M_cl)./((mat_v_esc).^2));
-figure(4)
-contour(M_cl/(10^5),r_h,f_N_bh_relative);
-colorbar
-%clabel(C,h)
-%title(['BH retention fraction. Z=' num2str(Z) '. ' num2str(kick)])
-title(['BH retention fraction. Z=' num2str(Z) '. W/ fallback , sigma 265'])
-xlabel('M_{cl} (10^5 M_{sun})')
-ylabel('r_h (in pc)')
-set(gca, 'YScale', 'log')
-%set(gca, 'XScale', 'log')
-%}
-
-%vec_M_cl=linspace(10^4,10^6,100);
-vec_M_cl=logspace(4,6,100);%100 logarithmically spaced numbers.
-vec_r_h=linspace(0.1,15,100);
-%vec_v_esc=linspace(5,60,100);
-%And with this we form the meshgrid:
-[M_cl,r_h]=meshgrid(vec_M_cl,vec_r_h);
-mat_v_esc=fun_v_esc(r_h,M_cl);
-fun_for_contour_fractions_all_plots(M_cl,Z,r_h)
-
+figure(2)
+subplot(4,1,1)
+h2(4)=plot(vec_ii,N_bh,'-k');
+%And now we define the legends for the plots defined above:
+legend(h1,{'Neutrino-driven NK','Hybrid model','Standard fallback-controlled NK','IMF (before kicks)'})
+legend(h2,{'Neutrino-driven NK','Hybrid model','Standard fallback-controlled NK','IMF (before kicks)'})
end
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/histogram_hybrid.fig b/results/cluster properties/Cluster properties/histogram_hybrid.fig
new file mode 100644
index 0000000000000000000000000000000000000000..a3ac9d87128c63b457e8a84306754223c7a1ecfa
Binary files /dev/null and b/results/cluster properties/Cluster properties/histogram_hybrid.fig differ
diff --git a/results/cluster properties/Cluster properties/histogram_hybrid.png b/results/cluster properties/Cluster properties/histogram_hybrid.png
new file mode 100644
index 0000000000000000000000000000000000000000..261cc727373672468152896e91a04879f6af1271
Binary files /dev/null and b/results/cluster properties/Cluster properties/histogram_hybrid.png differ
diff --git a/results/cluster properties/Cluster properties/histogram_hybrid_2.fig b/results/cluster properties/Cluster properties/histogram_hybrid_2.fig
new file mode 100644
index 0000000000000000000000000000000000000000..1bd4559f59e76993609f47541b3ed14b8e470c50
Binary files /dev/null and b/results/cluster properties/Cluster properties/histogram_hybrid_2.fig differ
diff --git a/results/cluster properties/Cluster properties/histogram_hybrid_2.png b/results/cluster properties/Cluster properties/histogram_hybrid_2.png
new file mode 100644
index 0000000000000000000000000000000000000000..bfd5bffcf817a279feb69a291c054e50efc0d574
Binary files /dev/null and b/results/cluster properties/Cluster properties/histogram_hybrid_2.png differ
diff --git a/results/cluster properties/Cluster properties/histogram_neutrino.png b/results/cluster properties/Cluster properties/histogram_neutrino.png
new file mode 100644
index 0000000000000000000000000000000000000000..6ea76e81f4ce185d1fa7be8537899f09a7222a42
Binary files /dev/null and b/results/cluster properties/Cluster properties/histogram_neutrino.png differ
diff --git a/results/cluster properties/Cluster properties/histogram_standard.fig b/results/cluster properties/Cluster properties/histogram_standard.fig
new file mode 100644
index 0000000000000000000000000000000000000000..b2e5fd9f5f4ea5c33efdccd7a4dacab6beaad892
Binary files /dev/null and b/results/cluster properties/Cluster properties/histogram_standard.fig differ
diff --git a/results/cluster properties/Cluster properties/histogram_standard.png b/results/cluster properties/Cluster properties/histogram_standard.png
new file mode 100644
index 0000000000000000000000000000000000000000..4ffb186e00e37c4544ea6730c5a2d42133ea2868
Binary files /dev/null and b/results/cluster properties/Cluster properties/histogram_standard.png differ
diff --git a/results/cluster properties/Cluster properties/plots_M_rem_and_probs_report.m b/results/cluster properties/Cluster properties/plots_M_rem_and_probs_report.m
index 9d59f52dc5174b2ad7e4b5374803ad87d7aa8007..9445ae43db8369de291f9b7ab88338fb3e22bfde 100644
--- a/results/cluster properties/Cluster properties/plots_M_rem_and_probs_report.m
+++ b/results/cluster properties/Cluster properties/plots_M_rem_and_probs_report.m
@@ -30,18 +30,23 @@ elseif k==1
elseif k==2
kick='W/fallback,sigma not reduced';
end
-
+%{
%We can also plot some information regarding M_rem, M_co, the probabilities
%of escape and the fallback fraction in one plot by using subplots:
vec_M_zams=linspace(0.07,150,1000);
vec_M_co=zeros(1,length(vec_M_zams));vec_M_rem=zeros(1,length(vec_M_zams));
-vec_f_fb=zeros(1,length(vec_M_zams));vec_M_zams=linspace(10,100,2000);
-for ii=1:length(vec_M_zams)
-M_zams=vec_M_zams(ii);
-vec_M_co(ii)=fun_M_co(Z,M_zams);vec_M_rem(ii)=fun_M_rem(Z,M_zams);
-vec_f_fb(ii)=fun_fb(Z,M_zams);
-end
-figure(3)
+vec_f_fb=zeros(1,length(vec_M_zams));vec_M_zams=linspace(10,150,2000);
+
+vec_Z=[0.0001 0.001 0.005 0.01];
+for kk=1:length(vec_Z)
+ Z=vec_Z(kk);
+ for ii=1:length(vec_M_zams)
+ M_zams=vec_M_zams(ii);
+ vec_M_co(ii)=fun_M_co(Z,M_zams);vec_M_rem(ii)=fun_M_rem(Z,M_zams);
+ vec_f_fb(ii)=fun_fb(Z,M_zams);
+ end
+figure(1)
+subplot(2,2,kk)
plot(vec_M_zams,vec_M_co,'-k')
title(['M_{rem} & M_{co}. Z=' num2str(Z)])
xlabel('M_{zams} (in M_{sun})')
@@ -62,9 +67,18 @@ line([M_zams_0_5 M_zams_0_5], get(gca, 'ylim'))
line([M_zams_0_75 M_zams_0_75], get(gca, 'ylim'))
line([M_zams_1 M_zams_1], get(gca, 'ylim'))
legend('M_{co}','M_{rem}','f_{fb}=0.5','f_{fb}=0.75','f_{fb}=1')
-
+end
+%}
+%{
%And now we look to plot the probabilities of escape for the three
%different kick prescriptions we have computed:
+
+vec_M_examples=linspace(10,150,1000);
+vec_M_rem_examples=fun_M_rem(Z,vec_M_examples);
+ll=find(abs(vec_M_bh_test-3)<1e-2);
+M_z_min_bh_test=vec_M_zams_test(ll(1))
+vec_M_zams=vec_M_examples(vec_M_examples>M_z_min_bh);
+vec_M_rem=fun_M_rem(Z,vec_M_zams);
for k=0:2
if k==0
kick='W/o fallback,reduced sigma';
@@ -73,8 +87,173 @@ for k=0:2
elseif k==2
kick='W/fallback,sigma not reduced';
end
- p=prob_escape(Z,M_rem,M_ns,sig_ns,v_esc,k)
+ probescape(k+1,:)=prob_escape(Z,M_zams,M_ns,sig_ns,v_esc,k);%We store it by rows
+end
+figure(2)
+plot(vec_M_rem,probescape(1,:),'-k')
+title('Probability
+%}
+%{
+%Now we are interested in plotting the fallback fraction according to the
+%BH mass:
+vec_M_zams=linspace(10,150,2000);
+vec_Z=[0.0001 0.001 0.005 0.01];
+for kk=1:length(vec_Z)
+ Z=vec_Z(kk);
+ for ii=1:length(vec_M_zams)
+ M_zams=vec_M_zams(ii);
+ vec_M_co(ii)=fun_M_co(Z,M_zams);vec_M_rem(ii)=fun_M_rem(Z,M_zams);
+ vec_f_fb(ii)=fun_fb(Z,M_zams);
+ end
+figure(3)
+h1(kk)=plot(vec_M_rem,vec_f_fb);
+title('Fallback fraction for different Z')
+xlabel('M_{rem} (M_{sun})')
+ylabel('Fallback fraction')
+xlim([0 40])
+hold on
+end
+%And now we define the legend:
+legend(h1,{['Z=' num2str(vec_Z(1))],['Z=' num2str(vec_Z(2))],['Z=' num2str(vec_Z(3))],['Z=' num2str(vec_Z(4))]})
+%}
+
+
+
+%And the last thing to do is to plot the probabilities themselves. As the
+%behaviour is difficult to capture in a single plot these will not appear
+%in the main work, but we will rathe put them in an Appendix with
+%additional plots.
+
+%We start by plotting the probability of escape as a function of BH mass
+%for the three different mechanisms for 3 different escape velocities of
+%10, 15 & 20 km/s for a metallicity of Z=0.0001.
+v_esc=10; Z=0.0001;
+vec_M_examples=linspace(10,150,10000);
+vec_M_rem_examples=fun_M_rem(Z,vec_M_examples);
+ll=find(abs(vec_M_rem_examples-3)<1e-2);
+M_z_min_bh_test=vec_M_examples(ll(1));
+vec_M_zams=vec_M_examples(vec_M_examples>M_z_min_bh_test);
+vec_M_rem=fun_M_rem(Z,vec_M_zams);
+prob_dif_vesc=[];
+for k=0:2
+ prob_dif_vesc(k+1,:)=prob_escape(Z,vec_M_zams,1.4,265,v_esc,k);
+end
+figure(1)
+subplot(3,1,1)
+h1(1)=plot(vec_M_rem,prob_dif_vesc(1,:));
+title(['Probability of escape. Z=' num2str(Z)])
+xlabel('M_{rem} (M_{sun})')
+ylabel('Probability of escape')
+xlim([3 100])
+hold on
+h1(2)=plot(vec_M_rem,prob_dif_vesc(2,:));
+hold on
+h1(3)=plot(vec_M_rem,prob_dif_vesc(3,:));
+hold on
+legend(h1,{'Neutrino-driven NK','Hybrid model','Standard fallback-controlled NK'})
+%And now for a different escape velocity:
+v_esc=15;
+prob_dif_vesc=[];
+for k=0:2
+ prob_dif_vesc(k+1,:)=prob_escape(Z,vec_M_zams,1.4,265,v_esc,k);
+end
+figure(1)
+subplot(3,1,2)
+plot(vec_M_rem,prob_dif_vesc(1,:))
+title(['Probability of escape. Z=' num2str(Z)])
+xlabel('M_{rem} (M_{sun})')
+ylabel('Probability of escape')
+xlim([3 100])
+hold on
+plot(vec_M_rem,prob_dif_vesc(2,:))
+hold on
+plot(vec_M_rem,prob_dif_vesc(3,:))
+hold on
+
+%And for another one:
+v_esc=20;
+prob_dif_vesc=[];
+for k=0:2
+ prob_dif_vesc(k+1,:)=prob_escape(Z,vec_M_zams,1.4,265,v_esc,k);
+end
+figure(1)
+subplot(3,1,3)
+plot(vec_M_rem,prob_dif_vesc(1,:))
+title(['Probability of escape. Z=' num2str(Z)])
+xlabel('M_{rem} (M_{sun})')
+ylabel('Probability of escape')
+xlim([3 100])
+hold on
+plot(vec_M_rem,prob_dif_vesc(2,:))
+hold on
+plot(vec_M_rem,prob_dif_vesc(3,:))
+hold on
+
+
+%And now would be interesting to plot this probability for a different
+%metallicity but for the same three velocities:
+
+v_esc=10; Z=0.001;
+vec_M_examples=linspace(10,150,10000);
+vec_M_rem_examples=fun_M_rem(Z,vec_M_examples);
+ll=find(abs(vec_M_rem_examples-3)<1e-2);
+M_z_min_bh_test=vec_M_examples(ll(1));
+vec_M_zams=vec_M_examples(vec_M_examples>M_z_min_bh_test);
+vec_M_rem=fun_M_rem(Z,vec_M_zams);
+prob_dif_vesc=[];
+for k=0:2
+ prob_dif_vesc(k+1,:)=prob_escape(Z,vec_M_zams,1.4,265,v_esc,k);
+end
+figure(2)
+subplot(3,1,1)
+h2(1)=plot(vec_M_rem,prob_dif_vesc(1,:));
+title(['Probability of escape. Z=' num2str(Z)])
+xlabel('M_{rem} (M_{sun})')
+ylabel('Probability of escape')
+xlim([3 100])
+hold on
+h2(2)=plot(vec_M_rem,prob_dif_vesc(2,:));
+hold on
+h2(3)=plot(vec_M_rem,prob_dif_vesc(3,:));
+hold on
+legend(h2,{'Neutrino-driven NK','Hybrid model','Standard fallback-controlled NK'})
+%And now for a different escape velocity:
+v_esc=15;
+prob_dif_vesc=[];
+for k=0:2
+ prob_dif_vesc(k+1,:)=prob_escape(Z,vec_M_zams,1.4,265,v_esc,k);
+end
+figure(2)
+subplot(3,1,2)
+plot(vec_M_rem,prob_dif_vesc(1,:))
+title(['Probability of escape. Z=' num2str(Z)])
+xlabel('M_{rem} (M_{sun})')
+ylabel('Probability of escape')
+xlim([3 100])
+hold on
+plot(vec_M_rem,prob_dif_vesc(2,:))
+hold on
+plot(vec_M_rem,prob_dif_vesc(3,:))
+hold on
+%And for another one:
+v_esc=20;
+prob_dif_vesc=[];
+for k=0:2
+ prob_dif_vesc(k+1,:)=prob_escape(Z,vec_M_zams,1.4,265,v_esc,k);
+end
+figure(2)
+subplot(3,1,3)
+plot(vec_M_rem,prob_dif_vesc(1,:))
+title(['Probability of escape. Z=' num2str(Z)])
+xlabel('M_{rem} (M_{sun})')
+ylabel('Probability of escape')
+xlim([3 100])
+hold on
+plot(vec_M_rem,prob_dif_vesc(2,:))
+hold on
+plot(vec_M_rem,prob_dif_vesc(3,:))
+hold on
-end
\ No newline at end of file
+end
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/NEUTRINO_DRIVEN.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/NEUTRINO_DRIVEN.fig
new file mode 100644
index 0000000000000000000000000000000000000000..7a06bcfc14ba467bb2c55c3fd5b6dc66310027a9
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/NEUTRINO_DRIVEN.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/NEUTRINO_DRIVEN.png b/results/cluster properties/Cluster properties/plots_after_big_changes/NEUTRINO_DRIVEN.png
new file mode 100644
index 0000000000000000000000000000000000000000..aa55fba64fafa06288bc81b2f11a2d09ad58a94b
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/NEUTRINO_DRIVEN.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/NEUTRINO_DRIVEN_low_radii.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/NEUTRINO_DRIVEN_low_radii.fig
new file mode 100644
index 0000000000000000000000000000000000000000..87f2f0ff2e937994949f176a7c7c63ce8ce5891f
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/NEUTRINO_DRIVEN_low_radii.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/NEUTRINO_DRIVEN_low_radii.png b/results/cluster properties/Cluster properties/plots_after_big_changes/NEUTRINO_DRIVEN_low_radii.png
new file mode 100644
index 0000000000000000000000000000000000000000..75fc7741da37ad17b46e70bb6301dd516775546d
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/NEUTRINO_DRIVEN_low_radii.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/NEUTRINO_DRIVEN_low_radii_2.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/NEUTRINO_DRIVEN_low_radii_2.fig
new file mode 100644
index 0000000000000000000000000000000000000000..f8c44183fcf48b71eb7f871134860adcd6ffc8b6
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/NEUTRINO_DRIVEN_low_radii_2.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/STANDARD.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/STANDARD.fig
new file mode 100644
index 0000000000000000000000000000000000000000..a2c92067ed51a6c60478bec418a01dae298f8b22
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/STANDARD.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/STANDARD.png b/results/cluster properties/Cluster properties/plots_after_big_changes/STANDARD.png
new file mode 100644
index 0000000000000000000000000000000000000000..9f543dde94b85567452f07dd27649eea79b71bff
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/STANDARD.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/STANDARD_low_radii.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/STANDARD_low_radii.fig
new file mode 100644
index 0000000000000000000000000000000000000000..579e9e02d00c51bbcbd88fe62577ef02af959552
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/STANDARD_low_radii.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/STANDARD_low_radii.png b/results/cluster properties/Cluster properties/plots_after_big_changes/STANDARD_low_radii.png
new file mode 100644
index 0000000000000000000000000000000000000000..d9da9daa70fe485c41cae22fe1e84c8d8a0b98ef
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/STANDARD_low_radii.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/STANDARD_low_radii_2.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/STANDARD_low_radii_2.fig
new file mode 100644
index 0000000000000000000000000000000000000000..d8f360f27cfb5581c0e4ad9a9df46eafaef56ad5
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/STANDARD_low_radii_2.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/STOCHASTICITY_RETENTION.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/STOCHASTICITY_RETENTION.fig
new file mode 100644
index 0000000000000000000000000000000000000000..ef2c25be2b9f48165913fd1e120f97cf6201a22c
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/STOCHASTICITY_RETENTION.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/STOCHASTICITY_RETENTION.png b/results/cluster properties/Cluster properties/plots_after_big_changes/STOCHASTICITY_RETENTION.png
new file mode 100644
index 0000000000000000000000000000000000000000..407c271031b6e5cd4f4ade8a9f7a5db3b91a403a
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/STOCHASTICITY_RETENTION.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/STOCHASTICITY_RETENTION_500000.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/STOCHASTICITY_RETENTION_500000.fig
new file mode 100644
index 0000000000000000000000000000000000000000..60682b8adb13a3f3afb5a0df1f0aa304c74e727f
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/STOCHASTICITY_RETENTION_500000.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/W_FB_REDUCED_SIGMA.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/W_FB_REDUCED_SIGMA.fig
new file mode 100644
index 0000000000000000000000000000000000000000..5059798b1a9e3a01894a26be329e5501929168d6
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/W_FB_REDUCED_SIGMA.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/W_FB_REDUCED_SIGMA.png b/results/cluster properties/Cluster properties/plots_after_big_changes/W_FB_REDUCED_SIGMA.png
new file mode 100644
index 0000000000000000000000000000000000000000..87a6d2de825609d114bb67c5a277da2f605163cb
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/W_FB_REDUCED_SIGMA.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/W_FB_REDUCED_SIGMA_low_radii.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/W_FB_REDUCED_SIGMA_low_radii.fig
new file mode 100644
index 0000000000000000000000000000000000000000..1d7de3d0f687935f17512247e3c95ae2ae332c93
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/W_FB_REDUCED_SIGMA_low_radii.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/W_FB_REDUCED_SIGMA_low_radii.png b/results/cluster properties/Cluster properties/plots_after_big_changes/W_FB_REDUCED_SIGMA_low_radii.png
new file mode 100644
index 0000000000000000000000000000000000000000..1a9944c4d9c81ca303756e83f462e57332a22dee
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/W_FB_REDUCED_SIGMA_low_radii.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/W_FB_REDUCED_SIGMA_low_radii_2.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/W_FB_REDUCED_SIGMA_low_radii_2.fig
new file mode 100644
index 0000000000000000000000000000000000000000..c672fd93d5955e6088b6c8d0892f0411c1c673fb
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/W_FB_REDUCED_SIGMA_low_radii_2.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/archive/NEUTRINO_DRIVEN.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/NEUTRINO_DRIVEN.fig
new file mode 100644
index 0000000000000000000000000000000000000000..58061db66c8eae490bd3beca7a4110719189af73
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/NEUTRINO_DRIVEN.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/archive/NEUTRINO_DRIVEN.png b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/NEUTRINO_DRIVEN.png
new file mode 100644
index 0000000000000000000000000000000000000000..9eeab7fe2899a009dcfedd17fce35fdb356e5f0a
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/NEUTRINO_DRIVEN.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/archive/NEUTRINO_DRIVEN_low_radii.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/NEUTRINO_DRIVEN_low_radii.fig
new file mode 100644
index 0000000000000000000000000000000000000000..4078ad4a369adbf3ab60d4244860e61faba501c6
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/NEUTRINO_DRIVEN_low_radii.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/archive/NEUTRINO_DRIVEN_low_radii.png b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/NEUTRINO_DRIVEN_low_radii.png
new file mode 100644
index 0000000000000000000000000000000000000000..e2f12fcc3d0b0f38c3f45ba157cc116c04f2c0c5
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/NEUTRINO_DRIVEN_low_radii.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STANDARD.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STANDARD.fig
new file mode 100644
index 0000000000000000000000000000000000000000..1e66fe46a3507c01b90c56d6f2968833320a6f28
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STANDARD.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STANDARD.png b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STANDARD.png
new file mode 100644
index 0000000000000000000000000000000000000000..1e00da9f8cbe80f31f199fa518d30b224d50a79f
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STANDARD.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STANDARD_low_radii.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STANDARD_low_radii.fig
new file mode 100644
index 0000000000000000000000000000000000000000..3bfa1a311a4e40b3203c32e0e4fb5a5e6f0c30a8
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STANDARD_low_radii.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STANDARD_low_radii.png b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STANDARD_low_radii.png
new file mode 100644
index 0000000000000000000000000000000000000000..93a201e3c782a31b4a810792c28f2a8f6b57263d
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STANDARD_low_radii.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STOCHASTICITY_RETENTION.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STOCHASTICITY_RETENTION.fig
new file mode 100644
index 0000000000000000000000000000000000000000..47dd989f49c03b324ca2160ceeabcf06f2b8b33e
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STOCHASTICITY_RETENTION.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STOCHASTICITY_RETENTION.png b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STOCHASTICITY_RETENTION.png
new file mode 100644
index 0000000000000000000000000000000000000000..15b011dbb1e9aa4470b3cdf6306bd5afc7ae7544
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/STOCHASTICITY_RETENTION.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/archive/W_FB_REDUCED_SIGMA.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/W_FB_REDUCED_SIGMA.fig
new file mode 100644
index 0000000000000000000000000000000000000000..d611fc54b97a0d2d1c41065a6b509f297564557f
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/W_FB_REDUCED_SIGMA.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/archive/W_FB_REDUCED_SIGMA.png b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/W_FB_REDUCED_SIGMA.png
new file mode 100644
index 0000000000000000000000000000000000000000..5c80364476dbc514118e79d42a4a281aee2d645b
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/W_FB_REDUCED_SIGMA.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/archive/W_FB_REDUCED_SIGMA_low_radii.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/W_FB_REDUCED_SIGMA_low_radii.fig
new file mode 100644
index 0000000000000000000000000000000000000000..844aa552cc990a9f49feaf036c49a13210937040
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/W_FB_REDUCED_SIGMA_low_radii.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/archive/W_FB_REDUCED_SIGMA_low_radii.png b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/W_FB_REDUCED_SIGMA_low_radii.png
new file mode 100644
index 0000000000000000000000000000000000000000..2262adc0654ace70c12be1dd21352cc4c16f8b7c
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/archive/W_FB_REDUCED_SIGMA_low_radii.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_hybrid_model.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_hybrid_model.fig
new file mode 100644
index 0000000000000000000000000000000000000000..3ccf879a31fc45cedf76685763a56084e22a8d8e
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_hybrid_model.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_hybrid_model.png b/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_hybrid_model.png
new file mode 100644
index 0000000000000000000000000000000000000000..67a3fa215ef6683aff3f109072a4158f40563626
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_hybrid_model.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_neutrino_driven.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_neutrino_driven.fig
new file mode 100644
index 0000000000000000000000000000000000000000..a50a961c7ece6639c1f1aa30a063ec654158dfbd
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_neutrino_driven.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_neutrino_driven.png b/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_neutrino_driven.png
new file mode 100644
index 0000000000000000000000000000000000000000..9993f717ec72f085eb84b8324dc06abd92a4800a
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_neutrino_driven.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_standard.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_standard.fig
new file mode 100644
index 0000000000000000000000000000000000000000..75645f12b885ee36e59d513677d0afcb531036c6
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_standard.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_standard.png b/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_standard.png
new file mode 100644
index 0000000000000000000000000000000000000000..7467ca60bd49f24d5be984cdedff461df8a0f4e6
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/histogram_standard.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_100000.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_100000.fig
new file mode 100644
index 0000000000000000000000000000000000000000..e1bce7d8dd96a5f6b5ec16b6df5e5ecc03fb86d7
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_100000.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_100000.png b/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_100000.png
new file mode 100644
index 0000000000000000000000000000000000000000..540865c116fd96edbdcebb7358b420f71ed82065
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_100000.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_1000000.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_1000000.fig
new file mode 100644
index 0000000000000000000000000000000000000000..f65f9f45073d0515410caa3c953912904b03598f
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_1000000.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_1000000.png b/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_1000000.png
new file mode 100644
index 0000000000000000000000000000000000000000..993bee6cf425718c98adf225d9be041c76f75365
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_1000000.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_50000.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_50000.fig
new file mode 100644
index 0000000000000000000000000000000000000000..880d9d313049f2ceff46b5269e3c34790add02bc
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_50000.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_50000.png b/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_50000.png
new file mode 100644
index 0000000000000000000000000000000000000000..63656a4f969de228f3318660bae93ee17ac956a7
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/mass_ret_dif_radii_50000.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_100000.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_100000.fig
new file mode 100644
index 0000000000000000000000000000000000000000..8a200e9d6fbcb4fbb50e208fc353c2d459635092
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_100000.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_100000.png b/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_100000.png
new file mode 100644
index 0000000000000000000000000000000000000000..5f68923a9c7f6d3a4dfcb3d7855a3cda50f23ba7
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_100000.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_1000000.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_1000000.fig
new file mode 100644
index 0000000000000000000000000000000000000000..65811703fe7d333a9a17e8b8510f970b84a4a3cb
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_1000000.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_1000000.png b/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_1000000.png
new file mode 100644
index 0000000000000000000000000000000000000000..5f61063bac594e5fb0d2f8cb9afa38ce0d9a93e1
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_1000000.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_50000.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_50000.fig
new file mode 100644
index 0000000000000000000000000000000000000000..4ff89756dc0b926f467ba50eaac352d66b9fe33c
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_50000.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_50000.png b/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_50000.png
new file mode 100644
index 0000000000000000000000000000000000000000..b77ddff341d378b3e17333623b09268ccb847b66
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/num_ret_dif_radii_50000.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/probs_dif_vesc.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/probs_dif_vesc.fig
new file mode 100644
index 0000000000000000000000000000000000000000..3d593889a6bfd8dac8d481875f19d7f6286e75a7
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/probs_dif_vesc.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/probs_dif_vesc.png b/results/cluster properties/Cluster properties/plots_after_big_changes/probs_dif_vesc.png
new file mode 100644
index 0000000000000000000000000000000000000000..9d122615b244df12d8aed9b11ac7c30047f87cc3
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/probs_dif_vesc.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/probs_dif_vesc_2.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/probs_dif_vesc_2.fig
new file mode 100644
index 0000000000000000000000000000000000000000..36c0d342b03d5d2732bfbd3f8a6837e0f1b259bc
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/probs_dif_vesc_2.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/probs_dif_vesc_2.png b/results/cluster properties/Cluster properties/plots_after_big_changes/probs_dif_vesc_2.png
new file mode 100644
index 0000000000000000000000000000000000000000..48a8ccb3d8eec09ae183bacb30b17ce9e42f1343
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/probs_dif_vesc_2.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/velocity_dispersion.fig b/results/cluster properties/Cluster properties/plots_after_big_changes/velocity_dispersion.fig
new file mode 100644
index 0000000000000000000000000000000000000000..57bfe8532fea6b5691f7eda8dbd4e48b072d9fdd
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/velocity_dispersion.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_big_changes/velocity_dispersion.png b/results/cluster properties/Cluster properties/plots_after_big_changes/velocity_dispersion.png
new file mode 100644
index 0000000000000000000000000000000000000000..360b217ca4b0fe0ff4483d621e3c7b0422d4210c
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_big_changes/velocity_dispersion.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_change_neutrino/NEUTRINO_DRIVEN.fig b/results/cluster properties/Cluster properties/plots_after_change_neutrino/NEUTRINO_DRIVEN.fig
new file mode 100644
index 0000000000000000000000000000000000000000..bba8640950b43d67658822647a1f3b6835524f8e
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_change_neutrino/NEUTRINO_DRIVEN.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_change_neutrino/NEUTRINO_DRIVEN.png b/results/cluster properties/Cluster properties/plots_after_change_neutrino/NEUTRINO_DRIVEN.png
new file mode 100644
index 0000000000000000000000000000000000000000..efb0d8c59bed40977de1d0b35f847bc51bf7a68e
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_change_neutrino/NEUTRINO_DRIVEN.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_change_neutrino/NEUTRINO_DRIVEN_0_001.fig b/results/cluster properties/Cluster properties/plots_after_change_neutrino/NEUTRINO_DRIVEN_0_001.fig
new file mode 100644
index 0000000000000000000000000000000000000000..90fa923102e81d2e4c50d04de28a83bd6e6e96be
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_change_neutrino/NEUTRINO_DRIVEN_0_001.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_change_neutrino/NEUTRINO_DRIVEN_low_radii.fig b/results/cluster properties/Cluster properties/plots_after_change_neutrino/NEUTRINO_DRIVEN_low_radii.fig
new file mode 100644
index 0000000000000000000000000000000000000000..8fd2cb302e84d3abb1312330c0825091bb8efc38
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_change_neutrino/NEUTRINO_DRIVEN_low_radii.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_change_neutrino/NEUTRINO_DRIVEN_low_radii.png b/results/cluster properties/Cluster properties/plots_after_change_neutrino/NEUTRINO_DRIVEN_low_radii.png
new file mode 100644
index 0000000000000000000000000000000000000000..2f7fd4c5f4ce0fe497f1c2b743a995bf886f6904
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_change_neutrino/NEUTRINO_DRIVEN_low_radii.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_change_neutrino/STANDARD.fig b/results/cluster properties/Cluster properties/plots_after_change_neutrino/STANDARD.fig
new file mode 100644
index 0000000000000000000000000000000000000000..1fac6aca836ef7e40fe799040560bd934e72a9bb
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_change_neutrino/STANDARD.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_change_neutrino/STANDARD_0_001.fig b/results/cluster properties/Cluster properties/plots_after_change_neutrino/STANDARD_0_001.fig
new file mode 100644
index 0000000000000000000000000000000000000000..4c1db2eb33b169b77c4229a4bb69c4c510959694
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_change_neutrino/STANDARD_0_001.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_change_neutrino/STANDARD_low_radii.fig b/results/cluster properties/Cluster properties/plots_after_change_neutrino/STANDARD_low_radii.fig
new file mode 100644
index 0000000000000000000000000000000000000000..72b6a7b7cb8b7fcb8dcecf57588823bd86ec4f73
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_change_neutrino/STANDARD_low_radii.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_change_neutrino/STOCHASTICITY_RETENTION.fig b/results/cluster properties/Cluster properties/plots_after_change_neutrino/STOCHASTICITY_RETENTION.fig
new file mode 100644
index 0000000000000000000000000000000000000000..a4c10e809808facd8c10b175494bf1f19afab579
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_change_neutrino/STOCHASTICITY_RETENTION.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_change_neutrino/STOCHASTICITY_RETENTION.png b/results/cluster properties/Cluster properties/plots_after_change_neutrino/STOCHASTICITY_RETENTION.png
new file mode 100644
index 0000000000000000000000000000000000000000..f53214d8d359cbece68b40ad234c2b2400911e84
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_change_neutrino/STOCHASTICITY_RETENTION.png differ
diff --git a/results/cluster properties/Cluster properties/plots_after_change_neutrino/W_FB_REDUCED_SIGMA.fig b/results/cluster properties/Cluster properties/plots_after_change_neutrino/W_FB_REDUCED_SIGMA.fig
new file mode 100644
index 0000000000000000000000000000000000000000..cf5c2226203257ff779bd51d70a7cf1789fca79d
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_change_neutrino/W_FB_REDUCED_SIGMA.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_change_neutrino/W_FB_REDUCED_SIGMA_0_001.fig b/results/cluster properties/Cluster properties/plots_after_change_neutrino/W_FB_REDUCED_SIGMA_0_001.fig
new file mode 100644
index 0000000000000000000000000000000000000000..a45fd76c290fc544e2697f39d227cd2c38cd55db
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_change_neutrino/W_FB_REDUCED_SIGMA_0_001.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_change_neutrino/W_FB_REDUCED_SIGMA_low_radii.fig b/results/cluster properties/Cluster properties/plots_after_change_neutrino/W_FB_REDUCED_SIGMA_low_radii.fig
new file mode 100644
index 0000000000000000000000000000000000000000..84813bc6520e72f54d7077469584fbddfbd9e988
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_change_neutrino/W_FB_REDUCED_SIGMA_low_radii.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_change_neutrino/probs_dif_vesc.fig b/results/cluster properties/Cluster properties/plots_after_change_neutrino/probs_dif_vesc.fig
new file mode 100644
index 0000000000000000000000000000000000000000..7b780028ecdded3db1d138b2678cb157a7eb5599
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_change_neutrino/probs_dif_vesc.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_after_change_neutrino/probs_dif_vesc.png b/results/cluster properties/Cluster properties/plots_after_change_neutrino/probs_dif_vesc.png
new file mode 100644
index 0000000000000000000000000000000000000000..0dd6433f4658980a13373a4e28ae144a8e2caf2e
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_after_change_neutrino/probs_dif_vesc.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/M_rem_co_dif_metap.fig b/results/cluster properties/Cluster properties/plots_report/M_rem_co_dif_metap.fig
new file mode 100644
index 0000000000000000000000000000000000000000..ba8ea088638fb7d10a819a8aac00318d6ef6bbb5
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/M_rem_co_dif_metap.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/M_rem_co_dif_metap.png b/results/cluster properties/Cluster properties/plots_report/M_rem_co_dif_metap.png
new file mode 100644
index 0000000000000000000000000000000000000000..c3158e82f60abc6d7655b461a0cb712a01f59d7a
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/M_rem_co_dif_metap.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_radii.fig b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_radii.fig
new file mode 100644
index 0000000000000000000000000000000000000000..91022ee492fbc4a1f56d287c00af045e751de814
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_radii.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_radii_5.fig b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_radii_5.fig
new file mode 100644
index 0000000000000000000000000000000000000000..a823a2bc081334dfddb38175b7f55c641fbb497d
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_radii_5.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_radii_5.png b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_radii_5.png
new file mode 100644
index 0000000000000000000000000000000000000000..34a89f711ebbf9a57a91a0f36e225d6b9eaf64ef
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_radii_5.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_radii_5_2.fig b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_radii_5_2.fig
new file mode 100644
index 0000000000000000000000000000000000000000..f7e51b300829f043fb4cfee5cd9ce5c58b39ee81
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_radii_5_2.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_with_real_clusters.fig b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_with_real_clusters.fig
new file mode 100644
index 0000000000000000000000000000000000000000..3629ece3f63aae6a98dc7c3f4c6f743ec658d0e3
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_with_real_clusters.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_with_real_clusters.png b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_with_real_clusters.png
new file mode 100644
index 0000000000000000000000000000000000000000..8cfc27bfa4eb8639781a60fdcbf2bbee145a5a3e
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_with_real_clusters.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_with_real_clusters_2.fig b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_with_real_clusters_2.fig
new file mode 100644
index 0000000000000000000000000000000000000000..beeacddc1c8ff1d88e4d311769c6c928a46c6b6e
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_with_real_clusters_2.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD_low_radii.fig b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_radii.fig
new file mode 100644
index 0000000000000000000000000000000000000000..a1ac42d05f00bd93e06763c9e42f985de6376592
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_radii.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD_low_radii_5.fig b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_radii_5.fig
new file mode 100644
index 0000000000000000000000000000000000000000..4641822ae8e18df6f8d4232a86b0df694a9c1a51
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_radii_5.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD_low_radii_5.png b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_radii_5.png
new file mode 100644
index 0000000000000000000000000000000000000000..df32a9d9412bc77598872d89a247f70000e593e7
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_radii_5.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD_with_real_clusters.fig b/results/cluster properties/Cluster properties/plots_report/STANDARD_with_real_clusters.fig
new file mode 100644
index 0000000000000000000000000000000000000000..8f7a09492e81af44fb47b9af044718e08ea0dc01
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD_with_real_clusters.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD_with_real_clusters.png b/results/cluster properties/Cluster properties/plots_report/STANDARD_with_real_clusters.png
new file mode 100644
index 0000000000000000000000000000000000000000..570e67bcd82fdc78583d4aec4fb540969c2bb5e7
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD_with_real_clusters.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD_with_real_clusters_2.fig b/results/cluster properties/Cluster properties/plots_report/STANDARD_with_real_clusters_2.fig
new file mode 100644
index 0000000000000000000000000000000000000000..f56fe162f9f92e32fa14c11f2c1b86764935084c
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD_with_real_clusters_2.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STOCHASTICITY_RETENTION.fig b/results/cluster properties/Cluster properties/plots_report/STOCHASTICITY_RETENTION.fig
new file mode 100644
index 0000000000000000000000000000000000000000..610fcc436dcee5083139a8f2a09497cfb3d3f235
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STOCHASTICITY_RETENTION.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STOCHASTICITY_RETENTION.png b/results/cluster properties/Cluster properties/plots_report/STOCHASTICITY_RETENTION.png
new file mode 100644
index 0000000000000000000000000000000000000000..f226f4285c78e71984c14e004c77c250c9b249f9
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STOCHASTICITY_RETENTION.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_radii.fig b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_radii.fig
new file mode 100644
index 0000000000000000000000000000000000000000..1cbec0a08195956ba73f82f0b2325ff145b173a0
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_radii.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_radii_5.fig b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_radii_5.fig
new file mode 100644
index 0000000000000000000000000000000000000000..c845803aafefc68c6743898e3c608566bb994c6b
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_radii_5.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_radii_5.png b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_radii_5.png
new file mode 100644
index 0000000000000000000000000000000000000000..789fa5d66683d02dedd19a2d57026c25c2b408ab
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_radii_5.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_with_real_clusters.fig b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_with_real_clusters.fig
new file mode 100644
index 0000000000000000000000000000000000000000..330075689aa6901fa1524c85b6dd654438de9070
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_with_real_clusters.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_with_real_clusters.png b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_with_real_clusters.png
new file mode 100644
index 0000000000000000000000000000000000000000..545f9f8cd1184a4534db970c5fc234bd2f5fd324
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_with_real_clusters.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_with_real_clusters_2.fig b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_with_real_clusters_2.fig
new file mode 100644
index 0000000000000000000000000000000000000000..a82f7246beeb67f366e45efef399a8d16a2028a1
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_with_real_clusters_2.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/fallback_fraction.fig b/results/cluster properties/Cluster properties/plots_report/fallback_fraction.fig
new file mode 100644
index 0000000000000000000000000000000000000000..b187f67c94e7935383d840b7b4796f607a3fcd22
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/fallback_fraction.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/fallback_fraction.png b/results/cluster properties/Cluster properties/plots_report/fallback_fraction.png
new file mode 100644
index 0000000000000000000000000000000000000000..16e1a9393194e98b523a783ad26ec83805e396d4
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/fallback_fraction.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/histogram_fb_reduced_sigma.fig b/results/cluster properties/Cluster properties/plots_report/histogram_fb_reduced_sigma.fig
new file mode 100644
index 0000000000000000000000000000000000000000..6ff60ff31192667dd105eb7ebcee42559676b5a6
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/histogram_fb_reduced_sigma.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/histogram_fb_reduced_sigma.png b/results/cluster properties/Cluster properties/plots_report/histogram_fb_reduced_sigma.png
new file mode 100644
index 0000000000000000000000000000000000000000..0b25149f5dd7230f363349d9b41165ae3edb6b38
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/histogram_fb_reduced_sigma.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/histogram_neutrino_driven.fig b/results/cluster properties/Cluster properties/plots_report/histogram_neutrino_driven.fig
new file mode 100644
index 0000000000000000000000000000000000000000..c0262398893f7b6a7bd04e8387ef28873f5faf71
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/histogram_neutrino_driven.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/histogram_neutrino_driven.png b/results/cluster properties/Cluster properties/plots_report/histogram_neutrino_driven.png
new file mode 100644
index 0000000000000000000000000000000000000000..6741b71274d9134a97d3293d86dcbdb7e7404530
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/histogram_neutrino_driven.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/histogram_standard.fig b/results/cluster properties/Cluster properties/plots_report/histogram_standard.fig
new file mode 100644
index 0000000000000000000000000000000000000000..7a7904af096171cf23329b349484fd90b77b1695
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/histogram_standard.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/histogram_standard.png b/results/cluster properties/Cluster properties/plots_report/histogram_standard.png
new file mode 100644
index 0000000000000000000000000000000000000000..2e6918e524e274903108afe538dc8ebb8154e81d
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/histogram_standard.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/probs_dif_vesc.fig b/results/cluster properties/Cluster properties/plots_report/probs_dif_vesc.fig
new file mode 100644
index 0000000000000000000000000000000000000000..be511e3f10826a3302f955b2038c9f69177b0217
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/probs_dif_vesc.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/probs_dif_vesc.png b/results/cluster properties/Cluster properties/plots_report/probs_dif_vesc.png
new file mode 100644
index 0000000000000000000000000000000000000000..e9a929b84882a6e7e074da1f2008e872ccd71828
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/probs_dif_vesc.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/probs_dif_vesc_2.fig b/results/cluster properties/Cluster properties/plots_report/probs_dif_vesc_2.fig
new file mode 100644
index 0000000000000000000000000000000000000000..5f4279e67ae624fe32fe247ba453d8f4bcfeb4b6
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/probs_dif_vesc_2.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/probs_dif_vesc_2.png b/results/cluster properties/Cluster properties/plots_report/probs_dif_vesc_2.png
new file mode 100644
index 0000000000000000000000000000000000000000..3163af1adb99229a1f6ce2306f1d7016c0108b18
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/probs_dif_vesc_2.png differ
diff --git a/results/cluster properties/Cluster properties/prob_escape.m b/results/cluster properties/Cluster properties/prob_escape.m
index cf1096301b5c8ff2d9b23ef44c0925787a5e5b50..807112d60c02e6af35bd4410af67b219b3ceb17c 100644
--- a/results/cluster properties/Cluster properties/prob_escape.m
+++ b/results/cluster properties/Cluster properties/prob_escape.m
@@ -7,16 +7,14 @@ Description:
%}
%Input:
% Z: Metallicity
-% M: Zero-age main sequence mass (in solar masses)
+% M_zams: Zero-age main sequence mass (in solar masses)
% M_ns: neutron star mass (in solar masses)
% sig_ns: neutron star velocity dispersion (in km/s)
% v_esc: escape velocity (in km/s)
% k: parameter that defines the kick prescription:
-% k=0 ==> no fallback and reduced velocity dispersion
-% k=1 ==> fallback considered and reduced velocity
-% dispersion
-% K=2 ==> fallback considered and not a reduced
-% velocity dispersion.
+% k=0 ==> neutrino-driven NK
+% k=1 ==> hybrid model
+% K=2 ==> standard kick
%Output:
% p: probability of v>v_esc
@@ -26,16 +24,42 @@ function p=prob_escape(Z,M_zams,M_ns,sig_ns,v_esc,k)
M_bh=fun_M_rem(Z,M_zams);
%And then we can obtain the BH's velocity dispersion and the fallback
%depending on the prescription:
-if k==1
- sig_bh=(M_ns./M_bh).*sig_ns;
- f_fb=fun_fb(Z,M_zams);
-elseif k==0
- sig_bh=(M_ns./M_bh).*sig_ns;
- f_fb=0;
-elseif k==2
- f_fb=fun_fb(Z,M_zams);
- sig_bh=sig_ns;
+f_fb=fun_fb(Z,M_zams);
+%And we have to take into account that the probability of escape is 0 if
+%sigma is 0, so we have to distinguish between two cases due to the
+%fallback fraction reaching one. To do that we define a M_zams_crit using
+%the vec_M_examples from cpdf.mat so that it does not depend on the cluster
+%generated.
+load('cpdf');
+M_bh_examples=fun_M_rem(Z,vec_M_examples);f_fb_examples=fun_fb(Z,vec_M_examples);
+ll=find(abs(f_fb_examples-0.99)<1e-2); M_zams_crit=vec_M_examples(ll(1));
+sig_bh=zeros(1,length(M_zams));
+for ii=1:length(M_zams)
+ if M_zams(ii)<M_zams_crit
+ if k==1
+ sig_bh(ii)=(M_ns./M_bh(ii)).*(1-f_fb(ii)).*sig_ns;
+ p(ii)=sqrt(2/pi)*(v_esc./sig_bh(ii)).*exp(-(v_esc.^2)./(2.*(sig_bh(ii)).^2))+(1-erf(v_esc./(sqrt(2).*sig_bh(ii))));
+ elseif k==0
+ sig_bh(ii)=(M_ns./M_bh(ii)).*sig_ns;
+ f_fb(ii)=0;
+ p(ii)=sqrt(2/pi)*(v_esc./sig_bh(ii)).*exp(-(v_esc.^2)./(2.*(sig_bh(ii)).^2))+(1-erf(v_esc./(sqrt(2).*sig_bh(ii))));
+ elseif k==2
+ sig_bh(ii)=(1-f_fb(ii)).*sig_ns;
+ p(ii)=sqrt(2/pi)*(v_esc./sig_bh(ii)).*exp(-(v_esc.^2)./(2.*(sig_bh(ii)).^2))+(1-erf(v_esc./(sqrt(2).*sig_bh(ii))));
+ end
+ else
+ if k==1
+ sig_bh(ii)=(M_ns./M_bh(ii)).*(1-f_fb(ii)).*sig_ns;
+ p(ii)=0;
+ elseif k==0
+ sig_bh(ii)=(M_ns./M_bh(ii)).*sig_ns;
+ f_fb(ii)=0;
+ p(ii)=sqrt(2/pi)*(v_esc./sig_bh(ii)).*exp(-(v_esc.^2)./(2.*(sig_bh(ii)).^2))+(1-erf(v_esc./(sqrt(2).*sig_bh(ii))));
+ elseif k==2
+ sig_bh(ii)=(1-f_fb(ii)).*sig_ns;
+ p(ii)=0;
+ end
+ end
end
-%And then we can compute the probability of escaping:
-p=sqrt(2/pi)*(v_esc./sig_bh).*exp(-(v_esc.^2)./(2*((1-f_fb).^2).*(sig_bh).^2))+(1-f_fb).*(1-erf(v_esc./(sqrt(2)*(1-f_fb).*sig_bh)));
end
+
diff --git a/results/cluster properties/maxwellian_python/kicks.py b/results/cluster properties/maxwellian_python/kicks.py
new file mode 100644
index 0000000000000000000000000000000000000000..c41d9c9f81255a54dc24b75912df145a0b44087c
--- /dev/null
+++ b/results/cluster properties/maxwellian_python/kicks.py
@@ -0,0 +1,143 @@
+# -*- coding: utf-8 -*-
+"""
+Created on Thu Aug 15 09:05:53 2019
+
+@author: m13239
+"""
+
+## KICK VELOCITIES ACCORDING TO A MAXWELLIAN DISTRIBUTION
+
+# We start by importing the necessary
+
+import math
+
+import matplotlib.pyplot as plt
+
+import numpy as np
+
+​
+
+from scipy.special import erf
+
+​
+
+import random
+
+​
+
+​
+
+# Hnn
+
+# http://pulsars.info/Thesis.pdf
+
+# https://arxiv.org/pdf/1807.11489.pdf
+
+# http://mathworld.wolfram.com/MaxwellDistribution.html
+
+# Give different formulae
+
+​
+
+# Especially the first factor is a bit unclear sometimes
+
+​#We define the probability density gunction and the cumulative probability density
+
+def maxwell_prob_density(vkick, dispersion):
+
+ return np.sqrt(2.0/np.pi)* (vkick**2)/(dispersion**3) * np.exp(-((vkick**2)/(2*(dispersion**2))))
+
+​
+
+def maxwell_cum_density(vkick, dispersion):
+
+ return erf((vkick)/(np.sqrt(2)*dispersion)) - ((vkick * np.exp(-((vkick**2)/(2*(dispersion**2)))))/(dispersion))*(np.sqrt(2.0/math.pi))
+
+​
+
+​
+
+# sampling from distribution:
+
+def sample_maxwell(sigma):
+
+ X = random.random()
+
+ Y = random.random()
+
+ W = random.random()
+
+ Z = random.random()
+
+​
+
+ s = sigma*math.sqrt(-2*math.log(1-X))
+
+ p = sigma*math.sqrt(-2*math.log(1-Y))
+
+​
+
+ theta = 2*math.pi*Y
+
+ phi = 2*math.pi*Z
+
+​
+
+ u = s * math.cos(theta)
+
+ v = s * math.sin(theta)
+
+ w = p * math.cos(phi)
+
+​
+
+ kick = math.sqrt((u**2) + (v**2) + (w**2))
+
+ return kick
+
+​
+
+​# We set the velocity dispersion of the maxwellian distribution:
+
+dispersion = 265
+
+​# We set the number of samplings:
+
+samples = 50000
+
+# And then we sample from that distribution according to the number of times set above (and obviously the defined velocity dispersion)
+
+sampled_array = [sample_maxwell(dispersion) for _ in range(samples)]
+
+
+
+v_arr = np.arange(0, 1500, 1) # This arranges values from 0 to 1500 with steps separated by 1 (0,1,2,3,...,1498,1499)
+# We then compute the probability density and the cumulative probability density for the v_arr defined above.
+mw_prob = maxwell_prob_density(v_arr, dispersion)
+
+mw_cum = maxwell_cum_density(v_arr, dispersion)
+
+​# The idea is that now we can compute on top of the sampling of velocities the maxwellian distribution itself.
+
+fig, axes = plt.subplots(nrows=1, ncols=2)
+
+axes[0].plot(v_arr, mw_prob)
+
+axes[0].hist(sampled_array, density=True, bins=100)
+
+axes[0].set_ylabel('Probability density')
+
+axes[0].set_xlabel('v (km s-1)')
+
+​
+
+axes[1].plot(v_arr, mw_cum)
+
+axes[1].set_ylabel('Cumulative density')
+
+axes[1].set_xlabel('v (km s-1)')
+
+​
+
+plt.show()
+