@@ -28,63 +28,73 @@ function ComparisonPlots(filename1, name1, filename2, name2)
2828 AU= 149597871e3 ; % m
2929 Rsun = 695500000 ; % m
3030
31- % Plot DCO mass distribution and BNS P-e distribution
32- figure(1 ); clf(1 ); figure(2 ); clf(2 );
33- [BNS ,NSBH ,BBH ,BNSCE ,NSBHCE ,BBHCE ]=DCOplot(filename1 , name1 , 1 , 2 , ' r' , 40 );
34- fprintf(' \n DCOs:\t\t #Merging DNS\t #Merging NSBH\t #Merging BBH\t%% BNS via CE\t%% NSBH via CE\t%% BBH via CE\n ' );
35- fprintf(' %s :\t%d\t\t%d\t\t%d\t\t%.0f\t\t%.0f\t\t%.0f\n ' , ...
36- name1 , BNS , NSBH , BBH , BNSCE / BNS * 100 , NSBHCE / NSBH * 100 , BBHCE / BBH * 100 );
31+ % Plot DCO mass distribution, BNS P-e distribution, chirp mass vs period
32+ % at DCO formation, BH mass vs secondary core mass for 2->1 CEs leading
33+ % to merging BBH formation
34+ figure(1 ); clf(1 ); figure(2 ); clf(2 ); figure(3 ); clf(3 ); figure(4 ); clf(4 ); figure(5 ); clf(5 );
35+ [BNS ,NSBH ,BBH ,BNSCE ,NSBHCE ,BBHCE ,CEBBH1 ]=DCOplot(filename1 , name1 , 1 , ' r' , 40 );
36+ fprintf(' \n DCOs:\t\t #Merging DNS\t #Merging NSBH\t #Merging BBH\t%% BNS via CE\t%% NSBH via CE\t%% BBH via CE\t%% BBH via CE with BH primary\n ' );
37+ fprintf(' %s :\t%d\t\t%d\t\t%d\t\t%.0f\t\t%.0f\t\t%.0f\t\t%.0f\n ' , ...
38+ name1 , BNS , NSBH , BBH , BNSCE / BNS * 100 , NSBHCE / NSBH * 100 , BBHCE / BBH * 100 , CEBBH1 / BBH * 100 );
3739 if (nargin == 4 ),
38- [BNS ,NSBH ,BBH ,BNSCE ,NSBHCE ,BBHCE ]=DCOplot(filename2 , name2 , 1 , 2 , ' b' , 20 );
39- fprintf(' %s :\t%d\t\t%d\t\t%d\t\t%.0f\t\t%.0f\t\t%.0f\n ' , ...
40- name2 , BNS , NSBH , BBH , BNSCE / BNS * 100 , NSBHCE / NSBH * 100 , BBHCE / BBH * 100 );
40+ [BNS ,NSBH ,BBH ,BNSCE ,NSBHCE ,BBHCE , CEBBH1 ]=DCOplot(filename2 , name2 , 1 , ' b' , 20 );
41+ fprintf(' %s :\t%d\t\t%d\t\t%d\t\t%.0f\t\t%.0f\t\t%.0f\t\t%.0f\ n ' , ...
42+ name2 , BNS , NSBH , BBH , BNSCE / BNS * 100 , NSBHCE / NSBH * 100 , BBHCE / BBH * 100 , CEBBH1 / BBH * 100 );
4143 end ;
42- figure(1 ), hold off ; figure(2 ), hold off ;
43- figure(1 ), set(gca ,' FontSize' ,20 ), xlabel(' Mass 1 (M$_\odot$)' , ' Interpreter' , ' latex' ),
44- ylabel(' Mass 2 (M$_\odot$)' , ' Interpreter' , ' latex' ), title(' Merging DCO masses' ); legend ;
44+ figure(1 ), hold off ; figure(2 ), hold off ; figure( 3 ), hold off ; figure( 4 ), hold off ;
45+ figure(1 ), set(gca ,' FontSize' ,20 ), xlabel(' $M_1$ (M$_\odot$)' , ' Interpreter' , ' latex' ),
46+ ylabel(' $M_2$ (M$_\odot$)' , ' Interpreter' , ' latex' ), title(' Merging DCO masses' ); legend ;
4547 figure(2 ), set(gca ,' FontSize' ,20 ), xlabel(' $\log_{10}$ (Orbital period/hr)' , ' Interpreter' , ' latex' );
4648 ylabel(' Eccentricity' ), title(' DNS at formation' ); legend ;
47-
49+ figure(3 ), set(gca ,' FontSize' ,20 ), xlabel(' Chirp Mass (M$_\odot$)' , ' Interpreter' , ' latex' ),
50+ ylabel(' $\log_{10} (P_\mathrm{orb}/\mathrm{d})$' , ' Interpreter' , ' latex' ), title(' Merging BBH at formation' ); legend ;
51+ figure(4 ), set(gca ,' FontSize' ,20 ); xlabel(' $M_1$ (M$_\odot$)' , ' Interpreter' , ' latex' ),
52+ ylabel(' $M_{\t extrm{core},2}$ (M$_\odot$)' , ' Interpreter' , ' latex' ), title(' CE from 2->1 en route to merging BBH' ), legend ;
53+ figure(5 ), set(gca ,' FontSize' ,20 ); xlabel(' $q \equiv M_2/M_1$' , ' Interpreter' , ' latex' ),
54+ ylabel(' CDF' ), title(' CDF of merging BBH' ), legend ;
55+ fignumber= 6 ;
4856
4957 % Plot BH HMXBs
50- figure(3 ); clf(3 );
51- HMXBplot(filename1 , name1 , 3 , ' r' , 40 );
58+ figure(fignumber ); clf(fignumber );
59+ HMXBplot(filename1 , name1 , fignumber , ' r' , 40 );
5260 if (nargin == 4 ),
53- HMXBplot(filename2 , name2 , 3 , ' b' , 20 );
61+ HMXBplot(filename2 , name2 , fignumber , ' b' , 20 );
5462 end ;
55- figure(3 ); hold off ; set(gca ,' FontSize' ,20 ), legend ; title(' HMXB masses' );
63+ figure(fignumber ); hold off ; set(gca ,' FontSize' ,20 ), legend ; title(' HMXB masses' );
5664 xlabel(' BH mass (M$_\odot$)' , ' Interpreter' , ' latex' ), ylabel(' Companion mass (M$_\odot$)' , ' Interpreter' , ' latex' );
57-
65+ fignumber = fignumber + 1 ;
5866
5967 % Plot BeXRBs
60- figure(4 ); clf(4 );
61- BeXRBplot(filename1 , name1 , 4 , ' r' , 40 );
68+ figure(fignumber ); clf(fignumber );
69+ BeXRBplot(filename1 , name1 , fignumber , ' r' , 40 );
6270 if (nargin == 4 ),
63- BeXRBplot(filename2 , name2 , 4 , ' b' , 20 );
71+ BeXRBplot(filename2 , name2 , fignumber , ' b' , 20 );
6472 end ;
65- figure(4 ); hold off ; set(gca ,' FontSize' ,20 ), legend ;
73+ figure(fignumber ); hold off ; set(gca ,' FontSize' ,20 ), legend ;
6674 xlabel(' Companion mass (M$_\odot$)' , ' Interpreter' , ' latex' );
6775 ylabel(' Formation time, Myr' ), title(' BeXRBs just after SN' );
76+ fignumber= fignumber + 1 ;
6877
6978 % Plot LMXBs/IMXBs
70- figure(5 ); clf(5 );
71- [LMXBcount , NSLMXBcount ]=LMXBplot(filename1 , name1 , 5 , ' r' , 40 );
79+ figure(fignumber ); clf(fignumber );
80+ [LMXBcount , NSLMXBcount ]=LMXBplot(filename1 , name1 , fignumber , ' r' , 40 );
7281 fprintf(' \n LMXBs:\t\t #LMXB\t\t #NS LMXB\n ' );
7382 fprintf(' %s :\t%d\t\t%d\n ' , name1 , LMXBcount , NSLMXBcount );
7483 if (nargin == 4 ),
75- [LMXBcount , NSLMXBcount ]=LMXBplot(filename2 , name2 , 5 , ' b' , 20 );
84+ [LMXBcount , NSLMXBcount ]=LMXBplot(filename2 , name2 , fignumber , ' b' , 20 );
7685 fprintf(' %s :\t%d\t\t%d\n ' , name2 , LMXBcount , NSLMXBcount );
7786 end ;
78- figure(5 ), hold off ; set(gca ,' FontSize' , 20 ); legend ; title(' LMXB on first MT onto CO' );
87+ figure(fignumber ), hold off ; set(gca ,' FontSize' , 20 ); legend ; title(' LMXB on first MT onto CO' );
7988 xlabel(' Compact object mass (M$_\odot$)' , ' Interpreter' , ' latex' ), ylabel(' Companion mass (M$_\odot$)' , ' Interpreter' , ' latex' );
89+ fignumber= fignumber + 1 ;
8090
8191 % Plot DWDs (just as a sanity check)
82- figure(6 ); clf(6 );
83- DWDplot(filename1 , name1 , 6 , ' r' , 40 );
92+ figure(fignumber ); clf(fignumber );
93+ DWDplot(filename1 , name1 , fignumber , ' r' , ' m ' , 10 );
8494 if (nargin == 4 ),
85- DWDplot(filename2 , name2 , 6 , ' b' , 20 );
95+ DWDplot(filename2 , name2 , fignumber , ' b' , ' g ' , 5 );
8696 end ;
87- figure(6 ); hold off ; set(gca ,' FontSize' ,20 ); title(' Double White Dwarfs' ); legend ;
97+ figure(fignumber ); hold off ; axis([- 3 5 - 3 5 ]); set(gca ,' FontSize' ,20 ); title(' Double White Dwarfs' ); legend ;
8898 xlabel(' $log_{10}(M*a) [M_\odot * R_\odot]$ @ ZAMS' , ' Interpreter' , ' latex' );
8999 ylabel(' $log_{10}(M*a) [M_\odot * R_\odot]$ @ end' , ' Interpreter' , ' latex' );
90100
@@ -151,8 +161,8 @@ function ComparisonPlots(filename1, name1, filename2, name2)
151161
152162
153163% Plot double compact objects; returns DCO counts
154- function [BNScount , NSBHcount , BBHcount , BNSCE , NSBHCE , BBHCE ] = ...
155- DCOplot(file , name , fignumberDCO , fignumberDNSPE , colour , point )
164+ function [BNScount , NSBHcount , BBHcount , BNSCE , NSBHCE , BBHCE , CEBBH1count ] = ...
165+ DCOplot(file , name , fignumber , colour , point )
156166 global Msunkg G AU
157167 type1= h5read(file ,' /BSE_Double_Compact_Objects/Stellar_Type(1)' );
158168 type2= h5read(file ,' /BSE_Double_Compact_Objects/Stellar_Type(2)' );
@@ -170,36 +180,45 @@ function ComparisonPlots(filename1, name1, filename2, name2)
170180 NSBH= (((type1 == 13 ) & (type2 == 14 )) | ((type1 == 14 ) & (type2 == 13 )));
171181 mergingDCO= mergingBNS | mergingNSBH | mergingBBH ;
172182 BNScount= sum(mergingBNS ); NSBHcount= sum(mergingNSBH ); BBHcount= sum(mergingBBH );
173- masschirp= mass1 .^ 0.6 .* mass2 .^ 0.6 ./(mass1 + mass2 ).^0.2 ;
183+ chirpmass= mass1 .^ 0.6 .* mass2 .^ 0.6 ./(mass1 + mass2 ).^0.2 ;
184+ q= mass2 ./ mass1 ;
174185 seedCE= h5read(file ,' /BSE_Common_Envelopes/SEED' );
175186 [isCE ,CEIndex ]=ismember(seedDCO ,seedCE );
176187 optCE= h5read(file ,' /BSE_Common_Envelopes/Optimistic_CE' );
177188 RLOFCE= h5read(file ,' /BSE_Common_Envelopes/Immediate_RLOF>CE' );
178189 OKCE= zeros(size(mergingDCO )); OKCE(CEIndex == 0 )=1 ; OKCE(CEIndex > 0)=(~optCE(CEIndex(CEIndex > 0))) & (~RLOFCE(CEIndex(CEIndex > 0)));
179190 BNSCE= sum(mergingBNS & isCE & OKCE ); NSBHCE= sum(mergingNSBH & isCE & OKCE ); BBHCE= sum(mergingBBH & isCE & OKCE );
191+ type1CE = h5read(file ,' /BSE_Common_Envelopes/Stellar_Type(1)<CE' );
192+ type2CE = h5read(file ,' /BSE_Common_Envelopes/Stellar_Type(2)<CE' );
193+ mass1CE = h5read(file ,' /BSE_Common_Envelopes/Mass(1)<CE' );
194+ mass2coreCE = h5read(file ,' /BSE_Common_Envelopes/Mass(2)<CE' )-h5read(file ,' /BSE_Common_Envelopes/Mass_Env(2)' );
195+ [CEtowardmergingBBH ,BBHindex ]=ismember(seedCE ,seedDCO(mergingBBH )); % pick CEs leading to merging BBHs
196+ % fprintf('Validation: %d valid CE episodes en route to merging BBH, including %d unique CEs, which should match %d merging BBHs produced via CE\n', ...
197+ % sum(CEtowardmergingBBH & ~optCE & ~RLOFCE), length(unique(seedCE(CEtowardmergingBBH & ~optCE & ~RLOFCE))), sum(mergingBBH & isCE & OKCE))
198+ CEBBH1 = CEtowardmergingBBH & ~optCE & ~RLOFCE & type1CE == 14 ; CEBBH1count= sum(CEBBH1 );
180199 % Merging DCO mass distribution
181- figure(fignumberDCO ), hold on ;
200+ figure(fignumber ), hold on ;
182201 scatter(mass1(mergingDCO & isCE & OKCE ), mass2(mergingDCO & isCE & OKCE ), point , ...
183202 ' filled' , colour , ' DisplayName' , [' CE, ' , name ]);
184203 scatter(mass1(mergingDCO & ~isCE ), mass2(mergingDCO & ~isCE ), point , colour , ' DisplayName' , [' Stable, ' , name ]);
185204 % P-e of Galactic DNS at formation
186205 P= 2 * pi *(a * AU ).^(3 / 2 )./sqrt(G * Msunkg *(mass1 + mass2 ))/3600 ; % orbital period at DCO, in hours
187- figure(fignumberDNSPE ), hold on ;
206+ figure(fignumber + 1 ), hold on ;
188207 scatter(log10(P(BNS & isCE & OKCE )), e(BNS & isCE & OKCE ), point , ' filled' , colour , ...
189208 ' DisplayName' , [' CE, ' , name ]); hold on ;
190209 scatter(log10(P(BNS & ~isCE )), e(BNS & ~isCE ), point , colour , ' DisplayName' , [' Stable, ' , name ]);
191-
192- % %% TO CLEAN
193- good = mergingBBH & isCE & OKCE ;
194- type1CE = h5read( file , ' /BSE_Common_Envelopes/Stellar_Type(1)<CE ' );
195- type2CE = h5read( file , ' /BSE_Common_Envelopes/Stellar_Type(2)<CE ' );
196- mass1CE = h5read( file , ' /BSE_Common_Envelopes/Mass(1)<CE ' );
197- mass2coreCE = h5read( file , ' /BSE_Common_Envelopes/Mass(2)<CE ' )-h5read( file , ' /BSE_Common_Envelopes/Mass_Env(2) ' ) ;
198- unique(type1CE(CEIndex( good ))), sum( good ), sum(type1CE(CEIndex( good ))== 14 )
199- figure( 7 ), scatter(mass1CE(CEIndex( good )), mass2coreCE(CEIndex( good )), 20 , ' filled ' ), set( gca , ' FontSize ' , 20 ); xlabel( ' M_1, Msun ' ); ylabel( ' M_{core,2}, Msun ' );
200- typesCE = [type1CE(CEIndex( good )) type2CE(CEIndex( good ))] ;
201- % good(find(type1CE(CEIndex(good))~=14))
202-
210+ % Chirp mass vs period at BBH formation
211+ figure( fignumber + 2 ), hold on ;
212+ scatter(chirpmass( mergingDCO & isCE & OKCE ), log10(P( mergingDCO & isCE & OKCE )/ 24 ), point , ...
213+ ' filled ' , colour , ' DisplayName ' , [ ' CE, ' , name ]);
214+ scatter(chirpmass( mergingDCO & ~ isCE ), log10(P( mergingDCO & ~ isCE )/ 24 ), point , colour , ' DisplayName ' , [ ' Stable, ' , name ] );
215+ % Masses at the time of CE
216+ figure( fignumber + 3 ), hold on ;
217+ scatter(mass1CE( CEBBH1 ), mass2coreCE( CEBBH1 ), point , ' filled ' , colour , ' DisplayName ' , name );
218+ % CSF of q distribution
219+ figure( fignumber + 4 ), hold on ;
220+ p = cdfplot(q( mergingBBH & isCE & OKCE )), p.Color = colour , p.LineStyle = ' - ' , p.LineWidth = point / 10 , p.DisplayName = [ ' CE, ' , name ];
221+ p = cdfplot(q( mergingBBH & ~ isCE )), p.Color = colour , p.LineStyle = ' : ' , p.LineWidth = point / 10 , p.DisplayName = [ ' Stable, ' , name ];
203222end % end of DCOplot
204223
205224% Plot BH HMXBs
@@ -310,7 +329,7 @@ function BeXRBplot(file, name, fignumberBeXRB, colour, point)
310329% Plot a*Mtot for white dwarfs
311330% Should be constant for non-mass-transferring WDs in the absence of tides
312331% and GR (and supernovae, which is why we focus on WDs as a sanity check
313- function DWDplot(file , name , fignumberDWD , colour , point )
332+ function DWDplot(file , name , fignumberDWD , colourNMT , colourMT , point )
314333 global AU Rsun
315334 M1ZAMS= h5read(file , ' /BSE_System_Parameters/Mass@ZAMS(1)' );
316335 M2ZAMS= h5read(file , ' /BSE_System_Parameters/Mass@ZAMS(2)' );
@@ -330,10 +349,12 @@ function DWDplot(file, name, fignumberDWD, colour, point)
330349 SeedRLOF= h5read(file , ' /BSE_RLOF/SEED' );
331350 [hadRLOF ,RLOFIndex ]=ismember(ind ,SeedRLOF );
332351 figure(fignumberDWD ), hold on ;
333- scatter(log10(MAZAMS(~hadRLOF )), log10(MAWD(~hadRLOF )), point , ' filled' , colour , ' DisplayName' , [' DWDs without mass transfer, ' , name ]);
334- scatter(log10(MAZAMS(hadRLOF )), log10(MAWD(hadRLOF )), point , colour , ' DisplayName' , [' DWDs after mass transfer, ' , name ]);
352+ scatter(log10(MAZAMS(~hadRLOF )), log10(MAWD(~hadRLOF )), point , ' filled' , colourNMT , ' DisplayName' , [' DWDs without mass transfer, ' , name ]);
353+ scatter(log10(MAZAMS(hadRLOF )), log10(MAWD(hadRLOF )), point , colourMT , ' DisplayName' , [' DWDs after mass transfer, ' , name ]);
335354end % end of DWDplot
336355
356+
357+
337358% SN varieties -- just count
338359function [binariescount , SNcount , BHcompletecount , SNbothcount , SNonecount , ...
339360 unboundcount , unbound1count , unbound2count , ...
0 commit comments