diff --git a/.command.txt.swp b/.command.txt.swp
new file mode 100644
index 0000000000000000000000000000000000000000..b5eff9f9c7071ccb54dff00dac9dda50e6b0b589
Binary files /dev/null and b/.command.txt.swp differ
diff --git a/mcluster_script.sh b/mcluster_script.sh
index 59492ca441367d615c0eb81d2a18a6d21b35b659..60b333e9dd01cf9ffdca7bc350fa15d7955b0017 100755
--- a/mcluster_script.sh
+++ b/mcluster_script.sh
@@ -19,21 +19,22 @@ echo "Name run: $NAME_RUN"
MCLUSTER_CONFIG_DIR=$PWD
#MCLUSTER_CONFIG_DIR='/vol/ph/astro_code/dhendriks/work_dir/mcluster'
-MCLUSTER_EXECUTABLE='/vol/ph/astro_code/dhendriks/mcluster/mcluster'
+MCLUSTER_EXECUTABLE='/user/HS103/m13239/Desktop/mcluster/mcluster'
# Change these settings
TARGET_DIR=$MCLUSTER_CONFIG_DIR
# Store settings for mcluster in here.
-MCLUSTER_COMMANDS="-C5 -N 1000 -R 0.8 -P 3 -W 3.0 -f 1 -b 0.5 -o"
+MCLUSTER_COMMANDS="-C5 -N 1000 -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"
echo "Writing files to $MCLUSTER_CONFIG_DIR"
# Execute mcluster cmd
-$MCLUSTER_EXECUTABLE $MCLUSTER_COMMANDS $MCLUSTER_CONFIG_DIR/$NAME_RUN | tee $MCLUSTER_CONFIG_DIR/$NAME_RUN.summary
+$MCLUSTER_EXECUTABLE $MCLUSTER_COMMANDS $NAME_RUN | tee $MCLUSTER_CONFIG_DIR/$NAME_RUN.summary
# Write command to file
-echo "$MCLUSTER_EXECUTABLE $MCLUSTER_COMMANDS $MCLUSTER_CONFIG_DIR/$NAME_RUN | tee $MCLUSTER_CONFIG_DIR/$NAME_RUN.summary" > $MCLUSTER_CONFIG_DIR/command.txt
+echo "$MCLUSTER_EXECUTABLE $MCLUSTER_COMMANDS $NAME_RUN | tee $MCLUSTER_CONFIG_DIR/$NAME_RUN.summary" > $MCLUSTER_CONFIG_DIR/command.txt
# Rename fïles:
if [ -f "$MCLUSTER_CONFIG_DIR/$NAME_RUN.input.fort.10" ]; then
@@ -46,4 +47,4 @@ fi
if [ -f "$MCLUSTER_CONFIG_DIR/$NAME_RUN.dat.10" ]; then
mv -f $MCLUSTER_CONFIG_DIR/$NAME_RUN.dat.10 $MCLUSTER_CONFIG_DIR/dat.10
-fi
\ No newline at end of file
+fi
diff --git a/requirements.txt b/requirements.txt
index 6973aea407e29071af5fbeccdc49e3eb73a23c3d..eab3860c771828dac5445c7ee788a12cb58108a6 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -7,8 +7,8 @@ h5py==2.9.0
kiwisolver==1.1.0
Mako==1.0.10
MarkupSafe==1.1.1
-matplotlib==3.1.0
-more-itertools==7.0.0
+#matplotlib==3.1.0
+#more-itertools==7.0.0
numpy==1.16.3
pandas==0.24.2
pluggy==0.11.0
diff --git a/results/cluster properties/Cluster properties/CLUSTER_PROPERTIES.m b/results/cluster properties/Cluster properties/CLUSTER_PROPERTIES.m
new file mode 100644
index 0000000000000000000000000000000000000000..1df557074ca57f839ffc7acb63dae05727686729
--- /dev/null
+++ b/results/cluster properties/Cluster properties/CLUSTER_PROPERTIES.m
@@ -0,0 +1,140 @@
+%% CLUSTER PROPERTIES
+% For now I think that the best way to do it is to only leave to choose the
+% metallicity and the kick prescription:
+function CLUSTER_PROPERTIES(Z,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
+
+
+%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);
+%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_zams=[10 20 30 40];
+p=prob_escape(Z,vec_M_zams(1),M_ns,sig_ns,v_esc,k);f_ret=1-p;
+figure(1)
+subplot(2,2,1)
+contour(M_cl./(10^5),r_h,p,'ShowText','on')
+title(['Retention fraction of stars. ' num2str(kick) '. Z=' num2str(Z)])
+xlabel('M_5 (M_{cl}/(10^5 M_{sun})')
+ylabel('r_h (pc)')
+legend(['M=' num2str(vec_M_zams(1))])
+
+p=prob_escape(Z,vec_M_zams(2),1.4,265,v_esc,k);f_ret=1-p;
+subplot(2,2,2)
+contour(M_cl./(10^5),r_h,f_ret,'ShowText','on')
+title(['Retention fraction of stars. ' num2str(kick) '. Z=' num2str(Z)])
+xlabel('M_5 (M_{cl}/(10^5 M_{sun})')
+ylabel('r_h (pc)')
+legend(['M=' num2str(vec_M_zams(2))])
+
+p=prob_escape(Z,vec_M_zams(3),1.4,265,v_esc,k);f_ret=1-p;
+subplot(2,2,3)
+contour(M_cl./(10^5),r_h,f_ret,'ShowText','on')
+title(['Retention fraction of stars. ' num2str(kick) '. Z=' num2str(Z)])
+xlabel('M_5 (M_{cl}/(10^5 M_{sun})')
+ylabel('r_h (pc)')
+legend(['M=' num2str(vec_M_zams(3))])
+
+p=prob_escape(Z,vec_M_zams(3),1.4,265,v_esc,k);f_ret=1-p;
+subplot(2,2,4)
+contour(M_cl./(10^5),r_h,f_ret,'ShowText','on')
+title(['Retention fraction of stars. ' num2str(kick) '. Z=' num2str(Z)])
+xlabel('M_5 (M_{cl}/(10^5 M_{sun})')
+ylabel('r_h (pc)')
+legend(['M=' num2str(vec_M_zams(4))])
+
+
+%We then want to plot the escape velocity as a function of the mass of the
+%cluster and the radius in parsecs:
+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 v escape 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:
+v_esc=linspace(20,200,1000);M_bh=linspace(3,50,1000);
+[M_bh,v_esc]=meshgrid(M_bh,v_esc); p=prob_escape(Z,M_bh,M_ns,sig_ns,v_esc,k);
+figure(3)
+contour(M_bh,v_esc,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title(['Probability of escape. ' num2str(kick) '. Z=' num2str(Z)])
+xlabel('M_{BH} (in M_{sun})')
+ylabel('Escape velocity (in km/s)')
+
+
+%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_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(4)
+subplot(2,2,[1,3])
+plot(vec_M_zams,vec_M_co,'-k')
+%text(15,25,'Z=2e-2')
+title(['M_{rem} & M_{co}. Z=' num2str(Z)])
+xlabel('M_{zams} (in M_{sun})')
+ylabel('M (in M_{sun})')
+hold on
+plot(vec_M_zams,vec_M_rem,'-r')
+hold on
+%We find the M_zams for which the f_fb has certain values.
+kk05=find(abs(vec_f_fb-0.5)<1e-2);M_zams_0_5=vec_M_zams(kk05(4));M_rem_0_5=fun_M_rem(Z,M_zams_0_5);
+kk075=find(abs(vec_f_fb-0.75)<1e-3);M_zams_0_75=vec_M_zams(kk075(1));M_rem_0_75=fun_M_rem(Z,M_zams_0_75);
+kk1=find(abs(vec_f_fb-1)<1e-5);M_zams_1=vec_M_zams(kk1(1));M_rem_1=fun_M_rem(Z,M_zams_1);
+plot(M_zams_0_5,M_rem_0_5,'k*','MarkerSize',10)
+hold on
+plot(M_zams_0_75,M_rem_0_75,'ko','MarkerSize',10)
+hold on
+plot(M_zams_1,M_rem_1,'k+','MarkerSize',10)
+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')
+%We compute probability of escape for different escape velocities with the kick
+%prescription chosen:
+vec_v_esc=[30 40 50];
+for ii=1:length(vec_v_esc)
+v_esc=vec_v_esc(ii);
+ for jj=1:length(vec_M_rem)
+ M_zams=vec_M_zams(jj);M_bh=fun_M_rem(Z,M_zams);
+ p(ii,jj)=prob_escape(Z,M_bh,M_ns,sig_ns,v_esc,k);
+ end
+end
+%We then plot the probability of escape as a function of BH mass:
+subplot(2,2,2)
+plot(vec_M_rem,p(1,:),'-k')
+title(['Probability of escape' num2str(kick), ' Z=' num2str(Z)])
+xlabel('M_{rem} (in M_{sun})')
+ylabel('Probability of escape')
+xlim([0 vec_M_rem(end)])
+%xticks([0 2.5 5 7.5 10 12.5 15 17.5 20])
+hold on
+plot(vec_M_rem,p(2,:),'-r')
+hold on
+plot(vec_M_rem,p(3,:),'-b')
+legend('v_{esc}=30','v_{esc}=40','v_{esc}=50')
+%And the next thing is to plot the fallback fraction:
+subplot(2,2,4)
+plot(vec_M_rem,vec_f_fb,'-k','LineWidth',1.2)
+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/SAMPLING.m b/results/cluster properties/Cluster properties/SAMPLING.m
new file mode 100644
index 0000000000000000000000000000000000000000..1bc25d3e9e5210b6b35476df458601eae59f5892
--- /dev/null
+++ b/results/cluster properties/Cluster properties/SAMPLING.m
@@ -0,0 +1,195 @@
+%% SAMPLING
+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 density
+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
+%}
+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);
+%So then the histogram is computed simply by doing (taking into account
+%that the vectors we obtain as an output for the function above are for the
+%ZAMS masses, so we need to transform that into BH masses.
+vec_m_bh_initial=fun_M_rem(Z,vec_m_bh_initial);
+vec_m_bh_final=fun_M_rem(Z,vec_m_bh_final);
+figure(1)
+h1=histogram(vec_m_bh_initial,45,'FaceColor','b');
+%h1.BinWidth = 0.25;
+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,45,'FaceColor','r');
+%h2.BinWidth = 0.25;
+legend('IMF','after kicks')
+
+
+%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.
+
+
+N_tot=[];
+vec_ii=[1:10];
+mat_m=[];
+for ii=1:length(vec_ii)
+ jj=1;
+ vec_m=[];%We will use this vector to obtain N_tot because if we try to do it via the matrix m_mat, as it is padded with 0's,we always get the same N_tot
+ M_tot=0;
+ while M_tot<M_cl
+ r=rand;
+ kk=find(abs(eval_fun-r)<1e-3);
+ m=vec_M_examples(kk(randi(length(kk))));%we use this randi command because for a given r, given the resolution we can ask for,
+ %there might be more than one mass that suits such requirements. This
+ %way, if there is more than one possible value, we asign that value
+ %randomly.
+ M_tot=M_tot+m;
+ mat_m(ii,jj)=m;
+ vec_m(jj)=m;
+ jj=jj+1;
+ end
+ N_tot(ii)=length(vec_m);
+ 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 functions presented by Spera, Mapelli & Bressan (2015), and Fryer
+%et al. (2012) for the delayed supernova mechanism because the mass is too low (minimum M_zams of 7 solar masses).
+%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 stars of each of them:
+N_bh(ii)=sum(vec_m>M_z_min_bh_test);
+N_stars(ii)=sum(vec_m<M_z_min_ns_test);
+N_ns(ii)=N_tot(ii)-N_bh(ii)-N_stars(ii);
+
+M_bh(ii)=sum(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-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:
+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
+%neutron star mass. Because this stars are retained and considering them in
+%an if statement only slows the computation:
+pos_stars=mat_m(ii,:)<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=mat_m(ii,:).*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_i=mat_m(ii,:)-vec_m_final_i;
+vec_bh_ns=nonzeros(vec_bh_ns_i)';%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));
+ for jj=1:length(vec_bh_ns)
+ if vec_bh_ns(ii,jj)>M_z_min_bh_test
+ %k=1 ==>with fallback
+ %k=0 ==>without fallback
+ %f_ret=(1-feval(@prob_escape,Z,m,1.4,265,vec_v_esc(dd,ff),1));
+ %We now generate a random number, if it is below f_ret it is
+ %retained, if not it is kicked out of the cluster.
+ rr=rand;
+ if rr<=f_ret(jj)
+ vec_m_final(ii,ll)=vec_bh_ns(ii,jj);
+ ll=ll+1;
+ else
+ ll=ll;
+ end
+ else
+ %Again we compute the probability of it being retained:
+ %f_ret=(1-feval(@prob_esc_ns,265,v_esc(dd,ff)));
+ rr=rand;
+ if rr<=(1-feval(@prob_esc_ns,265,v_esc))
+ vec_m_final(ii,ll)=vec_bh_ns(ii,jj);
+ ll=ll+1;
+ else
+ ll=ll;
+ end
+ end
+ end
+ N_tot_final(ii)=length(vec_m_final);
+ M_tot_final(ii)=sum(vec_m_final);
+ N_bh_final(ii)=sum(vec_m_final>M_z_min_bh_test);
+ N_stars_final(ii)=sum(vec_m_final<M_z_min_ns_test);
+ N_ns_final(ii)=N_tot_final(ii)-N_bh_final(ii)-N_stars_final(ii);
+ f_N_bh(ii)=N_bh_final(ii)/N_bh(ii);f_N_ns(ii)=N_ns_final(ii)/N_ns(ii);f_N_stars(ii)=N_stars_final(ii)/N_tot(ii);
+
+ M_bh_final(ii)=sum(vec_m_final(find(vec_m_final>M_z_min_bh_test)));
+ M_stars_final(ii)=sum(vec_m_final(find(vec_m_final<M_z_min_ns_test)));
+ M_ns_final(ii)=M_tot_final(ii)-M_bh_final(ii)-M_stars_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);
+end
+%And now we can plot this:
+figure(2)
+plot(vec_ii,N_tot,'-k')
+title(['Number of stars for 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')
+hold on
+plot(vec_ii,N_tot_final,'-r')
+legend(['IMF','after kicks,' num2str(kick)])
+figure(3)
+subplot(2,1,1)
+plot(vec_ii,N_bh,'-r')
+hold on
+plot(vec_ii,N_bh_final,'-b')
+
+
+end
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/SAMPLING.m~ b/results/cluster properties/Cluster properties/SAMPLING.m~
new file mode 100644
index 0000000000000000000000000000000000000000..01fbfe27a4981404e943907c07ba53a619737d39
--- /dev/null
+++ b/results/cluster properties/Cluster properties/SAMPLING.m~
@@ -0,0 +1,195 @@
+%% SAMPLING
+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 density
+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
+%}
+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);
+%So then the histogram is computed simply by doing (taking into account
+%that the vectors we obtain as an output for the function above are for the
+%ZAMS masses, so we need to transform that into BH masses.
+vec_m_bh_initial=fun_M_rem(Z,vec_m_bh_initial);
+vec_m_bh_final=fun_M_rem(Z,vec_m_bh_final);
+figure(1)
+h1=histogram(vec_m_bh_initial,45,'FaceColor','b');
+%h1.BinWidth = 0.25;
+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,45,'FaceColor','r');
+%h2.BinWidth = 0.25;
+legend('IMF','after kicks')
+
+
+%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.
+
+
+N_tot=[];
+vec_ii=[1:10];
+mat_m=[];
+for ii=1:length(vec_ii)
+ jj=1;
+ vec_m=[];%We will use this vector to obtain N_tot because if we try to do it via the matrix m_mat, as it is padded with 0's,we always get the same N_tot
+ M_tot=0;
+ while M_tot<M_cl
+ r=rand;
+ kk=find(abs(eval_fun-r)<1e-3);
+ m=vec_M_examples(kk(randi(length(kk))));%we use this randi command because for a given r, given the resolution we can ask for,
+ %there might be more than one mass that suits such requirements. This
+ %way, if there is more than one possible value, we asign that value
+ %randomly.
+ M_tot=M_tot+m;
+ mat_m(ii,jj)=m;
+ vec_m(jj)=m;
+ jj=jj+1;
+ end
+ N_tot(ii)=length(vec_m);
+ 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 functions presented by Spera, Mapelli & Bressan (2015), and Fryer
+%et al. (2012) for the delayed supernova mechanism because the mass is too low (minimum M_zams of 7 solar masses).
+%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 stars of each of them:
+N_bh(ii)=sum(vec_m>M_z_min_bh_test);
+N_stars(ii)=sum(vec_m<M_z_min_ns_test);
+N_ns(ii)=N_tot(ii)-N_bh(ii)-N_stars(ii);
+
+M_bh(ii)=sum(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-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:
+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
+%neutron star mass. Because this stars are retained and considering them in
+%an if statement only slows the computation:
+pos_stars=mat_m(ii,:)<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=mat_m(ii,:).*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_i=mat_m(ii,:)-vec_m_final_i;
+vec_bh_ns=nonzeros(vec_bh_ns_i)';%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));
+ for jj=1:length(vec_bh_ns)
+ if vec_bh_ns(ii,jj)>M_z_min_bh_test
+ %k=1 ==>with fallback
+ %k=0 ==>without fallback
+ %f_ret=(1-feval(@prob_escape,Z,m,1.4,265,vec_v_esc(dd,ff),1));
+ %We now generate a random number, if it is below f_ret it is
+ %retained, if not it is kicked out of the cluster.
+ rr=rand;
+ if rr<=f_ret(jj)
+ vec_m_final(ii,ll)=vec_bh_ns(ii,jj);
+ ll=ll+1;
+ else
+ ll=ll;
+ end
+ else
+ %Again we compute the probability of it being retained:
+ %f_ret=(1-feval(@prob_esc_ns,265,v_esc(dd,ff)));
+ rr=rand;
+ if rr<=(1-feval(@prob_esc_ns,265,v_esc))
+ vec_m_final(ii,ll)=vec_bh_ns(ii,jj);
+ ll=ll+1;
+ else
+ ll=ll;
+ end
+ end
+ end
+ N_tot_final(ii)=length(vec_m_final(ii,:));
+ M_tot_final(ii)=sum(vec_m_final(ii,:));
+ N_bh_final(ii)=sum(vec_m_final(ii,:)>M_z_min_bh_test);
+ N_stars_final(ii)=sum(vec_m_final(ii,:)<M_z_min_ns_test);
+ N_ns_final(ii)=N_tot_final(ii)-N_bh_final(ii)-N_stars_final(ii);
+ f_N_bh(ii)=N_bh_final(ii)/N_bh(ii);f_N_ns(ii)=N_ns_final(ii)/N_ns(ii);f_N_stars(ii)=N_stars_final(ii)/N_tot(ii);
+
+ M_bh_final(ii)=sum(vec_m_final(ii,(find(vec_m_final>M_z_min_bh_test))));
+ M_stars_final(ii)=sum(vec_m_final(ii,(find(vec_m_final<M_z_min_ns_test))));
+ M_ns_final(ii)=M_tot_final(ii)-M_bh_final(ii)-M_stars_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);
+end
+%And now we can plot this:
+figure(2)
+plot(vec_ii,N_tot,'-k')
+title(['Number of stars for 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')
+hold on
+plot(vec_ii,N_tot_final,'-r')
+legend(['IMF','after kicks,' num2str(kick)])
+figure(3)
+subplot(2,1,1)
+plot(vec_ii,N_bh,'-r')
+hold on
+plot(vec_ii,N_bh_final,'-b')
+
+
+end
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/bh_retention_fraction_fb_265.jpg b/results/cluster properties/Cluster properties/bh_retention_fraction_fb_265.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..d3fa7d12a87ce27556b8c16a1a634587ba43f5eb
Binary files /dev/null and b/results/cluster properties/Cluster properties/bh_retention_fraction_fb_265.jpg differ
diff --git a/results/cluster properties/Cluster properties/bh_retention_fraction_fb_265jpg.jpg b/results/cluster properties/Cluster properties/bh_retention_fraction_fb_265jpg.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..d3fa7d12a87ce27556b8c16a1a634587ba43f5eb
Binary files /dev/null and b/results/cluster properties/Cluster properties/bh_retention_fraction_fb_265jpg.jpg differ
diff --git a/results/cluster properties/Cluster properties/fun_M_co.m b/results/cluster properties/Cluster properties/fun_M_co.m
index 39da5df61da32b6252227d021a59dbfa85b192c9..a86bbe25db664b98c6a8673b320f7ed95a08ff5f 100644
--- a/results/cluster properties/Cluster properties/fun_M_co.m
+++ b/results/cluster properties/Cluster properties/fun_M_co.m
@@ -1,4 +1,10 @@
-%% M_co as a function of M_zams and Z
+%% MASS OF THE CO CORE AS A FUNCTON OF Z AND M_zams
+%{
+Description:
+ This function computes the mass of the carbon-oxygen core according to the
+ results presented by Spera, Mapelli & Bressan (2015), and Fryer
+ et al. (2012) for the delayed supernova mechanism.
+%}
%Input:
% Z:Metallicity
% M_zams:Zero age main sequence mass (in solar masses)
diff --git a/results/cluster properties/Cluster properties/fun_M_rem.m b/results/cluster properties/Cluster properties/fun_M_rem.m
index 0df6083132a694135e157fcba89d7eac82dfb547..afca185858ec53fc8fc3750562eefed87d2cbd77 100644
--- a/results/cluster properties/Cluster properties/fun_M_rem.m
+++ b/results/cluster properties/Cluster properties/fun_M_rem.m
@@ -1,4 +1,10 @@
-%% M_rem as a function of Z and M_zams
+%% REMNANT MASS (M_rem) AS A FUNCTION OF Z AND M_zams
+%{
+Description:
+ This function computes the mass of the remnant star according to the
+ results presented by Spera, Mapelli & Bressan (2015), and Fryer
+ et al. (2012) for the delayed supernova mechanism.
+%}
%Input:
% Z:Metallicity
% M_zams:Zero age main sequence mass (in solar masses)
@@ -6,7 +12,7 @@
%Output:
% M_rem:Mass of the stellar remnant (in solar masses)
function M_rem=fun_M_rem(Z,M_zams)
-%We start by computing the M_co via fun_M_co:
+%We start by computing the M_co via fun_M_co.m:
mat_M_co=fun_M_co(Z,M_zams);
for jj=1:length(mat_M_co(:,1))
for ii=1:length(mat_M_co(jj,:))
diff --git a/results/cluster properties/Cluster properties/fun_fb.m b/results/cluster properties/Cluster properties/fun_fb.m
index 22ae0f5d3c466039671810b8b279363b33a4af67..4c480561ff5415488812f567c8ff9115c1ec4530 100644
--- a/results/cluster properties/Cluster properties/fun_fb.m
+++ b/results/cluster properties/Cluster properties/fun_fb.m
@@ -1,4 +1,10 @@
-%% f_fb and M_fb as a function of Z and M_zams
+%% FALLBACK FRACTION (f_fb) AS A FUNCTION OF Z AND M_zams
+%{
+Description:
+ This function computes the fallback fraction according to the
+ results presented by Spera, Mapelli & Bressan (2015), and Fryer
+ et al. (2012) for the delayed supernova mechanism.
+%}
%Input:
% Z:Metallicity
% M_zams:Zero age main sequence mass (in solar masses)
@@ -7,10 +13,10 @@
% f_fb: Fallback fraction
% M_fb: Fallback mass (in solar masses)
function f_fb=fun_fb(Z,M_zams)
-%We start by computing the M_co via fun_M_co
+%We start by computing the M_co via fun_M_co.m
vec_M_co=fun_M_co(Z,M_zams);
%We then compute the value of M_proto:
-M_proto=zeros(1,length(vec_M_co));
+M_proto=zeros(1,length(vec_M_co));
for ii=1:length(vec_M_co)
M_co=vec_M_co(ii);
if M_co<3.5
@@ -25,7 +31,7 @@ end
%Now that we have M_proto we can compute a2 and b2:
a2=(0.133)-((0.093)./(M_zams(ii)-M_proto(ii)));
b2=-(11)*a2+1;
-%And then, in order to obtain M_fb and f_fb:
+%And then, in order to obtain f_fb:
if M_co<2.5
M_fb(ii)=0.2;
f_fb(ii)=0;%I don't know if this is what I have to do
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 5c0a02a45b832a1562f3cac2f5a06201785b7114..1fbc4b7c9d5594545873b2a4b4098d566a4d4225 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
@@ -79,7 +79,7 @@ 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),1));
+ f_ret=(1-feval(@prob_escape,Z,vec_bh_ns,1.4,265,mat_v_esc(dd,ff),2));
for jj=1:length(vec_bh_ns)
if vec_bh_ns(jj)>M_z_min_bh_test
%k=1 ==>with fallback
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 c04c5b7c3b1b8c064ceb61767ffbf88e97d67393..8ef9563af804ed4fed7a05699a3fbcd1912d3fc2 100644
--- a/results/cluster properties/Cluster properties/fun_for_mass_histograms.m
+++ b/results/cluster properties/Cluster properties/fun_for_mass_histograms.m
@@ -103,23 +103,4 @@ ll=length(vec_m_final); %We initialize ll so that we fill the vector after that
pos_bh_final=vec_m_final>M_z_min_bh_test;
vec_m_bh_final=vec_m_final.*pos_bh_final;
vec_m_bh_final=nonzeros(vec_m_bh_final)';
-end
-%{
-M_tot_initial=0;
-tic
-jj=1;
-while M_tot_initial<10^5
- 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;
- jj=jj+1;
-end
-toc
-
-vec=[1 2 3 4 5];
-for vec<<4.5
- vec_f=vec;
-end
-%}
\ No newline at end of file
+end
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/fun_v_esc.m b/results/cluster properties/Cluster properties/fun_v_esc.m
index bad101a92443c3b28252b39ed8d0136955345953..720cc0933eb2c2a1bc198b2409d1c5cadad80c7b 100644
--- a/results/cluster properties/Cluster properties/fun_v_esc.m
+++ b/results/cluster properties/Cluster properties/fun_v_esc.m
@@ -1,4 +1,9 @@
%% FUNCTION THAT GIVES THE ESCAPE VELOCITY OF THE CLUSTER AS A FUNCTION OF DENSITY AND MASS OF THE CLUSTER
+%{
+Description:
+ This function computes the escape velocity of a cluster
+ according to the results presented by Antonini & Gieles (2019)
+%}
%Input:
% r_h: half-mass radius (in pc)
% M_cl: cluster mass (in solar masses)
@@ -7,11 +12,7 @@
% v_esc: escape velocity (in km/s)
function v_esc=fun_v_esc(r_h,M_cl)
-%rho_h=3*M_cl./(8*pi*(r_h).^3);
-%M_5=M_cl./(10^5);
-%rho_5=rho_h./(10^5);
f_c=1;
-%v_esc=50*((M_5).^1/3).*((rho_5).^1/6)*f_c;
v_esc=50*(sqrt(M_cl./(1e5))).*(sqrt(1./r_h))*((3/(8*pi))^(1/6))*f_c;
end
diff --git a/results/cluster properties/Cluster properties/histogram_fb_265.jpg b/results/cluster properties/Cluster properties/histogram_fb_265.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..afe9fde7d4bb053d5b3305904cdb57eeca86c6d9
Binary files /dev/null and b/results/cluster properties/Cluster properties/histogram_fb_265.jpg differ
diff --git a/results/cluster properties/Cluster properties/html/CLUSTER_PROPERTIES.pdf b/results/cluster properties/Cluster properties/html/CLUSTER_PROPERTIES.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..bf66dbf72d14f960a36c51ec00b4f1b43037525f
Binary files /dev/null and b/results/cluster properties/Cluster properties/html/CLUSTER_PROPERTIES.pdf differ
diff --git a/results/cluster properties/Cluster properties/prob_esc_ns.m b/results/cluster properties/Cluster properties/prob_esc_ns.m
index 8b3de3656d0520719f021fba5d548fb5df239618..0c0234ec3988582c607a7fae6dd68784d1809b62 100644
--- a/results/cluster properties/Cluster properties/prob_esc_ns.m
+++ b/results/cluster properties/Cluster properties/prob_esc_ns.m
@@ -1,7 +1,13 @@
%% FUNCTION THAT GIVES THE PROBABILITY OF v>vesc FOR A NEUTRON STAR (NS)
+%{
+Description:
+ This function computes the probability of escape for a given
+ NS, according to the kick prescription, considering a
+ Maxwellian distribution.
+%}
%Input:
% sig_ns: neutron star velocity dispersion (in km/s)
-% v_esc: escape velocity
+% v_esc: escape velocity (in km/s)
%Output:
% p: probability of v>vesc for a neutron star.
diff --git a/results/cluster properties/Cluster properties/prob_escape.m b/results/cluster properties/Cluster properties/prob_escape.m
index 7744be78c4400365a61292b421761eab07ddd46c..30d79a655244b3a7e9cb2348c66494b7042eb6d7 100644
--- a/results/cluster properties/Cluster properties/prob_escape.m
+++ b/results/cluster properties/Cluster properties/prob_escape.m
@@ -1,13 +1,22 @@
-%% FUNCTION THAT GIVES THE PROBABILITY OF v>vesc TAKING INTO ACCOUNT FALLBACK f_fb
+%% FUNCTION THAT GIVES THE PROBABILITY OF ESCAPE FOR A BLACK HOLE (BH)
+%{
+Description:
+ This function computes the probability of escape for a given
+ BH, according to the kick prescription, considering a
+ Maxwellian distribution.
+%}
%Input:
% Z: Metallicity
% 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
-% k: parameter that defines whether there is fallback or not.
-% k=0 ==> no fallback
-% k=1 ==> fallback considered
+% 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:
% p: probability of v>v_esc
@@ -15,12 +24,13 @@ function p=prob_escape(Z,M_zams,M_ns,sig_ns,v_esc,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_zams.
M_bh=fun_M_rem(Z,M_zams);
-%And then we can obtain the BH's velocity dispersion.
-sig_bh=(M_ns./M_bh).*sig_ns;
-%We define the fallback fraction:
+%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);
diff --git a/results/cluster properties/Cluster properties/readout_hdf5/readout_hdf5.m b/results/cluster properties/Cluster properties/readout_hdf5/readout_hdf5.m
new file mode 100644
index 0000000000000000000000000000000000000000..32b52e7edb1a024afec4aa178f2ebbcffbeb10df
--- /dev/null
+++ b/results/cluster properties/Cluster properties/readout_hdf5/readout_hdf5.m
@@ -0,0 +1,14 @@
+hdf5_filename = '/user/HS103/m13239/Desktop/results_nbody6/snap.40_0.h5part'
+
+info = h5info(hdf5_filename, '/Step#0')
+data = h5read(hdf5_filename, '/Step#0/KW')
+
+%disp(data)
+
+disp(h5read(hdf5_filename, '/Step#0/M'))
+%h5disp(hdf5_filename, '/Step#0/')
+
+%h5disp(hdf5_filename)
+
+%level2 = info.Groups(1)
+%disp(level2)
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/readout_hdf5/test_test_3.m b/results/cluster properties/Cluster properties/readout_hdf5/test_test_3.m
new file mode 100644
index 0000000000000000000000000000000000000000..b4121ce846165f8b65a0af146d7382d5418b635c
--- /dev/null
+++ b/results/cluster properties/Cluster properties/readout_hdf5/test_test_3.m
@@ -0,0 +1,19 @@
+clear all; close all; format compact; format long g;
+%% SCRIPT TO TEST HOW TO EXTRACT INFORMATION FROM THE SNAP.40_0.h5part
+%We start by defining where the file is:
+hdf5_filename = '/user/HS103/m13239/Desktop/tests/test_3/snap.40_0.h5part';
+%Once we have this we can get some info on the file by doing:
+info = h5info(hdf5_filename, '/Step#0') %In this case we are obtaining information for the first step.
+%As we can see we have 15 data sets in this group.
+%We can be more specific, for example, to get the different masses by
+%doing:
+vec_m=h5read(hdf5_filename, '/Step#0/M');%As one can quickly see this command yields the different masses ordered in a decreasing manner.
+
+%And with this information we can easily plot a histogram:
+figure(1)
+h1=histogram(vec_m,45,'FaceColor','b');
+title('Number of stars for the cluster of test 3')
+xlabel('M (M_{sun})')
+ylabel('Number of stars')
+
+
diff --git a/results/cluster properties/Cluster properties/readout_hdf5/test_test_4.m b/results/cluster properties/Cluster properties/readout_hdf5/test_test_4.m
new file mode 100644
index 0000000000000000000000000000000000000000..8f34d09e8fbc9e1a8428bb55ed2009b81264c462
--- /dev/null
+++ b/results/cluster properties/Cluster properties/readout_hdf5/test_test_4.m
@@ -0,0 +1,17 @@
+clear all; close all; format compact; format long g;
+%% SCRIPT TO TEST HOW TO EXTRACT INFORMATION FROM THE SNAP.40_0.h5part
+%We start by defining where the file is:
+hdf5_filename = '/user/HS103/m13239/Desktop/tests/test_4/snap.40_0.h5part';
+%Once we have this we can get some info on the file by doing:
+info = h5info(hdf5_filename, '/Step#0') %In this case we are obtaining information for the first step.
+%As we can see we have 15 data sets in this group.
+%We can be more specific, for example, to get the different masses by
+%doing:
+vec_m=h5read(hdf5_filename, '/Step#0/M');%As one can quickly see this command yields the different masses ordered in a decreasing manner.
+
+%And with this information we can easily plot a histogram:
+figure(1)
+h1=histogram(vec_m,45,'FaceColor','b');
+title('Number of stars for the cluster of test 4')
+xlabel('M (M_{sun})')
+ylabel('Number of stars')
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/readout_hdf5/test_test_9.m b/results/cluster properties/Cluster properties/readout_hdf5/test_test_9.m
new file mode 100644
index 0000000000000000000000000000000000000000..3d7d8022bab5d90ec0537883d3009248f4c6dcfb
--- /dev/null
+++ b/results/cluster properties/Cluster properties/readout_hdf5/test_test_9.m
@@ -0,0 +1,111 @@
+clear all; close all; format compact; format long g;
+%% SCRIPT TO TEST HOW TO EXTRACT INFORMATION FROM THE SNAP.40_0.h5part
+%We start by defining where the file is:
+hdf5_filename = '/user/HS103/m13239/Desktop/tests/test_9/snap.40_0.h5part';
+%Once we have this we can get some info on the file by doing:
+info = h5info(hdf5_filename, '/Step#0');%In this case we are obtaining information for the first step.
+%As we can see we have 15 data sets in this group.
+%{
+%% HISTOGRAM FOR THE INITIAL MASSES
+%We can be more specific, for example, to get the different masses by
+%doing:
+vec_m_initial=h5read(hdf5_filename, '/Step#0/M');%As one can quickly see this command yields the different masses ordered in a decreasing manner.
+
+%And with this information we can easily plot a histogram:
+figure(1)
+h1=histogram(vec_m_initial,45,'FaceColor','b');
+title('Number of stars for N=5000,Z=0.0001 (IMF)')
+xlabel('M (M_{sun})')
+ylabel('Number of stars')
+set(gca, 'YScale', 'log')
+%}
+%% HISTOGRAM FOR THE MASSES AFTER 100 Myr (last Step)
+info = h5info(hdf5_filename);
+%We see that it has 274 groups. If we start by Step 0 that means that the
+%last step is Step 273:
+info = h5info(hdf5_filename, '/Step#273');
+%If we then want the dataset of the masses for that step:
+vec_m_final=h5read(hdf5_filename, '/Step#273/M');
+%And we can then plot the histogram:
+%{
+figure(2)
+h1=histogram(vec_m_final,45,'FaceColor','b');
+title('Number of stars for N=5000,Z=0.0001 (100 Myr after)')
+xlabel('M (M_{sun})')
+ylabel('Number of stars')
+set(gca, 'YScale', 'log')
+%}
+%% FINDING THE BLACK HOLES AND NEUTRON STARS
+%The next step is to see the number of black holes (BH) and neutron stars (NS) we have after 100
+%Myr and then plot a histogram for them.
+%To do that we need to find the dataset that stores the values for K*, that
+%is, the stellar type. In our case for now we care for BHs, NSs, and the
+%rest.
+%The values of interest of this K* are for us then:
+% 13 ==> NS
+% 14 ==> BH
+% the other stars
+
+%Then, to extract that information first we analyze what do we have for
+%this step:
+h5disp(hdf5_filename, '/Step#273');
+%After that we can see that this K* values are stored in KW. Then:
+vec_values_K=h5read(hdf5_filename, '/Step#273/KW');
+%As expected from the IMF (Kroupa 2001) the majority of stars are still
+%low-mass main sequence star (M<0.7), and that is why most K* values are
+%0.
+
+%As the stars are ordered in the same way for the different values we can
+%then, via these K* values, obtain the BH and NS masses.
+vec_m_bh_final=[];vec_m_ns_final=[];
+
+positions_bh_final=(vec_values_K==14); %1's for the stars that verify it, 0 for the others.
+vec_m_bh_final=positions_bh_final.*vec_m_final;
+%And now we eliminate the 0's:
+vec_m_bh_final=nonzeros(vec_m_bh_final)';
+N_bh_final=length(vec_m_bh_final);
+M_tot_bh_final=sum(vec_m_bh_final);
+
+%And we can do the same for the NSs:
+positions_ns_final=(vec_values_K==13);
+vec_m_ns_final=positions_ns_final.*vec_m_final;
+vec_m_ns_final=nonzeros(vec_m_ns_final)';
+N_ns_final=length(vec_m_ns_final);
+M_tot_ns_final=sum(vec_m_ns_final);
+
+%And if we want to find the total number of stars after 100 Myr we only
+%need to determine the length of the vector of final masses or values of
+%K*:
+N_tot_final=length(vec_m_final);
+%And the final total mass of the cluster:
+M_tot_final=sum(vec_m_final);
+
+%And if we then want the numbers and total mass of the stars that are
+%neither BHs nor NSs:
+N_tot_rest_final=N_tot_final-N_bh_final-N_ns_final;
+M_tot_rest_final=M_tot_final-M_tot_bh_final-M_tot_ns_final;
+
+
+%And we can then plot some pie charts:
+X(1)=M_tot_rest_final;X(2)=M_tot_bh_final;X(3)=M_tot_ns_final;
+%X=[M_tot_rest_final;M_tot_bh_final;M_tot_ns_final];
+Y=[N_tot_rest_final N_bh_final N_ns_final];
+%Before doing the pie charts it is important to notice that the values we
+%get from our datasets are 'single'. But to plot them in the pie chart we
+%need them to be 'double' (more precission but more memory required). To do
+%that for X and Y:
+X=double(X); Y=double(Y);
+figure(3)
+%labels = {'Others','BHs','NSs'};
+%ax1 = subplot(1,2,1);
+pie(X,[0 1 0])
+title('Mass of the cluster after 100 Myr');
+legend('Others','BHs','NSs')
+%{
+%It doesn't make a lot of sense to plot this numbers, because you cant
+really see anything of value due to the higher number of stars that do not
+become neither BHs nor NSs.
+ax2 = subplot(1,2,2);
+pie(ax2,Y,[0 1 0],labels)
+title(ax2,'Numbers of stars in the cluster after 100 Myr');
+%}
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/study_variation_of_parameters.m b/results/cluster properties/Cluster properties/study_variation_of_parameters.m
index 864129fb373489647a108f008f4dca09966e5e30..88d4cb3248fab2abd090635edd0cced6e6435f27 100644
--- a/results/cluster properties/Cluster properties/study_variation_of_parameters.m
+++ b/results/cluster properties/Cluster properties/study_variation_of_parameters.m
@@ -138,20 +138,20 @@ disp('And taking into account fallback:')
disp(T)
%}
-%% CONTOUR PLOT FOR THE MASS FRACTION OF BLACK HOLES AS A FUNCTION OF M_cl AND r_h:
-%{
+%% CONTOUR PLOT FOR THE NUMBER FRACTION OF BLACK HOLES AS A FUNCTION OF M_cl AND r_h:
+
Z=0.0001;
%We define the following vectors:
%vec_M_cl=[5*10^3 10^4 2*10^4 3*10^4 5*10^4 10^5 2*10^5 5*10^5 10^6];
-vec_M_cl=linspace(10^4,10^6,100);
+vec_M_cl=linspace(10^3,10^6,125);
%vec_M_cl=[10^3 2*10^3 3*10^3 4*10^3 5*10^3 10^4 2*10^4 5*10^4 10^5];
%vec_M_cl=[10^4 2*10^4 3*10^4 5*10^4 7.5*10^4 10^5 5*10^5];
%vec_r_h=[0.2 0.5 0.75 1 1.5 2 2.5 3 4];
%vec_r_h=[0.25 0.5 0.75 1 1.25 1.5 1.75 2];
%vec_v_esc=[10 15 20 25 30 35 40 45 50];
-vec_v_esc=linspace(5,60,100);
+vec_v_esc=linspace(5,60,125);
%And with this we form the meshgrid:
[M_cl,mat_v_esc]=meshgrid(vec_M_cl,vec_v_esc);
%And we need to precompute the probabilities so
@@ -164,7 +164,7 @@ r_h=(1/40)*((3/(8*pi))^(1/3))*((M_cl)./((mat_v_esc).^2));
figure(1)
[C,h]=contour(M_cl/(10^5),r_h,f_N_bh_relative,'ShowText','on');
clabel(C,h)
-title('BH retention fraction (Z=0.0001)')
+title('BH retention fraction (Z=0.0001) (with fb &