diff --git a/mcluster_script.sh b/mcluster_script.sh
index 60b333e9dd01cf9ffdca7bc350fa15d7955b0017..985606886bf4ee30bbf8f41e7362f6434b56f043 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 -N 1000 -R 2 -P 0 -f 1 -m 0.08 -m 150 -B 0 -Z 0.0001 -o"
+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 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"
 
 echo "Writing files to $MCLUSTER_CONFIG_DIR"
diff --git a/results/cluster properties/Cluster properties/CLUSTER_PROPERTIES.m b/results/cluster properties/Cluster properties/CLUSTER_PROPERTIES.m
index 1df557074ca57f839ffc7acb63dae05727686729..d770d6d6081ce16bd50465cf515f65eb7068d090 100644
--- a/results/cluster properties/Cluster properties/CLUSTER_PROPERTIES.m	
+++ b/results/cluster properties/Cluster properties/CLUSTER_PROPERTIES.m	
@@ -1,6 +1,32 @@
 %% 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:
+%{
+Description: 
+  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.
+%}
+%Input:
+%                Z: Metallicity
+%                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 CLUSTER_PROPERTIES(Z,k)
 if k==0
     kick='W/o fallback,reduced sigma';
@@ -21,62 +47,68 @@ v_esc=fun_v_esc(r_h,M_cl);
 %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;
+vec_M_bh=[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)
 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))])
+set(gca, 'YScale', 'log')
+legend(['M=' num2str(vec_M_bh(1))])
 
-p=prob_escape(Z,vec_M_zams(2),1.4,265,v_esc,k);f_ret=1-p;
+p=prob_escape(Z,vec_M_bh(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))])
+set(gca, 'YScale', 'log')
+legend(['M=' num2str(vec_M_bh(2))])
 
-p=prob_escape(Z,vec_M_zams(3),1.4,265,v_esc,k);f_ret=1-p;
+p=prob_escape(Z,vec_M_bh(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))])
+set(gca, 'YScale', 'log')
+legend(['M=' num2str(vec_M_bh(3))])
 
-p=prob_escape(Z,vec_M_zams(3),1.4,265,v_esc,k);f_ret=1-p;
+p=prob_escape(Z,vec_M_bh(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))])
+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:
 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')
+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:
 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')
+contour(M_bh,v_esc,p,'ShowText','on')
 title(['Probability of escape. ' num2str(kick) '. Z=' num2str(Z)])
 xlabel('M_{BH} (in M_{sun})')
 ylabel('Escape velocity (in km/s)')
-
-
+%set(gca, 'YScale', 'log')
+%}
 %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)
@@ -84,7 +116,7 @@ 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)
+figure(3)
 subplot(2,2,[1,3])
 plot(vec_M_zams,vec_M_co,'-k')
 %text(15,25,'Z=2e-2')
diff --git a/results/cluster properties/Cluster properties/CONTINUOUS_IMF.m b/results/cluster properties/Cluster properties/CONTINUOUS_IMF.m
new file mode 100644
index 0000000000000000000000000000000000000000..81d508a514d65f068a3d2b20799f4a9cd6814504
--- /dev/null
+++ b/results/cluster properties/Cluster properties/CONTINUOUS_IMF.m	
@@ -0,0 +1,102 @@
+%% STUDY OF NUMBER OF BHs VIA A CONTINUOUS IMF.
+%We rewrite this so it can be run from another script just selecting the
+%input variables:
+function [results_b_kicks,results_a_kicks]=CONTINUOUS_IMF(Z,v_esc,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 define in another script the values that we will use in order to stuy
+%this mass fractions:
+
+% MASS AND NUMBER FRACTION OF BHs BEFORE SUPERNOVA KICKS (AFTER FORMATION)
+%If we use a standard Kroupa IMF with M_min=0.1 solar masses and M_max=150
+%solar masses:
+%alfa1=0.3;alfa2=1.3;alfa3=2.3;
+%We then test it for different metallicites in order to plot the f_N and
+%f_M as a function of the metallicity:
+%vec_Z=[2e-4,2e-3,5e-3,2e-2];
+vec_f_N_b_kicks=zeros(1,length(Z));
+vec_f_m_b_kicks=zeros(1,length(Z));
+vec_M_z_min=zeros(1,length(Z));
+for ii=1:length(Z)
+    Zi=Z(ii);
+    [vec_M_z_min(ii),vec_f_N_b_kicks(ii),vec_f_m_b_kicks(ii)]=fraction_before_kicks(Zi,alfa1,alfa2,alfa3);
+end
+disp('FRACTIONS BEFORE SUPERNOVA KICKS')
+disp( '                      Z              f_N_b_kicks               f_m_b_kicks                 M_z,min')
+results_b_kicks=[Z',vec_f_N_b_kicks',vec_f_m_b_kicks',vec_M_z_min'];
+disp(results_b_kicks)
+
+    
+% MASS AND NUMBER FRACTION OF BHs AFTER SUPERNOVA KICKS  (AS A FUNCTION OF THE CURRENT MASS OF THE CLUSTER)
+%We choose an standard escape velocity of 40 km/s:
+%v_esc=20;
+%We then test it for different metallicites in order to plot the f_N and
+%f_M as a function of the metallicity:
+vec_f_N_a_kicks=zeros(1,length(Z));
+vec_f_m_a_kicks=zeros(1,length(Z));
+vec_M_z_min=zeros(1,length(Z));
+for ii=1:length(Z)
+    Zi=Z(ii);
+    [vec_M_z_min(ii),vec_f_N_a_kicks(ii),vec_f_m_a_kicks(ii)]=fraction_after_kicks(Zi,alfa1,alfa2,alfa3,v_esc,k);
+end
+disp('FRACTIONS (AS FUNCTIONS OF THE VALUES OF N AND M FOR THE CURRENT CLUSTER)')
+disp( '                      Z              f_N_a_kicks               f_m_a_kicks                 M_z,min')
+results_a_kicks=[Z',vec_f_N_a_kicks',vec_f_m_a_kicks',vec_M_z_min'];
+disp(results_a_kicks)
+
+%{
+% MASS AND NUMBER FRACTION OF BHs AFTER SUPERNOVA KICKS  (AS A FUNCTION OF THE CURRENT MASS OF THE CLUSTER)
+%We choose an standard escape velocity of 40 km/s:
+v_esc=20;
+%We then test it for different metallicites in order to plot the f_N and
+%f_M as a function of the metallicity:
+vec_Z=[2e-4,2e-3,5e-3,2e-2];
+vec_f_N_a_kicks=zeros(1,length(vec_Z));
+vec_f_m_a_kicks=zeros(1,length(vec_Z));
+vec_M_z_min=zeros(1,length(vec_Z));
+for ii=1:length(vec_Z)
+    Z=vec_Z(ii);
+    [vec_M_z_min(ii),vec_f_N_a_kicks(ii),vec_f_m_a_kicks(ii)]=fraction_after_kicks_initial(Z,alfa1,alfa2,alfa3,v_esc);
+end
+disp('FRACTIONS (AS FUNCTIONS OF THE VALUES OF N AND M FOR THE INITIAL CLUSTER)')
+disp( '                      Z              f_N_a_kicks               f_m_a_kicks                 M_z,min')
+results_a_kicks=[vec_Z',vec_f_N_a_kicks',vec_f_m_a_kicks',vec_M_z_min'];
+disp(results_a_kicks)
+%}
+
+% PLOTS OF FRACTIONS AS A FUNCTION OF METALLICITY FOR A KROUPA IMF
+vec_Z=linspace(1e-3,1e-2,200);
+for ii=1:length(vec_Z)
+    Z=vec_Z(ii);
+    [vec_M_z_min(ii),vec_f_N_b_kicks(ii),vec_f_m_b_kicks(ii)]=fraction_before_kicks(Z,alfa1,alfa2,alfa3);
+    [vec_M_z_min(ii),vec_f_N_a_kicks(ii),vec_f_m_a_kicks(ii)]=fraction_after_kicks(Z,alfa1,alfa2,alfa3,v_esc,k);
+end
+figure(1);
+loglog(vec_Z,vec_f_N_b_kicks,'-r')
+xlim([1e-3,1e-2])
+title('Number fraction of stars with M<5 M_{sun}')
+xlabel('Z')
+ylabel('Number fraction')
+hold on
+loglog(vec_Z,vec_f_N_a_kicks,'-b')
+legend('f_N for stars that will be BHs','f_N after kicks (for current GC)')
+
+figure(2)
+semilogx(vec_Z,vec_f_m_b_kicks,'-r')
+xlim([1e-3,1e-2])
+title('Mass fraction of stars with M<5 M_{sun}')
+xlabel('Z')
+ylabel('Mass fraction')
+hold on
+semilogx(vec_Z,vec_f_m_a_kicks,'-b')
+hold on
+legend('f_m for stars that will be BHs','f_m after kicks (for current GC)')
+end
+
+
+
diff --git a/results/cluster properties/Cluster properties/FB_REDUCED_SIGMA.fig b/results/cluster properties/Cluster properties/FB_REDUCED_SIGMA.fig
new file mode 100644
index 0000000000000000000000000000000000000000..9d3baecfeb7934078bf4eb06ff67edc95f52afb1
Binary files /dev/null and b/results/cluster properties/Cluster properties/FB_REDUCED_SIGMA.fig differ
diff --git a/results/cluster properties/Cluster properties/SAMPLING.m b/results/cluster properties/Cluster properties/SAMPLING.m
index 1bc25d3e9e5210b6b35476df458601eae59f5892..e8b0f994fc53d5a53098cd8845c8b6c127f602b4 100644
--- a/results/cluster properties/Cluster properties/SAMPLING.m	
+++ b/results/cluster properties/Cluster properties/SAMPLING.m	
@@ -1,14 +1,41 @@
-%% SAMPLING
+%% SCRIPT WITH EVERYTHING REGARDING THE RESULTS OBTAINED VIA SAMPLING
+%{
+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';
+    kick='W/ fallback,sigma not reduced';
 end
 %{
-We do the sampling by calculating the cumulative probability density
+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
@@ -43,20 +70,17 @@ for ii=1:length(vec_M_examples)
 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;
+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')
@@ -64,9 +88,18 @@ 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') 
+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.
@@ -75,21 +108,21 @@ legend('IMF','after kicks')
 %variations.
 
 
-N_tot=[];
-vec_ii=[1:10];
+vec_ii=[1:100];
 mat_m=[];
+N_tot=zeros(1,length(vec_ii));
+M_tot=zeros(1,length(vec_ii));
 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
+  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.
+  while M_tot(ii)<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;
+    M_tot(ii)=M_tot(ii)+m;
     mat_m(ii,jj)=m;
     vec_m(jj)=m;
     jj=jj+1;
@@ -112,9 +145,9 @@ 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_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-M_bh(ii)-M_stars(ii);
+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
@@ -124,7 +157,7 @@ for ii=1:length(vec_ii)
 %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. 
+pos_stars=mat_m(ii,:)<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=mat_m(ii,:).*pos_stars;
 %But we need to eliminate the 0's and only consider those that are
@@ -138,58 +171,153 @@ vec_bh_ns=nonzeros(vec_bh_ns_i)';%Vector that contains the initial stars that ar
 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)); 
+%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 after the explosion(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(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.
+         if vec_bh_ns(jj)>3
             rr=rand;
             if rr<=f_ret(jj)
-                vec_m_final(ii,ll)=vec_bh_ns(ii,jj);
-                ll=ll+1;
+                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;
+                ll=ll;nn=nn;
              end
-           else
-            %Again we compute the probability of it being retained:
-            %f_ret=(1-feval(@prob_esc_ns,265,v_esc(dd,ff)));
+         else
             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;
+                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;
+                ll=ll;nn=nn;
             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);
+     N_bh_final(ii)=sum(vec_bh_ns_final>3);
+     N_ns_final(ii)=sum(vec_bh_ns_final<3);
+     N_stars_final(ii)=N_tot_final(ii)-N_bh_final(ii)-N_ns_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);
+     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);
 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)])
+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')
 hold on
 plot(vec_ii,N_tot_final,'-r')
-legend(['IMF','after kicks,' num2str(kick)])
+legend('IMF (before kicks)',['after kicks,' num2str(kick)],'Location','southwest')
+legend('boxoff')
+
 figure(3)
 subplot(2,1,1)
 plot(vec_ii,N_bh,'-r')
+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')
 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')
+title('Fraction of BHs retained')
+xlabel('Different samplings')
+ylabel('Fraction of BHs retained')
+
+%And for these calculations we can study the stochastic nature, that is,
+%what is the mean value, standard deviations, etc:
+%IMF:
+disp('                       VARIATIONS IN THE NUMBER OF BHs')
+disp('Mean number of initial BHs:')
+mean_N_bh_initial=sum(N_bh)/length(vec_ii);
+disp(mean_N_bh_initial)
+var=(N_bh-mean_N_bh_initial).^2;
+sigma_N_bh_initial=sqrt(sum(var)/length(vec_ii));
+disp('Standard deviation in the initial number of BHs:')
+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_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(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('Mean number fraction of BHs retained:')
+mean_f_N_bh=sum(f_N_bh)/length(vec_ii);
+disp(mean_f_N_bh)
+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_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')
+%}
+%{
+%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_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.
+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_low_mass(M_cl,Z,r_h)
 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~
index 01fbfe27a4981404e943907c07ba53a619737d39..350e60b3806c56adcc03b4ac8dab30ad7a72f9b4 100644
--- a/results/cluster properties/Cluster properties/SAMPLING.m~	
+++ b/results/cluster properties/Cluster properties/SAMPLING.m~	
@@ -1,14 +1,41 @@
-%% SAMPLING
+%% SCRIPT WITH EVERYTHING REGARDING THE RESULTS OBTAINED VIA SAMPLING
+%{
+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';
+    kick='W/ fallback,sigma not reduced';
 end
 %{
-We do the sampling by calculating the cumulative probability density
+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
@@ -43,20 +70,17 @@ for ii=1:length(vec_M_examples)
 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;
+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')
@@ -64,9 +88,18 @@ 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') 
+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.
@@ -75,21 +108,21 @@ legend('IMF','after kicks')
 %variations.
 
 
-N_tot=[];
-vec_ii=[1:10];
+vec_ii=[1:100];
 mat_m=[];
+N_tot=zeros(1,length(vec_ii));
+M_tot=zeros(1,length(vec_ii));
 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
+  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.
+  while M_tot(ii)<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;
+    M_tot(ii)=M_tot(ii)+m;
     mat_m(ii,jj)=m;
     vec_m(jj)=m;
     jj=jj+1;
@@ -112,9 +145,9 @@ 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_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-M_bh(ii)-M_stars(ii);
+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
@@ -124,7 +157,7 @@ for ii=1:length(vec_ii)
 %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. 
+pos_stars=mat_m(ii,:)<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=mat_m(ii,:).*pos_stars;
 %But we need to eliminate the 0's and only consider those that are
@@ -138,58 +171,143 @@ vec_bh_ns=nonzeros(vec_bh_ns_i)';%Vector that contains the initial stars that ar
 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)); 
+%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 after the explosion(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(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.
+         if vec_bh_ns(jj)>3
             rr=rand;
             if rr<=f_ret(jj)
-                vec_m_final(ii,ll)=vec_bh_ns(ii,jj);
-                ll=ll+1;
+                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;
+                ll=ll;nn=nn;
              end
-           else
-            %Again we compute the probability of it being retained:
-            %f_ret=(1-feval(@prob_esc_ns,265,v_esc(dd,ff)));
+         else
             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;
+                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;
+                ll=ll;nn=nn;
             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);
+     N_tot_final(ii)=length(vec_m_final);
+     M_tot_final(ii)=sum(vec_m_final);
+     N_bh_final(ii)=sum(vec_bh_ns_final>3);
+     N_ns_final(ii)=sum(vec_bh_ns_final<3);
+     N_stars_final(ii)=N_tot_final(ii)-N_bh_final(ii)-N_ns_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);
+     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);
 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)])
+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')
 hold on
 plot(vec_ii,N_tot_final,'-r')
-legend(['IMF','after kicks,' num2str(kick)])
+legend('IMF (before kicks)',['after kicks,' num2str(kick)],'Location','southwest')
+legend('boxoff')
+
 figure(3)
 subplot(2,1,1)
 plot(vec_ii,N_bh,'-r')
+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')
 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')
+title('Fraction of BHs retained')
+xlabel('Different samplings')
+ylabel('Fraction of BHs retained')
+
+%And for these calculations we can study the stochastic nature, that is,
+%what is the mean value, standard deviations, etc:
+%IMF:
+disp('                       VARIATIONS IN THE NUMBER OF BHs')
+disp('Mean number of initial BHs:')
+mean_N_bh_initial=sum(N_bh)/length(vec_ii);
+disp(mean_N_bh_initial)
+var=(N_bh-mean_N_bh_initial).^2;
+sigma_N_bh_initial=sqrt(sum(var)/length(vec_ii));
+disp('Standard deviation in the initial number of BHs:')
+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_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(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('Mean number fraction of BHs retained:')
+mean_f_N_bh=sum(f_N_bh)/length(vec_ii);
+disp(mean_f_N_bh)
+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_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)
 
 end
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/WO_FB_REDUCED_SIGMA.fig b/results/cluster properties/Cluster properties/WO_FB_REDUCED_SIGMA.fig
new file mode 100644
index 0000000000000000000000000000000000000000..b6eb593fc5e8c40a3ae7a52d39ce4f6d84afacaf
Binary files /dev/null and b/results/cluster properties/Cluster properties/WO_FB_REDUCED_SIGMA.fig differ
diff --git a/results/cluster properties/Cluster properties/archive 2/FB_REDUCED_SIGMA.fig b/results/cluster properties/Cluster properties/archive 2/FB_REDUCED_SIGMA.fig
new file mode 100644
index 0000000000000000000000000000000000000000..980895a5553e2316a3359dcddcf790d0a8ddbab7
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/FB_REDUCED_SIGMA.fig differ
diff --git a/results/cluster properties/Cluster properties/archive 2/FB_REDUCED_SIGMA.jpg b/results/cluster properties/Cluster properties/archive 2/FB_REDUCED_SIGMA.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..b7f0500115489034e80178aa9c106c8a931d355d
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/FB_REDUCED_SIGMA.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/MASS_FRACTION_WO_FALLBACK_REDUCED_SIGMA.jpg b/results/cluster properties/Cluster properties/archive 2/MASS_FRACTION_WO_FALLBACK_REDUCED_SIGMA.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..d436cd1671e0d145ec091212f16bd3698e4e38fe
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/MASS_FRACTION_WO_FALLBACK_REDUCED_SIGMA.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/MASS_FRACTION_W_FALLBACK_REDUCED_SIGMA.jpg b/results/cluster properties/Cluster properties/archive 2/MASS_FRACTION_W_FALLBACK_REDUCED_SIGMA.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..e72d061c705dbfdfc0c1442fc744e3ce9c7a2967
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/MASS_FRACTION_W_FALLBACK_REDUCED_SIGMA.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/MASS_FRACTION_W_FALLBACK_SIGMA_265.jpg b/results/cluster properties/Cluster properties/archive 2/MASS_FRACTION_W_FALLBACK_SIGMA_265.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..c2840c0e31f40482d56a6700e70cc9c32ca4c3f8
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/MASS_FRACTION_W_FALLBACK_SIGMA_265.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/NO_FB_REDUCED_SIGMA.fig b/results/cluster properties/Cluster properties/archive 2/NO_FB_REDUCED_SIGMA.fig
new file mode 100644
index 0000000000000000000000000000000000000000..6d37787e589fabebb62f1ce614107d8375f619cf
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/NO_FB_REDUCED_SIGMA.fig differ
diff --git a/results/cluster properties/Cluster properties/archive 2/NO_FB_REDUCED_SIGMA.jpg b/results/cluster properties/Cluster properties/archive 2/NO_FB_REDUCED_SIGMA.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..6ba871e643cc111a03adf4d478f846d457a3d894
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/NO_FB_REDUCED_SIGMA.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/N_bh_fb_not_reduced_sigma.jpg b/results/cluster properties/Cluster properties/archive 2/N_bh_fb_not_reduced_sigma.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..61823852948c0c80be4de72afeee1248ab7aea63
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/N_bh_fb_not_reduced_sigma.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/N_bh_fb_reducedd_sigma.jpg b/results/cluster properties/Cluster properties/archive 2/N_bh_fb_reducedd_sigma.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..209a9e1bd70c04aca0caa3c4d05c169422be0c57
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/N_bh_fb_reducedd_sigma.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/N_stars_fb_not_reduced_sigma.jpg b/results/cluster properties/Cluster properties/archive 2/N_stars_fb_not_reduced_sigma.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..73ea59f136c4c59c1dcb0629df0a317dc9219af1
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/N_stars_fb_not_reduced_sigma.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/N_stars_fb_reduced_sigma.jpg b/results/cluster properties/Cluster properties/archive 2/N_stars_fb_reduced_sigma.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..fb8bc13656412db4e39d03aff0e36585961d74a4
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/N_stars_fb_reduced_sigma.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/RETENTION_WO_FALLBACK_REDUCED_SIGMA.jpg b/results/cluster properties/Cluster properties/archive 2/RETENTION_WO_FALLBACK_REDUCED_SIGMA.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..6e7d8a1bf904e818c09318b036c4bc52b88290a8
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/RETENTION_WO_FALLBACK_REDUCED_SIGMA.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/RETENTION_W_FALLBACK_REDUCED_SIGMA.jpg b/results/cluster properties/Cluster properties/archive 2/RETENTION_W_FALLBACK_REDUCED_SIGMA.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..377ef5c0ba20ffffb2426ea2a1b9afa862562aeb
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/RETENTION_W_FALLBACK_REDUCED_SIGMA.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/RETENTION_W_FALLBACK_SIGMA_265.jpg b/results/cluster properties/Cluster properties/archive 2/RETENTION_W_FALLBACK_SIGMA_265.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..f693c40c7fe865245e93dfac9f22df69f9024c74
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/RETENTION_W_FALLBACK_SIGMA_265.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/STANDARD.fig b/results/cluster properties/Cluster properties/archive 2/STANDARD.fig
new file mode 100644
index 0000000000000000000000000000000000000000..5b310cd18093a1f72168d2bf08a5e5ab37c921d4
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/STANDARD.fig differ
diff --git a/results/cluster properties/Cluster properties/archive 2/STANDARD.jpg b/results/cluster properties/Cluster properties/archive 2/STANDARD.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..e341403c10e931b5ea3a352e575f75e934e5f82d
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/STANDARD.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/histogram_fb_not_reduced_sigma.jpg b/results/cluster properties/Cluster properties/archive 2/histogram_fb_not_reduced_sigma.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..cb086688e247ff66df3e3ea0d7e1e971fb9f41a0
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/histogram_fb_not_reduced_sigma.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/histogram_fb_reduced_sigma.jpg b/results/cluster properties/Cluster properties/archive 2/histogram_fb_reduced_sigma.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..116cc8889fea6ead8d77f67ae27a54280f078549
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/histogram_fb_reduced_sigma.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/retention_fraction_fb_not_reduced_sigma.jpg b/results/cluster properties/Cluster properties/archive 2/retention_fraction_fb_not_reduced_sigma.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..19aaae3de6000aa9338d24bfb7154e42c300e41c
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/retention_fraction_fb_not_reduced_sigma.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/retention_fraction_fb_not_reduced_sigma_2.jpg b/results/cluster properties/Cluster properties/archive 2/retention_fraction_fb_not_reduced_sigma_2.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..582fff81eab587ec67b75817e4a2a45928d927d8
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/retention_fraction_fb_not_reduced_sigma_2.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/retention_fraction_fb_not_reduced_sigma_3.jpg b/results/cluster properties/Cluster properties/archive 2/retention_fraction_fb_not_reduced_sigma_3.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..351bcefd3e856adcb7988c8513f5706978a1e22c
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/retention_fraction_fb_not_reduced_sigma_3.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/retention_fraction_fb_reduced_sigma.jpg b/results/cluster properties/Cluster properties/archive 2/retention_fraction_fb_reduced_sigma.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..1cf2c571c04b1425bae49b39ccbf775ef151aed3
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/retention_fraction_fb_reduced_sigma.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive 2/retention_with_rh.jpg b/results/cluster properties/Cluster properties/archive 2/retention_with_rh.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..eae03b2601f8e5a8a36ca653eae8376fb0a9b905
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive 2/retention_with_rh.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive/N_bh_fb_reduced_sigma.jpg b/results/cluster properties/Cluster properties/archive/N_bh_fb_reduced_sigma.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..a1c2a4dc089155c739ba5051d3b1d5e1fc7630f2
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive/N_bh_fb_reduced_sigma.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive/N_bh_no_fb.jpg b/results/cluster properties/Cluster properties/archive/N_bh_no_fb.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..b6fe86096bd0f47fc646543ef5f88cf022ea4f9f
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive/N_bh_no_fb.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive/N_stars_bray_eldridge.jpg b/results/cluster properties/Cluster properties/archive/N_stars_bray_eldridge.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..25d2db0717106a36166449ea9bd9ed9c2a9ac091
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive/N_stars_bray_eldridge.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive/N_stars_fb_reduced_sigma.jpg b/results/cluster properties/Cluster properties/archive/N_stars_fb_reduced_sigma.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..47b55d19da1d534cd6aed57454150340f37af287
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive/N_stars_fb_reduced_sigma.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive/N_stars_no_fb.jpg b/results/cluster properties/Cluster properties/archive/N_stars_no_fb.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..f728361e7be37cbbbd886032a210301b323bdd8e
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive/N_stars_no_fb.jpg differ
diff --git a/results/cluster properties/Cluster properties/bh_retention_fraction.jpg b/results/cluster properties/Cluster properties/archive/bh_retention_fraction.jpg
similarity index 100%
rename from results/cluster properties/Cluster properties/bh_retention_fraction.jpg
rename to results/cluster properties/Cluster properties/archive/bh_retention_fraction.jpg
diff --git a/results/cluster properties/Cluster properties/bh_retention_fraction_Z_0001.jpg b/results/cluster properties/Cluster properties/archive/bh_retention_fraction_Z_0001.jpg
similarity index 100%
rename from results/cluster properties/Cluster properties/bh_retention_fraction_Z_0001.jpg
rename to results/cluster properties/Cluster properties/archive/bh_retention_fraction_Z_0001.jpg
diff --git a/results/cluster properties/Cluster properties/bh_retention_fraction_fb_265.jpg b/results/cluster properties/Cluster properties/archive/bh_retention_fraction_fb_265.jpg
similarity index 100%
rename from results/cluster properties/Cluster properties/bh_retention_fraction_fb_265.jpg
rename to results/cluster properties/Cluster properties/archive/bh_retention_fraction_fb_265.jpg
diff --git a/results/cluster properties/Cluster properties/bh_retention_fraction_fb_265jpg.jpg b/results/cluster properties/Cluster properties/archive/bh_retention_fraction_fb_265jpg.jpg
similarity index 100%
rename from results/cluster properties/Cluster properties/bh_retention_fraction_fb_265jpg.jpg
rename to results/cluster properties/Cluster properties/archive/bh_retention_fraction_fb_265jpg.jpg
diff --git a/results/cluster properties/Cluster properties/bh_retention_fraction_low_mass_with_fallback.jpg b/results/cluster properties/Cluster properties/archive/bh_retention_fraction_low_mass_with_fallback.jpg
similarity index 100%
rename from results/cluster properties/Cluster properties/bh_retention_fraction_low_mass_with_fallback.jpg
rename to results/cluster properties/Cluster properties/archive/bh_retention_fraction_low_mass_with_fallback.jpg
diff --git a/results/cluster properties/Cluster properties/bh_retention_fraction_low_mass_wo_fallback.jpg b/results/cluster properties/Cluster properties/archive/bh_retention_fraction_low_mass_wo_fallback.jpg
similarity index 100%
rename from results/cluster properties/Cluster properties/bh_retention_fraction_low_mass_wo_fallback.jpg
rename to results/cluster properties/Cluster properties/archive/bh_retention_fraction_low_mass_wo_fallback.jpg
diff --git a/results/cluster properties/Cluster properties/bh_retention_fraction_wo_fallback.jpg b/results/cluster properties/Cluster properties/archive/bh_retention_fraction_wo_fallback.jpg
similarity index 100%
rename from results/cluster properties/Cluster properties/bh_retention_fraction_wo_fallback.jpg
rename to results/cluster properties/Cluster properties/archive/bh_retention_fraction_wo_fallback.jpg
diff --git a/results/cluster properties/Cluster properties/archive/black_hole_retention_bray_eldridge.jpg b/results/cluster properties/Cluster properties/archive/black_hole_retention_bray_eldridge.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..16cef2e5deb47d16dd2f213cca0a6fb430663d2a
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive/black_hole_retention_bray_eldridge.jpg differ
diff --git a/results/cluster properties/Cluster properties/histogram_1.jpg b/results/cluster properties/Cluster properties/archive/histogram_1.jpg
similarity index 100%
rename from results/cluster properties/Cluster properties/histogram_1.jpg
rename to results/cluster properties/Cluster properties/archive/histogram_1.jpg
diff --git a/results/cluster properties/Cluster properties/histogram_2.jpg b/results/cluster properties/Cluster properties/archive/histogram_2.jpg
similarity index 100%
rename from results/cluster properties/Cluster properties/histogram_2.jpg
rename to results/cluster properties/Cluster properties/archive/histogram_2.jpg
diff --git a/results/cluster properties/Cluster properties/histogram_3.jpg b/results/cluster properties/Cluster properties/archive/histogram_3.jpg
similarity index 100%
rename from results/cluster properties/Cluster properties/histogram_3.jpg
rename to results/cluster properties/Cluster properties/archive/histogram_3.jpg
diff --git a/results/cluster properties/Cluster properties/histogram_4.jpg b/results/cluster properties/Cluster properties/archive/histogram_4.jpg
similarity index 100%
rename from results/cluster properties/Cluster properties/histogram_4.jpg
rename to results/cluster properties/Cluster properties/archive/histogram_4.jpg
diff --git a/results/cluster properties/Cluster properties/histogram_fb_265.jpg b/results/cluster properties/Cluster properties/archive/histogram_fb_265.jpg
similarity index 100%
rename from results/cluster properties/Cluster properties/histogram_fb_265.jpg
rename to results/cluster properties/Cluster properties/archive/histogram_fb_265.jpg
diff --git a/results/cluster properties/Cluster properties/archive/histogram_fb_reduced_sigma.jpg b/results/cluster properties/Cluster properties/archive/histogram_fb_reduced_sigma.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..cb9bb7799c1c8b62c9b80f236623ae50ae463a0e
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive/histogram_fb_reduced_sigma.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive/histogram_no_fb.jpg b/results/cluster properties/Cluster properties/archive/histogram_no_fb.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..7f8087c6e0728112e6ad6212925b9e99dc8077da
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive/histogram_no_fb.jpg differ
diff --git a/results/cluster properties/Cluster properties/mass_fraction_black_holes.fig b/results/cluster properties/Cluster properties/archive/mass_fraction_black_holes.fig
similarity index 100%
rename from results/cluster properties/Cluster properties/mass_fraction_black_holes.fig
rename to results/cluster properties/Cluster properties/archive/mass_fraction_black_holes.fig
diff --git a/results/cluster properties/Cluster properties/mass_fraction_black_holes.jpg b/results/cluster properties/Cluster properties/archive/mass_fraction_black_holes.jpg
similarity index 100%
rename from results/cluster properties/Cluster properties/mass_fraction_black_holes.jpg
rename to results/cluster properties/Cluster properties/archive/mass_fraction_black_holes.jpg
diff --git a/results/cluster properties/Cluster properties/number_of_black_holes.fig b/results/cluster properties/Cluster properties/archive/number_of_black_holes.fig
similarity index 100%
rename from results/cluster properties/Cluster properties/number_of_black_holes.fig
rename to results/cluster properties/Cluster properties/archive/number_of_black_holes.fig
diff --git a/results/cluster properties/Cluster properties/number_of_black_holes.jpg b/results/cluster properties/Cluster properties/archive/number_of_black_holes.jpg
similarity index 100%
rename from results/cluster properties/Cluster properties/number_of_black_holes.jpg
rename to results/cluster properties/Cluster properties/archive/number_of_black_holes.jpg
diff --git a/results/cluster properties/Cluster properties/number_of_stars.fig b/results/cluster properties/Cluster properties/archive/number_of_stars.fig
similarity index 100%
rename from results/cluster properties/Cluster properties/number_of_stars.fig
rename to results/cluster properties/Cluster properties/archive/number_of_stars.fig
diff --git a/results/cluster properties/Cluster properties/number_of_stars.jpg b/results/cluster properties/Cluster properties/archive/number_of_stars.jpg
similarity index 100%
rename from results/cluster properties/Cluster properties/number_of_stars.jpg
rename to results/cluster properties/Cluster properties/archive/number_of_stars.jpg
diff --git a/results/cluster properties/Cluster properties/archive/retention_fraction_fb_reduced_sigma.jpg b/results/cluster properties/Cluster properties/archive/retention_fraction_fb_reduced_sigma.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..a3cc820d3cdb988532fe074cd3328440f1f588d8
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive/retention_fraction_fb_reduced_sigma.jpg differ
diff --git a/results/cluster properties/Cluster properties/archive/retention_fraction_no_fb.jpg b/results/cluster properties/Cluster properties/archive/retention_fraction_no_fb.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..1b0f24a14bdf41a20836c18367ca025d0ec9e57a
Binary files /dev/null and b/results/cluster properties/Cluster properties/archive/retention_fraction_no_fb.jpg 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 1fbc4b7c9d5594545873b2a4b4098d566a4d4225..c8fb5a15a53eb077e09a06c58b4ab2a79ae3a0c8 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	
@@ -1,9 +1,18 @@
 %% FUNCTION THAT WILL GIVE US WHAT WE NEED IN ORDER TO FORM THE CONTOUR PLOT:
 %{
+Description: 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.
+%}
+%{
 Inputs:
        M_cl: Mass of the cluster (in solar masses)
        Z: Metallicity
        mat_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:
         f_N_bh_relative: number fraction of black holes retained compared
         to the initial number of black holes (as defined by the IMF)
@@ -14,9 +23,9 @@ DONE)
 For more information on the sampling and how we compute the retained stars see sampling.m
 %}
 
-function f_N_bh_relative=fun_for_contour_fractions_3(M_cl,Z,mat_v_esc)
+function f_N_bh_relative=fun_for_contour_fractions_3(M_cl,Z,r_h,k)
 %We start by sampling the IMF. And for that we need to load the cpdf
-%(cumulative probability density function):
+%(cumulative probability distribution function):
 
 load('cpdf')
 
@@ -50,16 +59,16 @@ M_z_min_bh_test=vec_M_zams_test(ll(1));
 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:
+%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
 %stars are retained and which ones are not:
-
-for dd=1:length(mat_v_esc(1,:))
-for ff=1:length(mat_v_esc(:,1))
+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  
 %In order to not use so many if statements we are going to find which
@@ -76,60 +85,48 @@ vec_m_final=nonzeros(vec_m_final_i)';
 %probabilities for them to be kicked out are:
 vec_bh_ns=vec_m_initial(ff,:)-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,mat_v_esc(dd,ff),2)); 
+    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
+    %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 after the explosion(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)>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.
+         if vec_bh_ns(jj)>3
             rr=rand;
             if rr<=f_ret(jj)
                 vec_m_final(ll)=vec_bh_ns(jj);
-                ll=ll+1;
+                vec_bh_ns_final(nn)=vec_bh_ns(jj);
+                ll=ll+1;nn=nn+1;
              else
-                ll=ll;
+                ll=ll;nn=nn;
              end
-           else
-            %Again we compute the probability of it being retained:
-            %f_ret=(1-feval(@prob_esc_ns,265,v_esc(dd,ff)));
+         else
             rr=rand;
             if rr<=(1-feval(@prob_esc_ns,265,mat_v_esc(dd,ff)))
                 vec_m_final(ll)=vec_bh_ns(jj);
-                ll=ll+1;
+                vec_bh_ns_final(nn)=vec_bh_ns(jj);
+                ll=ll+1;nn=nn+1;
             else
-                ll=ll;
+                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,ff)=sum(vec_m_final>M_z_min_bh_test);
+    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
-%{
-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
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
new file mode 100644
index 0000000000000000000000000000000000000000..52ad520762d8d9ba3d7ab5039eff9d1107e76363
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fun_for_contour_fractions_all_plots.m	
@@ -0,0 +1,168 @@
+%% FUNCTION THAT WILL GIVE US WHAT WE NEED IN ORDER TO FORM THE CONTOUR PLOT FOR ALL PRESCRIPTIONS AT ONCE
+%{
+Description: 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.
+%}
+%{
+Inputs:
+       M_cl: Mass of the cluster (in solar masses)
+       Z: Metallicity
+       mat_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:
+        f_N_bh_relative: number fraction of black holes retained compared
+        to the initial number of black holes (as defined by the IMF)
+
+IMPORTANT: THIS CODE ONLY WORKS IF THE MATRIX FOR THE POSSIBLE M_cl AND
+v_esc HAVE THE SAME SIZE. (SO IN THE END ONLY WORKS AFTER THE MESHGRID IS
+DONE)
+For more information on the sampling and how we compute the retained stars see sampling.m
+%}
+
+function f_N_bh_relative=fun_for_contour_fractions_all_plots(M_cl,Z,r_h)
+%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 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);
+        M_tot_initial=0;
+jj=1;%We initialize this variable before starting the while loop.
+  while M_tot_initial<vec_M_cl(s)
+    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(s,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(s,:)=length(vec_m_initial(s,:));%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(s)=sum(vec_m_initial(s,:)>M_z_min_bh_test);
+toc
+ end
+ 
+for k=0:2
+   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
+%Now, we take the values we obtained before and use them to determine which
+%stars are retained and which ones are not:
+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.
+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(ff,:)<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(ff,:).*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(ff,:)-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,mat_v_esc(dd,ff),k)); %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<=(1-feval(@prob_esc_ns,265,mat_v_esc(dd,ff)))
+                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,ff)=sum(vec_bh_ns_final>3);
+    f_N_bh_relative(dd,ff)=N_bh_final(dd,ff)./N_bh_initial(ff);
+    %And to obtain the mass fractions:
+    vec_m_bh_final=vec_bh_ns_final(find(vec_bh_ns_final>3));
+    f_M_bh(dd,ff)=sum(vec_m_bh_final)/(M_cl(dd,ff));
+   toc
+end
+end
+figure()
+subplot(1,2,1)
+contour(M_cl/(10^4),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^4 M_{sun})')
+ylabel('r_h (in pc)')
+set(gca, 'YScale', 'log')
+%set(gca, 'XScale', 'log')
+hold on
+subplot(1,2,2)
+contour(M_cl/(10^4),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})')
+ylabel('r_h (in pc)')
+set(gca, 'YScale', 'log')
+%set(gca, 'XScale', 'log')
+hold off
+
+end
+end
diff --git a/results/cluster properties/Cluster properties/fun_for_contour_fractions_low_mass.m b/results/cluster properties/Cluster properties/fun_for_contour_fractions_low_mass.m
new file mode 100644
index 0000000000000000000000000000000000000000..b04c72a1c1e340bfbdc01160ed16b6d40c57d6a9
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fun_for_contour_fractions_low_mass.m	
@@ -0,0 +1,215 @@
+%% FUNCTION THAT WILL GIVE US WHAT WE NEED IN ORDER TO FORM THE CONTOUR PLOT FOR ALL PRESCRIPTIONS AT ONCE
+%% FOR THE MEAN VALUES OF LOW-MASS CLUSTERS
+%{
+Description: 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. Concretely it focuses on low-mass clusters, such that it
+computes the mean value for different samplings, as for this lower masses
+the sampling can really affect the results.
+%}
+%{
+Inputs:
+       M_cl: Mass of the cluster (in solar masses)
+       Z: Metallicity
+       mat_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:
+        f_N_bh_relative: number fraction of black holes retained compared
+        to the initial number of black holes (as defined by the IMF)
+
+IMPORTANT: THIS CODE ONLY WORKS IF THE MATRIX FOR THE POSSIBLE M_cl AND
+v_esc HAVE THE SAME SIZE. (SO IN THE END ONLY WORKS AFTER THE MESHGRID IS
+DONE)
+For more information on the sampling and how we compute the retained stars see sampling.m
+%}
+
+function [f_N_bh_mean,f_M_bh_mean]=fun_for_contour_fractions_low_mass(M_cl,Z,r_h)
+%We start by sampling the IMF. And for that we need to load the cpdf
+%(cumulative probability distribution function):
+vec_ii=[1:5];
+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);
+        M_tot_initial=0;
+jj=1;%We initialize this variable before starting the while loop.
+  while M_tot_initial<vec_M_cl(s)
+    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(s,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(s,:)=length(vec_m_initial(s,:));%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(s)=sum(vec_m_initial(s,:)>M_z_min_bh_test);
+toc
+ end
+%}
+ 
+for k=0:2
+f_N_bh_total=zeros(size(M_cl));f_M_bh_total=zeros(size(M_cl));
+   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
+ for ii=length(vec_ii) 
+     
+     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);
+        M_tot_initial=0;
+jj=1;%We initialize this variable before starting the while loop.
+  while M_tot_initial<vec_M_cl(s)
+    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(s,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(s,:)=length(vec_m_initial(s,:));%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(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
+%stars are retained and which ones are not:
+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.
+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(ff,:)<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(ff,:).*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(ff,:)-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,mat_v_esc(dd,ff),k)); %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<=(1-feval(@prob_esc_ns,265,mat_v_esc(dd,ff)))
+                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,ff)=sum(vec_bh_ns_final>3);
+    f_N_bh_relative(dd,ff)=N_bh_final(dd,ff)./N_bh_initial(ff);
+    %And to obtain the mass fractions:
+    vec_m_bh_final=vec_bh_ns_final(find(vec_bh_ns_final>3));
+    f_M_bh(dd,ff)=sum(vec_m_bh_final)/(M_cl(dd,ff));
+   toc
+end
+end
+f_N_bh_total=f_N_bh_total+f_N_bh_relative;
+f_M_bh_total=f_M_bh_total+f_M_bh;
+end
+%And now we compute the mean values:
+f_N_bh_mean=f_N_bh_total./(length(vec_ii));
+f_M_bh_mean=f_M_bh_total./(length(vec_ii));
+figure()
+subplot(1,2,1)
+contour(M_cl/(10^4),r_h,f_N_bh_mean);
+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})')
+ylabel('r_h (in pc)')
+set(gca, 'YScale', 'log')
+%set(gca, 'XScale', 'log')
+hold on
+subplot(1,2,2)
+contour(M_cl/(10^4),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})')
+ylabel('r_h (in pc)')
+set(gca, 'YScale', 'log')
+%set(gca, 'XScale', 'log')
+hold off
+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 8ef9563af804ed4fed7a05699a3fbcd1912d3fc2..3c4ab6dd8246545b06ce1a35591a263012904135 100644
--- a/results/cluster properties/Cluster properties/fun_for_mass_histograms.m	
+++ b/results/cluster properties/Cluster properties/fun_for_mass_histograms.m	
@@ -15,7 +15,7 @@ For more information on the sampling and how we compute the retained stars see s
 
 function [vec_m_bh_initial,vec_m_bh_final]=fun_for_mass_histograms(M_cl,Z,v_esc,k)
 %We start by sampling the IMF. And for that we need to load the cpdf
-%(cumulative probability density function):
+%(cumulative probability distribution function):
 
 load('cpdf')
 
@@ -63,44 +63,47 @@ vec_m_final=nonzeros(vec_m_final_i)';
 %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:
-pos_bh=vec_bh_ns>M_z_min_bh_test;
-vec_m_bh_initial=vec_bh_ns.*pos_bh;
-vec_m_bh_initial=nonzeros(vec_m_bh_initial)';
 
+%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:
     f_ret=(1-feval(@prob_escape,Z,vec_bh_ns,1.4,265,v_esc,k)); 
+%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 after the explosion(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)>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.
+         if vec_bh_ns(jj)>3
             rr=rand;
             if rr<=f_ret(jj)
                 vec_m_final(ll)=vec_bh_ns(jj);
-                ll=ll+1;
+                vec_bh_ns_final(nn)=vec_bh_ns(jj);
+                ll=ll+1;nn=nn+1;
              else
-                ll=ll;
+                ll=ll;nn=nn;
              end
-           else
-            %Again we compute the probability of it being retained:
-            %f_ret=(1-feval(@prob_esc_ns,265,v_esc(dd,ff)));
+         else
             rr=rand;
             if rr<=(1-feval(@prob_esc_ns,265,v_esc))
                 vec_m_final(ll)=vec_bh_ns(jj);
-                ll=ll+1;
+                vec_bh_ns_final=vec_bh_ns(jj);
+                ll=ll+1;nn=nn+1;
             else
-                ll=ll;
+                ll=ll;nn=nn;
             end
          end
-     end
      %And now we store the values of the masses of black holes retained
      %after the kicks:
-     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)';
+     vec_m_bh_final=vec_bh_ns_final(find(vec_bh_ns_final>3));
 end
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/plot_4_prescriptions b/results/cluster properties/Cluster properties/plot_4_prescriptions
new file mode 100644
index 0000000000000000000000000000000000000000..56b97a8b684d6bb8d0e0160dcde07a8f5f2de7d1
--- /dev/null
+++ b/results/cluster properties/Cluster properties/plot_4_prescriptions	
@@ -0,0 +1,12 @@
+clear all; close all; format compact; format long g;
+%% SCRIPT TO PLOT THE RETENTION FRACTION OF BHs FOR THE 4 DIFFERENT KICK PRESCRIPTIONS THAT WE HAVE
+%We will load the results from the calculations already done separately:
+load('study_variation_final_Z_0001_without_fallback')
+clear all
+load('study_variation_final_Z_0001')
+clear all
+load('study_variation_final_Z_0001_fb_265')
+%{
+clear all
+load('study_variation_final_Z_0001_tugboat')
+%}
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/plot_4_prescriptions.m b/results/cluster properties/Cluster properties/plot_4_prescriptions.m
new file mode 100644
index 0000000000000000000000000000000000000000..56b97a8b684d6bb8d0e0160dcde07a8f5f2de7d1
--- /dev/null
+++ b/results/cluster properties/Cluster properties/plot_4_prescriptions.m	
@@ -0,0 +1,12 @@
+clear all; close all; format compact; format long g;
+%% SCRIPT TO PLOT THE RETENTION FRACTION OF BHs FOR THE 4 DIFFERENT KICK PRESCRIPTIONS THAT WE HAVE
+%We will load the results from the calculations already done separately:
+load('study_variation_final_Z_0001_without_fallback')
+clear all
+load('study_variation_final_Z_0001')
+clear all
+load('study_variation_final_Z_0001_fb_265')
+%{
+clear all
+load('study_variation_final_Z_0001_tugboat')
+%}
\ No newline at end of file
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
new file mode 100644
index 0000000000000000000000000000000000000000..9d59f52dc5174b2ad7e4b5374803ad87d7aa8007
--- /dev/null
+++ b/results/cluster properties/Cluster properties/plots_M_rem_and_probs_report.m	
@@ -0,0 +1,80 @@
+%% CLUSTER PROPERTIES
+%{
+Description: 
+  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.
+  It also shows the probability of escape for a cluster of a certain mass,
+  half-mass radius and metallicity for different BH masses, for the
+  different prescriptions we have computed
+%}
+%}
+%Input:
+%                Z: Metallicity
+%                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 plots_M_rem_and_probs_report(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 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)
+plot(vec_M_zams,vec_M_co,'-k')
+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')
+
+%And now we look to plot the probabilities of escape for the three
+%different kick prescriptions we have computed:
+for k=0:2
+    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
+   p=prob_escape(Z,M_rem,M_ns,sig_ns,v_esc,k)
+
+
+
+end
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN.fig b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN.fig
new file mode 100644
index 0000000000000000000000000000000000000000..7a4a79483dd6681cf8e2b80b86f45e7ea35944f7
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN.png b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN.png
new file mode 100644
index 0000000000000000000000000000000000000000..095c5e9068558e90c0348ebe7921b71dfd6d43d6
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_log_space.fig b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_log_space.fig
new file mode 100644
index 0000000000000000000000000000000000000000..81cd1192d38ac24e744485b04620424bbbc46ea1
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_log_space.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_log_space.png b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_log_space.png
new file mode 100644
index 0000000000000000000000000000000000000000..15b5276c255bb21229510e4451b01c6b17bc2bc8
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_log_space.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass.fig b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass.fig
new file mode 100644
index 0000000000000000000000000000000000000000..49714ac09d7edda5e715b09062d52d9a3331eb05
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass.png b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass.png
new file mode 100644
index 0000000000000000000000000000000000000000..3190e768c1593435162e464af9c916861ddd806c
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass_log_space.fig b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass_log_space.fig
new file mode 100644
index 0000000000000000000000000000000000000000..bede43e16501a2be0ad16a37ebd508cf38826393
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass_log_space.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass_log_space.png b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass_log_space.png
new file mode 100644
index 0000000000000000000000000000000000000000..c3947b8086ebb1def35f7573270dbd2327e31798
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass_log_space.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass_mean.fig b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass_mean.fig
new file mode 100644
index 0000000000000000000000000000000000000000..38126370be97168ba90f6059a65d61dcd2382463
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass_mean.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass_mean.png b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass_mean.png
new file mode 100644
index 0000000000000000000000000000000000000000..254e598f1c2fb5bfd4df4e2779b153f2113a32fd
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/NEUTRINO_DRIVEN_low_mass_mean.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD.fig b/results/cluster properties/Cluster properties/plots_report/STANDARD.fig
new file mode 100644
index 0000000000000000000000000000000000000000..4a6d43c1e4a84cc50f1e1b7e8ac6a2a0584ba866
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD.png b/results/cluster properties/Cluster properties/plots_report/STANDARD.png
new file mode 100644
index 0000000000000000000000000000000000000000..7c8a004a5b3915497274c1b616423b3086a00d02
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD_log_space.fig b/results/cluster properties/Cluster properties/plots_report/STANDARD_log_space.fig
new file mode 100644
index 0000000000000000000000000000000000000000..6527e0a18f95916ea11e3aa2daf4adcd801cdf14
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD_log_space.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD_log_space.png b/results/cluster properties/Cluster properties/plots_report/STANDARD_log_space.png
new file mode 100644
index 0000000000000000000000000000000000000000..50fc8a8262a69493585f1ae864ed86300b117fe6
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD_log_space.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mas.fig b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mas.fig
new file mode 100644
index 0000000000000000000000000000000000000000..b4d87cf539956278d57821006e60e35ce78f09c8
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mas.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mas.png b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mas.png
new file mode 100644
index 0000000000000000000000000000000000000000..b6d747b34a7595da1de79eb80cdf78d36bcdb8f0
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mas.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mas_log_space.fig b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mas_log_space.fig
new file mode 100644
index 0000000000000000000000000000000000000000..a8cee51bf79cc347d023fe70240fa635f02dc19a
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mas_log_space.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mas_log_space.png b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mas_log_space.png
new file mode 100644
index 0000000000000000000000000000000000000000..3ea8d69eee80ae31200a5838a8e3cf433666fd62
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mas_log_space.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mass_mean.fig b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mass_mean.fig
new file mode 100644
index 0000000000000000000000000000000000000000..af38f201fb87a4066187b99057bc83cd9cf63367
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mass_mean.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mass_mean.png b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mass_mean.png
new file mode 100644
index 0000000000000000000000000000000000000000..446a65a94d54892117953f006faed37ab53d7e10
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/STANDARD_low_mass_mean.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA.fig b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA.fig
new file mode 100644
index 0000000000000000000000000000000000000000..c553cbfe8ee26fdb96b6dfe7db2d8cd40f698e66
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA.png b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA.png
new file mode 100644
index 0000000000000000000000000000000000000000..a50dc98b89e3d6b329c75feee7f55f92a8ba3659
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_log_space.fig b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_log_space.fig
new file mode 100644
index 0000000000000000000000000000000000000000..3f31e60167b8779fd37848f9114386778ed0e02f
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_log_space.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_log_space.png b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_log_space.png
new file mode 100644
index 0000000000000000000000000000000000000000..b7c87463cd690b0bf44ce565530947d2558e7ee3
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_log_space.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass.fig b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass.fig
new file mode 100644
index 0000000000000000000000000000000000000000..30223ee7fc095f1ee600c9957bc4af1c0f7773af
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass.png b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass.png
new file mode 100644
index 0000000000000000000000000000000000000000..2dcf239dc64f04adfec5c5bf33bcd941864bfaa5
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass_log_space.fig b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass_log_space.fig
new file mode 100644
index 0000000000000000000000000000000000000000..6551ff8f711c050948c7c3db8ce00ea976bd6114
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass_log_space.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass_log_space.png b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass_log_space.png
new file mode 100644
index 0000000000000000000000000000000000000000..27f690cb71c29308d2844606272e2a062dd1321f
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass_log_space.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass_mean.fig b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass_mean.fig
new file mode 100644
index 0000000000000000000000000000000000000000..01af95ffc33f76e29170b18e921cf5dec085482b
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass_mean.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass_mean.png b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass_mean.png
new file mode 100644
index 0000000000000000000000000000000000000000..058d7980066ccff1f1734afb79cf7c31cd111cf2
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/W_FB_REDUCED_SIGMA_low_mass_mean.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/archive/NEUTRINO_DRIVEN.jpg b/results/cluster properties/Cluster properties/plots_report/archive/NEUTRINO_DRIVEN.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..c05c93b1383e766f123fe5de63ebf5fae5a919fc
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/archive/NEUTRINO_DRIVEN.jpg differ
diff --git a/results/cluster properties/Cluster properties/plots_report/archive/STANDARD.jpg b/results/cluster properties/Cluster properties/plots_report/archive/STANDARD.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..da9d3763dce9ad6f100a353f0e67bd34c264bf77
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/archive/STANDARD.jpg differ
diff --git a/results/cluster properties/Cluster properties/plots_report/archive/W_FB_REDUCED_SIGMA.jpg b/results/cluster properties/Cluster properties/plots_report/archive/W_FB_REDUCED_SIGMA.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..8af931490f7d15e05e6515e22fb11592e0a7d48b
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/archive/W_FB_REDUCED_SIGMA.jpg differ
diff --git a/results/cluster properties/Cluster properties/plots_report/archive/v_esc_as_function_Mcl_rh.jpg b/results/cluster properties/Cluster properties/plots_report/archive/v_esc_as_function_Mcl_rh.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..9c7745e123d093202b49f79659fe30d03e285389
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/archive/v_esc_as_function_Mcl_rh.jpg differ
diff --git a/results/cluster properties/Cluster properties/plots_report/archive/v_esc_as_function_Mcl_rh.png b/results/cluster properties/Cluster properties/plots_report/archive/v_esc_as_function_Mcl_rh.png
new file mode 100644
index 0000000000000000000000000000000000000000..503906b4f768b73c6160873e0a222d2f8b550ab0
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/archive/v_esc_as_function_Mcl_rh.png differ
diff --git a/results/cluster properties/Cluster properties/plots_report/v_esc_as_function_Mcl_rh.fig b/results/cluster properties/Cluster properties/plots_report/v_esc_as_function_Mcl_rh.fig
new file mode 100644
index 0000000000000000000000000000000000000000..11396ce453d24c6292179e107d67511e128ec0c5
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/v_esc_as_function_Mcl_rh.fig differ
diff --git a/results/cluster properties/Cluster properties/plots_report/v_esc_as_function_Mcl_rh.png b/results/cluster properties/Cluster properties/plots_report/v_esc_as_function_Mcl_rh.png
new file mode 100644
index 0000000000000000000000000000000000000000..6b8b2bfc80de22ff07b4caff058ecc5ad8a1cb91
Binary files /dev/null and b/results/cluster properties/Cluster properties/plots_report/v_esc_as_function_Mcl_rh.png differ
diff --git a/results/cluster properties/Cluster properties/prob_escape.m b/results/cluster properties/Cluster properties/prob_escape.m
index 30d79a655244b3a7e9cb2348c66494b7042eb6d7..cf1096301b5c8ff2d9b23ef44c0925787a5e5b50 100644
--- a/results/cluster properties/Cluster properties/prob_escape.m	
+++ b/results/cluster properties/Cluster properties/prob_escape.m	
@@ -7,10 +7,10 @@ Description:
 %}
 %Input:
 %                Z: Metallicity
-%                M_zams:Zero age main sequence mass (in solar masses)
+%                M: 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
+%                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
@@ -22,7 +22,7 @@ Description:
 %                p: probability of v>v_esc
 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.
+%mass of the stellar remnant as a function of Z and M.
 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:
diff --git a/results/cluster properties/Cluster properties/readout_hdf5/test_again.m b/results/cluster properties/Cluster properties/readout_hdf5/test_again.m
new file mode 100644
index 0000000000000000000000000000000000000000..608b7d21e3dab8f7b4b9e3eb23c22c796e7f5a27
--- /dev/null
+++ b/results/cluster properties/Cluster properties/readout_hdf5/test_again.m	
@@ -0,0 +1,122 @@
+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_again/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 154 groups. If we start by Step 0 that means that the
+%last step is Step 153:
+info = h5info(hdf5_filename, '/Step#153');
+%If we then want the dataset of the masses for that step:
+vec_m_final=h5read(hdf5_filename, '/Step#153/M');
+%And we can then plot the histogram:
+
+figure(2)
+h2=histogram(vec_m_final,45,'FaceColor','b');
+title('Number of stars for M_cl=10000 M_{sun},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#153');
+%After that we can see that this K* values are stored in KW. Then:
+vec_values_K=h5read(hdf5_filename, '/Step#153/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');
+%}
+
+%It would also be interesting to plot a histogram for the masses of the
+%BHs:
+figure(4)
+h4=histogram(vec_m_bh_final,45,'FaceColor','b');
+title('Number of BHs for M_cl=10000 M_{sun},Z=0.0001 (100 Myr after)')
+xlabel('M (M_{sun})')
+ylabel('Number of BHs')
+set(gca, 'YScale', 'log')
+
+
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 88d4cb3248fab2abd090635edd0cced6e6435f27..96cc71b956adb8533642bca06a20aca9532ad23d 100644
--- a/results/cluster properties/Cluster properties/study_variation_of_parameters.m	
+++ b/results/cluster properties/Cluster properties/study_variation_of_parameters.m	
@@ -140,7 +140,7 @@ disp(T)
 
 %% CONTOUR PLOT FOR THE NUMBER FRACTION OF BLACK HOLES AS A FUNCTION OF M_cl AND r_h:
 
-Z=0.0001;
+Z=0.0001;k=1;
 
 %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];
@@ -158,7 +158,7 @@ vec_v_esc=linspace(5,60,125);
 %study_variation_of_parameters.m does not need to compute them in every
 %iteration
 %So that:
-f_N_bh_relative=fun_for_contour_fractions_3(M_cl,Z,mat_v_esc);
+f_N_bh_relative=fun_for_contour_fractions_3(M_cl,Z,mat_v_esc,k);
 %But we want to express it as a function of r_h and nor mat_v_esc
 r_h=(1/40)*((3/(8*pi))^(1/3))*((M_cl)./((mat_v_esc).^2));
 figure(1)
diff --git a/results/cluster properties/Cluster properties/test_bray_eldridge.m b/results/cluster properties/Cluster properties/test_bray_eldridge.m
new file mode 100644
index 0000000000000000000000000000000000000000..e3b8bd926359b30589f44745a343e0926029e0e3
--- /dev/null
+++ b/results/cluster properties/Cluster properties/test_bray_eldridge.m	
@@ -0,0 +1,273 @@
+%% TEST BRAY & ELDRIDGE (2018)
+function [f_N_bh_relative,M_cl,r_h]= test_bray_eldridge(Z,M_cl,r_h)
+load('cpdf')
+%{
+v_esc=fun_v_esc(r_h,M_cl);
+vec_ii=[1:300];
+mat_m=[];
+N_tot=zeros(1,length(vec_ii));
+M_tot=zeros(1,length(vec_ii));
+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.
+  while M_tot(ii)<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(ii)=M_tot(ii)+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(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:
+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.
+%The thing now is that we do not have a probability, but a function for the
+%velocity of the star, so we can compare directly if the star is ejected or
+%not.
+M_rem=fun_M_rem(Z,vec_bh_ns); M_fin=0.9519.*vec_bh_ns +1.45;
+M_eject=M_fin-M_rem;
+v_bh_ns=abs(100.*(M_eject./M_rem) -170);
+%But what we want to store in vec_m_final are the masses of the remnants,
+and not the ZAMS mass that they had. That's why we do:
+vec_bh_ns=fun_M_rem(Z,vec_bh_ns);
+%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;
+     for jj=1:length(v_bh_ns)
+         if v_bh_ns(jj)>v_esc
+                ll=ll;nn=nn;
+         else
+                vec_m_final(ll)=vec_bh_ns(jj);
+                vec_bh_ns_final(nn)=vec_bh_ns(jj);
+                ll=ll+1;nn=nn+1;  
+         end
+     end
+N_tot_final(ii)=length(vec_m_final);
+     M_tot_final(ii)=sum(vec_m_final);
+     N_bh_final(ii)=sum(vec_bh_ns_final>3);
+     N_ns_final(ii)=sum(vec_bh_ns_final<3);
+     N_stars_final(ii)=N_tot_final(ii)-N_bh_final(ii)-N_ns_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_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);
+end
+%And now we can plot this:
+figure(2)
+plot(vec_ii,N_tot,'-k')
+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')
+hold on
+plot(vec_ii,N_tot_final,'-r')
+legend('IMF (before kicks)','after kicks (Bray & Eldridge 2018)','Location','southwest')
+legend('boxoff')
+
+figure(3)
+subplot(2,1,1)
+plot(vec_ii,N_bh,'-r')
+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')
+hold on
+plot(vec_ii,N_bh_final,'-b')
+legend('IMF (before kicks)','after kicks (Bray & Eldridge 2018)','Location','southwest')
+legend('boxoff')
+subplot(2,1,2)
+plot(vec_ii,f_N_bh,'-k')
+title('Fraction of BHs retained')
+xlabel('Different samplings')
+ylabel('Fraction of BHs retained')
+
+%And for these calculations we can study the stochastic nature, that is,
+%what is the mean value, standard deviations, etc:
+%IMF:
+disp('                       VARIATIONS IN THE NUMBER OF BHs')
+disp('Mean number of initial BHs:')
+mean_N_bh_initial=sum(N_bh)/length(vec_ii);
+disp(mean_N_bh_initial)
+var=(N_bh-mean_N_bh_initial).^2;
+sigma_N_bh_initial=sqrt(sum(var)/length(vec_ii));
+disp('Standard deviation in the initial number of BHs:')
+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_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(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('Mean number fraction of BHs retained:')
+mean_f_N_bh=sum(f_N_bh)/length(vec_ii);
+disp(mean_f_N_bh)
+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 we rewrite what we have in order to be able to plot the contour
+%plot:
+%We define the following vectors:
+vec_M_cl=linspace(10^3,10^6,100);
+vec_v_esc=linspace(5,60,100);
+%And with this we form the meshgrid:
+[M_cl,mat_v_esc]=meshgrid(vec_M_cl,vec_v_esc);
+
+
+
+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);
+        M_tot_initial=0;
+jj=1;%We initialize this variable before starting the while loop.
+  while M_tot_initial<vec_M_cl(s)
+    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(s,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(s,:)=length(vec_m_initial(s,:));%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:
+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
+%stars are retained and which ones are not:
+
+for dd=1:length(mat_v_esc(1,:))
+for ff=1:length(mat_v_esc(:,1))
+vec_m_final=[];%We initialize the vector that will store the retained masses.
+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(ff,:)<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(ff,:).*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(ff,:)-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 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);
+%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.
+%But now we are not dealing with probabilities:
+M_rem=fun_M_rem(Z,vec_bh_ns); M_fin=0.9519.*vec_bh_ns +1.45;
+M_eject=M_fin-M_rem;
+v_bh_ns=abs(100.*(M_eject./M_rem) -170);
+%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;
+for jj=1:length(v_bh_ns)
+         if v_bh_ns(jj)>mat_v_esc(dd,ff)
+                ll=ll;nn=nn;
+         else
+                vec_m_final(ll)=vec_bh_ns(jj);
+                vec_bh_ns_final(nn)=vec_bh_ns(jj);
+                ll=ll+1;nn=nn+1;
+         end
+     end
+
+     %And now we store the values for the fraction of BHs retained in the
+     %cluster after considering the supernova kick:
+     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
+%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)
+[C,h]=contour(M_cl/(10^5),r_h,f_N_bh_relative,'ShowText','on');
+clabel(C,h)
+title(['BH retention fraction. Z=' num2str(Z) 'Bray & Eldridge'])
+xlabel('M_{cl} (10^5 M_{sun})')
+ylabel('r_h (in pc)')
+set(gca, 'YScale', 'log')
+
+end
\ No newline at end of file