Skip to content

Commit bed99bc

Browse files
author
Ilya Mandel
committed
Done with cosmic integration
1 parent 4e2d2ee commit bed99bc

File tree

2 files changed

+79
-69
lines changed

2 files changed

+79
-69
lines changed

‎compas_matlab_utils/CosmicHistoryIntegrator.m‎

Lines changed: 62 additions & 48 deletions
Original file line numberDiff line numberDiff line change
@@ -1,12 +1,12 @@
1-
function [Zlist, MergerRateByRedshiftByZ, SFR, Rdetections, DetectableMergerRate, Mtzlist, etalist, zlistdetection]=...
1+
function [Zlist, MergerRateByRedshiftByZ, SFR, Rdetections, DetectableMergerRate, Mtzlistdetection, etalist, zlistdetection]=...
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)
55
% software package
66
%
77
% USAGE:
88
% [Zlist, MergerRateByRedshiftByZ, SFR, Zweight, Rdetections, DetectableMergerRate, Mtzlist, etalist, zlistdetection]=...
9-
% CosmicHistoryIntegrator(filename, zlistformation, zmaxdetection, Msimulated, [,makeplots])
9+
% CosmicHistoryIntegrator(filename, zlistformation, zmaxdetection, Msimulated [,makeplots])
1010
%
1111
% INPUTS:
1212
% filename: name of population synthesis input file
@@ -41,8 +41,12 @@
4141
% computed (a subset of zlistformation going up to zmaxdetection)
4242
%
4343
% EXAMPLE:
44-
% CosmicHistoryIntegrator('~/Work/COMPASresults/runs/Zdistalpha1-031803.h5', ...
45-
% zlist, zlist, 90e6, 1, '~/Work/COMPASresults/runs/Zdist2stage-031803.h5', 'Default', '2 Stage')
44+
% zlist=0:0.01:10;
45+
% [Zlist, MergerRateByRedshiftByZ, SFR, Rdetections, DetectableMergerRate,
46+
% Mtzlistdetection, etalist, zlistdetection] = ...
47+
% CosmicHistoryIntegrator('~/Work/COMPASresults/runs/Zdistalpha1-031803.h5', zlist, 1.5, 90e6, 1);
48+
%
49+
%
4650
% zlist=0:0.1:5;
4751
% filename=['~/Work/COMPASresults/popsynth/runs/',...
4852
% '20170628-Coen-Pessimistic/AllMergers.dat'];
@@ -91,18 +95,18 @@
9195
%metallicity-specific star formation rate
9296
dz=zlistformation(2)-zlistformation(1);
9397
tLmerge=tL(floor(zlistformation/dz)+1);
94-
etavec=0.01:0.01:0.25;
95-
Mtzvec=1:1:200; %ignore things on wrong side of mass gap, go up to z=1
98+
etalist=0.01:0.01:0.25;
99+
Mtzlist=1:1:ceil(max(M1+M2)*(1+max(zlistformation)));
96100
MergerRateByRedshiftByZ=zeros(length(zlistformation),length(Zlist));
97-
MergerRateByRedshiftByMtzByEta=zeros(length(zlistformation),length(etavec),length(Mtzvec));
101+
MergerRateByRedshiftByMtzByEta=zeros(length(zlistformation),length(Mtzlist),length(etalist));
98102
for(i=1:length(M1)),
99-
i, M1(i), M2(i)
100103
Zcounter=find(Zlist>=Z(i),1);
101-
zformindex=find(tL>=Tdelay(i),1);
102-
etaindex=M1(i)^0.6*M2(i)^0.6/(M1(i)+M2(i))^0.2;
103-
for(k=1:length(zlistformation)),
104-
zformindex=find(tL>=(Tdelay(i)+tLmerge(k)),1)
105-
Mtzindex=ceil((M1(i)+M2(i))*zlistformation(k))
104+
eta=M1(i)*M2(i)/(M1(i)+M2(i))^2;
105+
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)));
106110
MergerRateByRedshiftByZ(k,Zcounter)=...
107111
MergerRateByRedshiftByZ(k,Zcounter)+...
108112
SFR(zformindex)*Zweight(zformindex,Zcounter)/Msimulated;
@@ -116,12 +120,12 @@
116120
fin=load('~/Work/Rai/LIGOfuture_data/freqVector.txt');
117121
%noise=load('~/Work/Rai/LIGOfuture_data/dataNomaLIGO.txt');
118122
noise=load('~/Work/Rai/LIGOfuture_data/dataEarly_low.txt');
119-
[Rdetections,DetectableMergerRate]=...
120-
DetectionRate(zlistformation,Mtzlist,etalist,MergerRateByRedshiftByMtzByEta,zlistdetection,fin,noise,Dl,dVc)
123+
[Mtzlistdetection, Rdetections,DetectableMergerRate]=...
124+
DetectionRate(zlistformation,Mtzlist,etalist,MergerRateByRedshiftByMtzByEta,zlistdetection,fin,noise,Dl,dVc);
121125

122126
if(makeplots==1), %make a set of default plots
123127
MakePlots(M1,M2,Z,Tdelay,zlistformation,Zlist,SFR,Zweight,...
124-
MergerRateByRedshiftByZ, Rdetections, DetectableMergerRate, Mtzlist, etalist, 1);
128+
MergerRateByRedshiftByZ, Rdetections, DetectableMergerRate, zlistdetection, Mtzlistdetection, etalist, 1);
125129
end;
126130

127131
end %end of CosmicHistoryIntegrator
@@ -217,7 +221,7 @@
217221
dPdlogZ=dPdlogZ./(sum(dPdlogZ,1)*dlogZ); %normalise
218222
for(i=1:length(Zvec))
219223
index=find(exp(logZvec)>=Zvec(i), 1, 'first');
220-
Zweight(:,i)=dPdlogZ(index,:)*dlogZ;
224+
Zweight(:,i)=dPdlogZ(index,:)./PdrawZ; %weight is desired over sampled probabiluty
221225
end;
222226
else %relevant for single-metallicity runs -- just give all binaries the same unit weight
223227
Zvec=minZ;
@@ -226,10 +230,11 @@
226230
end %end of Metallicity
227231

228232

229-
%Compute detection rates per unit observer time and per unit source time
230-
%per unit comoving volume as a function of redshifted total mass and eta
231-
function [Rdetections,DetectableMergerRate]=...
233+
%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]=...
232236
DetectionRate(zlistformation,Mtzlist,etalist,MergerRateByRedshiftByMtzByEta,zlistdetection,freqfile,noisefile,Dl,dVc)
237+
233238
fin=load('~/Work/Rai/LIGOfuture_data/freqVector.txt');
234239
noise=load('~/Work/Rai/LIGOfuture_data/dataEarly_low.txt');
235240

@@ -251,50 +256,52 @@
251256
Thetas=sort(Theta);
252257

253258

254-
SNRat1Mpc=zeros(length(Mtzvec),length(etavec));
255-
for(i=1:length(Mtzvec)),
256-
for(j=1:length(etavec)),
257-
[h,Am,psi]=IMRSAWaveform(f, Mtzvec(i), etavec(j), 0, 0, 0, 1, flow);
259+
%save time by not doing calculations beyond maximum redshifted total
260+
%mass corresponding to detection redshift threshold
261+
Mtzlistdetection=Mtzlist(1:ceil(length(Mtzlist)*max(zlistdetection)/max(zlistformation)));
262+
SNRat1Mpc=zeros(length(Mtzlistdetection),length(etalist));
263+
for(i=1:length(Mtzlistdetection)),
264+
for(j=1:length(etalist)),
265+
[h,Am,psi]=IMRSAWaveform(f, Mtzlistdetection(i), etalist(j), 0, 0, 0, 1, flow);
258266
integral=sum(4*Am.^2./Sf*df);
259267
SNRat1Mpc(i,j)=sqrt(integral);
260268
end;
261269
end;
270+
SNRat1Mpc(1,25)
262271

263-
SNR=zeros(length(zlistdetection),length(Mtzlist),length(etalist));
264-
for(i=1:length(zlistdetection)), SNR(i,:,:)=1./(Dl(i)./Mpc); end;
265-
SNR=SNR.*SNRat1Mpc;
266-
272+
SNR=zeros(length(zlistdetection),length(Mtzlistdetection),length(etalist));
273+
for(i=1:length(zlistdetection)), SNR(i,:,:)=SNRat1Mpc./Dl(i); end;
267274

268-
SNR8pre=0.1:0.1:1000;
275+
SNR8pre=1:0.1:1000;
269276
theta=1./SNR8pre;
270277
pdetect=1-interp1([0,Thetas,1],[(0:Ntheta)/Ntheta,1],theta);
278+
pdetect(1)=0; %set of measure zero to exceed threshold, but enforce just in case
271279

272-
Rdetections=zeros(length(zlistdetection),length(Mtzlist),length(etalist)); %Detections per unit observer time
273-
DetectableMergerRate=zeros(length(zlistdetection),length(Mtzlist),length(etalist)); %Detections per unit source time per unit Vc
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
274282
SNR8=SNR/8;
275283
pdetection=zeros(size(Rdetections));
276-
pdetection(SNR8>1)=pdetect(floor(SNR8(SNR8>1 & SNR8<max(SNR8pre))*10));
277-
pdetection(SNR8>max(SNR8pre))=1;
278-
DetectableMergerRate=MergerRateByRedshiftByMtzByEta(i,1:length(zlistdetection)).*pdetection;
279-
Rdetections(i,:)=MergerRateByRedshiftByMtzByEta(i,1:length(zlistdetection)).*pdetection.*dVc(1:maxzdetindex)./(1+zlistdetection);
284+
pdetection=pdetect(max(min(floor(SNR8*10),length(pdetect)),1));
285+
DetectableMergerRate=MergerRateByRedshiftByMtzByEta(1:length(zlistdetection),1:length(Mtzlistdetection),:).*pdetection;
286+
Rdetections=DetectableMergerRate.*transpose(dVc(1:length(zlistdetection)))./(1+zlistdetection');
280287

281288
end %end of DetectionRate
282289

283290
%Make a set of default plots
284-
function MakePlots(M1,M2,Z,Tdelay,zformationlist,Zlist,SFR,Zweight,...
285-
MergerRateByRedshiftByZ, Rdetections, DetectableMergerRate, zdetectionlist, Mtzlist, etazlist, fignumber)
291+
function MakePlots(M1,M2,Z,Tdelay,zlistformation,Zlist,SFR,Zweight,...
292+
MergerRateByRedshiftByZ, Rdetections, DetectableMergerRate, zlistdetection, Mtzlistdetection, etazlist, fignumber)
286293

287-
zvec=zformationlist;
294+
zvec=zlistformation;
288295

289296
figure(fignumber), clf(fignumber); %,colormap jet;
290297
plot(zvec, sum(MergerRateByRedshiftByZ,2)*1e9, 'LineWidth', 3), hold on;
291298
plot(zvec, sum(MergerRateByRedshiftByZ(:,Zlist<=0.001),2)*1e9, 'LineWidth', 1);
292299
plot(zvec, sum(MergerRateByRedshiftByZ(:,Zlist>0.001 & Zlist<0.01),2)*1e9, 'LineWidth', 1);
293300
plot(zvec, sum(MergerRateByRedshiftByZ(:,Zlist>=0.01),2)*1e9, 'LineWidth', 1); hold off;
294-
legend('Total rate', 'From Z<=0.001', 'From 0.001<Z<0.01', 'From Z>=0.01'),
295301
set(gca, 'FontSize', 20); %for labels
296302
xlabel('z'),
297303
ylabel('DCO merger rate per Gpc^3 per yr')
304+
legend('Total rate', 'From Z<=0.001', 'From 0.001<Z<0.01', 'From Z>=0.01'),
298305
disp(['Total DCO merger rate at z=0: ', ...
299306
num2str(1e9*sum(MergerRateByRedshiftByZ(1,:))),...
300307
' per Gpc^3 per year']);
@@ -316,18 +323,25 @@ function MakePlots(M1,M2,Z,Tdelay,zformationlist,Zlist,SFR,Zweight,...
316323

317324
figure(fignumber+3), clf(fignumber+3);
318325
plot(zvec, SFR*1e9, 'LineWidth', 3); hold on;
319-
plot(zvec, SFR'.*sum(Zweight(:,Zlist<=0.001),2)*1e9, 'LineWidth', 1);
320-
plot(zvec, SFR'.*sum(Zweight(:,Zlist>0.001&Zlist<0.01),2)*1e9, 'LineWidth', 1);
321-
plot(zvec, SFR'.*sum(Zweight(:,Zlist>=0.01),2)*1e9, 'LineWidth', 1); hold off;
322-
legend('Total rate', 'From Z<=0.001', 'From 0.001<Z<0.01', 'From Z>=0.01'),
326+
plot(zvec, SFR'.*sum(Zweight(:,Zlist<=0.001),2)./sum(Zweight,2)*1e9, 'LineWidth', 1);
327+
plot(zvec, SFR'.*sum(Zweight(:,Zlist>0.001&Zlist<0.01),2)./sum(Zweight,2)*1e9, 'LineWidth', 1);
328+
plot(zvec, SFR'.*sum(Zweight(:,Zlist>=0.01),2)./sum(Zweight,2)*1e9, 'LineWidth', 1); hold off;
323329
set(gca, 'FontSize', 20); %for labels
324330
xlabel('z'), ylabel('Star-formation rate, M_o per Gpc^3 per yr');
325-
326-
figure(5), colormap jet;
327-
plot(zvec, Zweight, 'LineWidth', 3)
331+
legend('Total rate', 'From Z<=0.001', 'From 0.001<Z<0.01', 'From Z>=0.01'),
332+
333+
334+
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'),
328341
set(gca, 'FontSize', 20); %for labels
329-
legend(num2str(Zlist))
330-
xlabel('z'), ylabel('Z-specific SFR weight');
342+
xlabel('z'),
343+
ylabel('cumulative detection rate per observer yr')
344+
331345

332346
end %end of MakePlots
333347

‎compas_matlab_utils/SSEDetailedOutput.m‎

Lines changed: 17 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -1,16 +1,22 @@
1-
function SSEDetailedOutput(filename)
1+
function SSEDetailedOutput(filename, nfiles)
22
% Carries out some basic analysis and makes plots for a single SSE run
33
%
44
% USAGE:
5-
% ComparisonPlots(filename1, name1, filename2, name2)
5+
% SSEDetailedOutput(filename [, nfiles])
66
%
77
% INPUTS:
88
% filename: name of detailed output file in COMPAS h5 format
9+
% nfiles: optional if multiple files are being compared; in that case,
10+
% the files are all assumed to start with filename, followed by numbers
11+
% 0...nfiles-1, followed by .h5
912
%
10-
% example:
13+
% examples:
1114
% SSEDetailedOutput('~/Work/COMPAS/src/COMPAS_Output/Detailed_Output/SSE_Detailed_Output_0.h5')
12-
%
15+
% SSEDetailedOutput('~/Work/COMPAS/src/COMPAS_Output/Detailed_Output/SSE_Detailed_Output_', 5)
16+
%
17+
1318

19+
if(nargin==1)
1420

1521
time=h5read(file,'/Time');
1622
Z=h5read(file,'/Metallicity@ZAMS');
@@ -29,6 +35,8 @@ function SSEDetailedOutput(filename)
2935
binding=h5read(file,'/BE_ConvectiveEnvelope');
3036
rec=h5read(file,'/Record_Type');
3137

38+
39+
3240
figure(1),
3341
scatter(time(type==1),radius(type==1),20,'filled', 'b'); hold on;
3442
scatter(time(type==2),radius(type==2),20,'filled', 'g');
@@ -60,22 +68,6 @@ function SSEDetailedOutput(filename)
6068
%set(gca,'FontSize',20); xlabel('Radius, Rsun'); ylabel('\lambda')
6169

6270

63-
64-
file='~/Work/COMPAS/src/COMPAS_Output/COMPAS_Output.h5';
65-
radius=h5read(file,'/BSE_Common_Envelopes/Radius(1)<CE');
66-
apre=h5read(file,'/BSE_Common_Envelopes/SemiMajorAxis<CE');
67-
a1=h5read(file,'/BSE_Common_Envelopes/SemiMajorAxisStage1>CE');
68-
apost=h5read(file,'/BSE_Common_Envelopes/SemiMajorAxis>CE');
69-
primarydonor=(h5read(file,'/BSE_Common_Envelopes/RLOF(1)')==1);
70-
figure(3); semilogy(radius(primarydonor), apre(primarydonor), 'b-.', radius(primarydonor), a1(primarydonor), 'b--', radius(primarydonor), apost(primarydonor), 'b-', 'LineWidth', 3); hold on;
71-
file='~/Work/COMPAS/src/COMPAS_Output_1/COMPAS_Output.h5';
72-
radius=h5read(file,'/BSE_Common_Envelopes/Radius(1)<CE');
73-
apre=h5read(file,'/BSE_Common_Envelopes/SemiMajorAxis<CE');
74-
a1=h5read(file,'/BSE_Common_Envelopes/SemiMajorAxisStage1>CE');
75-
apost=h5read(file,'/BSE_Common_Envelopes/SemiMajorAxis>CE');
76-
primarydonor=(h5read(file,'/BSE_Common_Envelopes/RLOF(1)')==1);
77-
figure(3); semilogy(radius(primarydonor), apre(primarydonor), 'm-.', radius(primarydonor), a1(primarydonor), 'm--', radius(primarydonor), apost(primarydonor), 'm-', 'LineWidth', 3); hold on;
78-
file='~/Work/COMPAS/src/COMPAS_Output_2/COMPAS_Output.h5';
7971
radius=h5read(file,'/BSE_Common_Envelopes/Radius(1)<CE');
8072
apre=h5read(file,'/BSE_Common_Envelopes/SemiMajorAxis<CE');
8173
a1=h5read(file,'/BSE_Common_Envelopes/SemiMajorAxisStage1>CE');
@@ -84,9 +76,11 @@ function SSEDetailedOutput(filename)
8476
figure(3); semilogy(radius(primarydonor), apre(primarydonor), 'y-.', radius(primarydonor), a1(primarydonor), 'y--', radius(primarydonor), apost(primarydonor), 'y-', 'LineWidth', 3); hold off;
8577
set(gca,'FontSize',20); xlabel('Radius, Rsun'); ylabel('Semimajor axis, Rsun'), legend('Pre CE','After Stage1','After Stage 2')
8678

79+
end; % end of if(nargin==1)
8780

81+
if(nargin==2),
8882

89-
for k=0:4,
83+
for (k=0:nfiles-1),
9084
file=['~/Work/COMPAS/src/COMPAS_Output_3/Detailed_Output/SSE_Detailed_Output_',int2str(k),'.h5'];
9185
time=h5read(file,'/Time');
9286
Z=h5read(file,'/Metallicity@ZAMS');
@@ -113,3 +107,5 @@ function SSEDetailedOutput(filename)
113107
hold off;
114108
set(gca,'FontSize',20); xlabel('Radius, Rsun'); ylabel('log_{10}(Nanjing binding energy / erg)'), legend('MS','HG','FGB','CHeB','EAGB','TPAGB'), title(['ZAMS Mass in solar masses:', int2str(mass(1))])
115109
end;
110+
111+
end; %end of if(nargin==2)

0 commit comments

Comments
 (0)