Clean up random draw

Author: Ilya Mandel

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, zlistmerger, 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
 %   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,7 +54,7 @@
 Mpc=Mpcm/c;         %Gpc in seconds
 yr=3.15569e7;       %year in seconds
 
-if (nargin<2)
+if (nargin<3)
     error('Not enough input arguments.');
 end;
 if (nargin<4), makeplots=0; end;
@@ -95,17 +93,38 @@
 end %end of CosmicHistoryIntegrator
 
 
-%Load the data stored in COMPAS output format from a file
+%Load the data stored in COMPAS .h5 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), 
+function [M1,M2,Z,Tdelay]=DataRead(file, tHubble)
+    if(exist(file, 'file')~=2), 
         error('Input file does not exist');
-    end;
-    data = importdata(filename, '\t', 2);
-    
-    colnames=data.textdata(2,:);
+    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');
+    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(mergingDCO)); 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);
+
     Z1index=find(strcmp(colnames,'Metallicity1'));
     M1index=find(strcmp(colnames,'M1'));
     M2index=find(strcmp(colnames,'M2'));
@@ -133,11 +152,11 @@
 
 %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;
+    %zmax=10; dz=0.001; zvec=0:dz:zmax;
     Nz=length(zvec);
 
     %Planck cosmology

src/BaseStar.cpp

@@ -3907,12 +3907,10 @@ double BaseStar::DrawSNKickMagnitude(const double p_Sigma,
             // 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 = DISBERG_MANDEL_MAX_KICK + 1.0;
-                double rand   = p_Rand;                                                         // makes it possible to adjust if p_Rand is too high, to avoid getting stuck
-                while (utils::Compare(kickMagnitude, DISBERG_MANDEL_MAX_KICK) > 0) {            // maximum kick of 1000 km/s, following Disberg & Mandel (2025)
-                    kickMagnitude = gsl_cdf_lognormal_Pinv(rand, DISBERG_MANDEL_MU, DISBERG_MANDEL_SIGMA);
-                    rand          = 0.99 * rand;
-                }
+                // 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(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);