Merge pull request #1383 from TeamCOMPAS/BugFix

Ensure that Loveridge Lambdas are always positive; speed up Matlab post-processing

Author: jeffriley

compas_matlab_utils/ComparisonPlots.m

@@ -69,7 +69,7 @@ function ComparisonPlots(filename1, name1, filename2, name2)
     %Plot LMXBs/IMXBs
     figure(5); clf(5);
     [LMXBcount, NSLMXBcount]=LMXBplot(filename1, name1, 5, 'r', 40);
-    fprintf('\nLMXBs:\t\t#LMXB\t#NS LMXB\n');
+    fprintf('\nLMXBs:\t\t#LMXB\t\t#NS LMXB\n');
     fprintf('%s:\t%d\t\t%d\n', name1, LMXBcount, NSLMXBcount);
     if(nargin==4),
         [LMXBcount, NSLMXBcount]=LMXBplot(filename2, name2, 5, 'b', 20);
@@ -345,13 +345,18 @@ function DWDplot(file, name, fignumberDWD, colour, point)
     timeCE=h5read(file,'/BSE_Common_Envelopes/Time');
     seedAll=h5read(file, '/BSE_System_Parameters/SEED');
     [isCE,indexCE]=ismember(seedSN,seedCE); 
+    [isCE,lastindexCE]=ismember(seedSN,seedCE,'legacy'); %last index of CE seed matching given SN seed (for binaries with 2+ CE events)
     simultaneouswithCE = false(size(indexCE)); 
-    for(i=1:length(isCE)), if(isCE(i)), simultaneouswithCE(i) = ~isempty(find(timeCE==timeSN(i) & seedCE==seedSN(i))); end; end;
+    simultaneouswithCE(isCE) = (timeCE(indexCE(isCE)) == timeSN(isCE)) | (timeCE(lastindexCE(isCE)) == timeSN(isCE)); 
+    %Note: could (very rarely) miss a coincidence if >2 CE events and only intermediate CE matches SN in time
     precedingCE = false(size(indexCE)); 
-    for(i=1:length(isCE)), if(isCE(i)), precedingCE(i) = ~isempty(find(timeCE<timeSN(i) & seedCE==seedSN(i))); end; end;
+    precedingCE(isCE) = (timeCE(indexCE(isCE)) < timeSN(isCE)); 
     [experiencedRLOF,indexRLOF]=ismember(seedSN,seedRLOF);
+    [experiencedRLOF,lastindexRLOF]=ismember(seedSN,seedRLOF, 'legacy');
     simultaneouswithRLOF = false(size(indexRLOF)); 
-    for(i=1:length(indexRLOF)), if(experiencedRLOF(i)), simultaneouswithRLOF(i) = ~isempty(find(timeRLOF==timeSN(i) & seedRLOF==seedSN(i))); end; end;
+    simultaneouswithRLOF(experiencedRLOF) = (timeRLOF(indexRLOF(experiencedRLOF)) == timeSN(experiencedRLOF) ...
+                                            | timeRLOF(lastindexRLOF(experiencedRLOF)) == timeSN(experiencedRLOF)); 
+    %Note: could miss a coincidence if >2 RLOF events and only intermediate RLOF matches SN in time
     binariescount=length(seedAll);
     SNcount=length(seedSN);
     BHcompletecount=sum(starSNSN==14 & MpreSN==MSNSN);

src/HG.cpp

@@ -134,9 +134,10 @@ double HG::CalculateLambdaLoveridge(const double p_EnvMass, const bool p_IsMassL
     logBindingEnergy += 33.29866;                                                           // + logBE0
     double bindingEnergy = PPOW(10.0, logBindingEnergy);
 
-    return utils::Compare(bindingEnergy, 0.0) > 0 && utils::Compare(1.0 / bindingEnergy, 0.0) > 0 && utils::Compare(p_EnvMass, 0.0) > 0
+    double lambda = utils::Compare(bindingEnergy, 0.0) > 0 && utils::Compare(1.0 / bindingEnergy, 0.0) > 0 && utils::Compare(p_EnvMass, MAXIMUM_MASS_LOSS_FRACTION * m_Mass) > 0
             ? (G_CGS * m_Mass * MSOL_TO_G * p_EnvMass * MSOL_TO_G) / (m_Radius * RSOL_TO_AU * AU_TO_CM * bindingEnergy)
-            : 1.0;                                                                          // default to 1.0 (usual lambda default) if binding energy is not sensible [sometimes can be infinite if logBindingEnergy is too high] or if envelope mass is not positive [can be zero]
+            : 1.0;                                                                          // default to 1.0 (usual lambda default) if binding energy is not sensible [sometimes can be infinite if logBindingEnergy is too high] or if envelope mass is too low to reliably evaluate lambda [can be zero]
+    return utils::Compare(lambda, 0.0) > 0 ? lambda : 1.0;                                  // final check to avoid returning zero lambda
 }
 
 

src/changelog.h

@@ -1545,7 +1545,9 @@
 //                                      - Changed the default behaviour to use enhanced Nanjing lambdas (for common envelope calculations), interpolating in mass and metallicity
 //  03.18.03    IM - May 2, 2025    - Defect repair:
 //                                      - Fix for issue #1380, which appears when the Loveridge binding energy is so high that lambda is rounded off to zero
+//  03.18.04    IM - May 4, 2025    - Defect repair:
+//                                      - Added a check to avoid Loveridge lambda becoming zero when the envelope mass is positive but very small
 
-const std::string VERSION_STRING = "03.18.03";
+const std::string VERSION_STRING = "03.18.04";
 
 # endif // __changelog_h__