Skip to content

Commit 0007c69

Browse files
author
Ilya Mandel
committed
First complete Matlab upload
1 parent bed99bc commit 0007c69

File tree

3 files changed

+123
-110
lines changed

3 files changed

+123
-110
lines changed

‎compas_matlab_utils/BSEDetailedOutput.m‎

Lines changed: 15 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,16 @@
1-
file='~/Work/COMPAS/src/COMPAS_Output/Detailed_Output/BSE_Detailed_Output_0.h5';
1+
function BSEDetailedOutput(file)
2+
% Carries out some basic analysis and makes plots for a single BSE run
3+
%
4+
% USAGE:
5+
% BSEDetailedOutput(file)
6+
%
7+
% INPUTS:
8+
% file: name of detailed output file in COMPAS h5 format
9+
%
10+
% examples:
11+
% BSEDetailedOutput('~/Work/COMPAS/src/COMPAS_Output/Detailed_Output/BSE_Detailed_Output_0.h5')
12+
%
213

3-
tic
414
time=h5read(file,'/Time');
515
MThistory=h5read(file,'/MT_History');
616
Z=h5read(file,'/Metallicity@ZAMS(1)');
@@ -29,8 +39,6 @@
2939
omegab1=h5read(file, '/Omega_Break(1)');
3040
omegab2=h5read(file, '/Omega_Break(2)');
3141
MTtimescale=h5read(file, '/MassTransferTimescale');
32-
inf=h5info(file);
33-
%inf.Datasets.Name
3442

3543

3644
disp('Time (Myr), Event, M1 (M_o), type1, M2 (M_o), type2, a (R_o), e');
@@ -65,8 +73,7 @@
6573
disp([num2str(time(i)+Tdelay), ' GW merger in ', num2str(Tdelay,'%.0f'), ' Myr ', num2str(mass1(i)), ' ', num2str(type1(i)),...
6674
' ', num2str(mass2(i)), ' ', num2str(type2(i))]);
6775
end;
68-
toc
6976

70-
figure(81), plot(time, mass1, 'LineWidth', 3), set(gca,'FontSize',20), xlabel('Time, Myr'), ylabel('Total mass 1, M_o')
71-
%figure(81), scatter(time, radius1, 30, type1, 'filled'), set(gca,'FontSize',20), xlabel('Time, Myr'), ylabel('Radius 1, R_o')
72-
figure(82), plot(time, radius2, 'b', time, RL2, 'r'); set(gca,'FontSize',20), xlabel('Time, Myr'), ylabel('Radius, R_o'), legend('R_2', 'RL_2')
77+
figure(1), plot(time, mass1, 'LineWidth', 3), set(gca,'FontSize',20), xlabel('Time, Myr'), ylabel('Total mass 1, M_o')
78+
%figure(2), scatter(time, radius1, 30, type1, 'filled'), set(gca,'FontSize',20), xlabel('Time, Myr'), ylabel('Radius 1, R_o')
79+
figure(2), plot(time, radius2, 'b', time, RL2, 'r'); set(gca,'FontSize',20), xlabel('Time, Myr'), ylabel('Radius, R_o'), legend('R_2', 'RL_2')

‎compas_matlab_utils/CosmicHistoryIntegrator.m‎

Lines changed: 71 additions & 69 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,4 @@
1-
function [Zlist, MergerRateByRedshiftByZ, SFR, Rdetections, DetectableMergerRate, Mtzlistdetection, etalist, zlistdetection]=...
1+
function [Zlist, MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, Mtlist, etalist, SFR, Rdetections, DetectableMergerRate, zlistdetection, x]=...
22
CosmicHistoryIntegrator(filename, zlistformation, zmaxdetection, Msimulated, makeplots)
33
% Integrator for the binary black hole merger rate over cosmic history
44
% COMPAS (Compact Object Mergers: Population Astrophysics and Statistics)
@@ -20,51 +20,32 @@
2020
%
2121
% OUTPUTS:
2222
% Zlist is a vector of metallicities, taken from the COMPAS run input file
23-
% MergerRateByRedshiftByZ is a matrix of size length(Zlist) X length(zformationlist)
23+
% MergerRateByRedshiftByZ is a matrix of size length(zformationlist) X length(Zlist)
2424
% which contains a merger rate of merging compact objects in the given redshift
2525
% and metallicity bin, in units of mergers per Mpc^3 of comoving volume per
2626
% year of source time
27+
% MergerRateByRedshiftByMtByEta is a matrix of size llength(zformationlist)
28+
% X length(Mtlist) X length(etalist) which contains a merger rate of merging compact objects
29+
% in the given redshift, total mass and eta bin, in units of mergers per Mpc^3
30+
% of comoving volume per year of source time
31+
% Mtlist is a list of total mass bins
32+
% etalist is a list of symmetric mass ratio bins
2733
% SFR is a vector of size length(zlistformation) containing the star formation rate
2834
% (solar masses per Mpc^3 of comoving volume per year of source time)
29-
% Rdetection is a matrix of size length(zlistdetection) X length Mtzlist X
35+
% Rdetection is a matrix of size length(zlistdetection) X length(Mtlist) X
3036
% length(etalist) containing the detection rate per year of observer time
31-
% from a given redshift bin and redshifted total mass and symmetric mass
32-
% ratio pixel
37+
% from a given redshift bin and total mass and symmetric mass ratio pixel
3338
% DetectableMergerRate is a matrix of the same size as Rdetection but
3439
% containing the intrinsic rate of detectable mergers per Mpc^3 of comoving
3540
% volume per year of source time
36-
% Mtzlist is a list of redshifted total mass bins (for computational
37-
% efficiency, very rare sources with Mtz>200 Msun are folded into the
38-
% last bin)
39-
% etalist is a list of symmetric mass ratio bins
4041
% zlistdetection is a vector of redshifts at which detection rates are
4142
% computed (a subset of zlistformation going up to zmaxdetection)
4243
%
4344
% EXAMPLE:
4445
% zlist=0:0.01:10;
45-
% [Zlist, MergerRateByRedshiftByZ, SFR, Rdetections, DetectableMergerRate,
46-
% Mtzlistdetection, etalist, zlistdetection] = ...
46+
% [Zlist, MergerRateByRedshiftByZ, MergerRateByRedshiftByMtByEta, Mtlist, etalist, SFR, Rdetections, DetectableMergerRate, zlistdetection,] = ...
4747
% CosmicHistoryIntegrator('~/Work/COMPASresults/runs/Zdistalpha1-031803.h5', zlist, 1.5, 90e6, 1);
4848
%
49-
%
50-
% zlist=0:0.1:5;
51-
% filename=['~/Work/COMPASresults/popsynth/runs/',...
52-
% '20170628-Coen-Pessimistic/AllMergers.dat'];
53-
% [Zlist,MergerRateByRedshiftByZ]=CosmicHistoryIntegrator(filename,zlist);
54-
% figure(1),colormap jet;
55-
% plot(zlist, MergerRateByRedshiftByZ, 'LineWidth', 2),
56-
% legend(num2str(Zlist)),
57-
% set(gca, 'FontSize', 20); %for labels
58-
% xlabel('z'),
59-
% ylabel('Pessimistic non-RLOF BBH merger rate, per Gpc^3 per yr')
60-
% sum(MergerRateByRedshiftByZ(:,1)) %Total BBH merger rate at z=0
61-
% [Zlist,MergerRateByRedshiftByZLanger]=CosmicHistoryIntegrator(filename,zlist,1);
62-
% [Zlist,MergerRateByRedshiftByZLambert]=CosmicHistoryIntegrator(filename,zlist,2);
63-
% plot(zlist, sum(MergerRateByRedshiftByZLambert,2), 'r', ...
64-
% zlist, sum(MergerRateByRedshiftByZLanger,2), 'b', 'LineWidth', 2),
65-
% set(gca, 'FontSize', 20), xlabel('z'), legend('Langer & Norman','Lambert')
66-
% ylabel('Pessimistic non-RLOF BBH merger rate, per Gpc^3 per yr'),
67-
6849

6950

7051
%define constants
@@ -88,44 +69,58 @@
8869
%metallicity-specific SFR
8970
[SFR,Zlist,Zweight]=Metallicity(zlistformation,min(Z),max(Z));
9071

91-
9272
%Consider the contribution of every simulated binary to the merger rate
9373
%in every redshift bin by considering when it would have to be formed to
9474
%merge at that redshift and normalizing by the relevant
9575
%metallicity-specific star formation rate
9676
dz=zlistformation(2)-zlistformation(1);
97-
tLmerge=tL(floor(zlistformation/dz)+1);
9877
etalist=0.01:0.01:0.25;
99-
Mtzlist=1:1:ceil(max(M1+M2)*(1+max(zlistformation)));
78+
Mtlist=1:1:ceil(max(M1+M2));
10079
MergerRateByRedshiftByZ=zeros(length(zlistformation),length(Zlist));
101-
MergerRateByRedshiftByMtzByEta=zeros(length(zlistformation),length(Mtzlist),length(etalist));
80+
MergerRateByRedshiftByMtByEta=zeros(length(zlistformation),length(Mtlist),length(etalist));
81+
x=zeros(size(M1));
10282
for(i=1:length(M1)),
10383
Zcounter=find(Zlist>=Z(i),1);
10484
eta=M1(i)*M2(i)/(M1(i)+M2(i))^2;
10585
etaindex=ceil(eta*100);
106-
for(k=1:length(zlistformation)), %merger redshift index
107-
if((Tdelay(i)+tLmerge(k)) > max(tL)), continue; end; %binary can't merge this early
108-
zformindex=find(tL>=(Tdelay(i)+tLmerge(k)),1);
109-
Mtzindex=ceil((M1(i)+M2(i))*(1+zlistformation(k)));
110-
MergerRateByRedshiftByZ(k,Zcounter)=...
111-
MergerRateByRedshiftByZ(k,Zcounter)+...
112-
SFR(zformindex)*Zweight(zformindex,Zcounter)/Msimulated;
113-
MergerRateByRedshiftByMtzByEta(k,Mtzindex,etaindex) =...
114-
MergerRateByRedshiftByMtzByEta(k,Mtzindex,etaindex) + ...
115-
SFR(zformindex)*Zweight(zformindex,Zcounter)/Msimulated;
86+
tLform=tL+Tdelay(i); %lookback time of when binary would have to form in order to merge at lookback time tL
87+
firsttooearlyindex=find((tLform)>max(tL),1);
88+
if(isempty(firsttooearlyindex)), firsttooearlyindex=length(tL)+1; end;
89+
zForm=interp1(tL,zlistformation,tLform(1:firsttooearlyindex-1));
90+
zFormindex=ceil((zForm-zlistformation(1))./dz)+1;
91+
Mtindex=ceil(M1(i)+M2(i));
92+
if(~isempty(zFormindex))
93+
x(i)=SFR(zFormindex(1))*Zweight(zFormindex(1),Zcounter)/Msimulated;
94+
MergerRateByRedshiftByZ(1:firsttooearlyindex-1,Zcounter)=...
95+
MergerRateByRedshiftByZ(1:firsttooearlyindex-1,Zcounter)+...
96+
transpose(SFR(zFormindex)).*Zweight(zFormindex,Zcounter)/Msimulated;
97+
MergerRateByRedshiftByMtByEta(1:firsttooearlyindex-1,Mtindex,etaindex) =...
98+
MergerRateByRedshiftByMtByEta(1:firsttooearlyindex-1,Mtindex,etaindex) + ...
99+
transpose(SFR(zFormindex)).*Zweight(zFormindex,Zcounter)/Msimulated;
116100
end;
101+
%for(k=1:length(zlistformation)), %merger redshift index
102+
% if((Tdelay(i)+tL(k)) > max(tL)), continue; end; %binary can't merge this early
103+
% zformindex=find(tL>=(Tdelay(i)+tL(k)),1);
104+
% Mtzindex=ceil((M1(i)+M2(i))*(1+zlistformation(k)));
105+
% MergerRateByRedshiftByZ(k,Zcounter)=...
106+
% MergerRateByRedshiftByZ(k,Zcounter)+...
107+
% SFR(zformindex)*Zweight(zformindex,Zcounter)/Msimulated;
108+
% MergerRateByRedshiftByMtzByEta(k,Mtzindex,etaindex) =...
109+
% MergerRateByRedshiftByMtzByEta(k,Mtzindex,etaindex) + ...
110+
% SFR(zformindex)*Zweight(zformindex,Zcounter)/Msimulated;
111+
%end;
117112
end;
118113

119114
zlistdetection=zlistformation(1:find(zlistformation<=zmaxdetection,1,"last"));
120115
fin=load('~/Work/Rai/LIGOfuture_data/freqVector.txt');
121116
%noise=load('~/Work/Rai/LIGOfuture_data/dataNomaLIGO.txt');
122117
noise=load('~/Work/Rai/LIGOfuture_data/dataEarly_low.txt');
123-
[Mtzlistdetection, Rdetections,DetectableMergerRate]=...
124-
DetectionRate(zlistformation,Mtzlist,etalist,MergerRateByRedshiftByMtzByEta,zlistdetection,fin,noise,Dl,dVc);
118+
[Rdetections,DetectableMergerRate]=...
119+
DetectionRate(zlistformation,Mtlist,etalist,MergerRateByRedshiftByMtByEta,zlistdetection,fin,noise,Dl,dVc);
125120

126121
if(makeplots==1), %make a set of default plots
127122
MakePlots(M1,M2,Z,Tdelay,zlistformation,Zlist,SFR,Zweight,...
128-
MergerRateByRedshiftByZ, Rdetections, DetectableMergerRate, zlistdetection, Mtzlistdetection, etalist, 1);
123+
MergerRateByRedshiftByZ, Rdetections, DetectableMergerRate, zlistdetection, Mtlist, etalist, 1);
129124
end;
130125

131126
end %end of CosmicHistoryIntegrator
@@ -208,16 +203,17 @@
208203
Zmu=log(Zmean)-0.39^2/2;
209204
dlogZ=0.1;
210205
logZvec=-12:dlogZ:0; %natural log
211-
dPdlogZ=1/0.39/sqrt(2*pi)*exp(-(logZvec'-Zmu).^2/2/0.39^2);
206+
dPdlogZ=1/0.39/sqrt(2*pi)*exp(-(logZvec'-Zmu).^2/2/0.39^2); %size length(logZvec) x length(zvec)
212207
dPdlogZ=dPdlogZ./(sum(dPdlogZ,1)*dlogZ); %normalise
213208
minlogZindex=find(exp(logZvec)>=minZ,1, 'first');
214209
maxlogZindex=find(exp(logZvec)>=maxZ,1, 'first');
215210
Zrange=logZvec(maxlogZindex)-logZvec(minlogZindex); %ugly correction for not including tails
216211
PdrawZ=1/Zrange;
217212
Zvec=exp(logZvec(minlogZindex:maxlogZindex));
218-
dPdlogZ(minlogZindex,:)=dPdlogZ(minlogZindex,:)+sum(dPdlogZ(1:minlogZindex,:),1)*dlogZ/(sum(Zvec==min(Zvec))/length(Zvec))*PdrawZ;
219-
dPdlogZ(maxlogZindex,:)=dPdlogZ(maxlogZindex,:)+sum(dPdlogZ(maxlogZindex:end,:),1)*dlogZ/(sum(Zvec==max(Zvec))/length(Zvec))*PdrawZ;
220-
dPdlogZ(1:minlogZindex,:)=0; dPdlogZ(maxlogZindex:size(dPdlogZ,1),:)=0;
213+
Zweight=zeros(length(zvec),length(Zvec));
214+
dPdlogZ(minlogZindex,:)=dPdlogZ(minlogZindex,:)+sum(dPdlogZ(1:minlogZindex-1,:),1)*dlogZ/dlogZ;
215+
dPdlogZ(maxlogZindex,:)=dPdlogZ(maxlogZindex,:)+sum(dPdlogZ(maxlogZindex+1:end,:),1)*dlogZ/dlogZ;
216+
dPdlogZ(1:minlogZindex-1,:)=0; dPdlogZ(maxlogZindex+1:size(dPdlogZ,1),:)=0;
221217
dPdlogZ=dPdlogZ./(sum(dPdlogZ,1)*dlogZ); %normalise
222218
for(i=1:length(Zvec))
223219
index=find(exp(logZvec)>=Zvec(i), 1, 'first');
@@ -231,12 +227,12 @@
231227

232228

233229
%Compute detection rates per year of observer time and per year of source time
234-
%per Mpc^3 of comoving volume as a function of redshifted total mass and eta
235-
function [Mtzlistdetection, Rdetections, DetectableMergerRate]=...
236-
DetectionRate(zlistformation,Mtzlist,etalist,MergerRateByRedshiftByMtzByEta,zlistdetection,freqfile,noisefile,Dl,dVc)
230+
%per Mpc^3 of comoving volume as a function of total mass and eta
231+
function [Rdetections, DetectableMergerRate]=...
232+
DetectionRate(zlistformation,Mtlist,etalist,MergerRateByRedshiftByMtByEta,zlistdetection,freqfile,noisefile,Dl,dVc)
237233

238234
fin=load('~/Work/Rai/LIGOfuture_data/freqVector.txt');
239-
noise=load('~/Work/Rai/LIGOfuture_data/dataEarly_low.txt');
235+
noise=load('~/Work/Rai/LIGOfuture_data/dataMid_low.txt');
240236

241237
flow=10;
242238
df=1;
@@ -258,7 +254,7 @@
258254

259255
%save time by not doing calculations beyond maximum redshifted total
260256
%mass corresponding to detection redshift threshold
261-
Mtzlistdetection=Mtzlist(1:ceil(length(Mtzlist)*max(zlistdetection)/max(zlistformation)));
257+
Mtzlistdetection=1:1:ceil(max(Mtlist)*(1+max(zlistdetection)));
262258
SNRat1Mpc=zeros(length(Mtzlistdetection),length(etalist));
263259
for(i=1:length(Mtzlistdetection)),
264260
for(j=1:length(etalist)),
@@ -267,29 +263,35 @@
267263
SNRat1Mpc(i,j)=sqrt(integral);
268264
end;
269265
end;
270-
SNRat1Mpc(1,25)
271266

272-
SNR=zeros(length(zlistdetection),length(Mtzlistdetection),length(etalist));
273-
for(i=1:length(zlistdetection)), SNR(i,:,:)=SNRat1Mpc./Dl(i); end;
267+
SNR=zeros(length(zlistdetection),length(Mtlist),length(etalist));
268+
269+
for(i=1:length(zlistdetection)),
270+
for(j=1:length(Mtlist)),
271+
SNR(i,j,:)=SNRat1Mpc(ceil(j*(zlistdetection(i)+1)),:)./Dl(i);
272+
end;
273+
end;
274+
%for(i=1:length(zlistdetection)), SNR(i,:,:)=SNRat1Mpc./Dl(i); end;
274275

275276
SNR8pre=1:0.1:1000;
276277
theta=1./SNR8pre;
277278
pdetect=1-interp1([0,Thetas,1],[(0:Ntheta)/Ntheta,1],theta);
278279
pdetect(1)=0; %set of measure zero to exceed threshold, but enforce just in case
279280

280-
Rdetections=zeros(length(zlistdetection),length(Mtzlistdetection),length(etalist)); %Detections per unit observer time
281-
DetectableMergerRate=zeros(length(zlistdetection),length(Mtzlistdetection),length(etalist)); %Detections per unit source time per unit Vc
281+
Rdetections=zeros(length(zlistdetection),length(Mtlist),length(etalist)); %Detections per unit observer time
282+
DetectableMergerRate=zeros(length(zlistdetection),length(Mtlist),length(etalist)); %Detections per unit source time per unit Vc
282283
SNR8=SNR/8;
283284
pdetection=zeros(size(Rdetections));
284285
pdetection=pdetect(max(min(floor(SNR8*10),length(pdetect)),1));
285-
DetectableMergerRate=MergerRateByRedshiftByMtzByEta(1:length(zlistdetection),1:length(Mtzlistdetection),:).*pdetection;
286+
287+
DetectableMergerRate=MergerRateByRedshiftByMtByEta(1:length(zlistdetection),:,:).*pdetection;
286288
Rdetections=DetectableMergerRate.*transpose(dVc(1:length(zlistdetection)))./(1+zlistdetection');
287289

288290
end %end of DetectionRate
289291

290292
%Make a set of default plots
291293
function MakePlots(M1,M2,Z,Tdelay,zlistformation,Zlist,SFR,Zweight,...
292-
MergerRateByRedshiftByZ, Rdetections, DetectableMergerRate, zlistdetection, Mtzlistdetection, etazlist, fignumber)
294+
MergerRateByRedshiftByZ, Rdetections, DetectableMergerRate, zlistdetection, Mtlist, etazlist, fignumber)
293295

294296
zvec=zlistformation;
295297

@@ -332,12 +334,12 @@ function MakePlots(M1,M2,Z,Tdelay,zlistformation,Zlist,SFR,Zweight,...
332334

333335

334336
figure(fignumber+4), clf(fignumber+4);
335-
RdetectionsByzMtz=sum(Rdetections,3); %sum across eta
336-
plot(zlistdetection, cumsum(sum(RdetectionsByzMtz,2)), 'LineWidth', 3), hold on;
337-
plot(zlistdetection, cumsum(sum(RdetectionsByzMtz(:,Mtzlistdetection<=5),2)), 'LineWidth', 1);
338-
plot(zlistdetection, cumsum(sum(RdetectionsByzMtz(:,Mtzlistdetection>5 & Mtzlistdetection<20),2)), 'LineWidth', 1);
339-
plot(zlistdetection, cumsum(sum(RdetectionsByzMtz(:,Mtzlistdetection>=20),2)), 'LineWidth', 1); hold off;
340-
legend('Total rate', 'From M_t(1+z)<=5 M_o', 'From 5<M_t(1+z)/M_o<20', 'From M_t(1+z)>=20 M_o'),
337+
RdetectionsByzMt=sum(Rdetections,3); %sum across eta
338+
semilogy(zlistdetection, cumsum(sum(RdetectionsByzMt,2)), 'LineWidth', 3), hold on;
339+
semilogy(zlistdetection, cumsum(sum(RdetectionsByzMt(:,Mtlist<=5),2)), 'LineWidth', 1);
340+
semilogy(zlistdetection, cumsum(sum(RdetectionsByzMt(:,Mtlist>5 & Mtlist<20),2)), 'LineWidth', 1);
341+
semilogy(zlistdetection, cumsum(sum(RdetectionsByzMt(:,Mtlist>=20),2)), 'LineWidth', 1); hold off;
342+
legend('Total rate', 'From M_t<=5 M_o', 'From 5<M_t/M_o<20', 'From M_t>=20 M_o'),
341343
set(gca, 'FontSize', 20); %for labels
342344
xlabel('z'),
343345
ylabel('cumulative detection rate per observer yr')

0 commit comments

Comments
 (0)