Update changelog.h

Author: brcekadam

compas_matlab_utils/ComparisonPlots.m

@@ -160,7 +160,8 @@ function ComparisonPlots(filename1, name1, filename2, name2)
 end %end of ComparisonPlots
 
 
-%Plot double compact objects; returns DCO counts
+%Plot double compact objects; returns DCO counts,  and counts of CEs leading
+%to DCOs
 function [BNScount, NSBHcount, BBHcount, BNSCE, NSBHCE, BBHCE, CEBBH1count] = ...
         DCOplot(file, name, fignumber, colour, point)
     global Msunkg G AU

compas_matlab_utils/CosmicHistoryIntegrator.m

@@ -1,21 +1,19 @@
 function [Zlist, MergerRateByRedshiftByZ, SFRfractionZ]=...
-    CosmicHistoryIntegrator(filename,zlist,metallicitychoice, makeplots)
+    CosmicHistoryIntegrator(filename, zlistformation, zlistdetection, Msimulated, 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])
+%    CosmicHistoryIntegrator(filename, zlistformation, zlistmerger, [,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
+%   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)
@@ -56,89 +54,94 @@
 Mpc=Mpcm/c;         %Gpc in seconds
 yr=3.15569e7;       %year in seconds
 
-if (nargin<2)
+if (nargin<4)
     error('Not enough input arguments.');
 end;
-if (nargin<4), makeplots=0; end;
+if (nargin<5), makeplots=0; end;
 
 %cosmology calculator
-[zvec,tL]=Cosmology();           
+[tL]=Cosmology(zlistformation);           
 %load COMPAS data
-[M1,M2,Z,Tdelay]=DataRead(filename, max(tL));   
+[M1,M2,Z,Tdelay]=DataRead(filename); 
+Zlist=unique(Z);
 %metallicity-specific SFR
-[Zlist,SFR,SFRfractionZ]=Metallicity(Z,zvec,metallicitychoice);       
+[SFR,Zweight]=Metallicity(Zlist,zlistformation);       
 
 
 %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));
+dz=zlistformation(2)-zlistformation(1);
+tLmerge=tL(floor(zlistformation/dz)+1);
+MergerRateByRedshiftByZ=zeros(length(zlistformation),length(Zlist));
 for(i=1:length(M1)),
     Zcounter=find(Zlist==Z(i));
     zformindex=find(tL>=Tdelay(i),1);
-    for(k=1:length(zlist)),
+    for(k=1:length(zlistformation)),
         zformindex=find(tL>=(Tdelay(i)+tLmerge(k)),1);
         MergerRateByRedshiftByZ(k,Zcounter)=...
             MergerRateByRedshiftByZ(k,Zcounter)+...
-            SFR(zformindex)*SFRfractionZ(zformindex,Zcounter)/Msimulated;   
+            SFR(zformindex)*Zweight(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);
+    MakePlots(M1,M2,Z,Tdelay,zlistformation,Zlist,SFR,Zweight,...
+        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
+%Load the data stored in COMPAS .h5 output format from a file
+%Select only double compact object 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), 
+function [M1,M2,Z,Tdelay]=DataRead(file)
+    if(exist(file, '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<tHubble & nonRLOF & Pessimistic & ...
-        data.data(:,type1index)==14 & data.data(:,type2index)==14;
-
-    M1=data.data(select,M1index);
-    M2=data.data(select,M2index);
-    Z=data.data(select,Z1index);
-    Tdelay=Ttotal(select);
+    end;    
+    type1=h5read(file,'/BSE_Double_Compact_Objects/Stellar_Type(1)');
+    type2=h5read(file,'/BSE_Double_Compact_Objects/Stellar_Type(2)');
+    mass1=h5read(file,'/BSE_Double_Compact_Objects/Mass(1)');
+    mass2=h5read(file,'/BSE_Double_Compact_Objects/Mass(2)');
+    seedDCO=h5read(file,'/BSE_Double_Compact_Objects/SEED');
+    merges=h5read(file,'/BSE_Double_Compact_Objects/Merges_Hubble_Time');
+    a=h5read(file,'/BSE_Double_Compact_Objects/SemiMajorAxis@DCO');
+    e=h5read(file,'/BSE_Double_Compact_Objects/Eccentricity@DCO');
+    Ttotal=(h5read(file,'/BSE_Double_Compact_Objects/Time')+h5read(file,'/BSE_Double_Compact_Objects/Coalescence_Time'))*1e6; %to years
+    %mergingBBH=(type1==14) & (type2==14) & merges;
+    %BBH=(type1==14) & (type2==14);
+    %mergingBNS=(type1==13) & (type2==13) & merges;
+    %BNS=(type1==13) & (type2==13);
+    %mergingNSBH=(((type1==13) & (type2==14)) | ((type1==14) & (type2==13))) & merges;
+    %NSBH=(((type1==13) & (type2==14)) | ((type1==14) & (type2==13)));
+    %mergingDCO=mergingBNS | mergingNSBH | mergingBBH;
+    %BNScount=sum(mergingBNS); NSBHcount=sum(mergingNSBH); BBHcount=sum(mergingBBH);
+    chirpmass=mass1.^0.6.*mass2.^0.6./(mass1+mass2).^0.2;
+    q=mass2./mass1;
+    seedCE=h5read(file,'/BSE_Common_Envelopes/SEED');
+    [isCE,CEIndex]=ismember(seedDCO,seedCE);
+    optCE=h5read(file,'/BSE_Common_Envelopes/Optimistic_CE');
+    RLOFCE=h5read(file,'/BSE_Common_Envelopes/Immediate_RLOF>CE');
+    OKCE=zeros(size(seedDCO)); OKCE(CEIndex==0)=1; OKCE(CEIndex>0)=(~optCE(CEIndex(CEIndex>0))) & (~RLOFCE(CEIndex(CEIndex>0)));
+    %BNSCE=sum(mergingBNS & isCE & OKCE); NSBHCE=sum(mergingNSBH & isCE & OKCE); BBHCE=sum(mergingBBH & isCE & OKCE);
+    mergingDCO=merges & OKCE;
+    Zsys=h5read(file,'/BSE_System_Parameters/Metallicity@ZAMS(1)');
+    seedsys=h5read(file,'/BSE_System_Parameters/SEED');
+    [blah,sysIndex]=ismember(seedDCO,seedsys);
+    Zdco=Zsys(sysIndex);
+    M1=mass1(mergingDCO); M2=mass2(mergingDCO); Z=Zdco(mergingDCO); Tdelay=Ttotal(mergingDCO);
 end %end of DataRead
 
 %Compute the star formation rate and lookback time (in years) 
 %for an array of redshifts
-function [zvec,tL]=Cosmology()
+function [tL]=Cosmology(zvec)
     global Mpcm
     global Mpc
     global yr
-    zmax=10; dz=0.001; zvec=0:dz:zmax;
-    Nz=length(zvec);
+    %zmax=10; dz=0.001; zvec=0:dz:zmax;
+    Nz=length(zvec); dz=zvec(2)-zvec(1);
 
     %Planck cosmology
     OmegaM=0.236+0.046;  %2012arXiv1212.5226H
@@ -157,66 +160,72 @@
 end %end of Cosmology
 
 
-%Compute the fraction of star formation that happens in a given metallicity
-%bin as a function of redshift
-function [Zlist,SFR,SFRfractionZ]=Metallicity(zvec,Z)
+%Compute the weight of each star-forming metallicity as a function of redshift
+function [SFR,Zweight]=Metallicity(Zvec, zvec)
     %M_/odot per Mpc^3 per year -- Neijssel+ 2019 preferred model 
     %would be SFR=0.015*(1+zvec).^2.7./(1+((1+zvec)/2.9).^5.6) in Madau & Dickinson, 2014, (15)
     SFR=0.01*(1+zvec).^2.77./(1+((1+zvec)/2.9).^4.7); 
-    Zmean=0.035.*10.^(-0.23*zvec);
-    Zmu=log(Zmean)-0.39^2/2;
-    dlogZ=0.01;
-    logZvec=-12:dlogZ:0;  %natural log
-    dPdlogZ=1/0.39/sqrt(2*pi)*exp(-(logZvec'-Zmu).^2/2/0.39^2);
-    dPdlogZ=dPdlogZ./(sum(dPdlogZ,1)*dlogZ);    %normalise
-    Zrange=log(max(Z))-log(min(Z));   %ugly correction for not including tails
-    PdrawZ=1/Zrange;
-    minlogZindex=find(exp(logZvec)>=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
+    if(length(Zvec)>1),
+        Zmean=0.035.*10.^(-0.23*zvec);
+        Zmu=log(Zmean)-0.39^2/2;
+        dlogZ=0.01;
+        logZvec=-12:dlogZ:0;  %natural log
+        dPdlogZ=1/0.39/sqrt(2*pi)*exp(-(logZvec'-Zmu).^2/2/0.39^2);
+        dPdlogZ=dPdlogZ./(sum(dPdlogZ,1)*dlogZ);    %normalise
+        Zrange=log(max(Zvec))-log(min(Zvec));   %ugly correction for not including tails
+        PdrawZ=1/Zrange;
+        minlogZindex=find(exp(logZvec)>=min(Zvec),1, 'first');
+        maxlogZindex=find(exp(logZvec)>=max(Zvec),1, 'last');
+        dPdlogZ(minlogZindex,:)=dPdlogZ(minlogZindex,:)+sum(dPdlogZ(1:minlogZindex,:),1)*dlogZ/(sum(Zvec==min(Zvec))/length(Zvec))*PdrawZ;
+        dPdlogZ(maxlogZindex,:)=dPdlogZ(maxlogZindex,:)+sum(dPdlogZ(maxlogZindex:end,:),1)*dlogZ/(sum(Zvec==max(Zvec))/length(Zvec))*PdrawZ;
+        dPdlogZ(1:minlogZindex,:)=0; dPdlogZ(maxlogZindex:size(dPdlogZ,1),:)=0;
+        dPdlogZ=dPdlogZ./(sum(dPdlogZ,1)*dlogZ);    %normalise
+        for(i=1:length(Zvec))
+            index=find(exp(logZvec)>=Zvec(i), 1, 'first');
+            Zweight(:,i)=dPdlogZ(index,:)*dlogZ;
+        end;
+    else    %relevant for single-metallicity runs -- just give all binaries the same unit weight
+        Zweight(1:length(zvec),1:length(Zvec))=1;
+    end;
 end %end of Metallicity
 
 
 %Make a set of default plots
-function MakePlots(M1,M2,Z,Tdelay,zvec,SFR,SFRfractionZ,...
-    zlist,Zlist,MergerRateByRedshiftByZ)
+function MakePlots(M1,M2,Z,Tdelay,zvec,Zlist,SFR,Zweight,...
+        MergerRateByRedshiftByZ)
 
     figure(1),colormap jet;
-    plot(zlist, MergerRateByRedshiftByZ, 'LineWidth', 2), 
+    plot(zvec, MergerRateByRedshiftByZ*1e9, '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: ', ...
+    ylabel('DCO merger rate per Gpc^3 per yr')
+    disp(['Total DCO 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]');
+    xlabel('Primary mass [M_o]'), ylabel('Secondary mass [M_o]');
 
     figure(3), colormap jet;
-    scatter(M1+M2,log(Tdelay/1e6)/log(10),20,log(Z)/log(10),'filled');
+    scatter(M1+M2,log10(Tdelay/1e6),20,log10(Z),'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)');
+    xlabel('Total DCO mass [M_o]'), ylabel('log_{10}(Tdelay/Myr)');
 
     figure(4), colormap jet;
-    plot(zvec, SFR)
+    plot(zvec, SFR*1e9)
     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)
+    plot(zvec, Zweight)
     set(gca, 'FontSize', 20); %for labels
     legend(num2str(Zlist))
-    xlabel('z'), ylabel('Z-specific SFR fraction');
+    xlabel('z'), ylabel('Z-specific SFR weight');
 
 end %end of MakePlots
 

online-docs/pages/quick-links.rst

@@ -6,6 +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
+   ./User guide/COMPAS output/standard-logfiles-record-specification-stellar
+   ./User guide/COMPAS output/standard-logfiles-record-specification-binary
    ./how-to-cite

src/BaseStar.cpp

@@ -3904,10 +3904,14 @@ double BaseStar::DrawSNKickMagnitude(const double p_Sigma,
 
         case KICK_MAGNITUDE_DISTRIBUTION::LOGNORMAL: {                                          // LOGNORMAL
             // Only draw Disberg & Mandel (2025) kicks for CCSN or PPISN (if they receive a kick)
-	    // use Maxwellians with sigma set in CalculateSNKickMagnitude() for other SN types
+            // use Maxwellians with sigma set in CalculateSNKickMagnitude() for other SN types
             SN_EVENT thisSNevent = utils::SNEventType(m_SupernovaDetails.events.current);       // current SN event
-            if (thisSNevent == SN_EVENT::CCSN || (thisSNevent == SN_EVENT::PPISN && OPTIONS->NatalKickForPPISN()))
-                kickMagnitude = gsl_cdf_lognormal_Pinv(p_Rand, DISBERG_MANDEL_MU, DISBERG_MANDEL_SIGMA);
+            if (thisSNevent == SN_EVENT::CCSN || (thisSNevent == SN_EVENT::PPISN && OPTIONS->NatalKickForPPISN())) {
+                // maximum kick of DISBERG_MANDEL_MAX_KICK = 1000 km/s, following Disberg & Mandel (2025)
+                // normalise the draw so that the kick is uniformly drawn from the inverse CDF below this maximum
+                double cdfMax = gsl_cdf_lognormal_P(DISBERG_MANDEL_MAX_KICK, DISBERG_MANDEL_MU, DISBERG_MANDEL_SIGMA);
+                kickMagnitude = gsl_cdf_lognormal_Pinv(p_Rand*cdfMax, DISBERG_MANDEL_MU, DISBERG_MANDEL_SIGMA);
+            }
             else
                 kickMagnitude = DrawKickMagnitudeDistributionMaxwell(p_Sigma, p_Rand);
             } break;

src/changelog.h

@@ -1591,11 +1591,14 @@
 //  03.20.06   IM - June 25, 2025       - Enhancements:
 //                                          - The MAXWELLIAN NS CCSN kick changed from the Hobbs value of 265 km/s to 217 km/s based on 48 younger than 10 Myr pulsars with proper motions from Disberg & Mandel (2025) sample; corrects Hobbs+ 2005 missing Jacobian
 //                                          - Implemented a LOGNORMAL NS CCSN kick magnitude distribution based on Disberg & Mandel, 2025
-//  03.20.07   AB - Jun 25, 2025        - Defect repair:
+//  03.20.07   IM - June 25, 2025       - Enhancement:
+//                                          - Added a maximum threshold of 1000 km/s for Disberg & Mandel (2025) LOGNORMAL kicks, matching paper
+//  03.20.08   AB - Jun 26, 2025        - Defect repair:
 //                                          - Fix for issue #400; correct Zsol values are now used in stellar wind prescriptions
 //                                          - To avoid ambiguous ZSOL, we now use ZSOL_HURLEY = 0.02, ZSOL_ANDERS = 0.019, and ZSOL_ASPLUND = 0.0142
+//                                          - Fixed error in MainSequence::CalculateInitialMainSequenceCoreMass()
 
-const std::string VERSION_STRING = "03.20.07";
+const std::string VERSION_STRING = "03.20.08";
 
 
 # endif // __changelog_h__

src/constants.h

@@ -397,6 +397,7 @@ constexpr double MULLERMANDEL_SIGMAKICK                 = 0.3;
 // Constants for Disberg & Mandel (2025) SN kick prescription
 constexpr double DISBERG_MANDEL_MU                      = 5.60;
 constexpr double DISBERG_MANDEL_SIGMA                   = 0.68;
+constexpr double DISBERG_MANDEL_MAX_KICK                = 1000.0;
 
 constexpr double HOBBS_CORRECTED_SIGMA                  = 217.0;    // Best fit Maxwellian sigma for 48 younger than 10 Myr pulsars with proper motions from Disberg & Mandel (2025) sample; corrects Hobbs+ 2005 missing Jacobian