diff --git a/.gitignore b/.gitignore
index c18dd8d83ceed1806b50b0aaa46beb7e335fff13..61c47565a27707bf4febc98a49f34b12470662cb 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1 +1,5 @@
 __pycache__/
+results/cluster properties/reading/*
+
+*.mat
+!cpdf.mat
diff --git a/results/cluster properties/Cluster properties/Cluster_analysis_final_version_2.m b/results/cluster properties/Cluster properties/Cluster_analysis_final_version_2.m
new file mode 100644
index 0000000000000000000000000000000000000000..c973651eb72d2d36e07646b002f0edac9a17ebcb
--- /dev/null
+++ b/results/cluster properties/Cluster properties/Cluster_analysis_final_version_2.m	
@@ -0,0 +1,36 @@
+clear all; close all; format compact; format long g;
+%% CLUSTER PROPERTIES
+%% INFORMATION ON INPUT AND OUTPUT VARIABLES
+%Input: 
+%      alfa1,alfa2,alfa3: exponents of the IMF.
+%      Z: metallicity
+%      v_esc: escape velocity (in km/s)
+%      k: parameter that defines whether there is fallback or not. 
+%        k=0 ==> no fallback
+%        k=1 ==> fallback considered
+%Output:
+%      results_b_kicks: mass and number fraction of stars that will be BHs before supernova kicks
+%      results_a_kicks: mass and number fraction of BHs after supernova kicks
+%% DEFINITION OF INPUT VARIABLES
+%For now there are certain values fixed, those are:
+%       sig_ns:velocity dispersion for neutron stars, with a value of 265km/s.
+%       M_ns: mass of the neutron star, with a value of 1.5 solar masses
+%       M_bh_min: minimum mass of a black hole, fixed at 3 solar masses.
+%                (we have to follow the SEVN code)
+
+alfa1=0.3;alfa2=1.3;alfa3=2.3;
+Z=[2e-4,5e-3,2e-2];
+v_esc=30; %One value at a time.
+k=0;
+%% MASS AND NUMBER FRACTIONS OF BHs IN THE CLUSTER
+if k==0
+    disp('                                        Without fallback')
+end
+if k==1
+    disp('                                        With fallback')
+end
+[results_b_kicks,results_a_kicks]=fun_cluster_fractions(alfa1,alfa2,alfa3,Z,v_esc,k);
+%% PROBABILITIES OF ESCAPE WITH FALLBACK AND M_co AND M_rem AS A FUNCTION OF M_zams FOR Z=2E-2 AND Z=1E-3
+%Cluster_properties_final_version
+
+
diff --git a/results/cluster properties/Cluster properties/Cluster_conditions.m b/results/cluster properties/Cluster properties/Cluster_conditions.m
new file mode 100644
index 0000000000000000000000000000000000000000..b23614988df7835bbecd86ad9aac811657a20bf9
--- /dev/null
+++ b/results/cluster properties/Cluster properties/Cluster_conditions.m	
@@ -0,0 +1,458 @@
+clear all; close all; format compact; format long g;
+%% PROBABILITY OF ESCAPE AS A FUNCTION OF BH MASS
+%We calculate the probability of escape for a BH assuming a Maxwellian
+%distribution with velocity dispersion sigma. Assuming same momentum for
+%both NS and BH the velocity dispersion for the BH is reduced. The function
+%used is:
+%{
+%Input:
+%                M_ns: neutron star mass (in solar masses)
+%                M_bh: black hole mass (in solar masses)
+%                sig_ns: neutron star velocity dispersion (in km/s)
+%                v_esc: escape velocity
+
+%Output:
+%                p: probability of v>vesc
+function p=prob_esc(M_ns,M_bh,sig_ns,v_esc)
+sig_bh=(M_ns./M_bh).*sig_ns;
+p=sqrt(2/pi)*(v_esc./sig_bh).*exp(-(v_esc.^2)./(2*(sig_bh).^2))+(1-erf(v_esc./(sqrt(2)*sig_bh)));
+end
+As it can be seen with the function we still haven't implemented fallback,
+that is probably the next step in improving this approximation.
+
+
+We assume a velocity dispersion of 265km/s as argued by Morscher, Umbreit, Farr & Rasio 2002. We would, in this case,
+get higher probabilities of escape for a given mass compared to the
+results obtained via a dispersion of 190km/s. The observational evidence
+for this value is found in "A statistical study of 233 pulsar proper
+motions" by Hobbs, Lorimer,Lyne and Kramer 2005
+We could asssume a velocity dispersion of 190km/s as it is used by
+Hansen & Phinney in "The Pulsar Kick Velocity Distribution".
+%}
+M_ns=1.5;%In solar masses
+sig_ns=265;vec_M=linspace(0,40);
+%We try different values for the escape velocity:
+v_esc=50;
+vec_p=[];
+for ii=1:length(vec_M)
+    M_bh=vec_M(ii);
+    p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+    vec_p=[vec_p p];
+end
+figure (1)
+plot(vec_M,vec_p,'-b')
+title('Probability of escaping as a function of BH mass (w/o fallback')
+xlabel('BH mass (in solar masses)')
+ylabel('prob(v>vesc)')
+
+hold on
+v_esc=40;
+vec_p=[];
+for ii=1:length(vec_M)
+    M_bh=vec_M(ii);
+    p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+    vec_p=[vec_p p];
+end
+figure(1)
+plot(vec_M,vec_p,'-r')
+hold on
+
+v_esc=30;
+vec_p=[];
+for ii=1:length(vec_M)
+    M_bh=vec_M(ii);
+    p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+    vec_p=[vec_p p];
+end
+figure(1)
+plot(vec_M,vec_p,'-k')
+
+legend('vesc=50','vesc=40','vesc=30')
+
+%% CALCULATION OF THE ESCAPE VELOCITY AS A FUNCTION OF DENSITY AND MASS OF THE CLUSTER
+
+%In this section, in order to calculate the escape velocity we use the
+%results from Fabio Antonini & Mark Gieles from "Population synthesis of
+%black hole binary mergers from star clusters" 2019.
+
+%We start by defining some vectors for the possible masses and radius that
+%we are going to give to the cluster.
+vec_M=linspace(10^3,10^6,10^3); %In solar masses.
+vec_r=linspace(0.5,10,1000); %In parsecs.
+%We then define a meshgrid for this two vectors:
+[M,r]=meshgrid(vec_M,vec_r);
+%We then use the model to determine the escape velocity from the cluster:
+rho_h=3*M./(8*pi*(r).^3);
+M_5=M./(10^5);
+rho_5=rho_h./(10^5);
+f_c=1;
+v_esc=50*((M_5).^(1/3)).*((rho_5).^(1/6)).*f_c;
+%Then we plot it:
+figure(2)
+contour(M_5,r,v_esc,[10,20,30,40,50,75,100],'ShowText','on')
+title('Contour plot for the v escape as a function of M_5 and r_h')
+xlabel('M_5    (M_{cl}/(10^5 M_{sun})')
+ylabel('r_h (pc)')
+
+%% CONNECTION WITH THE PROBABILITY OF ESCAPE FOR DIFFERENT BLACK HOLE MASSES
+%We now try and plot the probability os escape as a function of the BH
+%mass. This will be useful if we combine it with the results from the
+%previous section in order to set the initial conditions for the globular
+%cluster.
+%Then:
+v_esc=linspace(20,200,1000); %In km/s
+M=linspace(3,40,1000); %In solar masses;
+[M,v_esc]=meshgrid(M,v_esc);
+M_ns=1.5; %Mass of the neutron star (in solar masses)
+sig_ns=265;%Velocity dispersion for the neutron star (in km/s);
+p=prob_esc(M_ns,M,sig_ns,v_esc); %Probability of escape
+%And then we plot it:
+figure(3)
+contour(M,v_esc,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Probability of escape for different BH masses and escape velocities (w/o fallback')
+xlabel('BH mass (in solar masses)')
+ylabel('Escape velocity (in km/s)')
+
+%% PROBABILITIES OF ESCAPING FOR  A CERTAIN BH MASS IN DIFFERENT CLUSTERS
+%We start with a BH of 10 solar masses:
+v_esc=linspace(0,200,1000); %In km/s;
+M_ns=1.5;M_bh=10; %In solar masses
+sig_ns=265; %In km/s;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc); %Probability of escape
+figure(4)
+plot(v_esc,p,'-k')
+title('Probabilities of escape (w/o fallback')
+xlabel('Escape velocity (km/s)')
+ylabel('Probability')
+hold on
+%20 solar masses:
+M_bh=20;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+plot(v_esc,p,'-b')
+hold on
+%30 solar masses:
+M_bh=30;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+plot(v_esc,p,'-r')
+%40 solar masses:
+M_bh=40;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+plot(v_esc,p,'-g')
+legend('M=10','M=20','M=30','M=40')
+
+%% PROBABILITIES OF ESCAPING FOR  A CERTAIN BH MASS IN DIFFERENT CLUSTERS BUT NOW AS A FUNCTION OF M_5 AND r_h IN THE FORM OF A CONTOUR PLOT
+vec_M=linspace(10^3,10^6,10^3); %In solar masses.
+vec_r=linspace(0.5,5,1000); %In parsecs.
+%We then define a meshgrid for this two vectors:
+[M,r]=meshgrid(vec_M,vec_r);
+%We then use the model to determine the escape velocity from the cluster:
+rho_h=3*M./(8*pi*(r).^3);
+M_5=M./(10^5);
+rho_5=rho_h./(10^5);
+f_c=1;
+v_esc=50*((M_5).^(1/3)).*((rho_5).^(1/6)).*f_c;
+M_ns=1.5;sig_ns=265;
+
+%10 solar mass BH:
+M_bh=10;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(5)
+subplot(2,2,1)
+contour(M_5,r,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Prob escape (w/o fb)')
+xlabel('M_5    (M_{cl}/(10^5 M_{sun})')
+ylabel('r_h (pc)')
+legend('M=10')
+
+%20 solar mass BH:
+M_bh=20;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(5)
+subplot(2,2,2)
+contour(M_5,r,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Prob escape (w/o fb)')
+xlabel('M_5    (M_{cl}/(10^5 M_{sun})')
+ylabel('r_h (pc)')
+legend('M=20')
+
+%30 solar mass BH:
+M_bh=30;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(5)
+subplot(2,2,3)
+contour(M_5,r,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Prob escape (w/o fb)')
+xlabel('M_5    (M_{cl}/(10^5 M_{sun})')
+ylabel('r_h (pc)')
+legend('M=30')
+
+%40 solar mass BH:
+M_bh=40;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(5)
+subplot(2,2,4)
+contour(M_5,r,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Prob escape (w/o fb)')
+xlabel('M_5    (M_{cl}/(10^5 M_{sun})')
+ylabel('r_h (pc)')
+legend('M=40')
+%{
+%% CONTOUR PLOT OF THE MASS OF THE STELLAR REMNANT AS A FUNCTION OF M_ZAMS AND Z
+%{
+We get all the fiting formulas for M_rem and M_co from "The mass spectrum
+of compact remnants from the PARSEC stellar evolution tracks" from Spera,
+Mapelli & Bressan 2015.
+We use the values they give for the outputs of SEVN with the delayed SN
+model and the PARSEC stellar evolution isochrones.
+%}
+
+vec_M_zams=linspace(0.1,100,1000);
+vec_Z=linspace(1e-5,1e-2,1000);
+%We then define a meshgrid for this two vectors:
+[M_zams,Z]=meshgrid(vec_M_zams,vec_Z);
+%And then using fun_M_rem we obtain:
+M_rem=fun_M_rem(Z,M_zams);
+%And then we plot it:
+figure(6)
+contour(M_zams,Z,M_rem,[20 30 40],'ShowText','on')
+title('M_{rem} as a function of Z and M_{zams}')
+xlabel('M_{zams} (In solar masses)')
+ylabel('Z')
+ylim([1e-5,4e-3])
+%}
+%% M_rem as a function of M_zams for different metallicities
+vec_M_zams=linspace(10,100,1000);
+vec_M_rem=zeros(1,length(vec_M_zams));
+%If we set Z=1e-4:
+Z=1e-4;
+for ii=1:length(vec_M_zams)
+ vec_M_rem(ii)=fun_M_rem(Z,vec_M_zams(ii));  
+end
+figure(7)
+plot(vec_M_zams,vec_M_rem,'-k')
+title('M_{rem} as a function of M_{zams}')
+xlabel('M_{zams} (In solar masses)')
+ylabel('M_{rem} (In solar masses)')
+hold on
+%Z=1e-3:
+Z=1e-3;
+for ii=1:length(vec_M_zams)
+ vec_M_rem(ii)=fun_M_rem(Z,vec_M_zams(ii));  
+end
+plot(vec_M_zams,vec_M_rem,'-b')
+hold on
+
+%Z=2e-2:
+Z=2e-2;
+for ii=1:length(vec_M_zams)
+ vec_M_rem(ii)=fun_M_rem(Z,vec_M_zams(ii));  
+end
+plot(vec_M_zams,vec_M_rem,'-r')
+hold on
+
+legend('Z=1e-4','Z=1e-3','Z=2e-2')
+
+%% Fallback fraction as a function of M_zams for given metallicities:
+vec_M_zams=linspace(10,100,1000);
+vec_f_fb=zeros(1,length(vec_M_zams));
+vec_M_fb=zeros(1,length(vec_M_zams));
+%If we set Z=1e-4:
+Z=1e-4;
+for ii=1:length(vec_M_zams)
+ vec_f_fb(ii)=fun_fb(Z,vec_M_zams(ii));  
+end
+figure(8)
+plot(vec_M_zams,vec_f_fb,'-k')
+title('Fallback fraction as a function of M_{zams}')
+xlabel('M_{zams} (In solar masses)')
+ylabel('f_{fb}')
+hold on
+%Z=1e-3:
+Z=1e-3;
+for ii=1:length(vec_M_zams)
+ vec_f_fb(ii)=fun_fb(Z,vec_M_zams(ii));   
+end
+plot(vec_M_zams,vec_f_fb,'-b')
+hold on
+
+%Z=2e-2:
+Z=2e-2;
+for ii=1:length(vec_M_zams)
+ vec_f_fb(ii)=fun_fb(Z,vec_M_zams(ii));  
+end
+plot(vec_M_zams,vec_f_fb,'-r')
+hold on
+
+legend('Z=1e-4','Z=1e-3','Z=2e-2')
+
+%% PROBABILITY OF ESCAPE AS A FUNCTION OF BH MASS WITH FALLBACK
+%We calculate the probability of escape for a BH assuming a Maxwellian
+%distribution with velocity dispersion sigma=265km/s. Assuming same momentum for
+%both NS and BH the velocity dispersion for the BH is reduced. We also take
+%into account fallback as described in the previous two sections and
+%fun_fb.m using the function prob_esc_with_fb.m
+
+%The values that define the distribution are:
+M_ns=1.5;%In solar masses
+sig_ns=265;
+%We start by defining different values for M_zams:
+vec_M_zams=linspace(10,200,2000);
+%We will assume a standard escape velocity of 30km/s:
+v_esc=30;
+%This way we compute for different metallicities:
+
+%Z=1e-4:
+Z=1e-4;
+prob_w_fb=zeros(1,length(vec_M_zams));
+vec_M_bh=zeros(1,length(vec_M_zams));
+for ii=1:length(vec_M_bh)
+    M_zams=vec_M_zams(ii);
+    prob_w_fb(ii)=prob_esc_with_fb(Z,M_zams,M_ns,sig_ns,v_esc);
+    vec_M_bh(ii)=fun_M_rem(Z,M_zams);
+end
+figure(9)
+plot(vec_M_bh,prob_w_fb,'-k')
+title('Probability as a function of M_{bh}')
+xlabel('M_{bh}  (in solar masses)')
+ylabel('Probability of escaping')
+hold on
+
+%Z=1e-3:
+Z=1e-3;
+prob_w_fb=zeros(1,length(vec_M_zams));
+vec_M_bh=zeros(1,length(vec_M_zams));
+for ii=1:length(vec_M_bh)
+    M_zams=vec_M_zams(ii);
+    prob_w_fb(ii)=prob_esc_with_fb(Z,M_zams,M_ns,sig_ns,v_esc);
+    vec_M_bh(ii)=fun_M_rem(Z,M_zams);
+end
+plot(vec_M_bh,prob_w_fb,'-b')
+hold on
+
+%Z=2e-2:
+Z=2e-2;
+prob_w_fb=zeros(1,length(vec_M_zams));
+vec_M_bh=zeros(1,length(vec_M_zams));
+for ii=1:length(vec_M_bh)
+    M_zams=vec_M_zams(ii);
+    prob_w_fb(ii)=prob_esc_with_fb(Z,M_zams,M_ns,sig_ns,v_esc);
+    vec_M_bh(ii)=fun_M_rem(Z,M_zams);
+end
+plot(vec_M_bh,prob_w_fb,'-r')
+hold on
+
+%And now we plot the probabilities without fallback to compare the results:
+vec_p_wo_fb=[];
+for ii=1:length(vec_M_zams)
+    M_zams=vec_M_zams(ii);
+    M_bh=fun_M_rem(Z,M_zams);
+    p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+    vec_p_wo_fb=[vec_p_wo_fb p];
+end
+plot(vec_M_bh,vec_p_wo_fb,'--')
+
+legend('Z=1e-4','Z=1e-3','Z=2e-2','w/o f_{fb}')
+
+
+%% PROBABILITY OF ESCAPE AS A FUNCTION OF M_zams
+%Z=1e-4:
+Z=1e-4;
+prob_w_fb=zeros(1,length(vec_M_zams));
+vec_M_bh=zeros(1,length(vec_M_zams));
+for ii=1:length(vec_M_bh)
+    M_zams=vec_M_zams(ii);
+    prob_w_fb(ii)=prob_esc_with_fb(Z,M_zams,M_ns,sig_ns,v_esc);
+    vec_M_bh(ii)=fun_M_rem(Z,M_zams);
+end
+figure(10)
+plot(vec_M_zams,prob_w_fb,'-k')
+title('Probability as a function of M_{zams}')
+xlabel('M_{zams}  (in solar masses)')
+ylabel('Probability of escaping')
+hold on
+
+%Z=1e-3:
+Z=1e-3;
+prob_w_fb=zeros(1,length(vec_M_zams));
+vec_M_bh=zeros(1,length(vec_M_zams));
+for ii=1:length(vec_M_bh)
+    M_zams=vec_M_zams(ii);
+    prob_w_fb(ii)=prob_esc_with_fb(Z,M_zams,M_ns,sig_ns,v_esc);
+    vec_M_bh(ii)=fun_M_rem(Z,M_zams);
+end
+plot(vec_M_zams,prob_w_fb,'-b')
+hold on
+
+%Z=2e-2:
+Z=2e-2;
+prob_w_fb=zeros(1,length(vec_M_zams));
+vec_M_bh=zeros(1,length(vec_M_zams));
+for ii=1:length(vec_M_bh)
+    M_zams=vec_M_zams(ii);
+    prob_w_fb(ii)=prob_esc_with_fb(Z,M_zams,M_ns,sig_ns,v_esc);
+    vec_M_bh(ii)=fun_M_rem(Z,M_zams);
+end
+plot(vec_M_zams,prob_w_fb,'-r')
+hold on
+
+%And now we plot the probabilities without fallback to compare the results:
+vec_p_wo_fb=[];
+for ii=1:length(vec_M_zams)
+    M_zams=vec_M_zams(ii);
+    M_bh=fun_M_rem(Z,M_zams);
+    p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+    vec_p_wo_fb=[vec_p_wo_fb p];
+end
+plot(vec_M_zams,vec_p_wo_fb,'--')
+
+legend('Z=1e-4','Z=1e-3','Z=2e-2','w/o f_{fb}')
+%% CONNECTION WITH THE PROBABILITY OF ESCAPE FOR DIFFERENT BLACK HOLE MASSES (W/ FALLBACK)
+%We now try and plot the probability os escape as a function of the BH
+%mass. This will be useful if we combine it with the results from the
+%previous section in order to set the initial conditions for the globular
+%cluster.
+sig_ns=265;M_ns=1.5;Z=2e-2;
+%Then:
+v_esc=linspace(20,50,1000); %In km/s
+M=linspace(3,60,1000); %In solar masses;
+[M,v_esc]=meshgrid(M,v_esc);
+M_ns=1.5; %Mass of the neutron star (in solar masses)
+sig_ns=265;%Velocity dispersion for the neutron star (in km/s);
+p=prob_esc_with_fb(Z,M,M_ns,sig_ns,v_esc); %Probability of escape
+%And then we plot it:
+figure(11)
+contour(M,v_esc,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Probability of escape for different BH masses and escape velocities (w fallback')
+xlabel('BH mass (in solar masses)')
+ylabel('Escape velocity (in km/s)')
+
+%% CONNECTION WITH THE PROBABILITY OF ESCAPE FOR DIFFERENT BLACK HOLE MASSES
+%We now try and plot the probability os escape as a function of the BH
+%mass. This will be useful if we combine it with the results from the
+%previous section in order to set the initial conditions for the globular
+%cluster.
+%For solar metallicity:
+Z=2e-2;
+%Then:
+v_esc=linspace(20,200,1000); %In km/s
+M=linspace(3,80,1000); %In solar masses;
+[M,v_esc]=meshgrid(M,v_esc);
+M_bh=fun_M_rem(Z,M);
+M_ns=1.5; %Mass of the neutron star (in solar masses)
+sig_ns=265;%Velocity dispersion for the neutron star (in km/s);
+p=prob_esc_with_fb(Z,M,M_ns,sig_ns,v_esc); %Probability of escape
+%And then we plot it:
+figure(12)
+contour(M_bh,v_esc,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Probability of escape for different BH masses and escape velocities (w/ fallback')
+xlabel('BH mass (in solar masses)')
+ylabel('Escape velocity (in km/s)')
+
+
+
+
+
+
+
+
+
diff --git a/results/cluster properties/Cluster properties/Cluster_properties_final_version.m b/results/cluster properties/Cluster properties/Cluster_properties_final_version.m
new file mode 100644
index 0000000000000000000000000000000000000000..2f4b74d487e1befa0bb4fe6c6e95f6ed679235a1
--- /dev/null
+++ b/results/cluster properties/Cluster properties/Cluster_properties_final_version.m	
@@ -0,0 +1,319 @@
+clear all; format compact; format long g
+%% CLUSTER PROPERTIES FOR DIFFERENT METALLICITIES
+vec_M_zams=linspace(10,100,2000);
+M_ns=1.5; %in solar masses
+sig_ns=265; %in km/s
+%% Z=2E-2:
+Z=2e-2;
+%We compute M_co and M_rem as a function of M_zams
+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));
+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)
+
+subplot(2,2,[1,3])
+plot(vec_M_zams,vec_M_co,'-k')
+text(15,25,'Z=2e-2')
+title('M_{rem} and M_{co} (M_{zams})')
+xlabel('M_{zams} (in M_{sun})')
+ylabel('M in solar masses')
+hold on
+plot(vec_M_zams,vec_M_rem,'-r')
+hold on
+%We find the M_zams for which the f_fb has certain values.
+kk05=find(abs(vec_f_fb-0.5)<1e-2);
+M_zams_0_5=vec_M_zams(kk05(4));
+M_rem_0_5=fun_M_rem(Z,M_zams_0_5);
+kk075=find(abs(vec_f_fb-0.75)<1e-3);
+M_zams_0_75=vec_M_zams(kk075(1));
+M_rem_0_75=fun_M_rem(Z,M_zams_0_75);
+kk1=find(abs(vec_f_fb-1)<1e-5);
+M_zams_1=vec_M_zams(kk1(1));
+M_rem_1=fun_M_rem(Z,M_zams_1);
+plot(M_zams_0_5,M_rem_0_5,'k*','MarkerSize',10)
+hold on
+plot(M_zams_0_75,M_rem_0_75,'ko','MarkerSize',10)
+hold on
+plot(M_zams_1,M_rem_1,'k+','MarkerSize',10)
+line([M_zams_0_5 M_zams_0_5], get(gca, 'ylim'))
+line([M_zams_0_75 M_zams_0_75], get(gca, 'ylim'))
+line([M_zams_1 M_zams_1], get(gca, 'ylim'))
+legend('M_{co}','M_{rem}','f_{fb}=0.5','f_{fb}=0.75','f_{fb}=1')
+%We compute probability of escape for v_esc=30,40,50 with and without
+%fallback:
+vec_v_esc=[30 40 50];
+prob_w_fb=zeros(3,length(vec_v_esc));%We save the probabilites for each v_esc in a given row
+prob_wo_fb=zeros(3,length(vec_v_esc));
+for ii=1:length(vec_v_esc)
+    v_esc=vec_v_esc(ii);
+    for jj=1:length(vec_M_rem)
+        M_zams=vec_M_zams(jj);
+        M_bh=fun_M_rem(Z,M_zams);
+        prob_w_fb(ii,jj)=prob_esc_with_fb(Z,M_zams,M_ns,sig_ns,v_esc);
+        prob_wo_fb(ii,jj)=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+    end
+end
+%We then plot the probability of escape as a function of BH mass:
+subplot(2,2,2)
+plot(vec_M_rem,prob_w_fb(1,:),'-b')
+title('Probability of escape w/ fb')
+xlabel('M_{rem} (in M_{sun})')
+ylabel('Probability of escape')
+xlim([0 15])
+xticks([0 2.5 5 7.5 10 12.5 15 17.5 20])
+hold on
+plot(vec_M_rem,prob_w_fb(2,:),'-r')
+hold on
+plot(vec_M_rem,prob_w_fb(3,:),'-k')
+legend('v_{esc}=30','v_{esc}=40','v_{esc}=50')
+
+
+%And now we compute the differences in probability that we get with and
+%without fallback:
+dif_prob1=abs(prob_w_fb(1,:)-prob_wo_fb(1,:));
+dif_prob2=abs(prob_w_fb(2,:)-prob_wo_fb(2,:));
+dif_prob3=abs(prob_w_fb(3,:)-prob_wo_fb(3,:));
+subplot(2,2,4)
+plot(vec_M_rem,dif_prob1,'-.b')
+title('Difference in probability w/ and w/o fallback')
+xlabel('M_{rem} (in M_{sun})')
+ylabel('Difference in probability')
+xlim([0 15])%Before 3.5 there is no fallback
+xticks([0 2.5 5 7.5 10 12.5 15 17.5 20])
+hold on
+plot(vec_M_rem,dif_prob2,'-.r')
+hold on
+plot(vec_M_rem,dif_prob3,'-.k')
+legend('v_{esc}=30','v_{esc}=40','v_{esc}=50')
+%{
+%% Z=1E-3:
+Z=1e-3;
+%We compute M_co and M_rem as a function of M_zams
+vec_M_co=zeros(1,length(vec_M_zams));
+vec_M_rem=zeros(1,length(vec_M_zams));
+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);
+end
+figure(2)
+
+subplot(2,2,1)
+plot(vec_M_zams,vec_M_co)
+text(15,40,'Z=1e-3')
+title('M_{co} (M_{zams})')
+xlabel('M_{zams} (in M_{sun})')
+ylabel('M_{co} (in M_{sun})')
+subplot(2,2,2)
+plot(vec_M_zams,vec_M_rem)
+title('M_{rem} (M_{zams})')
+xlabel('M_{zams} (in M_{sun})')
+ylabel('M_{rem} (in M_{sun})')
+%We compute probability of escape for v_esc=30,40,50 with and without
+%fallback:
+vec_v_esc=[30 40 50];
+prob_w_fb=zeros(3,length(vec_v_esc));%We save the probabilites for each v_esc in a given row
+prob_wo_fb=zeros(3,length(vec_v_esc));
+for ii=1:length(vec_v_esc)
+    v_esc=vec_v_esc(ii);
+    for jj=1:length(vec_M_rem)
+        M_zams=vec_M_zams(jj);
+        M_bh=fun_M_rem(Z,M_zams);
+        prob_w_fb(ii,jj)=prob_esc_with_fb(Z,M_zams,M_ns,sig_ns,v_esc);
+        prob_wo_fb(ii,jj)=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+    end
+end
+%We then plot the probability of escape as a function of BH mass:
+subplot(2,2,3)
+plot(vec_M_rem,prob_w_fb(1,:),'-b')
+title('Probability of escape')
+xlabel('M_{rem} (in M_{sun})')
+ylabel('Probability of escape')
+xlim([0 20])
+hold on
+plot(vec_M_rem,prob_w_fb(2,:),'-r')
+hold on
+plot(vec_M_rem,prob_w_fb(3,:),'-k')
+legend('v_{esc}=30','v_{esc}=40','v_{esc}=50')
+
+%And now we compute the differences in probability that we get with and
+%without fallback:
+dif_prob1=abs(prob_w_fb(1,:)-prob_wo_fb(1,:));
+dif_prob2=abs(prob_w_fb(2,:)-prob_wo_fb(2,:));
+dif_prob3=abs(prob_w_fb(3,:)-prob_wo_fb(3,:));
+subplot(2,2,4)
+plot(vec_M_rem,dif_prob1,'-.b')
+title('Difference in probability w/ and w/o fallback')
+xlabel('M_{rem} (in M_{sun})')
+ylabel('Difference in probability')
+xlim([3.5 25])%Because befor 3.5 there is no fallback
+hold on
+plot(vec_M_rem,dif_prob2,'-.r')
+hold on
+plot(vec_M_rem,dif_prob3,'-.k')
+legend('v_{esc}=30','v_{esc}=40','v_{esc}=50')
+
+%% Z=1E-4:
+Z=1e-4;
+%We compute M_co and M_rem as a function of M_zams
+vec_M_co=zeros(1,length(vec_M_zams));
+vec_M_rem=zeros(1,length(vec_M_zams));
+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);
+end
+figure(3)
+
+subplot(2,2,1)
+plot(vec_M_zams,vec_M_co)
+text(15,40,'Z=1e-4')
+title('M_{co} (M_{zams})')
+xlabel('M_{zams} (in M_{sun})')
+ylabel('M_{co} (in M_{sun})')
+subplot(2,2,2)
+plot(vec_M_zams,vec_M_rem)
+title('M_{rem} (M_{zams})')
+xlabel('M_{zams} (in M_{sun})')
+ylabel('M_{rem} (in M_{sun})')
+%We compute probability of escape for v_esc=30,40,50 with and without
+%fallback:
+vec_v_esc=[30 40 50];
+prob_w_fb=zeros(3,length(vec_v_esc));%We save the probabilites for each v_esc in a given row
+prob_wo_fb=zeros(3,length(vec_v_esc));
+for ii=1:length(vec_v_esc)
+    v_esc=vec_v_esc(ii);
+    for jj=1:length(vec_M_rem)
+        M_zams=vec_M_zams(jj);
+        M_bh=fun_M_rem(Z,M_zams);
+        prob_w_fb(ii,jj)=prob_esc_with_fb(Z,M_zams,M_ns,sig_ns,v_esc);
+        prob_wo_fb(ii,jj)=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+    end
+end
+%We then plot the probability of escape as a function of BH mass:
+subplot(2,2,3)
+plot(vec_M_rem,prob_w_fb(1,:),'-b')
+title('Probability of escape')
+xlabel('M_{rem} (in M_{sun})')
+ylabel('Probability of escape')
+xlim([0 20])
+hold on
+plot(vec_M_rem,prob_w_fb(2,:),'-r')
+hold on
+plot(vec_M_rem,prob_w_fb(3,:),'-k')
+legend('v_{esc}=30','v_{esc}=40','v_{esc}=50')
+
+%And now we compute the differences in probability that we get with and
+%without fallback:
+dif_prob1=abs(prob_w_fb(1,:)-prob_wo_fb(1,:));
+dif_prob2=abs(prob_w_fb(2,:)-prob_wo_fb(2,:));
+dif_prob3=abs(prob_w_fb(3,:)-prob_wo_fb(3,:));
+subplot(2,2,4)
+plot(vec_M_rem,dif_prob1,'-.b')
+title('Difference in probability w/ and w/o fallback')
+xlabel('M_{rem} (in M_{sun})')
+ylabel('Difference in probability')
+xlim([3.5 25])%Because befor 3.5 there is no fallback
+hold on
+plot(vec_M_rem,dif_prob2,'-.r')
+hold on
+plot(vec_M_rem,dif_prob3,'-.k')
+legend('v_{esc}=30','v_{esc}=40','v_{esc}=50')
+
+%}
+
+%% Z=1E-3:
+Z=2e-3;
+%We compute M_co and M_rem as a function of M_zams
+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));
+for ii=1:length(vec_M_zams)
+    M_zams=vec_M_zams(ii);
+    vec_M_co(ii)=fun_M_co(Z,M_zams);
+    vec_M_rem(ii)=fun_M_rem(Z,M_zams);
+    vec_f_fb(ii)=fun_fb(Z,M_zams);
+end
+figure(4)
+
+subplot(2,2,[1,3])
+plot(vec_M_zams,vec_M_co,'-k')
+text(15,65,'Z=1e-3')
+title('M_{rem} and M_{co} (M_{zams})')
+xlabel('M_{zams} (in M_{sun})')
+ylabel('M in solar masses')
+hold on
+plot(vec_M_zams,vec_M_rem,'-r')
+hold on
+%We find the M_zams for which the f_fb has certain values.
+kk05=find(abs(vec_f_fb-0.5)<1e-2);
+M_zams_0_5=vec_M_zams(kk05(4));
+M_rem_0_5=fun_M_rem(Z,M_zams_0_5);
+kk075=find(abs(vec_f_fb-0.75)<1e-3);
+M_zams_0_75=vec_M_zams(kk075(1));
+M_rem_0_75=fun_M_rem(Z,M_zams_0_75);
+kk1=find(abs(vec_f_fb-1)<1e-5);
+M_zams_1=vec_M_zams(kk1(1));
+M_rem_1=fun_M_rem(Z,M_zams_1);
+plot(M_zams_0_5,M_rem_0_5,'k*','MarkerSize',10)
+hold on
+plot(M_zams_0_75,M_rem_0_75,'ko','MarkerSize',10)
+hold on
+plot(M_zams_1,M_rem_1,'k+','MarkerSize',10)
+line([M_zams_0_5 M_zams_0_5], get(gca, 'ylim'))
+line([M_zams_0_75 M_zams_0_75], get(gca, 'ylim'))
+line([M_zams_1 M_zams_1], get(gca, 'ylim'))
+legend('M_{co}','M_{rem}','f_{fb}=0.5','f_{fb}=0.75','f_{fb}=1')
+%We compute probability of escape for v_esc=30,40,50 with and without
+%fallback:
+vec_v_esc=[30 40 50];
+prob_w_fb=zeros(3,length(vec_v_esc));%We save the probabilites for each v_esc in a given row
+prob_wo_fb=zeros(3,length(vec_v_esc));
+for ii=1:length(vec_v_esc)
+    v_esc=vec_v_esc(ii);
+    for jj=1:length(vec_M_rem)
+        M_zams=vec_M_zams(jj);
+        M_bh=fun_M_rem(Z,M_zams);
+        prob_w_fb(ii,jj)=prob_esc_with_fb(Z,M_zams,M_ns,sig_ns,v_esc);
+        prob_wo_fb(ii,jj)=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+    end
+end
+%We then plot the probability of escape as a function of BH mass:
+subplot(2,2,2)
+plot(vec_M_rem,prob_w_fb(1,:),'-b')
+title('Probability of escape w/ fb')
+xlabel('M_{rem} (in M_{sun})')
+ylabel('Probability of escape')
+xlim([0 20])
+xticks([0 2.5 5 7.5 10 12.5 15 17.5 20])
+hold on
+plot(vec_M_rem,prob_w_fb(2,:),'-r')
+hold on
+plot(vec_M_rem,prob_w_fb(3,:),'-k')
+legend('v_{esc}=30','v_{esc}=40','v_{esc}=50')
+
+
+%And now we compute the differences in probability that we get with and
+%without fallback:
+dif_prob1=abs(prob_w_fb(1,:)-prob_wo_fb(1,:));
+dif_prob2=abs(prob_w_fb(2,:)-prob_wo_fb(2,:));
+dif_prob3=abs(prob_w_fb(3,:)-prob_wo_fb(3,:));
+subplot(2,2,4)
+plot(vec_M_rem,dif_prob1,'-.b')
+title('Difference in probability w/ and w/o fallback')
+xlabel('M_{rem} (in M_{sun})')
+ylabel('Difference in probability')
+xlim([0 20])%Before 3.5 there is no fallback
+xticks([0 2.5 5 7.5 10 12.5 15 17.5 20])
+hold on
+plot(vec_M_rem,dif_prob2,'-.r')
+hold on
+plot(vec_M_rem,dif_prob3,'-.k')
+legend('v_{esc}=30','v_{esc}=40','v_{esc}=50')
+
diff --git a/results/cluster properties/Cluster properties/bh_fractions.asv b/results/cluster properties/Cluster properties/bh_fractions.asv
new file mode 100644
index 0000000000000000000000000000000000000000..88510470dea6eb36204eaf31ae665bc9a75ce4e2
--- /dev/null
+++ b/results/cluster properties/Cluster properties/bh_fractions.asv	
@@ -0,0 +1,96 @@
+%% STUDY OF CLUSTER PROPERTIES
+%We rewrite this so it can be run from another script just selecting the
+%input variables:
+function [results_b_kicks,results_a_kicks
+
+%We define in another script the values that we will use in order to study
+%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(vec_Z));
+vec_f_m_b_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_b_kicks(ii),vec_f_m_b_kicks(ii)]=fraction_before_kicks(Z,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=[vec_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_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(Z,alfa1,alfa2,alfa3,v_esc);
+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=[vec_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);
+end
+number_fraction=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)')
+
+mass_fraction=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)')
+
+
+
diff --git a/results/cluster properties/Cluster properties/bh_fractions.m b/results/cluster properties/Cluster properties/bh_fractions.m
new file mode 100644
index 0000000000000000000000000000000000000000..fb7de4869d3c66f555e4b54a2d705a3576f6796c
--- /dev/null
+++ b/results/cluster properties/Cluster properties/bh_fractions.m	
@@ -0,0 +1,96 @@
+%% CLUSTER'S MASS AND NUMBER FRACTIONS
+%We rewrite this so it can be run from another script just selecting the
+%input variables:
+function [results_b_kicks,results_a_kicks,number_fraction,mass_fraction]=fun_cluster_fractions(alfa1,alfa2,alfa3,Z,v_esc)
+
+%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(vec_Z));
+vec_f_m_b_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_b_kicks(ii),vec_f_m_b_kicks(ii)]=fraction_before_kicks(Z,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=[vec_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_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(Z,alfa1,alfa2,alfa3,v_esc);
+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=[vec_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);
+end
+number_fraction=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)')
+
+mass_fraction=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)')
+
+
+
diff --git a/results/cluster properties/Cluster properties/bh_retention_fraction.jpg b/results/cluster properties/Cluster properties/bh_retention_fraction.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..bc065ef985f26bc6057244f15fe54e347ff5f5a9
Binary files /dev/null and b/results/cluster properties/Cluster properties/bh_retention_fraction.jpg differ
diff --git a/results/cluster properties/Cluster properties/bh_retention_fraction_Z_0001.jpg b/results/cluster properties/Cluster properties/bh_retention_fraction_Z_0001.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..99d3756bbfb83d67397427737bebbba4e94a0188
Binary files /dev/null and b/results/cluster properties/Cluster properties/bh_retention_fraction_Z_0001.jpg differ
diff --git a/results/cluster properties/Cluster properties/bh_retention_fraction_low_mass_with_fallback.jpg b/results/cluster properties/Cluster properties/bh_retention_fraction_low_mass_with_fallback.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..56b90ea8bffd6fa6aa00d0fe41fa698ec6925aac
Binary files /dev/null and b/results/cluster properties/Cluster properties/bh_retention_fraction_low_mass_with_fallback.jpg differ
diff --git a/results/cluster properties/Cluster properties/bh_retention_fraction_low_mass_wo_fallback.jpg b/results/cluster properties/Cluster properties/bh_retention_fraction_low_mass_wo_fallback.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..11f73d034a1011809c45ac41f83c951f62bd4189
Binary files /dev/null and b/results/cluster properties/Cluster properties/bh_retention_fraction_low_mass_wo_fallback.jpg differ
diff --git a/results/cluster properties/Cluster properties/bh_retention_fraction_wo_fallback.jpg b/results/cluster properties/Cluster properties/bh_retention_fraction_wo_fallback.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..f48154db001fb5b9b4ac4a57a9a5c5cafae77a34
Binary files /dev/null and b/results/cluster properties/Cluster properties/bh_retention_fraction_wo_fallback.jpg differ
diff --git a/results/cluster properties/Cluster properties/cpdf.mat b/results/cluster properties/Cluster properties/cpdf.mat
new file mode 100644
index 0000000000000000000000000000000000000000..a23b6c3dd0e23f9c6fc88b526b9c965d47ec6d29
Binary files /dev/null and b/results/cluster properties/Cluster properties/cpdf.mat differ
diff --git a/results/cluster properties/Cluster properties/expulsion.m b/results/cluster properties/Cluster properties/expulsion.m
new file mode 100644
index 0000000000000000000000000000000000000000..c88812497c366b83c3a05c4400aa96c7d2d08834
--- /dev/null
+++ b/results/cluster properties/Cluster properties/expulsion.m	
@@ -0,0 +1,59 @@
+clear all; close all; format compact; format long g;
+%% EXPULSION OF BLACK HOLES
+M_ns=1.5;sig_ns=190;vec_M=linspace(0,40);
+%We try different values for the escape velocity:
+v_esc=50;
+vec_p=[];
+for ii=1:length(vec_M)
+    M_bh=vec_M(ii);
+    p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+    vec_p=[vec_p p];
+end
+figure (1)
+plot(vec_M,vec_p,'-b')
+title('Probability of escaping as a function of BH mass')
+xlabel('BH mass (in solar masses)')
+ylabel('prob(v>vesc)')
+
+hold on
+
+v_esc=40;
+vec_p=[];
+for ii=1:length(vec_M)
+    M_bh=vec_M(ii);
+    p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+    vec_p=[vec_p p];
+end
+figure(1)
+plot(vec_M,vec_p,'-r')
+hold on
+
+v_esc=30;
+vec_p=[];
+for ii=1:length(vec_M)
+    M_bh=vec_M(ii);
+    p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+    vec_p=[vec_p p];
+end
+figure(1)
+plot(vec_M,vec_p,'-k')
+
+%If we want to try and fit it to a Plummer Model:
+M=(10^5)*2*10^30;%Total mass of the cluster;
+b=5*3.086*10^16;
+G=6.67408*10^-11;
+v_esc=(1/1000)*sqrt(2*G*M/b);
+vec_p=[];
+for ii=1:length(vec_M)
+    M_bh=vec_M(ii);
+    p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+    vec_p=[vec_p p];
+end
+figure(1)
+%plot(vec_M,vec_p,'-g')
+
+
+legend('vesc=50','vesc=40','vesc=30','Plummer model,M=10^5')
+
+
+    
diff --git a/results/cluster properties/Cluster properties/fraction_after_kicks.m b/results/cluster properties/Cluster properties/fraction_after_kicks.m
new file mode 100644
index 0000000000000000000000000000000000000000..2ea0bba4c1426704367e991ecfb37191d4a51ff0
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fraction_after_kicks.m	
@@ -0,0 +1,61 @@
+%% FUNCTION THAT GIVES THE MASS AND NUMBER FRACTION OF BH AFTER SUPERNOVA KICKS
+%Input:
+%{
+We coud set:
+      Z:The metallicity will let us obtain M_zams_min if we consider M_bh_min=4.4
+        We get this value from Farr et al. 2011. (They say between 4.3 and
+        4.5 with 90% confidence.
+But the SEVN code cites M_bh_min as 3, and we have to be consistent with
+that.
+%}
+%      alfa1 and alfa2: in order to use the IMF.
+%      k: parameter that defines whether there is fallback or not. 
+%        k=0 ==> no fallback
+%        k=1 ==> fallback considered
+%Output:
+%      f_N: number fraction of BHs that remain after supernova kicks.
+%      f_M: mass fraction of BHs that remain after supernova kicks.
+function [M_z_min,f_N,f_m]=fraction_after_kicks(Z,alfa1,alfa2,alfa3,v_esc,k)
+vec_M_zams=linspace(7,40,2000);
+vec_M_bh=fun_M_rem(Z,vec_M_zams);
+%Then we find the value for M_zams that gives M_bh=3.5:
+kk=find(abs(vec_M_bh-3)<1e-2);
+M_z_min=vec_M_zams(kk(1));
+%But we need to take into account the retention fraction, that is, we need
+%to compute the probability of escaping, so p_ret=1-p_esc.
+sig_ns=265;M_ns=1.4;
+fun= @(M) (M<0.08).*((M.^(-alfa1)).*(1-feval(@prob_escape,Z,M,M_ns,sig_ns,v_esc,k))) + (M<0.5 & M>=0.08).*((M.^(-alfa2)).*(1-feval(@prob_escape,Z,M,M_ns,sig_ns,v_esc,k))*(0.08)) + (M>=0.5).*((M.^(-alfa3)).*(1-feval(@prob_escape,Z,M,M_ns,sig_ns,v_esc,k))*(0.5));
+%This way we can compute the number fraction:
+%fun is the imf times the retention probability.
+q1=integral(fun,M_z_min,150,'RelTol',1e-6,'AbsTol',1e-6);
+%q1=integral(feval(@p_times_imf,Z,alfa1,alfa2,M,M_ns,sig_ns,v_esc),M_z_min,150,'RelTol',1e-6,'AbsTol',1e-6);
+fun2= @(M) (M<0.08).*(M.^(-alfa1)) + (M<0.5 & M>=0.08).*((M.^(-alfa2))*(0.08)) + (M>=0.5).*((M.^(-alfa3))*(0.5));
+%But we also need to take into account the retention fraction for neutron
+%stars(NS). So we need to find M_z_min_ns and M_z_max_ns corresponding to a
+%minimum mass given by M_zams=7 and a maximum mass given by 3.
+M_z_min_ns=vec_M_zams(1);
+f_ret_ns=1-feval(@prob_esc_ns,sig_ns,v_esc);
+fun_ns= @(M) (M<0.08).*((M.^(-alfa1))*f_ret_ns) + (M<0.5 & M>=0.08).*((M.^(-alfa2))*f_ret_ns*(0.08)) + (M>=0.5).*((M.^(-alfa3))*f_ret_ns*(0.5));
+q2_1=integral(fun2,0.07,M_z_min_ns,'RelTol',1e-6,'AbsTol',1e-6);
+q2_2=integral(fun_ns,M_z_min_ns,M_z_min,'RelTol',1e-6,'AbsTol',1e-6);
+q2_3=integral(fun,M_z_min,150,'RelTol',1e-6,'AbsTol',1e-6);
+f_N=q1/(q2_1+q2_2+q2_3);
+
+%Defining fun2 we obtain the BH  mass as a function of the initial total
+%mass (or total number of stars) of the cluster.
+
+%Now, in order to compute the mass fraction we have to multiply the previous function by
+%the mass. Meaning:
+alfa1=alfa1-1;alfa2=alfa2-1;alfa3=alfa3-1;
+%Here fun refers to the imf times the mass times the retention probability(as we are looking for the mass
+%fraction).
+fun= @(M) (M<0.08).*((M.^(-alfa1)).*(1-feval(@prob_escape,Z,M,M_ns,sig_ns,v_esc,k))) + (M<0.5 & M>=0.08).*((M.^(-alfa2)).*(1-feval(@prob_escape,Z,M,M_ns,sig_ns,v_esc,k))*(0.08)) + (M>=0.5).*((M.^(-alfa3)).*(1-feval(@prob_escape,Z,M,M_ns,sig_ns,v_esc,k))*(0.5));
+q1=integral(fun,M_z_min,150,'RelTol',1e-6,'AbsTol',1e-6);
+fun2= @(M) (M<0.08).*(M.^(-alfa1)) + (M<0.5 & M>=0.08).*((M.^(-alfa2))*(0.08)) + (M>=0.5).*((M.^(-alfa3))*(0.5));
+fun_ns= @(M) (M<0.08).*((M.^(-alfa1))*f_ret_ns) + (M<0.5 & M>=0.08).*((M.^(-alfa2))*f_ret_ns*(0.08)) + (M>=0.5).*((M.^(-alfa3))*f_ret_ns*(0.5));q2_1=integral(fun2,0.07,M_z_min_ns,'RelTol',1e-6,'AbsTol',1e-6);
+q2_2=integral(fun_ns,M_z_min_ns,M_z_min,'RelTol',1e-6,'AbsTol',1e-6);
+q2_3=integral(fun,M_z_min,150,'RelTol',1e-6,'AbsTol',1e-6);
+f_m=q1/(q2_1+q2_2+q2_3);
+%}
+end
+
diff --git a/results/cluster properties/Cluster properties/fraction_after_kicks_initial.m b/results/cluster properties/Cluster properties/fraction_after_kicks_initial.m
new file mode 100644
index 0000000000000000000000000000000000000000..41d1230e0706a836487f3772486face4287ee7d0
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fraction_after_kicks_initial.m	
@@ -0,0 +1,44 @@
+%% FUNCTION THAT GIVES THE MASS AND NUMBER FRACTION OF BH AFTER SUPERNOVA KICKS (AS A FUNCTION OF THE INITIAL CLUSTER)
+%Input:
+%{
+      Z:The metallicity will let us obtain M_zams_min if we consider M_bh_min=4.4
+        We get this value from Farr et al. 2011. (They say between 4.3 and
+        4.5 with 90% confidence.
+%}
+%      alfa1 and alfa2: in order to use the IMF.
+%Output:
+%      f_N: number fraction of BHs that remain after supernova kicks.
+%      f_M: mass fraction of BHs that remain after supernova kicks.
+function [M_z_min,f_N,f_m]=fraction_after_kicks_initial(Z,alfa1,alfa2,alfa3,v_esc)
+vec_M_zams=linspace(10,40,2000);
+vec_M_bh=fun_M_rem(Z,vec_M_zams);
+%Then we find the value for M_zams that gives M_bh=3.5:
+kk=find(abs(vec_M_bh-4.4)<1e-2);
+M_z_min=vec_M_zams(kk(1));
+%But we need to take into account the retention fraction, that is, we need
+%to compute the probability of escaping, so p_ret=1-p_esc
+sig_ns=265;M_ns=1.5;
+fun= @(M) (M<0.08).*((M.^(-alfa1)).*(1-feval(@prob_esc_with_fb,Z,M,M_ns,sig_ns,v_esc))) + (M<0.5 & M>=0.08).*((M.^(-alfa2)).*(1-feval(@prob_esc_with_fb,Z,M,M_ns,sig_ns,v_esc))) + (M>=0.5).*((M.^(-alfa3)).*(1-feval(@prob_esc_with_fb,Z,M,M_ns,sig_ns,v_esc)));
+%This way we can compute the number fraction:
+%fun is the imf times the retention probability.
+q1=integral(fun,M_z_min,150,'RelTol',1e-6,'AbsTol',1e-6);
+%q1=integral(feval(@p_times_imf,Z,alfa1,alfa2,M,M_ns,sig_ns,v_esc),M_z_min,150,'RelTol',1e-6,'AbsTol',1e-6);
+fun2= @(M) (M<0.08).*(M.^(-alfa1)) + (M<0.5 & M>=0.08).*(M.^(-alfa2)) + (M>=0.5).*(M.^(-alfa3));
+q2=integral(fun2,0.07,150,'RelTol',1e-6,'AbsTol',1e-6);
+f_N=q1/q2;
+
+%Defining fun2 we obtain the BH  mass as a function of the initial total
+%mass (or total number of stars) of the cluster.
+
+%Now, in order to compute the mass fraction we have to multiply the previous function by
+%the mass. Meaning:
+alfa1=alfa1-1;alfa2=alfa2-1;alfa3=alfa3-1;
+%Here fun refers to the imf times the mass times the retention probability(as we are looking for the mass
+%fraction).
+fun= @(M) (M<0.08).*((M.^(-alfa1)).*(1-feval(@prob_esc_with_fb,Z,M,M_ns,sig_ns,v_esc))) + (M<0.5 &M>=0.08).*((M.^(-alfa2)).*(1-feval(@prob_esc_with_fb,Z,M,M_ns,sig_ns,v_esc))) + (M>=0.5).*((M.^(-alfa3)).*(1-feval(@prob_esc_with_fb,Z,M,M_ns,sig_ns,v_esc)));
+q1=integral(fun,M_z_min,150,'RelTol',1e-6,'AbsTol',1e-6);
+fun2= @(M) (M<0.08).*(M.^(-alfa1)) + (M<0.5 & M>=0.08).*(M.^(-alfa2)) + (M>=0.5).*(M.^(-alfa3));
+q2=integral(fun2,0.07,150,'RelTol',1e-6,'AbsTol',1e-6);
+f_m=q1/q2;
+%}
+end
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/fraction_after_kicks_ns.m b/results/cluster properties/Cluster properties/fraction_after_kicks_ns.m
new file mode 100644
index 0000000000000000000000000000000000000000..ecc57b05c9faeef80b26ca4e6e1ec414afef8c02
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fraction_after_kicks_ns.m	
@@ -0,0 +1,64 @@
+%% FUNCTION THAT GIVES THE MASS AND NUMBER FRACTION OF BH AFTER SUPERNOVA KICKS
+%Input:
+%{
+We coud set:
+      Z:The metallicity will let us obtain M_zams_min if we consider M_bh_min=4.4
+        We get this value from Farr et al. 2011. (They say between 4.3 and
+        4.5 with 90% confidence.
+But the SEVN code cites M_bh_min as 3, and we have to be consistent with
+that.
+%}
+%      alfa1 and alfa2: in order to use the IMF.
+%      k: parameter that defines whether there is fallback or not. 
+%        k=0 ==> no fallback
+%        k=1 ==> fallback considered
+%Output:
+%      f_N: number fraction of BHs that remain after supernova kicks.
+%      f_M: mass fraction of BHs that remain after supernova kicks.
+function [f_N,f_m]=fraction_after_kicks_ns(Z,alfa1,alfa2,alfa3,v_esc,k)
+vec_M_zams=linspace(7,40,2000);
+vec_M_bh=fun_M_rem(Z,vec_M_zams);
+%Then we find the value for M_zams that gives M_bh=3.5:
+kk=find(abs(vec_M_bh-3)<1e-2);
+M_z_min=vec_M_zams(kk(1));
+%And the minimum NS mass is:
+M_z_min_ns=vec_M_zams(1);
+%But we need to take into account the retention fraction, that is, we need
+%to compute the probability of escaping, so p_ret=1-p_esc.
+sig_ns=265;M_ns=1.4;
+fun= @(M) (M<0.08).*((M.^(-alfa1)).*(1-feval(@prob_escape,Z,M,M_ns,sig_ns,v_esc,k))) + (M<0.5 & M>=0.08).*((M.^(-alfa2)).*(1-feval(@prob_escape,Z,M,M_ns,sig_ns,v_esc,k))*(0.08)) + (M>=0.5).*((M.^(-alfa3)).*(1-feval(@prob_escape,Z,M,M_ns,sig_ns,v_esc,k))*(0.5));
+%This way we can compute the number fraction:
+%fun is the imf times the retention probability.
+
+fun2= @(M) (M<0.08).*(M.^(-alfa1)) + (M<0.5 &M>=0.08).*(M.^(-alfa2)) + (M>=0.5).*(M.^(-alfa3));
+%But we also need to take into account the retention fraction for neutron
+%stars(NS). So we need to find M_z_min_ns and M_z_max_ns corresponding to a
+%minimum mass given by M_zams=7 and a maximum mass given by 3.
+M_z_min_ns=vec_M_zams(1);
+f_ret_ns=1-feval(@prob_esc_ns,sig_ns,v_esc);
+fun_ns= @(M) (M<0.08).*((M.^(-alfa1))*f_ret_ns) + (M<0.5 & M>=0.08).*((M.^(-alfa2))*f_ret_ns*(0.08)) + (M>=0.5).*((M.^(-alfa3))*f_ret_ns*(0.5));q2_1=integral(fun2,0.07,M_z_min_ns,'RelTol',1e-6,'AbsTol',1e-6);
+q1=integral(fun_ns,M_z_min_ns,M_z_min,'RelTol',1e-6,'AbsTol',1e-6);
+q2_1=integral(fun2,0.07,M_z_min_ns,'RelTol',1e-6,'AbsTol',1e-6);
+q2_2=integral(fun_ns,M_z_min_ns,M_z_min,'RelTol',1e-6,'AbsTol',1e-6);
+q2_3=integral(fun,M_z_min,150,'RelTol',1e-6,'AbsTol',1e-6);
+f_N=q1/(q2_1+q2_2+q2_3);
+
+%Defining fun2 we obtain the BH  mass as a function of the initial total
+%mass (or total number of stars) of the cluster.
+
+%Now, in order to compute the mass fraction we have to multiply the previous function by
+%the mass. Meaning:
+alfa1=alfa1-1;alfa2=alfa2-1;alfa3=alfa3-1;
+%Here fun refers to the imf times the mass times the retention probability(as we are looking for the mass
+%fraction).
+fun= @(M) (M<0.08).*((M.^(-alfa1)).*(1-feval(@prob_escape,Z,M,M_ns,sig_ns,v_esc,k))) + (M<0.5 & M>=0.08).*((M.^(-alfa2)).*(1-feval(@prob_escape,Z,M,M_ns,sig_ns,v_esc,k))*(0.08)) + (M>=0.5).*((M.^(-alfa3)).*(1-feval(@prob_escape,Z,M,M_ns,sig_ns,v_esc,k))*(0.5));
+fun2= @(M) (M<0.08).*(M.^(-alfa1)) + (M<0.5 & M>=0.08).*((M.^(-alfa2))*(0.08)) + (M>=0.5).*((M.^(-alfa3))*(0.5));
+fun_ns= @(M) (M<0.08).*((M.^(-alfa1))*f_ret_ns) + (M<0.5 & M>=0.08).*((M.^(-alfa2))*f_ret_ns*(0.08)) + (M>=0.5).*((M.^(-alfa3))*f_ret_ns*(0.5));q2_1=integral(fun2,0.07,M_z_min_ns,'RelTol',1e-6,'AbsTol',1e-6);
+q1=integral(fun_ns,M_z_min_ns,M_z_min,'RelTol',1e-6,'AbsTol',1e-6);
+q2_1=integral(fun2,0.07,M_z_min_ns,'RelTol',1e-6,'AbsTol',1e-6);
+q2_2=integral(fun_ns,M_z_min_ns,M_z_min,'RelTol',1e-6,'AbsTol',1e-6);
+q2_3=integral(fun,M_z_min,150,'RelTol',1e-6,'AbsTol',1e-6);
+f_m=q1/(q2_1+q2_2+q2_3);
+%}
+end
+
diff --git a/results/cluster properties/Cluster properties/fraction_before_kicks.m b/results/cluster properties/Cluster properties/fraction_before_kicks.m
new file mode 100644
index 0000000000000000000000000000000000000000..f9466a022bce0433f0cd458968abcce5ff6cb78c
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fraction_before_kicks.m	
@@ -0,0 +1,36 @@
+%% FUNCTION THAT GIVES THE MASS AND NUMBER FRACTION OF BH BEFORE SUPERNOVA KICKS (AFTER FORMATION)
+%Input:
+%{
+      Z:The metallicity will let us obtain M_zams_min if we consider M_bh_min=4.4
+        We get this value from Farr et al.2011. (They say between 4.3 and
+        4.5 with 90% confidence.
+But the SEVN code cites M_bh_min as 3, and we have to be consistent with
+that.
+%}
+%      alfa1,alfa2,alfa3: in order to use the IMF.
+%Output:
+%      f_N: number fraction of stars that will become black holes.
+%      f_M: mass fraction of stars that will become black holes.
+function [M_z_min,f_N,f_m]=fraction_before_kicks(Z,alfa1,alfa2,alfa3)
+vec_M_zams=linspace(10,40,2000);
+vec_M_bh=fun_M_rem(Z,vec_M_zams);
+%Then we find the value for M_zams that gives M_bh=3:
+kk=find(abs(vec_M_bh-3)<1e-2);
+M_z_min=vec_M_zams(kk(1));
+%This way we can compute the number fraction:
+%fun is the imf basically.
+fun= @(M) (M<0.08).*(M.^(-alfa1)) + (M<0.5 & M>=0.08).*((M.^(-alfa2))*(0.08)) + (M>=0.5).*((M.^(-alfa3))*(0.5));
+q1=integral(fun,M_z_min,150,'RelTol',1e-6,'AbsTol',1e-6);
+q2=integral(fun,0.07,150,'RelTol',1e-6,'AbsTol',1e-6);
+f_N=q1/q2;
+%Now, in order to compute the mass fraction we have to multiply the IMF by
+%the mass. Meaning:
+alfa1=alfa1-1;alfa2=alfa2-1;alfa3=alfa3-1;
+%Here fun refers to the imf times the mass (as we are looking for the mass
+%fraction).
+fun= @(M) (M<0.08).*(M.^(-alfa1)) + (M<0.5 & M>=0.08).*((M.^(-alfa2))*(0.08)) + (M>=0.5).*((M.^(-alfa3))*(0.5));
+q1=integral(fun,M_z_min,150,'RelTol',1e-6,'AbsTol',1e-6);
+q2=integral(fun,0.07,150,'RelTol',1e-6,'AbsTol',1e-6);
+f_m=q1/q2;
+end
+
diff --git a/results/cluster properties/Cluster properties/fraction_before_kicks_ns.m b/results/cluster properties/Cluster properties/fraction_before_kicks_ns.m
new file mode 100644
index 0000000000000000000000000000000000000000..321b26194be364eec76c32f5384744379317aa5d
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fraction_before_kicks_ns.m	
@@ -0,0 +1,42 @@
+%% FUNCTION THAT GIVES THE MASS AND NUMBER FRACTION OF NS BEFORE SUPERNOVA KICKS (AFTER FORMATION)
+%Considerations:
+%{
+The function we have for M_rem and M_co are only valid for M_zams>7, so wi
+we will set that as the lower limit for NS, M=3 being the upper limit.
+%}
+%Input:
+%{
+      Z:The metallicity will let us obtain M_zams_min if we consider M_bh_min=4.4
+        We get this value from Farr et al. 2011. (They say between 4.3 and
+        4.5 with 90% confidence.
+But the SEVN code cites M_bh_min as 3, and we have to be consistent with
+that.
+%}
+%      alfa1,alfa2,alfa3: in order to use the IMF.
+%Output:
+%      f_N: number fraction of stars that will become black holes.
+%      f_M: mass fraction of stars that will become black holes.
+function [f_N,f_m]=fraction_before_kicks_ns(Z,alfa1,alfa2,alfa3)
+vec_M_zams=linspace(7,40,2000);
+vec_M_bh=fun_M_rem(Z,vec_M_zams);
+%The minimum NS mass is:
+M_z_min_ns=vec_M_zams(1);
+%Then we find the value for M_zams that gives M_bh=3:
+kk=find(abs(vec_M_bh-3)<1e-2);
+M_z_min=vec_M_zams(kk(1));
+%This way we can compute the number fraction:
+%fun is the imf basically.
+fun= @(M) (M<0.08).*(M.^(-alfa1)) + (M<0.5 & M>=0.08).*((M.^(-alfa2))*(0.08)) + (M>=0.5).*((M.^(-alfa3))*(0.5));
+q1=integral(fun,M_z_min_ns,M_z_min,'RelTol',1e-6,'AbsTol',1e-6);
+q2=integral(fun,0.07,150,'RelTol',1e-6,'AbsTol',1e-6);
+f_N=q1/q2;
+%Now, in order to compute the mass fraction we have to multiply the IMF by
+%the mass. Meaning:
+alfa1=alfa1-1;alfa2=alfa2-1;alfa3=alfa3-1;
+%Here fun refers to the imf times the mass (as we are looking for the mass
+%fraction).
+fun= @(M) (M<0.08).*(M.^(-alfa1)) + (M<0.5 & M>=0.08).*((M.^(-alfa2))*(0.08)) + (M>=0.5).*((M.^(-alfa3))*(0.5));
+q1=integral(fun,M_z_min_ns,M_z_min,'RelTol',1e-6,'AbsTol',1e-6);
+q2=integral(fun,0.07,150,'RelTol',1e-6,'AbsTol',1e-6);
+f_m=q1/q2;
+end
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/fraction_pie_charta.m b/results/cluster properties/Cluster properties/fraction_pie_charta.m
new file mode 100644
index 0000000000000000000000000000000000000000..fe12dd155f4d0d3e5780643a5a9e272b8aa15cb7
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fraction_pie_charta.m	
@@ -0,0 +1,46 @@
+clear all; close all; format compact; format long g;
+%% PIE CHARTS OF MASS AND NUMBER FRACTIONS FOR CONTINUOUS CASE (NO SAMPLING)
+%If we use a Kroupa IMF:
+alfa1=0.3;alfa2=1.3;alfa3=2.3;
+%We define the metallicity:
+Z=2e-2;
+%We want fallback:
+k=1;
+%We define the escape velocity: (in km/s)
+v_esc=30;
+%% UPON FORMATION (JUST AS A FUNCTION OF THE IMF
+%We then get:
+[M_z_min,f_N_bh_bef_kicks,f_m_bh_bef_kicks]=fraction_before_kicks(Z,alfa1,alfa2,alfa3);
+[f_N_ns_bef_kicks,f_m_ns_bef_kicks]=fraction_before_kicks_ns(Z,alfa1,alfa2,alfa3);
+f_N_stars_bef_kicks=1-f_N_bh_bef_kicks-f_N_ns_bef_kicks;
+f_m_stars_bef_kicks=1-f_m_bh_bef_kicks-f_m_ns_bef_kicks;
+%We plot the number fraction:
+figure(1)
+labels={'BHs','NSs','stars'};
+results_N_bef=[f_N_bh_bef_kicks f_N_ns_bef_kicks f_N_stars_bef_kicks];
+pie(results_N_bef,labels)
+title('Number fraction (just the IMF)')
+%We plot the mass fraction:
+figure(2)
+labels={'BHs','NSs','stars'};
+results_m_bef=[f_m_bh_bef_kicks f_m_ns_bef_kicks f_m_stars_bef_kicks];
+pie(results_m_bef,labels)
+title('Mass fraction (just the IMF)')
+%% AFTER SUPERNOVA KICKS
+[M_z_min,f_N_bh_after_kicks,f_m_bh_after_kicks]=fraction_after_kicks(Z,alfa1,alfa2,alfa3,v_esc,k);
+[f_N_ns_after_kicks,f_m_ns_after_kicks]=fraction_after_kicks_ns(Z,alfa1,alfa2,alfa3,v_esc,k);
+f_N_stars_after_kicks=1-f_N_bh_after_kicks-f_N_ns_after_kicks;
+f_m_stars_after_kicks=1-f_m_bh_after_kicks-f_m_ns_after_kicks;
+%We plot the number fraction:
+figure(3)
+labels={'BHs','NSs','stars'};
+results_N_aft=[f_N_bh_after_kicks f_N_ns_after_kicks f_N_stars_after_kicks];
+pie(results_N_aft,labels)
+title('Number fraction (after supernova kicks)')
+%We plot the mass fraction:
+figure(4)
+labels={'BHs','NSs','stars'};
+results_m_aft=[f_m_bh_after_kicks f_m_ns_after_kicks f_m_stars_after_kicks];
+pie(results_m_aft,labels)
+title('Mass fraction (after supernova kicks)')
+
diff --git a/results/cluster properties/Cluster properties/fun_M_co.m b/results/cluster properties/Cluster properties/fun_M_co.m
new file mode 100644
index 0000000000000000000000000000000000000000..39da5df61da32b6252227d021a59dbfa85b192c9
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fun_M_co.m	
@@ -0,0 +1,25 @@
+%% M_co as a function of M_zams and Z
+%Input:
+%      Z:Metallicity
+%      M_zams:Zero age main sequence mass (in solar masses)
+
+%Output:
+%      M_co:Mass of the CO core (in solar masses)
+function M_co=fun_M_co(Z,M_zams)
+if Z>(4e-3)
+    B1=(59.63)-(2.969e3)*Z+(4.988e4)*Z.^2;
+    K1=(45.04)-(2.176e3)*Z+(3.806e4)*Z.^2;
+    K2=(1.389e2)-(4.664e3)*Z+(5.106e4)*Z.^2;
+    delta1=(2.790e-2)-(1.78e-2)*Z+(77.05)*Z.^2;
+    delta2=(6.730e-3)+(2.690)*Z-(52.39)*Z.^2;
+elseif Z<(1e-3)
+    B1=67.07; K1=46.89; K2=1.138e2; delta1=2.199e-2; delta2=2.602e-2;
+else
+    B1=(40.98)+(3.415e4)*Z-(8.064e6)*Z.^2;
+    K1=(35.17)+(1.548e4)*Z-(3.759e6)*Z.^2;
+    K2=(20.36)+(1.162e5)*Z-(2.276e7)*Z.^2;
+    delta1=(2.500e-2)-(4.346)*Z+(1.340e3)*Z.^2;
+    delta2=(1.750e-2)+(11.39)*Z-(2.902e3)*Z.^2;
+end
+M_co=-2.0+((B1+2).*(((0.5)./(1+10.^((K1-M_zams).*delta1)))+((0.5)./(1+10.^((K2-M_zams).*delta2)))));
+end
diff --git a/results/cluster properties/Cluster properties/fun_M_rem.m b/results/cluster properties/Cluster properties/fun_M_rem.m
new file mode 100644
index 0000000000000000000000000000000000000000..0df6083132a694135e157fcba89d7eac82dfb547
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fun_M_rem.m	
@@ -0,0 +1,62 @@
+%% M_rem as a function of Z and M_zams
+%Input:
+%      Z:Metallicity
+%      M_zams:Zero age main sequence mass (in solar masses)
+
+%Output:
+%      M_rem:Mass of the stellar remnant (in solar masses)
+function M_rem=fun_M_rem(Z,M_zams)
+%We start by computing the M_co via fun_M_co:
+mat_M_co=fun_M_co(Z,M_zams);
+for jj=1:length(mat_M_co(:,1))
+for ii=1:length(mat_M_co(jj,:))
+    M_co=mat_M_co(jj,ii);
+%We then compute m and q as a function of the metallicity Z
+if Z>=(2e-3)
+    m=1.217;q=1.061;
+elseif Z<(1e-3)
+    m=-(6.476e2)*Z+(1.911);
+    q=(2.300e3)*Z+(11.67);
+else
+    m=-(43.82)*Z+(1.304);
+    q=-(1.296e4)*Z+26.98;
+end
+%We then compute A1,A2,L,eta:
+if Z>=(1e-3)
+    A1=(1.340)-((29.46)./(1+(Z./(1.110e-3)).^(2.361)));
+    A2=(80.22)-(74.73)*((Z.^(0.965))./((2.720e-3)+(Z.^(0.965))));
+    L=(5.683)+((3.533)./(1+(Z./(7.430e-3)).^(1.993)));
+    eta=(1.066)-((1.121)./(1+(Z./(2.558e-2)).^(0.609)));
+else
+    A1=(1.105e5)*Z-(1.258e2);
+    A2=(91.56)-(1.957e4)*Z-(1.558e7)*Z.^2;
+    L=(1.134e4)*Z-(2.143);
+    eta=(3.090e-2)-(22.30)*Z+(7.363e4)*Z.^2;
+end
+%And for the M_rem:
+if Z<=(5e-4)
+ %Then we compute p and f:
+ p=-2.333+(0.1559)*M_co+(0.2700)*(M_co).^2;
+ f=m.*M_co+q;
+    if M_co<=5
+        M_rem(jj,ii)=max(p,1.27);
+    elseif M_co>=10
+        M_rem(jj,ii)=min(p,f);
+    else
+        M_rem(jj,ii)=p;
+    end
+else
+ %We calculate h and f so that:
+ h=A1+((A2-A1)./(1+10.^((L-M_co).*eta)));
+ f=m.*M_co+q;
+    if M_co<=5
+        M_rem(jj,ii)=max(h,1.27);
+    elseif M_co>=10
+        M_rem(jj,ii)=max(h,f);
+    else 
+        M_rem(jj,ii)=h;
+    end
+end
+end
+end
+end
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/fun_cluster_fractions.m b/results/cluster properties/Cluster properties/fun_cluster_fractions.m
new file mode 100644
index 0000000000000000000000000000000000000000..64528988ead4d3e0b1d9aadf5c73ee0fdec9bafb
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fun_cluster_fractions.m	
@@ -0,0 +1,96 @@
+%% CLUSTER'S MASS AND NUMBER FRACTIONS
+%We rewrite this so it can be run from another script just selecting the
+%input variables:
+function [results_b_kicks,results_a_kicks]=fun_cluster_fractions(alfa1,alfa2,alfa3,Z,v_esc,k)
+
+%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/fun_fb.m b/results/cluster properties/Cluster properties/fun_fb.m
new file mode 100644
index 0000000000000000000000000000000000000000..22ae0f5d3c466039671810b8b279363b33a4af67
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fun_fb.m	
@@ -0,0 +1,45 @@
+%% f_fb and M_fb as a function of Z and M_zams
+%Input:
+%      Z:Metallicity
+%      M_zams:Zero age main sequence mass (in solar masses)
+
+%Output:
+%      f_fb: Fallback fraction
+%      M_fb: Fallback mass (in solar masses)
+function f_fb=fun_fb(Z,M_zams)
+%We start by computing the M_co via fun_M_co
+vec_M_co=fun_M_co(Z,M_zams);
+%We then compute the value of M_proto:
+M_proto=zeros(1,length(vec_M_co));
+for ii=1:length(vec_M_co)
+    M_co=vec_M_co(ii);
+if M_co<3.5 
+    M_proto(ii)=1.2;
+elseif (3.5<=M_co)&(M_co<6.0)
+    M_proto(ii)=1.3;
+elseif (6.0<=M_co)&(M_co<11.0)
+    M_proto(ii)=1.4;
+else
+    M_proto(ii)=1.6;
+end
+%Now that we have M_proto we can compute a2 and b2:
+a2=(0.133)-((0.093)./(M_zams(ii)-M_proto(ii)));
+b2=-(11)*a2+1;
+%And then, in order to obtain M_fb and f_fb:
+if M_co<2.5
+    M_fb(ii)=0.2;
+    f_fb(ii)=0;%I don't know if this is what I have to do
+elseif (2.5<=M_co)&(M_co<3.5)
+    M_fb(ii)=0.5*M_co-1.05;
+    f_fb(ii)=0;%I don't know if this is what I have to do
+elseif (3.5<=M_co)&(M_co<11.0)
+    f_fb(ii)=a2.*M_co+b2;
+    M_fb(ii)=f_fb(ii).*(M_zams(ii)-M_proto(ii));
+else
+    f_fb(ii)=1.0;
+    M_fb(ii)=f_fb(ii).*(M_zams(ii)-M_proto(ii));
+end
+end
+end
+    
+    
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/fun_for_contour_fractions.m b/results/cluster properties/Cluster properties/fun_for_contour_fractions.m
new file mode 100644
index 0000000000000000000000000000000000000000..c2c65965c89501dc2a0fb5196dc659a8f87fbab8
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fun_for_contour_fractions.m	
@@ -0,0 +1,94 @@
+%% ONLY USABLE AFTER THE USE OF A MESHGRID
+%% FUNCTION THAT WILL GIVE US WHAT WE NEED IN ORDER TO FORM THE CONTOUR PLOT:
+%In this case we only want to consider the number fraction of black holes,
+%but not compared to the current cluster, but compared to the initial
+%number of black holes in such cluster. That is, compare the number of
+%black holes after kicks (w/ or w/o fallbkack) to the numbers obtained via
+%the sampling of the IMF.
+%
+
+%For more information on the processes followed see sampling.m
+function f_N_bh_relative=fun_for_contour_fractions(M_cl,Z,vec_v_esc)
+%We start by sampling the IMF. And for that we need to load the cpdf
+%(cumulative probability density function):
+load('cpdf')
+load('precomp_probs')
+%And then we can sample the IMF:
+jj=1;%We initialize this variable before starting the while loop.
+vec_m_initial=[];%We will use this vector to save the different masses.
+for dd=1:length(vec_v_esc(1,:))
+for ff=1:length(vec_v_esc(:,1))
+        %M_cl_it=M_cl(dd,ff);
+        M_tot_initial=0;
+tic
+while M_tot_initial<M_cl(dd,ff)
+    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
+N_tot_initial(dd,ff)=length(vec_m_initial);
+%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 stars of each of them:
+N_bh_initial(dd,ff)=sum(vec_m_initial>M_z_min_bh_test);
+
+%The important thing to consider in this case is that we are not interested
+%in storing the different masses, only the fraction of BHs. This helps a
+%lot because it means that do not have to worry about storing all those
+%variables.
+f_N_bh=[];
+vec_m_final=[];
+    ll=1;%We initialize this and then change the value by one if the star is retained.
+    tic
+    %We calculate only once the probability for the initial vector of masses:
+    f_ret=(1-feval(@prob_escape,Z,vec_m_initial,1.4,265,vec_v_esc(dd,ff),1));
+    for jj=1:length(vec_m_initial)        
+        if m<M_z_min_ns_test
+            vec_m_final(ll)=vec_m_initial(jj);
+            ll=ll+1;%Because it is always retained.
+        elseif m>M_z_min_bh_test
+            %k=1 ==>with fallback
+            %f_ret=(1-feval(@prob_escape,Z,m,1.4,265,vec_v_esc(dd,ff),1));
+            %We now generate a random number, if it is below f_ret it is
+            %retained, if not it is kicked out of the cluster.
+            rr=rand;
+            if rr<=f_ret(jj)
+                vec_m_final(ll)=vec_m_initial(jj);
+                ll=ll+1;
+            else
+                ll=ll;
+            end
+        else
+            %Again we compute the probability of it being retained:
+            %f_ret=(1-feval(@prob_esc_ns,265,v_esc(dd,ff)));
+            rr=rand;
+            if rr<=(1-feval(@prob_esc_ns,265,vec_v_esc(dd,ff)))
+                vec_m_final(ll)=vec_m_initial(jj);
+                ll=ll+1;
+            else
+                ll=ll;
+            end
+        end
+    N_bh_final(dd,ff)=sum(vec_m_final>M_z_min_bh_test);
+    f_N_bh_relative(dd,ff)=N_bh_final(dd,ff)./N_bh_initial(dd,ff);
+   
+    end
+    toc
+end
+end
+end
diff --git a/results/cluster properties/Cluster properties/fun_for_contour_fractions_2.m b/results/cluster properties/Cluster properties/fun_for_contour_fractions_2.m
new file mode 100644
index 0000000000000000000000000000000000000000..d8ad52e87fd3a74d308fe03c7d0428cd011e6f36
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fun_for_contour_fractions_2.m	
@@ -0,0 +1,122 @@
+%% ONLY USABLE AFTER THE USE OF A MESHGRID
+%% FUNCTION THAT WILL GIVE US WHAT WE NEED IN ORDER TO FORM THE CONTOUR PLOT:
+%In this case we only want to consider the number fraction of black holes,
+%but not compared to the current cluster, but compared to the initial
+%number of black holes in such cluster. That is, compare the number of
+%black holes after kicks (w/ or w/o fallbkack) to the numbers obtained via
+%the sampling of the IMF.
+%
+
+%For more information on the processes followed see sampling.m
+function f_N_bh_relative=fun_for_contour_fractions_2(M_cl,Z,vec_v_esc)
+%We start by sampling the IMF. And for that we need to load the cpdf
+%(cumulative probability density function):
+load('cpdf')
+load('precomp_probs')
+%And then we can sample the IMF:
+vec_m_initial=[];%We will use this vector to save the different masses.
+for dd=1:length(vec_v_esc(1,:))
+for ff=1:length(vec_v_esc(:,1))
+        %M_cl_it=M_cl(dd,ff);
+        M_tot_initial=0;
+tic
+jj=1;%We initialize this variable before starting the while loop.
+while M_tot_initial<M_cl(dd,ff)
+    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
+N_tot_initial(dd,ff)=length(vec_m_initial);
+%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 stars of each of them:
+N_bh_initial(dd,ff)=sum(vec_m_initial>M_z_min_bh_test);
+
+%The important thing to consider in this case is that we are not interested
+%in storing the different masses, only the fraction of BHs. This helps a
+%lot because it means that do not have to worry about storing all those
+%variables.
+f_N_bh=[];
+vec_m_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:
+pos_stars=vec_m_initial<M_z_min_ns_test;
+vec_m_final_i=vec_m_initial.*pos_stars;
+%But we need to eliminate the 0's and only consider those that are
+%different  from 0.
+vec_m_final=nonzeros(vec_m_final_i)';
+%And the vector of masses for the NS and BHs before considering the
+%probabilities for them to be kicked out are:
+vec_bh_ns=vec_m_initial-vec_m_final_i;
+vec_bh_ns=nonzeros(vec_bh_ns)';
+%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,vec_v_esc(dd,ff),1)); 
+     for jj=1:length(vec_bh_ns)
+         if vec_bh_ns(jj)>M_z_min_bh_test
+            %k=1 ==>with fallback
+            %f_ret=(1-feval(@prob_escape,Z,m,1.4,265,vec_v_esc(dd,ff),1));
+            %We now generate a random number, if it is below f_ret it is
+            %retained, if not it is kicked out of the cluster.
+            rr=rand;
+            if rr<=f_ret(jj)
+                vec_m_final(ll)=vec_bh_ns(jj);
+                ll=ll+1;
+             else
+                ll=ll;
+             end
+           else
+            %Again we compute the probability of it being retained:
+            %f_ret=(1-feval(@prob_esc_ns,265,v_esc(dd,ff)));
+            rr=rand;
+            if rr<=(1-feval(@prob_esc_ns,265,vec_v_esc(dd,ff)))
+                vec_m_final(ll)=vec_bh_ns(jj);
+                ll=ll+1;
+            else
+                ll=ll;
+            end
+         end
+     end
+    N_bh_final(dd,ff)=sum(vec_m_final>M_z_min_bh_test);
+    f_N_bh_relative(dd,ff)=N_bh_final(dd,ff)./N_bh_initial(dd,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_3.m b/results/cluster properties/Cluster properties/fun_for_contour_fractions_3.m
new file mode 100644
index 0000000000000000000000000000000000000000..5c0a02a45b832a1562f3cac2f5a06201785b7114
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fun_for_contour_fractions_3.m	
@@ -0,0 +1,135 @@
+%% FUNCTION THAT WILL GIVE US WHAT WE NEED IN ORDER TO FORM THE CONTOUR PLOT:
+%{
+Inputs:
+       M_cl: Mass of the cluster (in solar masses)
+       Z: Metallicity
+       mat_v_esc: escape velocity (in km/s)
+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_3(M_cl,Z,mat_v_esc)
+%We start by sampling the IMF. And for that we need to load the cpdf
+%(cumulative probability density 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:
+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.
+%And now we do an if statement for this stars:
+ll=length(vec_m_final); %We initialize ll so that we fill the vector after that position.
+%We calculate only once the probability for the initial vector of masses:
+    f_ret=(1-feval(@prob_escape,Z,vec_bh_ns,1.4,265,mat_v_esc(dd,ff),1)); 
+     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.
+            rr=rand;
+            if rr<=f_ret(jj)
+                vec_m_final(ll)=vec_bh_ns(jj);
+                ll=ll+1;
+             else
+                ll=ll;
+             end
+           else
+            %Again we compute the probability of it being retained:
+            %f_ret=(1-feval(@prob_esc_ns,265,v_esc(dd,ff)));
+            rr=rand;
+            if rr<=(1-feval(@prob_esc_ns,265,mat_v_esc(dd,ff)))
+                vec_m_final(ll)=vec_bh_ns(jj);
+                ll=ll+1;
+            else
+                ll=ll;
+            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);
+    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_mass_histograms.m b/results/cluster properties/Cluster properties/fun_for_mass_histograms.m
new file mode 100644
index 0000000000000000000000000000000000000000..c04c5b7c3b1b8c064ceb61767ffbf88e97d67393
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fun_for_mass_histograms.m	
@@ -0,0 +1,125 @@
+%% FUNCTION THAT WILL GIVE US WHAT WE NEED IN ORDER TO FORM THE CONTOUR PLOT:
+%{
+Inputs:
+       M_cl: Mass of the cluster (in solar masses)
+       Z: Metallicity
+       mat_v_esc: escape velocity (in km/s)
+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 FOR A GIVEN MASS AND ESCAPE VELOCITY (OR
+HALF-MASS RADIUS BASICALLY)
+For more information on the sampling and how we compute the retained stars see sampling.m
+%}
+
+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):
+
+load('cpdf')
+
+vec_m_initial=[];%We will use this vector to save the different masses. We initialize this variable so that we store every mass.
+%And then we can sample the IMF:
+ M_tot_initial=0;
+jj=1;%We initialize this variable before starting the while loop.
+  while M_tot_initial<M_cl
+    r=rand;
+    kk=find(abs(eval_fun-r)<1e-3);
+    m=vec_M_examples(kk(randi(length(kk))));
+    M_tot_initial=M_tot_initial+m;
+    vec_m_initial(jj)=m;% As it can be seen we store in each row the masses of a given sampling (for a given cluster mass)
+    jj=jj+1;
+  end  
+N_tot_initial=length(vec_m_initial);%In this vector we save the total number of stars.
+%We now have to define the minimum M_zams for NS and BHs:
+vec_M_zams_test=linspace(7,40,2000);
+vec_M_bh_test=fun_M_rem(Z,vec_M_zams_test);
+ll=find(abs(vec_M_bh_test-3)<1e-2);
+M_z_min_bh_test=vec_M_zams_test(ll(1));
+%And to find the minimum NS mass: The minimum mass of 1.17 solar masses is
+%given by Yudai Suwa, Takashi Yoshida, Masaru Shibata, Hideyuki Umeda & Koh
+%Takahashi in "On the minimum mass of neutron stars" in 2009. The problem
+%is that we cannot get M_rem=1.17 with the function we obtained via the
+%PARSEC paper.That is why we set a minimum value of 1.4:
+ll=find(abs(vec_M_bh_test-1.4)<1e-2);
+M_z_min_ns_test=vec_M_zams_test(ll(1));
+
+ 
+%Now, we take the values we obtained before and use them to determine which
+%stars are retained and which ones are not:
+vec_m_final=[];%We initialize the vector that will store the retained masses. 
+%In order to not use so many if statements we are going to find which
+%positions in the vector vec_m_initial have a mass smaller than the minimum
+%neutron star mass. Because this stars are retained and considering them in
+%an if statement only slows the computation:
+pos_stars=vec_m_initial<M_z_min_ns_test;%Position of such stars in the vector vec_m_initial. 
+%(we get 1's in the positions of the stars that verify the requirement, and 0 for the others)
+vec_m_final_i=vec_m_initial.*pos_stars;
+%But we need to eliminate the 0's and only consider those that are
+%different  from 0.
+vec_m_final=nonzeros(vec_m_final_i)';
+%And the vector of masses for the NS and BHs before considering the
+%probabilities for them to be kicked out are:
+vec_bh_ns=vec_m_initial-vec_m_final_i;
+vec_bh_ns=nonzeros(vec_bh_ns)';%Vector that contains the initial stars that are (or will be, basically, as these are ZAMS masses) BHs or NSs.
+%But we want a vector that contains the initial masses of BHs:
+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)';
+
+%And now we do an if statement for this stars:
+ll=length(vec_m_final); %We initialize ll so that we fill the vector after that position.
+%We calculate only once the probability for the initial vector of masses:
+    f_ret=(1-feval(@prob_escape,Z,vec_bh_ns,1.4,265,v_esc,k)); 
+     for jj=1:length(vec_bh_ns)
+         if vec_bh_ns(jj)>M_z_min_bh_test
+            %k=1 ==>with fallback
+            %k=0 ==>without fallback
+            %f_ret=(1-feval(@prob_escape,Z,m,1.4,265,vec_v_esc(dd,ff),1));
+            %We now generate a random number, if it is below f_ret it is
+            %retained, if not it is kicked out of the cluster.
+            rr=rand;
+            if rr<=f_ret(jj)
+                vec_m_final(ll)=vec_bh_ns(jj);
+                ll=ll+1;
+             else
+                ll=ll;
+             end
+           else
+            %Again we compute the probability of it being retained:
+            %f_ret=(1-feval(@prob_esc_ns,265,v_esc(dd,ff)));
+            rr=rand;
+            if rr<=(1-feval(@prob_esc_ns,265,v_esc))
+                vec_m_final(ll)=vec_bh_ns(jj);
+                ll=ll+1;
+            else
+                ll=ll;
+            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)';
+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_v_esc.m b/results/cluster properties/Cluster properties/fun_v_esc.m
new file mode 100644
index 0000000000000000000000000000000000000000..bad101a92443c3b28252b39ed8d0136955345953
--- /dev/null
+++ b/results/cluster properties/Cluster properties/fun_v_esc.m	
@@ -0,0 +1,17 @@
+%% FUNCTION THAT GIVES THE ESCAPE VELOCITY OF THE CLUSTER AS A FUNCTION OF DENSITY AND MASS OF THE CLUSTER
+%Input:
+%      r_h: half-mass radius (in pc)
+%      M_cl: cluster mass (in solar masses)
+
+%Output: 
+%      v_esc: escape velocity (in km/s)
+
+function v_esc=fun_v_esc(r_h,M_cl)
+%rho_h=3*M_cl./(8*pi*(r_h).^3);
+%M_5=M_cl./(10^5);
+%rho_5=rho_h./(10^5);
+f_c=1;
+%v_esc=50*((M_5).^1/3).*((rho_5).^1/6)*f_c;
+v_esc=50*(sqrt(M_cl./(1e5))).*(sqrt(1./r_h))*((3/(8*pi))^(1/6))*f_c;
+end
+
diff --git a/results/cluster properties/Cluster properties/histogram_1.jpg b/results/cluster properties/Cluster properties/histogram_1.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..fa21e587e8c21373de329da41ac42908fb87d061
Binary files /dev/null and b/results/cluster properties/Cluster properties/histogram_1.jpg differ
diff --git a/results/cluster properties/Cluster properties/histogram_2.jpg b/results/cluster properties/Cluster properties/histogram_2.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..0de7bfed223beabfc4633f0b931ea4fd1b538305
Binary files /dev/null and b/results/cluster properties/Cluster properties/histogram_2.jpg differ
diff --git a/results/cluster properties/Cluster properties/histogram_3.jpg b/results/cluster properties/Cluster properties/histogram_3.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..8b40651e596a09bee40747230b5b4714c1bd9afd
Binary files /dev/null and b/results/cluster properties/Cluster properties/histogram_3.jpg differ
diff --git a/results/cluster properties/Cluster properties/histogram_4.jpg b/results/cluster properties/Cluster properties/histogram_4.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..9a8eb40dc12ef92fe9b30d7ea4adfbac6db6ab3c
Binary files /dev/null and b/results/cluster properties/Cluster properties/histogram_4.jpg differ
diff --git a/results/cluster properties/Cluster properties/html/Cluster_conditions.pdf b/results/cluster properties/Cluster properties/html/Cluster_conditions.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..25b0d1f12d0b6a1f40fde3e7cfd5afaadff2ed83
Binary files /dev/null and b/results/cluster properties/Cluster properties/html/Cluster_conditions.pdf differ
diff --git a/results/cluster properties/Cluster properties/html/Cluster_properties_final_version.pdf b/results/cluster properties/Cluster properties/html/Cluster_properties_final_version.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..545de7f7f29c58df021e1d648cf4a9d350a552a0
Binary files /dev/null and b/results/cluster properties/Cluster properties/html/Cluster_properties_final_version.pdf differ
diff --git a/results/cluster properties/Cluster properties/html/intial_conditions.pdf b/results/cluster properties/Cluster properties/html/intial_conditions.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..5e9bc1ea0b7e86686c83d7efaecd40e9280fb4f8
Binary files /dev/null and b/results/cluster properties/Cluster properties/html/intial_conditions.pdf differ
diff --git a/results/cluster properties/Cluster properties/html/sampling.pdf b/results/cluster properties/Cluster properties/html/sampling.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..b66b86549b7d952d4744ce8114c031e9e924b840
Binary files /dev/null and b/results/cluster properties/Cluster properties/html/sampling.pdf differ
diff --git a/results/cluster properties/Cluster properties/imf.m b/results/cluster properties/Cluster properties/imf.m
new file mode 100644
index 0000000000000000000000000000000000000000..a0bef89eac4effaea2862127b7f4ca70ede09067
--- /dev/null
+++ b/results/cluster properties/Cluster properties/imf.m	
@@ -0,0 +1,23 @@
+%% IMF
+%This will be a function that will give us the IMF if we determine the
+%values of alfa in dN/dm proportional to m^(-alfa).
+
+%Input:
+%      alfa1:exponent for M<0.08 solar masses
+%      alfa2:exponent for 0.08<M<0.5 solar masses
+%      alfa3:exponent for M>0.5 solar masses
+%      M_zams: zero age main sequence mass (in solar masses)
+%Output:
+%      phi: dN/dm
+function phi=imf(alfa1,alfa2,alfa3,M_zams)
+for ii=1:length(M_zams)
+    M_zams_i=M_zams(ii);
+if M_zams_i<0.08
+    phi(ii)=(M_zams_i).^(-alfa1);
+elseif M_zams_i<0.5 && M_zams_i>=0.08
+    phi(ii)=((M_zams_i).^(-alfa2))*0.08;
+elseif M_zams_i>=0.5
+    phi(ii)=((M_zams_i).^(-alfa3))*0.5;
+end
+end
+end
diff --git a/results/cluster properties/Cluster properties/intial_condition_2.asv b/results/cluster properties/Cluster properties/intial_condition_2.asv
new file mode 100644
index 0000000000000000000000000000000000000000..352087c30281c9e17f2f5db12012fd98629253a5
--- /dev/null
+++ b/results/cluster properties/Cluster properties/intial_condition_2.asv	
@@ -0,0 +1,128 @@
+clear all; close all; format compact; format long g;
+%% CALCULATION OF THE ESCAPE VELOCITY AS A FUNCTION OF DENSITY AND MASS OF THE CLUSTER
+%We start by defining some vectors for the possible masses and radius that
+%we are going to give to the cluster.
+vec_M=linspace(10^3,10^6,10^3); %In solar masses.
+vec_r=linspace(0.5,10,1000); %In parsecs.
+%We then define a meshgrid for this two vectors:
+[M,r]=meshgrid(vec_M,vec_r);
+%We then use the model to determine the escape velocity from the cluster:
+rho_h=3*M./(8*pi*(r).^3);
+M_5=M./(10^5);
+rho_5=rho_h./(10^5);
+f_c=1;
+v_esc=50*((M_5).^1/3).*((rho_5).^1/6).*f_c;
+%Then we plot it:
+figure(1)
+contour(M_5,r,v_esc,[25,50,100,150],'ShowText','on')
+title('Contour plot for the v escape as a function of M_5 and r_h')
+xlabel('M_5    (M_{cl}/(10^5 M_{sun})')
+ylabel('r (pc)')
+
+%% CONNECTION WITH THE PROBABILITY OF ESCAPE FOR DIFFERENT BLACK HOLE MASSES
+v_esc=linspace(20,200,1000); %In km/s
+M=linspace(3,40,1000); %In solar masses;
+[M,v_esc]=meshgrid(M,v_esc);
+M_ns=1.5; %Mass of the neutron star (in solar masses)
+sig_ns=190;%Velocity dispersion for the neutron star (in km/s);
+p=prob_esc(M_ns,M,sig_ns,v_esc); %Probability of escape
+%And then we plot it:
+figure(2)
+contour(M,v_esc,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Probability of escape for different BH masses and escape velocities')
+xlabel('BH mass (in solar masses)')
+ylabel('Escape velocity (in km/s)')
+
+%% PROBABILITIES OF ESCAPING FOR  A CERTAIN BH MASS IN DIFFERENT CLUSTERS
+%We start with a BH of 10 solar masses:
+v_esc=linspace(0,200,1000); %In km/s;
+M_ns=1.5;M_bh=10; %In solar masses
+sig_ns=190; %In km/s;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc); %Probability of escape
+figure(3)
+plot(v_esc,p,'-k')
+title('Probabilities of escape')
+xlabel('Escape velocity (km/s)')
+ylabel('Probability')
+hold on
+%20 solar masses:
+M_bh=20;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+plot(v_esc,p,'-b')
+hold on
+%30 solar masses:
+M_bh=30;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+plot(v_esc,p,'-r')
+%40 solar masses:
+M_bh=40;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+plot(v_esc,p,'-g')
+legend('M=10','M=20','M=30','M=40')
+
+%% PROBABILITIES OF ESCAPING FOR  A CERTAIN BH MASS IN DIFFERENT CLUSTERS 
+%% BUT NOW AS A FUNCTION OF M_5 AND RHO_5 IN THE FORM OF A CONTOUR PLOT
+vec_M=linspace(10^3,10^6,10^3); %In solar masses.
+vec_r=linspace(0.5,5,1000); %In parsecs.
+%We then define a meshgrid for this two vectors:
+[M,r]=meshgrid(vec_M,vec_r);
+%We then use the model to determine the escape velocity from the cluster:
+rho_h=3*M./(8*pi*(r).^3);
+M_5=M./(10^5);
+rho_5=rho_h./(10^5);
+f_c=1;
+v_esc=50*((M_5).^1/3).*((rho_5).^1/6).*f_c;
+M_ns=1.5;sig_ns=190;
+
+%10 solar mass BH:
+M_bh=10;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(4)
+subplot(2,2,1)
+contourf(M_5,rho_5,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+colorbar('eastoutside')
+title('Prob escape')
+xlabel('M_5')
+ylabel('rho_5')
+legend('M=10')
+
+%20 solar mass BH:
+M_bh=20;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(4)
+subplot(2,2,2)
+contour(M_5,rho_5,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Prob escape')
+xlabel('M_5')
+ylabel('rho_5')
+legend('M=20')
+
+%30 solar mass BH:
+M_bh=30;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(4)
+subplot(2,2,3)
+contour(M_5,rho_5,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Prob escape')
+xlabel('M_5')
+ylabel('rho_5')
+legend('M=30')
+
+%40 solar mass BH:
+M_bh=40;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(4)
+subplot(2,2,4)
+contour(M_5,rho_5,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Prob escape')
+xlabel('M_5')
+ylabel('rho_5')
+legend('M=40')
+
+
+
+
+
+
+
+
diff --git a/results/cluster properties/Cluster properties/intial_condition_2.m b/results/cluster properties/Cluster properties/intial_condition_2.m
new file mode 100644
index 0000000000000000000000000000000000000000..eafc27af69513f2c9c62e866e175eadfd42f3046
--- /dev/null
+++ b/results/cluster properties/Cluster properties/intial_condition_2.m	
@@ -0,0 +1,127 @@
+clear all; close all; format compact; format long g;
+%% CALCULATION OF THE ESCAPE VELOCITY AS A FUNCTION OF DENSITY AND MASS OF THE CLUSTER
+%We start by defining some vectors for the possible masses and radius that
+%we are going to give to the cluster.
+vec_M=linspace(10^3,10^6,10^3); %In solar masses.
+vec_r=linspace(0.5,10,1000); %In parsecs.
+%We then define a meshgrid for this two vectors:
+[M,r]=meshgrid(vec_M,vec_r);
+%We then use the model to determine the escape velocity from the cluster:
+rho_h=3*M./(8*pi*(r).^3);
+M_5=M./(10^5);
+rho_5=rho_h./(10^5);
+f_c=1;
+v_esc=50*((M_5).^(1/3)).*((rho_5).^(1/6)).*f_c;
+%Then we plot it:
+figure(1)
+contour(M_5,r,v_esc,[10,20,30,40,50,75,100],'ShowText','on')
+title('Contour plot for the v escape as a function of M_5 and r_h')
+xlabel('M_5    (M_{cl}/(10^5 M_{sun})')
+ylabel('r_h (pc)')
+
+%% CONNECTION WITH THE PROBABILITY OF ESCAPE FOR DIFFERENT BLACK HOLE MASSES
+v_esc=linspace(20,200,1000); %In km/s
+M=linspace(3,40,1000); %In solar masses;
+[M,v_esc]=meshgrid(M,v_esc);
+M_ns=1.5; %Mass of the neutron star (in solar masses)
+sig_ns=190;%Velocity dispersion for the neutron star (in km/s);
+p=prob_esc(M_ns,M,sig_ns,v_esc); %Probability of escape
+%And then we plot it:
+figure(2)
+contour(M,v_esc,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Probability of escape for different BH masses and escape velocities')
+xlabel('BH mass (in solar masses)')
+ylabel('Escape velocity (in km/s)')
+
+%% PROBABILITIES OF ESCAPING FOR  A CERTAIN BH MASS IN DIFFERENT CLUSTERS
+%We start with a BH of 10 solar masses:
+v_esc=linspace(0,200,1000); %In km/s;
+M_ns=1.5;M_bh=10; %In solar masses
+sig_ns=190; %In km/s;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc); %Probability of escape
+figure(3)
+plot(v_esc,p,'-k')
+title('Probabilities of escape')
+xlabel('Escape velocity (km/s)')
+ylabel('Probability')
+hold on
+%20 solar masses:
+M_bh=20;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+plot(v_esc,p,'-b')
+hold on
+%30 solar masses:
+M_bh=30;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+plot(v_esc,p,'-r')
+%40 solar masses:
+M_bh=40;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+plot(v_esc,p,'-g')
+legend('M=10','M=20','M=30','M=40')
+
+%% PROBABILITIES OF ESCAPING FOR  A CERTAIN BH MASS IN DIFFERENT CLUSTERS 
+%% BUT NOW AS A FUNCTION OF M_5 AND RHO_5 IN THE FORM OF A CONTOUR PLOT
+vec_M=linspace(10^3,10^6,10^3); %In solar masses.
+vec_r=linspace(0.5,5,1000); %In parsecs.
+%We then define a meshgrid for this two vectors:
+[M,r]=meshgrid(vec_M,vec_r);
+%We then use the model to determine the escape velocity from the cluster:
+rho_h=3*M./(8*pi*(r).^3);
+M_5=M./(10^5);
+rho_5=rho_h./(10^5);
+f_c=1;
+v_esc=50*((M_5).^(1/3)).*((rho_5).^(1/6)).*f_c;
+M_ns=1.5;sig_ns=190;
+
+%10 solar mass BH:
+M_bh=10;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(4)
+subplot(2,2,1)
+contour(M_5,r,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Prob escape')
+xlabel('M_5')
+ylabel('r_h (pc)')
+legend('M=10')
+
+%20 solar mass BH:
+M_bh=20;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(4)
+subplot(2,2,2)
+contour(M_5,r,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Prob escape')
+xlabel('M_5')
+ylabel('r_h (pc)')
+legend('M=20')
+
+%30 solar mass BH:
+M_bh=30;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(4)
+subplot(2,2,3)
+contour(M_5,r,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Prob escape')
+xlabel('M_5')
+ylabel('r_h (pc)')
+legend('M=30')
+
+%40 solar mass BH:
+M_bh=40;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(4)
+subplot(2,2,4)
+contour(M_5,r,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Prob escape')
+xlabel('M_5')
+ylabel('r_h (pc)')
+legend('M=40')
+
+
+
+
+
+
+
+
diff --git a/results/cluster properties/Cluster properties/intial_conditions.m b/results/cluster properties/Cluster properties/intial_conditions.m
new file mode 100644
index 0000000000000000000000000000000000000000..f91fc92f2463ff0586d9453599747bdebd7080f9
--- /dev/null
+++ b/results/cluster properties/Cluster properties/intial_conditions.m	
@@ -0,0 +1,131 @@
+clear all; close all; format compact; format long g;
+%% CALCULATION OF THE ESCAPE VELOCITY AS A FUNCTION OF DENSITY AND MASS OF THE CLUSTER
+%We start by defining some vectors for the possible masses and radius that
+%we are going to give to the cluster.
+vec_M=linspace(10^3,10^6,10^3); %In solar masses.
+vec_r=linspace(0.5,5,1000); %In parsecs.
+%We then define a meshgrid for this two vectors:
+[M,r]=meshgrid(vec_M,vec_r);
+%We then use the model to determine the escape velocity from the cluster:
+rho_h=3*M./(8*pi*(r).^3);
+M_5=M./(10^5);
+rho_5=rho_h./(10^5);
+f_c=1;
+v_esc=50*((M_5).^(1/3)).*((rho_5).^(1/6)).*f_c;
+%Then we plot it:
+figure(1)
+contour(M_5,rho_5,v_esc,[25,50,100,150],'ShowText','on')
+title('Contour plot for the v escape as a function of M_5 and r_h')
+xlabel('M_5    (M_{cl}/(10^5 M_{sun})')
+ylabel('rho_5')
+
+%% CONNECTION WITH THE PROBABILITY OF ESCAPE FOR DIFFERENT BLACK HOLE MASSES
+v_esc=linspace(20,200,1000); %In km/s
+M=linspace(3,40,1000); %In solar masses;
+[M,v_esc]=meshgrid(M,v_esc);
+M_ns=1.5; %Mass of the neutron star (in solar masses)
+sig_ns=190;%Velocity dispersion for the neutron star (in km/s);
+p=prob_esc(M_ns,M,sig_ns,v_esc); %Probability of escape
+%And then we plot it:
+figure(2)
+contour(M,v_esc,p,[0.01, 0.25, 0.50, 0.75, 0.9],'ShowText','on')
+title('Probability of escape for different BH masses and escape velocities')
+xlabel('BH mass (in solar masses)')
+ylabel('Escape velocity (in km/s)')
+
+%% PROBABILITIES OF ESCAPING FOR  A CERTAIN BH MASS IN DIFFERENT CLUSTERS
+%We start with a BH of 10 solar masses:
+v_esc=linspace(0,200,1000); %In km/s;
+M_ns=1.5;M_bh=10; %In solar masses
+sig_ns=190; %In km/s;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc); %Probability of escape
+figure(3)
+plot(v_esc,p,'-k')
+title('Probabilities of escape')
+xlabel('Escape velocity (km/s)')
+ylabel('Probability')
+hold on
+%20 solar masses:
+M_bh=20;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+plot(v_esc,p,'-b')
+hold on
+%30 solar masses:
+M_bh=30;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+plot(v_esc,p,'-r')
+%40 solar masses:
+M_bh=40;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+plot(v_esc,p,'-g')
+legend('M=10','M=20','M=30','M=40')
+
+%% PROBABILITIES OF ESCAPING FOR  A CERTAIN BH MASS IN DIFFERENT CLUSTERS 
+%% BUT NOW AS A FUNCTION OF M_5 AND RHO_5 IN THE FORM OF A CONTOUR PLOT
+vec_M=linspace(10^3,10^6,10^3); %In solar masses.
+vec_r=linspace(0.5,5,1000); %In parsecs.
+%We then define a meshgrid for this two vectors:
+[M,r]=meshgrid(vec_M,vec_r);
+%We then use the model to determine the escape velocity from the cluster:
+rho_h=3*M./(8*pi*(r).^3);
+M_5=M./(10^5);
+rho_5=rho_h./(10^5);
+f_c=1;
+v_esc=50*((M_5).^(1/3)).*((rho_5).^(1/6)).*f_c;
+M_ns=1.5;sig_ns=190;
+
+%10 solar mass BH:
+M_bh=10;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(4)
+subplot(2,2,1)
+contourf(M_5,rho_5,p,'ShowText','on')
+colorbar('eastoutside')
+title('Prob escape')
+xlabel('M_5')
+ylabel('rho_5')
+legend('M=10')
+
+%20 solar mass BH:
+M_bh=20;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(4)
+subplot(2,2,2)
+contourf(M_5,rho_5,p,'ShowText','on')
+colorbar('eastoutside')
+title('Prob escape')
+xlabel('M_5')
+ylabel('rho_5')
+legend('M=20')
+
+%30 solar mass BH:
+M_bh=30;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(4)
+subplot(2,2,3)
+contourf(M_5,rho_5,p)
+colorbar('eastoutside')
+title('Prob escape')
+xlabel('M_5')
+ylabel('rho_5')
+legend('M=30')
+
+%40 solar mass BH:
+M_bh=40;
+p=prob_esc(M_ns,M_bh,sig_ns,v_esc);
+figure(4)
+subplot(2,2,4)
+contourf(M_5,rho_5,p)
+colorbar('eastoutside')
+title('Prob escape')
+xlabel('M_5')
+ylabel('rho_5')
+legend('M=40')
+
+
+
+
+
+
+
+
diff --git a/results/cluster properties/Cluster properties/mass_fraction_black_holes.fig b/results/cluster properties/Cluster properties/mass_fraction_black_holes.fig
new file mode 100644
index 0000000000000000000000000000000000000000..5487f010bd8600b7f990673bd6632221008bbf6c
Binary files /dev/null and b/results/cluster properties/Cluster properties/mass_fraction_black_holes.fig differ
diff --git a/results/cluster properties/Cluster properties/mass_fraction_black_holes.jpg b/results/cluster properties/Cluster properties/mass_fraction_black_holes.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..ac71dc8344d69c9865bc1489e4c06a70f4ab0514
Binary files /dev/null and b/results/cluster properties/Cluster properties/mass_fraction_black_holes.jpg differ
diff --git a/results/cluster properties/Cluster properties/number_of_black_holes.fig b/results/cluster properties/Cluster properties/number_of_black_holes.fig
new file mode 100644
index 0000000000000000000000000000000000000000..89f8c4c232bfcf9e6913493aacf6fea06365bf95
Binary files /dev/null and b/results/cluster properties/Cluster properties/number_of_black_holes.fig differ
diff --git a/results/cluster properties/Cluster properties/number_of_black_holes.jpg b/results/cluster properties/Cluster properties/number_of_black_holes.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..ac22fe2cd572ece03648de1df9a84f70e420a26f
Binary files /dev/null and b/results/cluster properties/Cluster properties/number_of_black_holes.jpg differ
diff --git a/results/cluster properties/Cluster properties/number_of_stars.fig b/results/cluster properties/Cluster properties/number_of_stars.fig
new file mode 100644
index 0000000000000000000000000000000000000000..55d46aaf1f5203265883141ca9f1c23bbf3f3e7b
Binary files /dev/null and b/results/cluster properties/Cluster properties/number_of_stars.fig differ
diff --git a/results/cluster properties/Cluster properties/number_of_stars.jpg b/results/cluster properties/Cluster properties/number_of_stars.jpg
new file mode 100644
index 0000000000000000000000000000000000000000..42d3f45c2dcd509b0953e0997d0d700da3588a14
Binary files /dev/null and b/results/cluster properties/Cluster properties/number_of_stars.jpg differ
diff --git a/results/cluster properties/Cluster properties/p_times_imf.m b/results/cluster properties/Cluster properties/p_times_imf.m
new file mode 100644
index 0000000000000000000000000000000000000000..b8be1b7e7ed080869e75e140252d7cf8ff24113e
--- /dev/null
+++ b/results/cluster properties/Cluster properties/p_times_imf.m	
@@ -0,0 +1,5 @@
+function p_imf=p_times_imf(Z,alfa1,alfa2,M_zams,M_ns,sig_ns,v_esc)
+p=feval(@prob_esc_with_fb,Z,M_zams,M_ns,sig_ns,v_esc);
+imf= @(M) (M<0.5).*(M.^(-alfa1)) + (M>=0.5).*(M.^(-alfa1));
+p_imf=p*imf(M_zams);
+end
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/precompute_probabilities.m b/results/cluster properties/Cluster properties/precompute_probabilities.m
new file mode 100644
index 0000000000000000000000000000000000000000..85a3b784861f3ce1346f182d9b09f893983ca092
--- /dev/null
+++ b/results/cluster properties/Cluster properties/precompute_probabilities.m	
@@ -0,0 +1,14 @@
+clear all;close all; format compact; format long g;
+%% SCRIPT THAT PRECOMPUTES PROBABILITIES.
+load('cpdf')
+%% Z=2e-2:
+Z=2e-2;
+vec_v_esc=[10 15 20 25 30 35 40 45 50];
+%We need to decide which possible escape velocities we will give:
+for ii=1:length(vec_v_esc)
+    v_esc=vec_v_esc(ii);
+    prob_esc(ii,:)=feval(@prob_escape,2e-2,vec_M_examples,1.4,265,v_esc,1);
+    vec_f_ret(ii,:)=1-prob_esc(ii,:);
+    vec_f_ret_ns(ii)=1-feval(@prob_esc_ns,265,v_esc);
+end
+    
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/prob_esc.m b/results/cluster properties/Cluster properties/prob_esc.m
new file mode 100644
index 0000000000000000000000000000000000000000..9a4962c15260dce8104ad22aef9f35049a8807a8
--- /dev/null
+++ b/results/cluster properties/Cluster properties/prob_esc.m	
@@ -0,0 +1,27 @@
+%% FUNCTION THAT GIVES THE PROBABILITY OF v>vesc
+%Input:
+%                M_ns: neutron star mass (in solar masses)
+%                M_bh: black hole mass (in solar masses)
+%                sig_ns: neutron star velocity dispersion (in km/s)
+%                v_esc: escape velocity
+
+%Output:
+%                p: probability of v>vesc
+function p=prob_esc(M_ns,M_bh,sig_ns,v_esc)
+sig_bh=(M_ns./M_bh).*sig_ns;
+%{
+%We could maybe try to implement fallback by doing:
+if M_bh<15 
+    f=0;
+else
+    f=(1/10)*atan((M_bh-15)/15);
+end
+p=sqrt(2/pi)*(v_esc./sig_bh).*exp(-(v_esc.^2)./(2*(sig_bh).^2))+(1-erf(v_esc./(sqrt(2)*sig_bh)))-f;
+if p>0
+    p=p;
+else 
+    p=0;
+end
+%}
+p=sqrt(2/pi)*(v_esc./sig_bh).*exp(-(v_esc.^2)./(2*(sig_bh).^2))+(1-erf(v_esc./(sqrt(2)*sig_bh)));
+end
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/prob_esc_ns.m b/results/cluster properties/Cluster properties/prob_esc_ns.m
new file mode 100644
index 0000000000000000000000000000000000000000..8b3de3656d0520719f021fba5d548fb5df239618
--- /dev/null
+++ b/results/cluster properties/Cluster properties/prob_esc_ns.m	
@@ -0,0 +1,10 @@
+%% FUNCTION THAT GIVES THE PROBABILITY OF v>vesc FOR A NEUTRON STAR (NS)
+%Input:
+%                sig_ns: neutron star velocity dispersion (in km/s)
+%                v_esc: escape velocity
+
+%Output:
+%                p: probability of v>vesc for a neutron star.
+function p=prob_esc_ns(sig_ns,v_esc)
+p=sqrt(2/pi)*(v_esc./sig_ns).*exp(-(v_esc.^2)./(2*(sig_ns).^2))+(1-erf(v_esc./(sqrt(2)*sig_ns)));
+end
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/prob_esc_with_fb.m b/results/cluster properties/Cluster properties/prob_esc_with_fb.m
new file mode 100644
index 0000000000000000000000000000000000000000..0c0c8246cb1aca6231c90bacdfd9d0f75023b104
--- /dev/null
+++ b/results/cluster properties/Cluster properties/prob_esc_with_fb.m	
@@ -0,0 +1,22 @@
+%% FUNCTION THAT GIVES THE PROBABILITY OF v>vesc TAKING INTO ACCOUNT FALLBACK f_fb
+%Input:
+%                Z: Metallicity
+%                M_zams:Zero age main sequence mass (in solar masses)
+%                M_ns: neutron star mass (in solar masses)
+%                sig_ns: neutron star velocity dispersion (in km/s)
+%                v_esc: escape velocity
+
+%Output:
+%                p: probability of v>v_esc
+function p=prob_esc_with_fb(Z,M_zams,M_ns,sig_ns,v_esc)
+%In order to compute the BH mass we make use of fun_M_rem.m that gives the
+%mass of the stellar remnant as a function of Z and M_zams.
+M_bh=fun_M_rem(Z,M_zams);
+%And then we can obtain the BH's velocity dispersion.
+sig_bh=(M_ns./M_bh).*sig_ns;
+%We define the fallback fraction:
+    f_fb=fun_fb(Z,M_zams);
+%And then we can compute the probability of escaping:
+p=sqrt(2/pi)*(v_esc./sig_bh).*exp(-(v_esc.^2)./(2*((1-f_fb).^2).*(sig_bh).^2))+(1-f_fb).*(1-erf(v_esc./(sqrt(2)*(1-f_fb).*sig_bh)));
+end
+
diff --git a/results/cluster properties/Cluster properties/prob_escape.m b/results/cluster properties/Cluster properties/prob_escape.m
new file mode 100644
index 0000000000000000000000000000000000000000..7744be78c4400365a61292b421761eab07ddd46c
--- /dev/null
+++ b/results/cluster properties/Cluster properties/prob_escape.m	
@@ -0,0 +1,31 @@
+%% FUNCTION THAT GIVES THE PROBABILITY OF v>vesc TAKING INTO ACCOUNT FALLBACK f_fb
+%Input:
+%                Z: Metallicity
+%                M_zams:Zero age main sequence mass (in solar masses)
+%                M_ns: neutron star mass (in solar masses)
+%                sig_ns: neutron star velocity dispersion (in km/s)
+%                v_esc: escape velocity
+%                k: parameter that defines whether there is fallback or not. 
+%                       k=0 ==> no fallback
+%                       k=1 ==> fallback considered
+
+%Output:
+%                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.
+M_bh=fun_M_rem(Z,M_zams);
+%And then we can obtain the BH's velocity dispersion.
+sig_bh=(M_ns./M_bh).*sig_ns;
+%We define the fallback fraction:
+if k==1
+    f_fb=fun_fb(Z,M_zams);
+elseif k==0
+    f_fb=0;
+elseif k==2
+    f_fb=fun_fb(Z,M_zams);
+    sig_bh=sig_ns;
+end
+%And then we can compute the probability of escaping:
+p=sqrt(2/pi)*(v_esc./sig_bh).*exp(-(v_esc.^2)./(2*((1-f_fb).^2).*(sig_bh).^2))+(1-f_fb).*(1-erf(v_esc./(sqrt(2)*(1-f_fb).*sig_bh)));
+end
diff --git a/results/cluster properties/Cluster properties/prueba.asv b/results/cluster properties/Cluster properties/prueba.asv
new file mode 100644
index 0000000000000000000000000000000000000000..5d94da56820c119b63c806877554d28d38c36046
--- /dev/null
+++ b/results/cluster properties/Cluster properties/prueba.asv	
@@ -0,0 +1,4 @@
+clear all; close all; format compact; format long g;
+%% Prueba
+x=linspace(0,100,10000);
+plot(x,1-exp(-x)/(40+exp(-x)))
\ No newline at end of file
diff --git a/results/cluster properties/Cluster properties/prueba.m b/results/cluster properties/Cluster properties/prueba.m
new file mode 100644
index 0000000000000000000000000000000000000000..fd980b71066d4a0ee8621e3a4532f2530f0eac7d
--- /dev/null
+++ b/results/cluster properties/Cluster properties/prueba.m	
@@ -0,0 +1,2 @@
+clear all; close all; format compact; format long g;
+%% Prueba
diff --git a/results/cluster properties/Cluster properties/prueba_fit_poly.m b/results/cluster properties/Cluster properties/prueba_fit_poly.m
new file mode 100644
index 0000000000000000000000000000000000000000..24777a43633ae2d1878ce4f30ae1cdbde80e1cae
--- /dev/null
+++ b/results/cluster properties/Cluster properties/prueba_fit_poly.m	
@@ -0,0 +1,10 @@
+clear all; close all; format compact; format long g;
+Z=2e-2;
+%But we also need dMzams/dMbh so:
+vec_M_zams=linspace(10,100,2000);
+vec_M_bh=fun_M_rem(Z,vec_M_zams);
+%And then we fit this into a polynomial:
+p = polyfit(vec_M_bh,vec_M_zams,15);%In decreasing power
+plot(vec_M_bh,vec_M_zams,'-r')
+hold on
+plot(vec_M_bh,polyval(p,vec_M_bh),'-b')
diff --git a/results/cluster properties/Cluster properties/sampling.m b/results/cluster properties/Cluster properties/sampling.m
new file mode 100644
index 0000000000000000000000000000000000000000..2d039a23aa24319333675168d60509f0fc3559c8
--- /dev/null
+++ b/results/cluster properties/Cluster properties/sampling.m	
@@ -0,0 +1,397 @@
+clearvars -except eval_fun vec_M_examples
+%% INFORMATION ABOUT THE SCRIPT
+%INPUT:
+%      Z: metallicity
+%      v_esc: escape velocity (in km/s)
+%OUTPUT:
+%{
+We get the number of stars, NSs and BHs for 100 different samplings, as
+well as the mass and number fractions of NSs and BHs.
+We obtain this results considering only the IMF, but also for the cluster
+after the supernova kicks. We do this by considering the scenario with
+fallback, but also without it.
+%}
+%We then define the metallicity:
+Z=2e-2;
+%And the escape velocity:
+v_esc=30; %km/s.
+%We also need the cumulative probability distribution function (cpdf) so we
+%load the results obtained from another script:
+load('cpdf')
+%% SAMPLING OF THE IMF
+%{
+%We start by determining the cumulative probability densitiy function (cpdf) for the masses in the vector vec_M_examples.
+This is a costly process, so we load theresults from a previous
+calculation by writing load('cpdf'). The idea of the process is then to
+asign a random value between 0 and 1 for every star we "generate" so that
+then we find which mass yields a value for the cpdf most similar to the
+random value generated.
+%}
+%Below we show what that process is:
+%{
+%We start by defining the IMF:
+alfa1=0.3;alfa2=1.3;alfa3=2.3;
+imf= @(M) (M<0.08).*(M.^(-alfa1)) + (M<0.5 & M>=0.08).*((M.^(-alfa2))*(0.08)) + (M>=0.5).*((M.^(-alfa3))*(0.5));
+%We then define the maximum value for the cumulative probability density
+%function so that we can normalise the function itself so it can be treated
+%as a probability.
+cpdf_max=integral(imf,0.07,150,'RelTol',1e-6,'AbsTol',1e-6);
+%Once we have this function defined we need to invert it. That it, we need
+%to have it such that we give a random number from 0 to 1 and it returns a
+%mass. So that we asign different masses until we reach the maximum value
+%we want it to have.
+%For example:
+%We initialize the total mass to be 0 before we start "creating" stars:
+M_tot=0;
+vec_M_examples=[linspace(0.07,2,15000) linspace(2,150,25000)];
+eval_fun=[];
+for ii=1:length(vec_M_examples)
+    M_example=vec_M_examples(ii);
+    fun=@(M) (integral(imf,0.07,M,'RelTol',1e-6,'AbsTol',1e-6)/(cpdf_max));
+    eval_fun(ii)=feval(fun,M_example);
+end
+%}
+%Below we have the calculations followed in order to plot the histogram for
+%the stars "generated" in a single sampling of the IMF. This single result
+%does not show the stochastic nature of this process, so we will focus on
+%doing 100 samplings in order to obtain the mean value and the sigma for
+%the results we obtain. However, the process is described below anyway.
+%{
+M_tot=0;
+vec_m=[];
+jj=1;
+while M_tot<10^5
+    r=rand;
+    kk=find(abs(eval_fun-r)<1e-3);
+    m=vec_M_examples(kk(randi(length(kk))));
+    vec_m(jj)=m;
+    M_tot=M_tot+m;
+    jj=jj+1;
+end
+
+%For stars less massive than 5 solar masses:
+less=find(vec_m<25);
+vec_m_stars=vec_m(less);
+figure(1)
+histogram(vec_m_stars)
+title('Number of stars for M_{cl}=10^5')
+xlabel('M_{zams}')
+ylabel('Number of stars')
+set(gca,'YScale','log')
+massive=find(vec_m>25);
+vec_m_massive=vec_m(massive);
+figure(2)
+title('Number of stars for M_{cl}=10^5')
+xlabel('M_{zams}')
+ylabel('Number of stars')
+histogram(vec_m_massive)
+
+%% FRACTIONS
+
+%Total number of stars:
+N_tot=length(vec_m);
+%We define the metallicity so we find M_z_min that gives M_rem=3 solar
+%masses:
+Z=2e-2;
+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:
+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 need to transform this test values into the most similar ones
+%for the given masses we have:
+M_z_min_bh=min(vec_m-M_z_min_bh_test);
+M_z_min_ns=min(vec_m-M_z_min_ns_test);
+%}
+%And then we compute the total number of stars of each of them:
+N_bh=sum(vec_m>M_z_min_bh_test);
+N_stars=sum(vec_m<M_z_min_ns_test);
+N_ns=N_tot-N_bh-N_stars;
+f_N_bh=N_bh/N_tot;f_N_ns=N_ns/N_tot;f_N_stars=N_stars/N_tot;
+results_N_before=[N_bh/N_tot N_stars/N_tot N_ns/N_tot];
+figure(3)
+label={'BH','stars','NS'};
+pie(results_N_before,label)
+title('Number fraction of stars (just IMF)')
+%And then we do the same for the mass:
+M_bh=sum(vec_m(find(vec_m>M_z_min_bh_test)));
+M_stars=sum(vec_m(find(vec_m<M_z_min_ns_test)));
+M_ns=M_tot-M_bh-M_stars;
+f_m_bh=M_bh/M_tot;f_m_ns=M_ns/M_tot;f_m_stars=M_stars/M_tot;
+results_M_before=[M_bh/M_tot M_stars/M_tot M_ns/M_tot];
+figure(4)
+label={'BH','stars','NS'};
+pie(results_M_before,label)
+title('Number fraction of stars (just IMF)')
+%}
+%% FLUCTUATIONS IN THE TOTAL NUMBER OF STARS AND DIFFERENT FRACTIONS (JUST THE IMF)
+%We start by computing the number of stars, number of BHs and NS, and their
+%mass and number fractions for different samplings in order to observe the
+%stochastic nature of the sampling itself.
+%We first start by studying this upon formation, that is, just focusing on
+%the IMF:
+N_tot=[];
+vec_ii=[1:100];
+mat_m=[];%We save by rows the different samplings taken.
+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<10^5
+    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.
+    mat_m(ii,jj)=m;
+    M_tot=M_tot+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 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 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);
+f_N_bh(ii)=N_bh(ii)/N_tot(ii);f_N_ns(ii)=N_ns(ii)/N_tot(ii);f_N_stars(ii)=N_stars(ii)/N_tot(ii);
+
+M_bh(ii)=sum(vec_m(find(vec_m>M_z_min_bh_test)));
+M_stars(ii)=sum(vec_m(find(vec_m<M_z_min_ns_test)));
+M_ns(ii)=M_tot-M_bh(ii)-M_stars(ii);
+f_m_bh(ii)=M_bh(ii)/M_tot;f_m_ns(ii)=M_ns(ii)/M_tot;f_m_stars=M_stars(ii)/M_tot;
+end
+figure(5)
+semilogy(vec_ii,N_tot,'-k')
+title('Number of stars in a cluster of M_{cl}=10^5 M_{sun} (just IMF)')
+xlabel('Different samplings')
+ylabel('Number of stars')
+figure(6)
+semilogy(vec_ii,N_bh,'-b')
+title('Number of stars in a cluster of M_{cl}=10^5 M_{sun} (just IMF)')
+xlabel('Different samplings')
+ylabel('Number of stars')
+hold on
+semilogy(vec_ii,N_ns,'-r')
+legend('N_{bh}','N_{ns}')
+figure(7)
+plot(vec_ii,f_m_bh,'-b')
+title('Mass fractions in a cluster of M_{cl}=10^5 M_{sun} (just IMF)')
+xlabel('Different samplings')
+ylabel('Mass fractions')
+hold on
+plot(vec_ii,f_m_ns,'-r')
+legend('f_{bh}','f_{ns}')
+
+%% FLUCTUATIONS IN THE TOTAL NUMBER OF STARS AND DIFFERENT FRACTIONS (AFTER SUPERNOVA KICKS) (WITH FALLBACK)
+%We now have to take every mass one by one and see, given the probabilities
+%that we have computed, if they are retained or not. For this we make use
+%of the masses that we have computed in the previous section. It is
+%important to remember that we were storing in each row of mat_m the masses
+%for a given sampling. It is also important to remember that we do not have
+%the same number of stars in each sampling, and that matlab (given we have
+%formed a matrix) fills the remaining positions with 0 so that they all
+%have the same length. That is why the first thing to do, when studying one
+%sampling, is to eliminate those 0's, and we can do that easily because we
+%know N_tot in every iteration.
+for ii=1:length(vec_ii)
+    vec_m=mat_m(ii,1:N_tot(ii));
+    vec_m_final=[];
+%Now that we have defined the mass vector we make use of the probabilities:
+    ll=1;
+    for jj=1:length(vec_m)
+        m=vec_m(jj);
+        if m<M_z_min_ns_test
+            vec_m_final(ll)=m;
+            ll=ll+1;%Because it is always retained.
+        elseif m>M_z_min_bh_test
+            %k=1 ==>with fallback
+            f_ret=(1-feval(@prob_escape,Z,m,1.4,265,v_esc,1));
+            %We now generate a random number, if it is below f_ret it is
+            %retained, if not it is kicked out of the cluster.
+            rr=rand;
+            if rr<=f_ret
+                vec_m_final(ll)=m;
+                ll=ll+1;
+            else
+                ll=ll;
+            end
+        else
+            %Again we compute the probability of it being retained:
+            f_ret=(1-feval(@prob_esc_ns,265,v_esc));
+            rr=rand;
+            if rr<=f_ret
+                vec_m_final(ll)=m;
+                ll=ll+1;
+            else
+                ll=ll;
+            end
+        end
+    end
+N_tot_final(ii)=length(vec_m_final);
+M_tot_final(ii)=sum(vec_m_final);
+N_bh_final(ii)=sum(vec_m_final>M_z_min_bh_test);
+N_stars_final(ii)=sum(vec_m_final<M_z_min_ns_test);
+N_ns_final(ii)=N_tot_final(ii)-N_bh_final(ii)-N_stars_final(ii);
+f_N_bh_final(ii)=N_bh_final(ii)/N_tot_final(ii);f_N_ns_final(ii)=N_ns_final(ii)/N_tot_final(ii);f_N_stars_final(ii)=N_stars_final(ii)/N_tot_final(ii);
+
+M_bh_final(ii)=sum(vec_m_final(find(vec_m_final>M_z_min_bh_test)));
+M_stars_final(ii)=sum(vec_m_final(find(vec_m_final<M_z_min_ns_test)));
+M_ns_final(ii)=M_tot_final(ii)-M_bh_final(ii)-M_stars_final(ii);
+f_m_bh_final(ii)=M_bh_final(ii)/M_tot_final(ii);f_m_ns_final(ii)=M_ns_final(ii)/M_tot_final(ii);f_m_stars_final=M_stars_final(ii)/M_tot_final(ii);
+end
+figure(8)
+semilogy(vec_ii,N_tot_final,'-k')
+title('Number of stars in a cluster of M_{cl}=10^5 M_{sun} (after supernova kicks)')
+xlabel('Different samplings')
+ylabel('Number of stars')
+figure(9)
+semilogy(vec_ii,N_bh_final,'-b')
+title('Number of stars in a cluster of M_{cl}=10^5 M_{sun} (after supernova kicks)')
+xlabel('Different samplings')
+ylabel('Number of stars')
+hold on
+semilogy(vec_ii,N_ns_final,'-r')
+legend('N_{bh}','N_{ns}')
+figure(10)
+plot(vec_ii,f_m_bh_final,'-b')
+title('Mass fractions in a cluster of M_{cl}=10^5 M_{sun} (after supernova kicks)')
+xlabel('Different samplings')
+ylabel('Mass fractions')
+hold on
+plot(vec_ii,f_m_ns_final,'-r')
+legend('f_{bh}','f_{ns}')
+
+
+%% MEAN VALUES AND DEVIATIONS WITH FALLBACK
+%IMF:
+mean_f_m_bh_initial=sum(f_m_bh)/100;
+disp('Mean BH mass fraction (IMF)')
+disp(mean_f_m_bh_initial)
+var=(f_m_bh-mean_f_m_bh_initial).^2;
+sigma_f_m_bh_initial=sqrt(sum(var)/100);
+disp('Mean deviation in the BH mass fraction (IMF)')
+disp(sigma_f_m_bh_initial)
+
+%After supernova kicks:
+mean_f_m_bh_final=sum(f_m_bh_final)/100;
+disp('Mean BH mass fraction (after kicks)')
+disp(mean_f_m_bh_final)
+var=(f_m_bh_final-mean_f_m_bh_final).^2;
+sigma_f_m_bh_final=sqrt(sum(var)/100);
+disp('Mean deviation in the BH mass fraction (after kicks)')
+disp(sigma_f_m_bh_final)
+
+%% FLUCTUATIONS IN THE TOTAL NUMBER OF STARS AND DIFFERENT FRACTIONS (AFTER SUPERNOVA KICKS) (WITHOUT FALLBACK)
+%We do the same again, after supernova kicks, but in this case without
+%considering fallback. In this case we would expect a higher number (and
+%fraction as a consequence) of BHs lost compared to the results obtained in
+%the previous section (with fallback implemented):
+for ii=1:length(vec_ii)
+    vec_m=mat_m(ii,1:N_tot(ii));
+    vec_m_final=[];
+%Now that we have defined the mass vector we make use of the probabilities:
+    ll=1;
+    for jj=1:length(vec_m)
+        m=vec_m(jj);
+        if m<M_z_min_ns_test
+            vec_m_final(ll)=m;
+            ll=ll+1;%Because it is always retained.
+        elseif m>M_z_min_bh_test
+            %k=0 ==>without fallback
+            f_ret=(1-feval(@prob_escape,Z,m,1.4,265,v_esc,0));
+            %We now generate a random number, if it is below f_ret it is
+            %retained, if not it is kicked out of the cluster.
+            rr=rand;
+            if rr<=f_ret
+                vec_m_final(ll)=m;
+                ll=ll+1;
+            else
+                ll=ll;
+            end
+        else
+            %Again we compute the probability of it being retained:
+            f_ret=(1-feval(@prob_esc_ns,265,v_esc));
+            rr=rand;
+            if rr<=f_ret
+                vec_m_final(ll)=m;
+                ll=ll+1;
+            else
+                ll=ll;
+            end
+        end
+    end
+N_tot_final_wo_fb(ii)=length(vec_m_final);
+M_tot_final_wo_fb(ii)=sum(vec_m_final);
+N_bh_final_wo_fb(ii)=sum(vec_m_final>M_z_min_bh_test);
+N_stars_final_wo_fb(ii)=sum(vec_m_final<M_z_min_ns_test);
+N_ns_final_wo_fb(ii)=N_tot_final_wo_fb(ii)-N_bh_final_wo_fb(ii)-N_stars_final_wo_fb(ii);
+f_N_bh_final_wo_fb(ii)=N_bh_final_wo_fb(ii)/N_tot_final_wo_fb(ii);f_N_ns_final_wo_fb(ii)=N_ns_final_wo_fb(ii)/N_tot_final_wo_fb(ii);f_N_stars_final_wo_fb(ii)=N_stars_final_wo_fb(ii)/N_tot_final_wo_fb(ii);
+
+M_bh_final_wo_fb(ii)=sum(vec_m_final(find(vec_m_final>M_z_min_bh_test)));
+M_stars_final_wo_fb(ii)=sum(vec_m_final(find(vec_m_final<M_z_min_ns_test)));
+M_ns_final_wo_fb(ii)=M_tot_final_wo_fb(ii)-M_bh_final_wo_fb(ii)-M_stars_final_wo_fb(ii);
+f_m_bh_final_wo_fb(ii)=M_bh_final_wo_fb(ii)/M_tot_final_wo_fb(ii);f_m_ns_final_wo_fb(ii)=M_ns_final_wo_fb(ii)/M_tot_final_wo_fb(ii);f_m_stars_final_wo_fb=M_stars_final_wo_fb(ii)/M_tot_final_wo_fb(ii);
+end
+%Now we compare the results without fallback with the ones obtained for the
+%IMF and the case with fallback considered:
+%And now we compare the initial and final state
+figure(11)
+plot(vec_ii,N_tot,'-k')
+title(['Number of stars in a cluster of M_{cl}=10^5 M_{sun}  (Z=' num2str(Z) ' & v_{esc}=' num2str(v_esc) 'km/s)'])
+xlabel('Different samplings')
+ylabel('Number of stars')
+hold on
+plot(vec_ii,N_tot_final,'-r')
+hold on
+plot(vec_ii,N_tot_final_wo_fb,'-b')
+legend('IMF','w/ fb','w/o fb')
+
+figure(12)
+plot(vec_ii,N_bh,'-k')
+title(['Number of BHs in a cluster of M_{cl}=10^5 M_{sun}  (Z=' num2str(Z) ' & v_{esc}=' num2str(v_esc) 'km/s)'])
+xlabel('Different samplings')
+ylabel('Number of BHs')
+hold on
+plot(N_bh_final,'-r')
+hold on
+plot(N_bh_final_wo_fb,'-b')
+legend('IMF','w/ fb','w/o fb')
+
+figure(13)
+plot(vec_ii,f_m_bh,'-k')
+title(['Mass fraction of BHs in a cluster of M_{cl}=10^5 M_{sun}  (Z=' num2str(Z) ' & v_{esc}=' num2str(v_esc) 'km/s)'])
+xlabel('Different samplings')
+ylabel('Mass fraction of BHs')
+hold on
+plot(vec_ii,f_m_bh_final,'-r')
+hold on
+plot(vec_ii,f_m_bh_final_wo_fb,'-b')
+legend('IMF','w/ fb','w/o fb')
+
+%% MEAN VALUES AND DEVIATIONS WITHOUT FALLBACK
+
+%After supernova kicks: (without fallback)
+mean_f_m_bh_final_wo_fb=sum(f_m_bh_final_wo_fb)/100;
+disp('Mean BH mass fraction (after kicks) (w/o fallback)')
+disp(mean_f_m_bh_final_wo_fb)
+var=(f_m_bh_final_wo_fb-mean_f_m_bh_final).^2;
+sigma_f_m_bh_final_wo_fb=sqrt(sum(var)/100);
+disp('Mean deviation in the BH mass fraction (after kicks) (w/o fallback)')
+disp(sigma_f_m_bh_final_wo_fb)
+
+
diff --git a/results/cluster properties/Cluster properties/study_variation_of_parameters.m b/results/cluster properties/Cluster properties/study_variation_of_parameters.m
new file mode 100644
index 0000000000000000000000000000000000000000..864129fb373489647a108f008f4dca09966e5e30
--- /dev/null
+++ b/results/cluster properties/Cluster properties/study_variation_of_parameters.m	
@@ -0,0 +1,239 @@
+clear all
+%% ANALYSIS OF THE CONSEQUENCES OF THE VARIATION OF PARAMETERS IN THE RESULTS OBTAINED VIA IMF SAMPLING:
+%First of all we need to load the values for the cumulative probability
+%distribution function (cpdf) corresponding to the vector of possible
+%masses vec_M_examples
+tic
+
+load('cpdf')
+load('precomp_probs')
+%For doubts regarding the nature of these calculations see sampling.m
+%{
+tic
+%% VARATIONS IN THE ESCAPE VELOCITY  (M_cl=10^5 M_sun)
+%We start by considering how does a variation in the escape velocity affect
+%the retention fraction of BHs. Later we will try to connect this to the
+%half-mass radius r_h so that we get a sense of the cluster itself.
+
+%First of all we need to sample the IMF. Once that is done we can start
+%considering what are the effects of a variation in the escape velocity.
+
+
+%We start by defining a vector for the values of the escape velocity we
+%want to consider so that then we can treat 
+%We consider the following metallicity:
+vec_v_esc=[10 15 20 25 30 35 40 45 50];
+Z=2e-2;
+
+jj=1;%We initialize this variable before starting the while loop.
+vec_m_initial=[];%We will use this vector to save the different masses.
+M_tot_initial=0;%We initialize M_tot before "filling" it with stars as we go along.
+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  
+N_tot_initial=length(vec_m_initial);
+%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 stars of each of them:
+N_bh_initial=sum(vec_m_initial>M_z_min_bh_test);
+N_stars_initial=sum(vec_m_initial<M_z_min_ns_test);
+N_ns_initial=N_tot_initial-N_bh_initial-N_stars_initial;
+f_N_bh_initial=N_bh_initial/N_tot_initial;f_N_ns_initial=N_ns_initial/N_tot_initial;f_N_stars_initial=N_stars_initial/N_tot_initial;
+
+M_bh_initial=sum(vec_m_initial(find(vec_m_initial>M_z_min_bh_test)));
+M_stars_initial=sum(vec_m_initial(find(vec_m_initial<M_z_min_ns_test)));
+M_ns_initial=M_tot_initial-M_bh_initial-M_stars_initial;
+f_m_bh_initial=M_bh_initial/M_tot_initial;f_m_ns_initial=M_ns_initial/M_tot_initial;f_m_stars_initial=M_stars_initial/M_tot_initial;
+
+toc
+%This way we have the initial values for everything of interest for us, so
+%we can study what are the effects of a variation in the escape velocity:
+
+%Now that we have defined the mass vector we make use of the probabilities:
+%But we need to consider that for each value of v_esc we have to save a
+%vector of different masses. For that we form a matrix, but it is important
+%to consider that we will not have the same number of stars for every
+%escape velocity, and matlab fills the ones that are "missing" with 0's. In
+%order, then, to not consider these zeros we have to take only till the
+%N_tot_final position for every escape velocity.
+tic
+mat_m=[];
+for ii=1:length(vec_v_esc)
+    v_esc=vec_v_esc(ii);
+    vec_m_final=[];
+    ll=1;%We initialize this and then change the value by one if the star is retained.
+    for jj=1:length(vec_m_initial)
+        m=vec_m_initial(jj);
+        if m<M_z_min_ns_test
+            vec_m_final(ll)=m;
+            ll=ll+1;%Because it is always retained.
+            
+        elseif m>M_z_min_bh_test
+            %k=1 ==>with fallback
+            %f_ret=(1-feval(@prob_escape,Z,m,1.4,265,v_esc,1));
+            %We now generate a random number, if it is below f_ret it is
+            %retained, if not it is kicked out of the cluster.
+            rr=rand;
+            if rr<=vec_f_ret(ii,find(vec_M_examples==m))
+                vec_m_final(ll)=m;
+                ll=ll+1;
+            else
+                ll=ll;
+            end
+        else
+            %Again we compute the probability of it being retained:
+            %f_ret=vec_f_ret_ns(ii);
+            %f_ret=(1-feval(@prob_esc_ns,265,v_esc));
+            rr=rand;
+            if rr<=vec_f_ret_ns(ii)
+                vec_m_final(ll)=m;
+                ll=ll+1;
+            else
+                ll=ll;
+            end
+        end
+%We save in each row of the matrix mat_m every mass retained for a given
+%escape velocity:
+    mat_m(ii,1:length(vec_m_final))=vec_m_final;
+    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);
+    
+    
+    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);
+    %Fractions compared to the current cluster:
+    f_N_bh_final(ii)=N_bh_final(ii)/N_tot_final(ii);f_N_ns_final(ii)=N_ns_final(ii)/N_tot_final(ii);f_N_stars_final(ii)=N_stars_final(ii)/N_tot_final(ii);
+    f_m_bh_final(ii)=M_bh_final(ii)/M_tot_final(ii);f_m_ns_final(ii)=M_ns_final(ii)/M_tot_final(ii);f_m_stars_final=M_stars_final(ii)/M_tot_final(ii);
+    %(BHs now)/(BH initially):
+    f_N_bh_relative(ii)=N_bh_final(ii)/N_bh_initial;
+    f_m_bh_relative(ii)=M_bh_final(ii)/M_bh_initial;
+    end
+end
+toc
+%And now with this we can obtain a table:
+T=table(vec_v_esc',N_bh_final',N_ns_final',f_m_bh_final');
+T.Properties.VariableNames = {'v_esc','N_BH','N_NS','f_m_BH'};
+disp(['Initially we had: Z=' num2str(Z) ' & N_tot=' num2str(N_tot_initial)])
+disp('And taking into account fallback:')
+disp(T)
+%}
+
+%% CONTOUR PLOT FOR THE MASS FRACTION OF BLACK HOLES AS A FUNCTION OF M_cl AND r_h:
+%{
+Z=0.0001;
+
+%We define the following vectors:
+%vec_M_cl=[5*10^3 10^4 2*10^4 3*10^4 5*10^4 10^5 2*10^5 5*10^5 10^6];
+vec_M_cl=linspace(10^4,10^6,100);
+%vec_M_cl=[10^3 2*10^3 3*10^3 4*10^3 5*10^3 10^4 2*10^4 5*10^4 10^5];
+
+%vec_M_cl=[10^4 2*10^4 3*10^4 5*10^4 7.5*10^4 10^5 5*10^5];
+%vec_r_h=[0.2 0.5 0.75 1 1.5 2 2.5 3 4]; 
+%vec_r_h=[0.25 0.5 0.75 1 1.25 1.5 1.75 2];
+%vec_v_esc=[10 15 20 25 30 35 40 45 50];
+vec_v_esc=linspace(5,60,100);
+%And with this we form the meshgrid:
+[M_cl,mat_v_esc]=meshgrid(vec_M_cl,vec_v_esc);
+%And we need to precompute the probabilities so
+%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);
+%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)
+[C,h]=contour(M_cl/(10^5),r_h,f_N_bh_relative,'ShowText','on');
+clabel(C,h)
+title('BH retention fraction (Z=0.0001)')
+xlabel('M_{cl} (10^5 M_{sun})')
+ylabel('r_h (in pc)')
+set(gca, 'YScale', 'log')
+%}
+%% CONTOUR PLOT FOR LOW MASS CLUSTERS
+%{
+% When working with low mass globular clusters we have to consider that one
+% sampling is not really representative of a cluster with that mass, as
+% one can get really different results depending on the sampling. This is
+% something that is not that important in more massive clusters with masses
+% of 10^5 or 10^6 solar masses.
+%That is why in this case we will do more than one sampling and then
+%compute and plot the mean values:
+
+%We define the following vectors:
+%vec_M_cl=[5*10^3 10^4 2*10^4 3*10^4 5*10^4 10^5 2*10^5 5*10^5 10^6];
+vec_M_cl=linspace(10^3,10^4,100);
+%vec_M_cl=[10^3 2*10^3 3*10^3 4*10^3 5*10^3 10^4 2*10^4 5*10^4 10^5];
+
+%vec_M_cl=[10^4 2*10^4 3*10^4 5*10^4 7.5*10^4 10^5 5*10^5];
+%vec_r_h=[0.2 0.5 0.75 1 1.5 2 2.5 3 4]; 
+%vec_r_h=[0.25 0.5 0.75 1 1.25 1.5 1.75 2];
+%vec_v_esc=[10 15 20 25 30 35 40 45 50];
+vec_v_esc=linspace(5,40,100);
+%And with this we form the meshgrid:
+[M_cl,mat_v_esc]=meshgrid(vec_M_cl,vec_v_esc);
+%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));
+
+%So that:
+n_times=20;
+ii=2;
+f_N_bh_relative=fun_for_contour_fractions_3(M_cl,2e-2,mat_v_esc);
+vec_N_bh_initial=[];
+while ii<n_times
+f_N_bh_relative=f_N_bh_relative + fun_for_contour_fractions_3(M_cl,2e-2,mat_v_esc);
+ii=ii+1;
+end 
+%And now we determine the mean value by dividing the sum we have done in
+%the while loop by the number of different samplings done for each mass.
+mean_f_bh_relative=f_N_bh_relative/(n_times);
+
+figure(1)
+[C,h]=contour(M_cl/(10^3),r_h,f_bh_relative,'ShowText','on');
+clabel(C,h)
+title('BH retention fraction (Z=2e-2) (w/ fallback)')
+xlabel('M_{cl} (10^3 M_{sun})')
+ylabel('r_h (in pc)')
+%set(gca, 'YScale', 'log')
+%}
+%% HISTOGRAM WITH THE MASSES RETAINED AND THE FINAL MASSES
+M_cl=5*10^4;r_h=2;Z=0.0001;
+v_esc=fun_v_esc(r_h,M_cl);
+
+[vec_m_bh_initial,vec_m_bh_final]=fun_for_mass_histograms(M_cl,Z,v_esc,1);
+%But we have to take into account that we obtain the vector of masses, but
+%as M_zams, we net to transform that into black hole masses:
+vec_m_bh_initial=fun_M_rem(Z,vec_m_bh_initial);
+vec_m_bh_final=fun_M_rem(Z,vec_m_bh_final);
+figure(1)
+h1=histogram(vec_m_bh_initial,45,'FaceColor','b');
+h1.BinWidth = 0.25;
+title('Number of BHs for [M_{cl}=5*10^4 M_{sun}, r_h=2 pc, Z=0.02]')
+xlabel('M_{BH} (M_{sun})')
+ylabel('Number of BHs')
+xlim([0 30])
+xticks([0 5 10 20 25 30])
+%set(gca,'YScale','log')
+hold on
+h2=histogram(vec_m_bh_final,45,'FaceColor','r');
+h2.BinWidth = 0.25;
+legend('IMF','after kicks')