From bf16afa6c12c2512bc984c888959b6265d847ff4 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Sun, 15 Jun 2025 19:12:15 +1000 Subject: [PATCH 1/4] Starting work on public cosmic history integrator --- compas_matlab_utils/BSEDetailedOutput.m | 72 ++++++ compas_matlab_utils/CosmicHistoryIntegrator.m | 224 ++++++++++++++++++ compas_matlab_utils/SSEDetailedOutput.m | 115 +++++++++ 3 files changed, 411 insertions(+) create mode 100644 compas_matlab_utils/BSEDetailedOutput.m create mode 100644 compas_matlab_utils/CosmicHistoryIntegrator.m create mode 100644 compas_matlab_utils/SSEDetailedOutput.m diff --git a/compas_matlab_utils/BSEDetailedOutput.m b/compas_matlab_utils/BSEDetailedOutput.m new file mode 100644 index 000000000..e04301c7e --- /dev/null +++ b/compas_matlab_utils/BSEDetailedOutput.m @@ -0,0 +1,72 @@ +file='~/Work/COMPAS/src/COMPAS_Output/Detailed_Output/BSE_Detailed_Output_0.h5'; + +tic +time=h5read(file,'/Time'); +MThistory=h5read(file,'/MT_History'); +Z=h5read(file,'/Metallicity@ZAMS(1)'); +mass1=h5read(file,'/Mass(1)'); +massCO1=h5read(file,'/Mass_CO_Core(1)'); +mass2=h5read(file,'/Mass(2)'); +massCO2=h5read(file,'/Mass_CO_Core(2)'); +massHe1=h5read(file,'/Mass_He_Core(1)'); +massHe2=h5read(file,'/Mass_He_Core(2)'); +luminosity1=h5read(file,'/Luminosity(1)'); +luminosity2=h5read(file,'/Luminosity(2)'); +type1=h5read(file,'/Stellar_Type(1)'); +type2=h5read(file,'/Stellar_Type(2)'); +radius1=h5read(file,'/Radius(1)'); +radius2=h5read(file,'/Radius(2)'); +RL1=h5read(file,'/RocheLobe(1)'); +RL2=h5read(file,'/RocheLobe(2)'); +a=h5read(file,'/SemiMajorAxis'); +e=h5read(file,'/Eccentricity'); +rec=h5read(file,'/Record_Type'); +AM1=h5read(file, '/Ang_Momentum(1)'); +AM2=h5read(file, '/Ang_Momentum(2)'); +AMtot=h5read(file, '/Ang_Momentum_Total'); +omega1=h5read(file, '/Omega(1)'); +omega2=h5read(file, '/Omega(2)'); +omegab1=h5read(file, '/Omega_Break(1)'); +omegab2=h5read(file, '/Omega_Break(2)'); +MTtimescale=h5read(file, '/MassTransferTimescale'); +inf=h5info(file); +%inf.Datasets.Name + + +disp('Time (Myr), Event, M1 (M_o), type1, M2 (M_o), type2, a (R_o), e'); +disp([num2str(time(1)), ' start: Z=', num2str(Z(1)), ' ', num2str(mass1(1)), ' ', num2str(type1(1)),... + ' ', num2str(mass2(1)), ' ', num2str(type2(1)), ' ', num2str(a(1)), ' ', num2str(e(1))]); +for(i=2:length(time)), + if(rec(i)==4), %only print records at end of timestep + prev=find(rec(1:i-1)==4, 1, 'last'); %previous timestep index + if(MThistory(i)~=0),% & MThistory(i)~=MThistory(i-1)), + switch (MThistory(i)), case 1, MTstring='Stable MT from 1 TO 2'; case 2, MTstring='Stable MT from 2 TO 1';... + case 3, MTstring='CE from 1 to 2'; case 4, MTstring='CE from 2 to 1'; case 5, MTstring='CE Double Core'; ... + case 6, MTstring='CE merger'; end; + switch(MTtimescale(i)), case 1, MTtimescalestring='Nuclear'; case 2, MTtimescalestring='Thermal'; case 3, ... + MTtimescalestring='Dynamical'; end; + disp([num2str(time(i)), ' ', MTstring, ', @', MTtimescalestring, ' ', num2str(mass1(i)), ' ', num2str(type1(i)),... + ' ', num2str(mass2(i)), ' ', num2str(type2(i)), ' ', num2str(a(i)), ' ', num2str(e(i))]); + end; + if(type1(i)~=type1(prev)) + disp([num2str(time(i)), ' star 1 ', num2str(type1(prev)), '->', num2str(type1(i)), ' ', num2str(mass1(i)), ' ', num2str(type1(i)),... + ' ', num2str(mass2(i)), ' ', num2str(type2(i)), ' ', num2str(a(i)), ' ', num2str(e(i))]); + end; + if(type2(i)~=type2(prev)), + disp([num2str(time(i)), ' star 2 ', num2str(type2(prev)), '->', num2str(type2(i)), ' ', num2str(mass1(i)), ' ', num2str(type1(i)),... + ' ', num2str(mass2(i)), ' ', num2str(type2(i)), ' ', num2str(a(i)), ' ', num2str(e(i))]); + end; + end; +end; +if((type1(i)==13 || type1(i)==14) && (type2(i)==13 || type2(i)==14) ), + Msunkg=1.98892e30; c=299792458; G=6.67428e-11; Rsun = 695500000; + beta=64/5*G^3*mass1(i)*mass2(i)*(mass1(i)+mass2(i))*Msunkg^3/c^5; T0=(a(i)*Rsun)^4/4/beta; + Tdelay=T0*(1-e(i)^2).^(7/2).*(1+0.27*e(i)^10 + 0.33*e(i)^20 + 0.2*e(i)^1000)/3.15e7/1e6; + disp([num2str(time(i)+Tdelay), ' GW merger in ', num2str(Tdelay,'%.0f'), ' Myr ', num2str(mass1(i)), ' ', num2str(type1(i)),... + ' ', num2str(mass2(i)), ' ', num2str(type2(i))]); +end; +toc + +figure(81), plot(time, mass1, 'LineWidth', 3), set(gca,'FontSize',20), xlabel('Time, Myr'), ylabel('Total mass 1, M_o') +%figure(81), scatter(time, radius1, 30, type1, 'filled'), set(gca,'FontSize',20), xlabel('Time, Myr'), ylabel('Radius 1, R_o') +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') diff --git a/compas_matlab_utils/CosmicHistoryIntegrator.m b/compas_matlab_utils/CosmicHistoryIntegrator.m new file mode 100644 index 000000000..03a8444bf --- /dev/null +++ b/compas_matlab_utils/CosmicHistoryIntegrator.m @@ -0,0 +1,224 @@ +function [Zlist, MergerRateByRedshiftByZ, SFRfractionZ]=... + CosmicHistoryIntegrator(filename,zlist,metallicitychoice, makeplots) +% Integrator for the binary black hole merger rate over cosmic history +% COMPAS (Compact Object Mergers: Population Astrophysics and Statistics) +% software package +% +% USAGE: +% [Zlist, MergerRateByRedshiftByZ]=... +% CosmicHistoryIntegrator(filename, zlistformation, zlistmerger, zlistdetection, Msimulated, [,makeplots]) +% +% INPUTS: +% filename: name of population synthesis input file +% should be in COMPAS output h5 format +% zlistformation: vector of redshifts at which the formation rate is +% computed +% zlistmerger: vector of redshifts at which the merger rate is computed +% zlistdetection: vector of redshifts at which the detection rate is +% computed +% Msimulated: total star forming mass represented by the simulation (for +% normalisation) +% makeplots: if set to 1, generates a set of useful plots (default = 0) +% +% Zlist is a vector of metallicities, taken from input file +% MergerRateByRedshiftByZ is a matrix of size length(Zlist) X length(zlist) +% which contains a merger rate of binary black holes in the given redshift +% and metallicity bin, in units of mergers per Gpc^3 of comoving volume per +% year of source time +% +% EXAMPLE: +% zlist=0:0.1:5; +% filename=['~/Work/COMPASresults/popsynth/runs/',... +% '20170628-Coen-Pessimistic/AllMergers.dat']; +% [Zlist,MergerRateByRedshiftByZ]=CosmicHistoryIntegrator(filename,zlist); +% figure(1),colormap jet; +% plot(zlist, MergerRateByRedshiftByZ, 'LineWidth', 2), +% legend(num2str(Zlist)), +% set(gca, 'FontSize', 20); %for labels +% xlabel('z'), +% ylabel('Pessimistic non-RLOF BBH merger rate, per Gpc^3 per yr') +% sum(MergerRateByRedshiftByZ(:,1)) %Total BBH merger rate at z=0 +% [Zlist,MergerRateByRedshiftByZLanger]=CosmicHistoryIntegrator(filename,zlist,1); +% [Zlist,MergerRateByRedshiftByZLambert]=CosmicHistoryIntegrator(filename,zlist,2); +% plot(zlist, sum(MergerRateByRedshiftByZLambert,2), 'r', ... +% zlist, sum(MergerRateByRedshiftByZLanger,2), 'b', 'LineWidth', 2), +% set(gca, 'FontSize', 20), xlabel('z'), legend('Langer & Norman','Lambert') +% ylabel('Pessimistic non-RLOF BBH merger rate, per Gpc^3 per yr'), + + + +%define constants +global Mpcm; +global Mpc; +global yr; +Mpcm=1*10^6 * 3.0856775807e16; %Mpc in meters +c=299792458; %speed of light, m/s +Mpc=Mpcm/c; %Gpc in seconds +yr=3.15569e7; %year in seconds + +if (nargin<2) + error('Not enough input arguments.'); +end; +if (nargin<4), makeplots=0; end; + +%cosmology calculator +[zvec,tL]=Cosmology(); +%load COMPAS data +[M1,M2,Z,Tdelay]=DataRead(filename, max(tL)); +%metallicity-specific SFR +[Zlist,SFR,SFRfractionZ]=Metallicity(Z,zvec,metallicitychoice); + + +%Consider the contribution of every simulated binary to the merger rate +%in every redshift bin by considering when it would have to be formed to +%merge at that redshift and normalizing by the relevant +%metallicity-specific star formation rate +dz=zvec(2)-zvec(1); +tLmerge=tL(floor(zlist/(zvec(2)-zvec(1)))+1); +MergerRateByRedshiftByZ=zeros(length(zlist),length(Zlist)); +for(i=1:length(M1)), + Zcounter=find(Zlist==Z(i)); + zformindex=find(tL>=Tdelay(i),1); + for(k=1:length(zlist)), + zformindex=find(tL>=(Tdelay(i)+tLmerge(k)),1); + MergerRateByRedshiftByZ(k,Zcounter)=... + MergerRateByRedshiftByZ(k,Zcounter)+... + SFR(zformindex)*SFRfractionZ(zformindex,Zcounter)/Msimulated; + end; +end; + +if(makeplots==1), %make a set of default plots + MakePlots(M1,M2,Z,Tdelay,zvec,SFR,SFRfractionZ,... + zlist,Zlist,MergerRateByRedshiftByZ); +end; + +end %end of CosmicHistoryIntegrator + + +%Load the data stored in COMPAS output format from a file +%Select only binary black hole mergers of interest, and return the +%component masses, metallicities, and star formation to merger delay times +function [M1,M2,Z,Tdelay]=DataRead(filename, tHubble) + global BH + if(exist(filename, 'file')~=2), + error('Input file does not exist'); + end; + data = importdata(filename, '\t', 2); + + colnames=data.textdata(2,:); + Z1index=find(strcmp(colnames,'Metallicity1')); + M1index=find(strcmp(colnames,'M1')); + M2index=find(strcmp(colnames,'M2')); + tcindex=find(strcmp(colnames,'tc')); + tformindex=find(strcmp(colnames,'tform')); + type1index=find(strcmp(colnames,'stellarType1')); + type2index=find(strcmp(colnames,'stellarType2')); + nonRLOF=data.data(:,find(strcmp(colnames,'RLOFSecondaryAfterCEE')))==0; + Pessimistic=data.data(:,find(strcmp(colnames,'optimisticCEFlag')))==0; + tc=data.data(:,tcindex); + tform=data.data(:,tformindex); + Ttotal=(tc+tform)*1e6; %convert Megayears to years + %select only binaries where both members are black holes at the last + %step, the total delay time is less than the age of the Universe, there + %is no Roche-lobe overflow immediately after the common-envelope phase, + %and the Pessimistic common-envelope conditions are met + select=Ttotal=min(Z),1, 'first'); + maxlogZindex=find(exp(logZvec)<=max(Z),1, 'last'); + dPdlogZ(minlogZindex,:)=dPdlogZ(minlogZindex,:)+sum(dPdlogZ(1:minlogZindex,:),1)*dlogZ/(sum(Z==min(Z))/length(Z))*PdrawZ; + dPdlogZ(maxlogZindex,:)=dPdlogZ(maxlogZindex,:)+sum(dPdlogZ(maxlogZindex:end,:),1)*dlogZ/(sum(Z==max(Z))/length(Z))*PdrawZ; + dPdlogZ=dPdlogZ./(sum(dPdlogZ,1)*dlogZ); %normalise +end %end of Metallicity + + +%Make a set of default plots +function MakePlots(M1,M2,Z,Tdelay,zvec,SFR,SFRfractionZ,... + zlist,Zlist,MergerRateByRedshiftByZ) + + figure(1),colormap jet; + plot(zlist, MergerRateByRedshiftByZ, 'LineWidth', 2), + legend(num2str(Zlist)), + set(gca, 'FontSize', 20); %for labels + xlabel('z'), + ylabel('Pessimistic non-RLOF BBH merger rate, per Gpc^3 per yr') + disp(['Total BBH merger rate at z=0: ', ... + num2str(sum(MergerRateByRedshiftByZ(:,1))),... + ' per Gpc^3 per year']); + + + figure(2), colormap jet; + scatter(M1,M2,20,log(Z)/log(10),'filled'); + set(gca, 'FontSize', 20); %for labels + H=colorbar; H.Label.String='log_{10} metallicity'; + xlabel('Primary BH mass [M_o]'), ylabel('Secondary BH mass [M_o]'); + + figure(3), colormap jet; + scatter(M1+M2,log(Tdelay/1e6)/log(10),20,log(Z)/log(10),'filled'); + set(gca, 'FontSize', 20); %for labels + H=colorbar; H.Label.String='log_{10} metallicity'; + xlabel('Total BBH mass [M_o]'), ylabel('log(Tdelay/Myr)'); + + figure(4), colormap jet; + plot(zvec, SFR) + set(gca, 'FontSize', 20); %for labels + xlabel('z'), ylabel('Star-formation rate, M_o per Gpc^3 per yr'); + + + figure(5), colormap jet; + plot(zvec, SFRfractionZ) + set(gca, 'FontSize', 20); %for labels + legend(num2str(Zlist)) + xlabel('z'), ylabel('Z-specific SFR fraction'); + +end %end of MakePlots + + + diff --git a/compas_matlab_utils/SSEDetailedOutput.m b/compas_matlab_utils/SSEDetailedOutput.m new file mode 100644 index 000000000..e58a023a9 --- /dev/null +++ b/compas_matlab_utils/SSEDetailedOutput.m @@ -0,0 +1,115 @@ +function SSEDetailedOutput(filename) +% Carries out some basic analysis and makes plots for a single SSE run +% +% USAGE: +% ComparisonPlots(filename1, name1, filename2, name2) +% +% INPUTS: +% filename: name of detailed output file in COMPAS h5 format +% +% example: +% SSEDetailedOutput('~/Work/COMPAS/src/COMPAS_Output/Detailed_Output/SSE_Detailed_Output_0.h5') +% + + +time=h5read(file,'/Time'); +Z=h5read(file,'/Metallicity@ZAMS'); +mass=h5read(file,'/Mass'); +massCO=h5read(file,'/Mass_CO_Core'); +massHe=h5read(file,'/Mass_He_Core'); +luminosity=h5read(file,'/Luminosity'); +type=h5read(file,'/Stellar_Type'); +radius=h5read(file,'/Radius'); +temperature=h5read(file,'/Teff'); +mdot=h5read(file,'/Mdot'); +ml=h5read(file,'/Dominant_Mass_Loss_Rate'); +envMass=h5read(file,'/Mass_Env'); +convEnvMass=h5read(file,'/Mass_Convective_Env'); +lambda=h5read(file,'/Lambda_Convective'); +binding=h5read(file,'/BE_ConvectiveEnvelope'); +rec=h5read(file,'/Record_Type'); + +figure(1), +scatter(time(type==1),radius(type==1),20,'filled', 'b'); hold on; +scatter(time(type==2),radius(type==2),20,'filled', 'g'); +scatter(time(type==3),radius(type==3),20,'filled', 'y'); +scatter(time(type==4),radius(type==4),20,'filled', 'r'); +scatter(time(type==5),radius(type==5),20,'filled', 'm'); +scatter(time(type==6),radius(type==6),20,'filled', 'k'); +hold off; +set(gca,'FontSize',20); xlabel('Time, Myr'); ylabel('Radius, Rsun'), legend('MS','HG','FGB','CHeB','EAGB','TPAGB') + +for(i=2:length(radius)), %remove points where radius is less than the radius at previous time + if(~isempty(find(radius(1:i-1)>=radius(i)))), + envMass(i)=0; convEnvMass(i)=0; radius(i)=0; + end; +end; + + +goodindices=type>1 & type<=6 & radius>0; +figure(2); +subplot(2,1,1), scatter(radius(goodindices), envMass(goodindices)-convEnvMass(goodindices), 20, 'b', 'filled'); hold on; +scatter(radius(goodindices), convEnvMass(goodindices), 20, 'm', 'filled'); hold on; +scatter(radius(goodindices), mass(goodindices)-envMass(goodindices), 20, 'y', 'filled'); +scatter(radius(goodindices), mass(goodindices), 20, 'k', 'filled'); +hold off; +set(gca,'FontSize',20); xlabel('Radius, Rsun'); ylabel('Mass, Msun'), legend('Radiative intershell','Convective Envelope','Core', 'Total') +subplot(2,1,2); scatter(radius(goodindices), binding(goodindices), 20, 'm', 'filled') +set(gca,'FontSize',20); xlabel('Radius, Rsun'); ylabel('Binding E conv env, erg') +%subplot(2,1,2); scatter(radius(goodindices), lambda(goodindices), 20, 'm', 'filled') +%set(gca,'FontSize',20); xlabel('Radius, Rsun'); ylabel('\lambda') + + + +file='~/Work/COMPAS/src/COMPAS_Output/COMPAS_Output.h5'; +radius=h5read(file,'/BSE_Common_Envelopes/Radius(1)CE'); +apost=h5read(file,'/BSE_Common_Envelopes/SemiMajorAxis>CE'); +primarydonor=(h5read(file,'/BSE_Common_Envelopes/RLOF(1)')==1); +figure(3); semilogy(radius(primarydonor), apre(primarydonor), 'b-.', radius(primarydonor), a1(primarydonor), 'b--', radius(primarydonor), apost(primarydonor), 'b-', 'LineWidth', 3); hold on; +file='~/Work/COMPAS/src/COMPAS_Output_1/COMPAS_Output.h5'; +radius=h5read(file,'/BSE_Common_Envelopes/Radius(1)CE'); +apost=h5read(file,'/BSE_Common_Envelopes/SemiMajorAxis>CE'); +primarydonor=(h5read(file,'/BSE_Common_Envelopes/RLOF(1)')==1); +figure(3); semilogy(radius(primarydonor), apre(primarydonor), 'm-.', radius(primarydonor), a1(primarydonor), 'm--', radius(primarydonor), apost(primarydonor), 'm-', 'LineWidth', 3); hold on; +file='~/Work/COMPAS/src/COMPAS_Output_2/COMPAS_Output.h5'; +radius=h5read(file,'/BSE_Common_Envelopes/Radius(1)CE'); +apost=h5read(file,'/BSE_Common_Envelopes/SemiMajorAxis>CE'); +primarydonor=(h5read(file,'/BSE_Common_Envelopes/RLOF(1)')==1); +figure(3); semilogy(radius(primarydonor), apre(primarydonor), 'y-.', radius(primarydonor), a1(primarydonor), 'y--', radius(primarydonor), apost(primarydonor), 'y-', 'LineWidth', 3); hold off; +set(gca,'FontSize',20); xlabel('Radius, Rsun'); ylabel('Semimajor axis, Rsun'), legend('Pre CE','After Stage1','After Stage 2') + + + +for k=0:4, +file=['~/Work/COMPAS/src/COMPAS_Output_3/Detailed_Output/SSE_Detailed_Output_',int2str(k),'.h5']; +time=h5read(file,'/Time'); +Z=h5read(file,'/Metallicity@ZAMS'); +mass=h5read(file,'/Mass'); +massCO=h5read(file,'/Mass_CO_Core'); +massHe=h5read(file,'/Mass_He_Core'); +luminosity=h5read(file,'/Luminosity'); +type=h5read(file,'/Stellar_Type'); +radius=h5read(file,'/Radius'); +binding=h5read(file,'/BE_Nanjing'); +rec=h5read(file,'/Record_Type'); +for(i=2:length(radius)), %remove points where radius is less than one of previously achieved radii + if(~isempty(find(radius(1:i-1)>=radius(i)))), + binding(i)=0; radius(i)=0; + end; +end; +figure(k+1), +scatter(radius(type==1),log10(binding(type==1)),20,'filled'); hold on; +scatter(radius(type==2),log10(binding(type==2)),20,'filled'); +scatter(radius(type==3),log10(binding(type==3)),20,'filled'); +scatter(radius(type==4),log10(binding(type==4)),20,'filled'); +scatter(radius(type==5),log10(binding(type==5)),20,'filled'); +scatter(radius(type==6),log10(binding(type==6)),20,'filled'); +hold off; +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))]) +end; From a5794c11eb8ca524270b9a1ad33cbc96757f7f49 Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Wed, 18 Jun 2025 18:40:02 +1000 Subject: [PATCH 2/4] Fix to behaviour for stars that do not qualify for ECSN, see #1305 --- online-docs/COMPAS-2025methodsPaper.bib | 16 ++++++ online-docs/pages/how-to-cite.rst | 69 +++++++++++++------------ online-docs/pages/quick-links.rst | 2 + src/EAGB.cpp | 8 +-- src/GiantBranch.cpp | 14 ++--- src/GiantBranch.h | 2 +- src/HeHG.cpp | 12 ++--- src/TPAGB.cpp | 11 ++-- src/TPAGB.h | 4 +- src/changelog.h | 11 ++-- src/typedefs.h | 1 - 11 files changed, 81 insertions(+), 69 deletions(-) create mode 100644 online-docs/COMPAS-2025methodsPaper.bib diff --git a/online-docs/COMPAS-2025methodsPaper.bib b/online-docs/COMPAS-2025methodsPaper.bib new file mode 100644 index 000000000..a9a9a9579 --- /dev/null +++ b/online-docs/COMPAS-2025methodsPaper.bib @@ -0,0 +1,16 @@ +@ARTICLE{COMPAS:2025, + author = {{Team COMPAS: Mandel}, Ilya and {Riley}, Jeff and {Boesky}, Adam and {Brcek}, Adam and {Hirai}, Ryosuke and {Kapil}, Veome and {Lau}, Mike Y.~M. and {Merritt}, JD and {Rodr{\'\i}guez-Segovia}, Nicol{\'a}s and {Romero-Shaw}, Isobel and {Song}, Yuzhe and {Stevenson}, Simon and {Vajpeyi}, Avi and {van Son}, L.~A.~C. and {Vigna-G{\'o}mez}, Alejandro and {Willcox}, Reinhold}, + title = "{Rapid stellar and binary population synthesis with COMPAS: methods paper II}", + journal = {arXiv e-prints}, + keywords = {Solar and Stellar Astrophysics, High Energy Astrophysical Phenomena, Instrumentation and Methods for Astrophysics}, + year = 2025, + month = jun, + eid = {arXiv:2506.02316}, + pages = {arXiv:2506.02316}, + doi = {10.48550/arXiv.2506.02316}, +archivePrefix = {arXiv}, + eprint = {2506.02316}, + primaryClass = {astro-ph.SR}, + adsurl = {https://ui.adsabs.harvard.edu/abs/2025arXiv250602316M}, + adsnote = {Provided by the SAO/NASA Astrophysics Data System} +} diff --git a/online-docs/pages/how-to-cite.rst b/online-docs/pages/how-to-cite.rst index 6688a2434..30cef9ed4 100644 --- a/online-docs/pages/how-to-cite.rst +++ b/online-docs/pages/how-to-cite.rst @@ -1,33 +1,36 @@ -Citing COMPAS -------------- - -If you use this code or parts of this code for results presented in a scientific publication, we would greatly appreciate you sending -us your paper reference and making your input settings and output data publicly available by uploading it to the COMPAS Zenodo community. - -Please also cite: - -.. _cite-compas: - - Team COMPAS: J. Riley `et al.` [:cite:year:`compas2021`]. |_| |_| |_| |_| |_| |_| :download:`Bibtex citation <../COMPAS-2021methodsPaper.bib>` - -|br| -We would also greatly appreciate an acknowledgement of the form: - -"Simulations in this paper made use of the COMPAS rapid binary population synthesis code (version x.y.z), which is freely available at -http://github.com/TeamCOMPAS/COMPAS." - -|br| -Furthermore, - -If using the COMPAS model of gravitational wave selection effects, please cite :cite:t:`Barrett2018`. - -If you use COMPAS's importance sampling algorithm STROOPWAFEL, please cite :cite:t:`Broekgaarden2019`. - -If using COMPAS's integration over cosmic star formation history, please cite :cite:t:`Neijssel2019`. - -If using the COMPAS model of (pulsational) pair instability supernova, please cite :cite:t:`Stevenson2019`. - -If evolving pulsar spins and magnetic fields with COMPAS, please cite :cite:t:`Chattopadhyay2020`. - -If you use the COMPAS model of chemically homogeneous evolution, please cite :cite:t:`Riley2021`. - +Citing COMPAS +------------- + +If you use this code or parts of this code for results presented in a scientific publication, we would greatly appreciate you sending +us your paper reference and making your input settings and output data publicly available by uploading it to the COMPAS Zenodo community. + +Please also cite: + +.. _cite-compas: + + Team COMPAS: J. Riley `et al.` [:cite:year:`compas2021`]. |_| |_| |_| |_| |_| |_| :download:`Bibtex citation <../COMPAS-2021methodsPaper.bib>` + + Team COMPAS: I. Mandel `et al.` [:cite:year:`compas2025`]. |_| |_| |_| |_| |_| |_| :download:`Bibtex citation <../COMPAS-2025methodsPaper.bib>` + + +|br| +We would also greatly appreciate an acknowledgement of the form: + +"Simulations in this paper made use of the COMPAS rapid binary population synthesis code (version x.y.z), which is freely available at +http://github.com/TeamCOMPAS/COMPAS." + +|br| +Furthermore, + +If using the COMPAS model of gravitational wave selection effects, please cite :cite:t:`Barrett2018`. + +If you use COMPAS's importance sampling algorithm STROOPWAFEL, please cite :cite:t:`Broekgaarden2019`. + +If using COMPAS's integration over cosmic star formation history, please cite :cite:t:`Neijssel2019`. + +If using the COMPAS model of (pulsational) pair instability supernova, please cite :cite:t:`Stevenson2019`. + +If evolving pulsar spins and magnetic fields with COMPAS, please cite :cite:t:`Chattopadhyay2020`. + +If you use the COMPAS model of chemically homogeneous evolution, please cite :cite:t:`Riley2021`. + diff --git a/online-docs/pages/quick-links.rst b/online-docs/pages/quick-links.rst index 371a13d1e..11e5e431a 100644 --- a/online-docs/pages/quick-links.rst +++ b/online-docs/pages/quick-links.rst @@ -6,4 +6,6 @@ Quick Links ./User guide/Program options/program-options-list-defaults ./User guide/COMPAS output/standard-logfiles-record-types + ./User guide/COMPAS output/standard-logfiles-record-specification-stellar.html + ./User guide/COMPAS output/standard-logfiles-record-specification-binary.html ./how-to-cite diff --git a/src/EAGB.cpp b/src/EAGB.cpp index e4845aa5f..97ddc6eeb 100755 --- a/src/EAGB.cpp +++ b/src/EAGB.cpp @@ -1078,14 +1078,12 @@ STELLAR_TYPE EAGB::ResolveEnvelopeLoss(bool p_Force) { */ bool EAGB::ShouldSkipPhase() const { #define timescales(x) m_Timescales[static_cast(TIMESCALE::x)] // for convenience and readability - undefined at end of function -#define gbParams(x) m_GBParams[static_cast(GBP::x)] // for convenience and readability - undefined at end of function double McCOBAGB = CalculateCOCoreMassOnPhase(timescales(tHeI) + timescales(tHe)); - double McSN = std::max(gbParams(McSN), 1.05 * McCOBAGB); // hack from Hurley fortran code, doesn't seem to be in the paper JR: do we know why? **Ilya** + double McSN = std::max(CalculateCoreMassAtSupernova_Static(MCH, m_GBParams[static_cast(GBP::McBAGB)]), 1.05 * McCOBAGB); // hack from Hurley fortran code, doesn't seem to be in the paper return (utils::Compare(McSN, m_COCoreMass) < 0); // skip phase if core is heavy enough to go supernova -#undef gbParams #undef timescales } @@ -1121,14 +1119,12 @@ bool EAGB::ShouldEvolveOnPhase() const { */ bool EAGB::IsSupernova() const { #define timescales(x) m_Timescales[static_cast(TIMESCALE::x)] // for convenience and readability - undefined at end of function -#define gbParams(x) m_GBParams[static_cast(GBP::x)] // for convenience and readability - undefined at end of function double McCOBAGB = CalculateCOCoreMassOnPhase(timescales(tHeI) + timescales(tHe)); - double McSN = std::max(gbParams(McSN), 1.05 * McCOBAGB); // hack from Hurley fortran code, doesn't seem to be in the paper JR: do we know why? **Ilya** + double McSN = std::max(CalculateCoreMassAtSupernova_Static(MCH, m_GBParams[static_cast(GBP::McBAGB)]), 1.05 * McCOBAGB); // hack from Hurley fortran code, doesn't seem to be in the paper JR: do we know why? **Ilya** return (utils::Compare(McSN, m_COCoreMass) <= 0); // core is heavy enough to go Supernova -#undef gbParams #undef timescales } diff --git a/src/GiantBranch.cpp b/src/GiantBranch.cpp index abe010436..2f057549f 100644 --- a/src/GiantBranch.cpp +++ b/src/GiantBranch.cpp @@ -276,9 +276,6 @@ void GiantBranch::CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams) gbParams(McBAGB) = CalculateCoreMassAtBAGB(p_Mass); gbParams(McDU) = CalculateCoreMassAt2ndDredgeUp_Static(gbParams(McBAGB)); gbParams(McBGB) = CalculateCoreMassAtBGB(p_Mass, p_GBParams); - - gbParams(McSN) = CalculateCoreMassAtSupernova_Static(gbParams(McBAGB)); - #undef gbParams } @@ -341,8 +338,6 @@ void GiantBranch::CalculateGBParams_Static(const double p_Mass, gbParams(McBAGB) = CalculateCoreMassAtBAGB_Static(p_Mass, p_BnCoefficients); gbParams(McBGB) = CalculateCoreMassAtBGB_Static(p_Mass, p_MassCutoffs, p_AnCoefficients, p_GBParams); - gbParams(McSN) = CalculateCoreMassAtSupernova_Static(gbParams(McBAGB)); - #undef gbParams } @@ -813,16 +808,17 @@ double GiantBranch::CalculateCoreMassAtBGB_Static(const double p_Mass, /* * Calculate the core mass at which the Asymptotic Giant Branch phase is terminated in a SN/loss of envelope * - * Hurley et al. 2000, eq 75 -- but note we use MECS rather than MCH + * Hurley et al. 2000, eq 75 -- but note we may use MECS rather than MCH as the CO core mass threshold * * - * double CalculateCoreMassAtSupernova_Static(const double p_McBAGB) + * double CalculateCoreMassAtSupernova_Static(const double p_Mthreshold, const double p_McBAGB) * + * @param [IN] p_Mthreshold Threshold mass, typically either MECS for stars that may explode in ECSNe or MCH otherwise * @param [IN] p_McBAGB Core mass at the Base of the Asymptotic Giant Branch in Msol * @return Maximum core mass before supernova on the Asymptotic Giant Branch */ -double GiantBranch::CalculateCoreMassAtSupernova_Static(const double p_McBAGB) { - return std::max(MECS, (0.773 * p_McBAGB) - 0.35); +double GiantBranch::CalculateCoreMassAtSupernova_Static(const double p_Mthreshold, const double p_McBAGB) { + return std::max(p_Mthreshold, (0.773 * p_McBAGB) - 0.35); } diff --git a/src/GiantBranch.h b/src/GiantBranch.h index 230c49eb1..640b13bab 100755 --- a/src/GiantBranch.h +++ b/src/GiantBranch.h @@ -36,7 +36,7 @@ class GiantBranch: virtual public BaseStar, public MainSequence { double CalculateCoreMassAtBGB(const double p_Mass, const DBL_VECTOR &p_GBParams); static double CalculateCoreMassAtBGB_Static(const double p_Mass, const DBL_VECTOR &p_MassCutoffs, const DBL_VECTOR &p_AnCoefficients, const DBL_VECTOR &p_GBParams); double CalculateCoreMassAtHeIgnition(const double p_Mass) const; - static double CalculateCoreMassAtSupernova_Static(const double p_McBAGB); + static double CalculateCoreMassAtSupernova_Static(const double p_Mthreshold, const double p_McBAGB); static double CalculateCoreMass_Luminosity_B_Static(const double p_Mass); static double CalculateCoreMass_Luminosity_D_Static(const double p_Mass, const double p_LogMetallicityXi, const DBL_VECTOR &p_MassCutoffs); diff --git a/src/HeHG.cpp b/src/HeHG.cpp index ef4102ef8..3a6235eae 100755 --- a/src/HeHG.cpp +++ b/src/HeHG.cpp @@ -65,8 +65,6 @@ void HeHG::CalculateGBParams(const double p_Mass, DBL_VECTOR &p_GBParams) { gbParams(McBAGB) = CalculateCoreMassAtBAGB(); gbParams(McBGB) = CalculateCoreMassAtBGB(p_Mass, p_GBParams); - gbParams(McSN) = CalculateCoreMassAtSupernova_Static(gbParams(McBAGB)); - #undef gbParams } @@ -134,8 +132,6 @@ void HeHG::CalculateGBParams_Static(const double p_Mass0, gbParams(McBAGB) = p_Mass0; gbParams(McBGB) = GiantBranch::CalculateCoreMassAtBGB_Static(p_Mass, p_MassCutoffs, p_AnCoefficients, p_GBParams); - gbParams(McSN) = CalculateCoreMassAtSupernova_Static(gbParams(McBAGB)); - #undef gbParams } @@ -412,12 +408,10 @@ double HeHG::ChooseTimestep(const double p_Time) const { * @return Boolean flag: true if evolution on this phase should continue, false if not */ bool HeHG::ShouldEvolveOnPhase() const { -#define gbParams(x) m_GBParams[static_cast(GBP::x)] // for convenience and readability - undefined at end of function double McMax = CalculateMaximumCoreMass(m_Mass); - return ((utils::Compare(m_COCoreMass, McMax) <= 0 || utils::Compare(McMax, gbParams(McSN)) >= 0) && !ShouldEnvelopeBeExpelledByPulsations()); // Evolve on HeHG phase if McCO <= McMax or McMax >= McSN and envelope is not ejected by pulsations - -#undef gbParams + double McSN = CalculateCoreMassAtSupernova_Static(MECS, m_GBParams[static_cast(GBP::McBAGB)]); + return ((utils::Compare(m_COCoreMass, McMax) <= 0 || utils::Compare(McMax, McSN) >= 0) && !ShouldEnvelopeBeExpelledByPulsations()); // Evolve on HeHG phase if McCO <= McMax or McMax >= McSN and envelope is not ejected by pulsations } @@ -484,7 +478,7 @@ bool HeHG::IsSupernova() const { return (utils::Compare(m_Mass, MECS) > 0); } - return (utils::Compare(m_COCoreMass, CalculateCoreMassAtSupernova_Static(m_GBParams[static_cast(GBP::McBAGB)])) >= 0); // Go supernova if CO core mass large enough + return (utils::Compare(m_COCoreMass, CalculateCoreMassAtSupernova_Static(MECS, m_GBParams[static_cast(GBP::McBAGB)])) >= 0); // Go supernova if CO core mass large enough } /* diff --git a/src/TPAGB.cpp b/src/TPAGB.cpp index f3fde9bca..a322b0332 100755 --- a/src/TPAGB.cpp +++ b/src/TPAGB.cpp @@ -978,8 +978,11 @@ double TPAGB::ChooseTimestep(const double p_Time) const { * @return Boolean flag: true if star has gone Supernova, false if not */ bool TPAGB::IsSupernova() const { - // no supernova if CO core mass is too low or helium core mass is too low at base of AGB or the envelope has already been removed - return utils::Compare(m_COCoreMass, m_GBParams[static_cast(GBP::McSN)]) >= 0 && - utils::Compare(CalculateInitialSupernovaMass(), OPTIONS->MCBUR1()) >= 0 && - utils::Compare(m_COCoreMass, m_Mass) < 0; + double snMass = CalculateInitialSupernovaMass(); + bool isCCSN = utils::Compare(m_COCoreMass, CalculateCoreMassAtSupernova_Static(MCH, m_GBParams[static_cast(GBP::McBAGB)])) >= 0 && + utils::Compare(snMass, OPTIONS->MCBUR1()) >= 0 && utils::Compare(m_COCoreMass, m_Mass) < 0; + bool isECSN = utils::Compare(snMass, MCBUR2) < 0 && (!m_MassTransferDonorHistory.empty() || OPTIONS->AllowNonStrippedECSN()) && + utils::Compare(m_COCoreMass, CalculateCoreMassAtSupernova_Static(MECS, m_GBParams[static_cast(GBP::McBAGB)])) >= 0 && + utils::Compare(snMass, OPTIONS->MCBUR1()) >= 0 && utils::Compare(m_COCoreMass, m_Mass) < 0; + return isCCSN || isECSN; } diff --git a/src/TPAGB.h b/src/TPAGB.h index d979ff9fc..5906770b5 100755 --- a/src/TPAGB.h +++ b/src/TPAGB.h @@ -48,7 +48,7 @@ class TPAGB: virtual public BaseStar, public EAGB { // member functions - alphabetically DBL_DBL CalculateConvectiveEnvelopeMass() const { return std::tuple (m_Mass-m_CoreMass, m_Mass-m_CoreMass); } // assume entire envelope is convective for TPAGB stars - double CalculateCOCoreMassAtPhaseEnd() const { return (utils::Compare(m_COCoreMass, m_GBParams[static_cast(GBP::McSN)]) >= 0 && utils::Compare(m_COCoreMass, m_Mass) < 0) ? m_COCoreMass : m_Mass; } + double CalculateCOCoreMassAtPhaseEnd() const { return (utils::Compare(m_COCoreMass, m_Mass) < 0) ? m_COCoreMass : m_Mass; } double CalculateCOCoreMassOnPhase() const { return CalculateCoreMassOnPhase(m_Mass0, m_Age); } // McCO(TPAGB) = Mc(TPAGB)Same as on phase double CalculateConvectiveCoreRadius() const { return std::min(5.0 * CalculateRemnantRadius(), m_Radius); } // Last paragraph of section 6 of Hurley+ 2000 @@ -103,7 +103,7 @@ class TPAGB: virtual public BaseStar, public EAGB { void ResolveHeliumFlash() { } // NO-OP STELLAR_TYPE ResolveSkippedPhase() { return m_StellarType; } // NO-OP - bool ShouldEvolveOnPhase() const { return ((utils::Compare(m_COCoreMass, std::min(m_GBParams[static_cast(GBP::McSN)], m_Mass)) < 0) && !ShouldEnvelopeBeExpelledByPulsations()); } // Evolve on TPAGB phase if envelope is not lost and not going supernova + bool ShouldEvolveOnPhase() const { return (utils::Compare(m_COCoreMass, m_Mass) < 0 && !IsSupernova() && !ShouldEnvelopeBeExpelledByPulsations()); } // Evolve on TPAGB phase if envelope is not lost and not going supernova bool ShouldSkipPhase() const { return false; } // Never skip TPAGB phase }; diff --git a/src/changelog.h b/src/changelog.h index dcce42312..c978e8127 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -944,10 +944,10 @@ // - Cleaned up stability check functions in BaseBinaryStar.cpp for clarity, and to allow for critical mass ratios to be checked correctly // 02.33.01 RTW - Sep 26, 2022 - Defect repair: // - Fixed interpolation of MACLEOD_LINEAR gamma for specific angular momentum. Previously interpolated on the gamma value, now interpolates in orbital separation -// 02.33.02 IM - Nov 27, 2022 - Defect repair: +// 02.33.02 IM - Nov 27, 2022 - Defect repair: // - Fixed ignored value of input radius when computing the thermal timescale, relevant if using Roche lobe radius instead (issue #853) // - Cleaned code and comments around the use of MT_THERMALLY_LIMITED_VARIATION::RADIUS_TO_ROCHELOBE vs. C_FACTOR (issue #850) -// 02.34.00 IM - Nov 28, 2022 - Enhancement: +// 02.34.00 IM - Nov 28, 2022 - Enhancement: // - Adding framework for Hirai & Mandel 2-stage common envelope formalism // (placeholders for now -- will have identical results to default version) // - Placed Dewi CE prescription on parity with others @@ -1579,9 +1579,12 @@ // 03.20.02 IM - May 30, 2026 - Defect repair, enhancement: // - Included unit conversion in WhiteDwarfs::CalculateEtaPTY() // - All critical mass ratios now return the HURLEY_HJELLMING_WEBBINK_QCRIT_WD for white dwarfs and 0 for other remnant donors (only stable mass transfer) as fix for issue #1385 +// 03.20.03 IM - June 18, 2026 - Defect repair, enhancement: +// - TPAGB stars should no longer experience supernovae if SN conditions are not satisfied, rather than defaulting to CCSN (corrects the partial fix in 03.10.02) +// - Added new parameter (threshold mass, generally expected to be MCH or MECS) to CalculateCoreMassAtSupernova_Static() +// - Removed McSN from GBParams, instead computed on the fly when needed - -const std::string VERSION_STRING = "03.20.02"; +const std::string VERSION_STRING = "03.20.03"; # endif // __changelog_h__ diff --git a/src/typedefs.h b/src/typedefs.h index 1180fd2a1..5dc107c64 100755 --- a/src/typedefs.h +++ b/src/typedefs.h @@ -522,7 +522,6 @@ enum class GBP: int { McBGB, // Core mass at BGB (Base of Giant Branch) McBAGB, // Core mass at BAGB (Base of Asymptotic Giant Branch). Hurley et al. 2000, eq 66 (also see eq 75 and discussion) McDU, // Core mass at second dredge up. Hurley et al. 2000, eq 69 - McSN, // Core mass at which the Asymptotic Giant Branch phase is terminated in a SN/loss of envelope COUNT // Sentinel for entry count }; From b48aa8d67ab480076655cdb207063350f5dc454a Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Wed, 18 Jun 2025 18:52:09 +1000 Subject: [PATCH 3/4] Whatsnew mod, changelog year fix --- online-docs/pages/whats-new.rst | 4 ++++ src/changelog.h | 4 ++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/online-docs/pages/whats-new.rst b/online-docs/pages/whats-new.rst index 684da3662..3ca033980 100644 --- a/online-docs/pages/whats-new.rst +++ b/online-docs/pages/whats-new.rst @@ -3,6 +3,10 @@ What's new Following is a brief list of important updates to the COMPAS code. A complete record of changes can be found in the file ``changelog.h``. +**03.20.03 June 18, 2025** + +* Resolved an issue that appeared in 03.10.02 with some TPAGB stars in the ECSN range (but not satisfying ECSN conditions) exploding as core-collapse supernovae + **03.20.01 May 26, 2025** * Updates to mass accretion for massive ONe WD diff --git a/src/changelog.h b/src/changelog.h index c978e8127..8940c8ab8 100644 --- a/src/changelog.h +++ b/src/changelog.h @@ -1576,10 +1576,10 @@ // - Update white dwarf mass-radius relation (WhiteDwarfs::CalculateRadiusOnPhase_Static) // - Moved white dwarf related constants to constants.h (resolves issue #1351) // - Set merger on unstable RLOF from WD -// 03.20.02 IM - May 30, 2026 - Defect repair, enhancement: +// 03.20.02 IM - May 30, 2025 - Defect repair, enhancement: // - Included unit conversion in WhiteDwarfs::CalculateEtaPTY() // - All critical mass ratios now return the HURLEY_HJELLMING_WEBBINK_QCRIT_WD for white dwarfs and 0 for other remnant donors (only stable mass transfer) as fix for issue #1385 -// 03.20.03 IM - June 18, 2026 - Defect repair, enhancement: +// 03.20.03 IM - June 18, 2025 - Defect repair, enhancement: // - TPAGB stars should no longer experience supernovae if SN conditions are not satisfied, rather than defaulting to CCSN (corrects the partial fix in 03.10.02) // - Added new parameter (threshold mass, generally expected to be MCH or MECS) to CalculateCoreMassAtSupernova_Static() // - Removed McSN from GBParams, instead computed on the fly when needed From 87806fa199285930663f45c585dfbc8e43d4e80b Mon Sep 17 00:00:00 2001 From: Ilya Mandel Date: Thu, 19 Jun 2025 08:42:35 +1000 Subject: [PATCH 4/4] Plot mod --- compas_matlab_utils/ComparisonPlots.m | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/compas_matlab_utils/ComparisonPlots.m b/compas_matlab_utils/ComparisonPlots.m index 30af1e03c..9a09af8f2 100644 --- a/compas_matlab_utils/ComparisonPlots.m +++ b/compas_matlab_utils/ComparisonPlots.m @@ -85,8 +85,8 @@ function ComparisonPlots(filename1, name1, filename2, name2) DWDplot(filename2, name2, 6, 'b', 20); end; figure(6); hold off; set(gca,'FontSize',20); title('Double White Dwarfs'); legend; - xlabel('$M*a [M_\odot * R_\odot]$ @ ZAMS', 'Interpreter', 'latex'); - ylabel('$M*a [M_\odot * R_\odot]$ @ end', 'Interpreter', 'latex'); + xlabel('$log_{10}(M*a) [M_\odot * R_\odot]$ @ ZAMS', 'Interpreter', 'latex'); + ylabel('$log_{10}(M*a) [M_\odot * R_\odot]$ @ end', 'Interpreter', 'latex'); [binariescount1, SNcount1, BHcompletecount1, SNbothcount1, SNonecount1, ... @@ -189,7 +189,7 @@ function ComparisonPlots(filename1, name1, filename2, name2) 'DisplayName', ['CE, ', name]); hold on; scatter(log10(P(BNS & ~isCE)), e(BNS & ~isCE), point, colour, 'DisplayName', ['Stable, ', name]); - %%% + %%% TO CLEAN good = mergingBBH & isCE & OKCE; type1CE = h5read(file,'/BSE_Common_Envelopes/Stellar_Type(1)