From 136d6ea1cc8bb65348b263ca3d7a82ef8f4f3be2 Mon Sep 17 00:00:00 2001 From: TristanBruel Date: Wed, 17 Jun 2026 17:09:29 +0200 Subject: [PATCH 1/4] Implement two-stage CE evolution from Hirai & Mandel 2022 --- examples/Params.ini | 9 ++ src/cosmic/consts.py | 4 +- src/cosmic/data/cosmic-settings.json | 17 +++ src/cosmic/src/comenv.f | 124 +++++++++++++++++- src/cosmic/src/const_bse.h | 2 + src/cosmic/src/test_bse.f | 1 + src/cosmic/tests/data/Params.ini | 5 + .../data/initial_conditions_for_testing.hdf5 | Bin 3193691 -> 1099583 bytes .../tests/data/kick_initial_conditions.h5 | Bin 3226940 -> 1132192 bytes 9 files changed, 154 insertions(+), 8 deletions(-) diff --git a/examples/Params.ini b/examples/Params.ini index 16f178042..1bfb135da 100644 --- a/examples/Params.ini +++ b/examples/Params.ini @@ -419,6 +419,15 @@ alpha1 = [1.0, 1.0] ; Default: 0.0 lambdaf = 0.0 +; ce2stageflag +; Selects whether to use the two-stage CE formalism from Hirai & Mandel 2022 +; +; Options: +; 0 - Use the standard energy formalism +; 1 - Use the two-stage CE formalism from Hirai & Mandel 2022 +; Default: 0 +ce2stageflag = 0 + ; ceflag ; Selects the model to set the initial orbital energy (choose whether to use total mass of the stars instead of the core masses) ; diff --git a/src/cosmic/consts.py b/src/cosmic/consts.py index f3d94c627..d2967c9c8 100644 --- a/src/cosmic/consts.py +++ b/src/cosmic/consts.py @@ -42,7 +42,7 @@ "epsnov", "eddfac", "gamma", "LBV_flag", ], "cevars": ["alpha1", "lambdaf", "qcrit_array"], - "ceflags": ["ceflag", "cekickflag", "cemergeflag", "cehestarflag", "ussn"], + "ceflags": ["ce2stageflag", "ceflag", "cekickflag", "cemergeflag", "cehestarflag", "ussn"], "flags": [ "tflag", "ifflag", "wdflag", "rtmsflag", "bhflag", "remnantflag", "maltsev_mode", "grflag", "bhms_coll_flag", "wd_mass_lim", "aic", @@ -63,4 +63,4 @@ "mixvars": ["rejuv_fac"], "metvars": ["zsun"], "rand1": ["randomseed"] -} \ No newline at end of file +} diff --git a/src/cosmic/data/cosmic-settings.json b/src/cosmic/data/cosmic-settings.json index 6f7936db8..a60cbaccb 100644 --- a/src/cosmic/data/cosmic-settings.json +++ b/src/cosmic/data/cosmic-settings.json @@ -865,6 +865,23 @@ } ] }, + { + "name": "ce2stageflag", + "description": "Sets whether to use the standard energy formalism or the two-stage CE formalism from Hirai & Mandel 2022", + "type": "dropdown", + "options-preface": "", + "options": [ + { + "name": 0, + "description": "Use the standard energy formalism", + "default": true + }, + { + "name": 1, + "description": "Use the two-stage CE formalism from Hirai & Mandel 2022" + } + ] + }, { "name": "ceflag", "description": "Selects the model to set the initial orbital energy (choose whether to use total mass of the stars instead of the core masses)", diff --git a/src/cosmic/src/comenv.f b/src/cosmic/src/comenv.f index 849b0fe06..d0a57c52b 100644 --- a/src/cosmic/src/comenv.f +++ b/src/cosmic/src/comenv.f @@ -75,6 +75,13 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, REAL*8 rad(2),tms(2),lumin(2),B_0(2),bacc(2),tacc(2),epoch(2) REAL*8 menv_bpp(2),renv_bpp(2) REAL*8 ALPHA_CE + + REAL*8 met,reagb,tmin,tset,mconvmax + REAL*8 ragbf,menvmax,mconvenv,CELAMHE,tonset + EXTERNAL ragbf,menvmax,mconvenv,CELAMHE,tonset + REAL*8 mcgbtf,lmcgbf + EXTERNAL mcgbtf,lmcgbf + REAL*8 m1endstage1,m2endstage1 * * Initialize * @@ -120,6 +127,44 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, RZAMS = RZAMSF(M01) LAMB1 = CELAMF(KW,M01,L1,R1,RZAMS,MENVD,LAMBDAF) ENDIF +* +* Two-stage energy formalism +* + IF(CE2STAGEFLAG.EQ.1)THEN + IF (using_SSE.eq.1) THEN +* Estimate the mass of the convective envelope +* from the fits in Picker, Hirai, & Mandel 2024 + met = 10**(LOG10(zpars(14))/0.4) + reagb = ragbf(M1,LUMS(7),zpars(2)) + tmin = 1000.d0*((1130.d0*LUMS(7)/ + & (reagb**2.d0))**(1.d0/4.d0)) + tset = tonset(tmin,met) + teff1 = 1000.d0*((1130.d0*lumin(1)/ + & (rad1_bpp**2.d0))**(1.d0/4.d0)) + MENV = mconvenv(KW,M1,met,teff1,tmin,tset,AJ1,TM1) + mconvmax = menvmax(KW,M1,met) + LAMB1 = CELAMHE(M1,met,MENV,mconvmax) + ENDIF +* if > 8 Msun, the envelope participating in the first stage is only the convective one +* if < 2 Msun, the entire envelope participates in the first stage +* linear interpolation in between + IF(M1.GE.8.0d0)THEN + m1endstage1 = M1-MENV + ELSEIF(M1.LT.2.0d0)THEN + m1endstage1 = MC1 + ELSE + m1endstage1 = MC1 + (M1-MC1 - MENV) * (M1-2.0d0) /6.0d0 + ENDIF +* +* Calculate the binding energy of the convective envelope (multiplied by lambda). +* + EBINDI = M1*(M1-m1endstage1)/(LAMB1*R1) + ELSE +* +* Calculate the binding energy of the giant envelope (multiplied by lambda). +* + EBINDI = M1*(M1-MC1)/(LAMB1*R1) + ENDIF KW = KW2 CALL star(KW2,M02,M2,TM2,TN,TSCLS2,LUMS,GB,ZPARS,dtm,star2) CALL hrdiag(M02,AJ2,M2,TM2,TN,TSCLS2,LUMS,GB,ZPARS, @@ -127,10 +172,6 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, & bhspin2,star2) OSPIN2 = JSPIN2/(K22*R2*R2*(M2-MC2)+K3*RC2*RC2*MC2) * -* Calculate the binding energy of the giant envelope (multiplied by lambda). -* - EBINDI = M1*(M1-MC1)/(LAMB1*R1) -* * If the secondary star is also giant-like add its envelopes energy. * Determine EORBI based on CEFLAG (CEFLAG=1 for de Kool prescription) * @@ -142,7 +183,44 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, RZAMS = RZAMSF(M02) LAMB2 = CELAMF(KW,M02,L2,R2,RZAMS,MENVD,LAMBDAF) ENDIF - EBINDI = EBINDI + M2*(M2-MC2)/(LAMB2*R2) +* +* Two-stage energy formalism +* + IF(CE2STAGEFLAG.EQ.1)THEN + IF (using_SSE.eq.1) THEN +* Estimate the mass of the convective envelope +* from the fits in Picker, Hirai, & Mandel 2024 + met = 10**(LOG10(zpars(14))/0.4) + reagb = ragbf(M2,LUMS(7),zpars(2)) + tmin = 1000.d0*((1130.d0*LUMS(7)/ + & (reagb**2.d0))**(1.d0/4.d0)) + tset = tonset(tmin,met) + teff2 = 1000.d0*((1130.d0*lumin(2)/ + & (rad2_bpp**2.d0))**(1.d0/4.d0)) + MENV = mconvenv(KW,M2,met,teff2,tmin,tset,AJ2,TM2) + mconvmax = menvmax(KW,M2,met) + LAMB2 = CELAMHE(M2,met,MENV,mconvmax) + ENDIF +* if > 8 Msun, the envelope participating in the first stage is only the convective one +* if < 2 Msun, the entire envelope participates in the first stage +* linear interpolation in between + IF(M2.GE.8.0d0)THEN + m2endstage1 = M2-MENV + ELSEIF(M1.LT.2.0d0)THEN + m2endstage1 = MC2 + ELSE + m2endstage1 = MC2 + (M2-MC2 - MENV) * (M2-2.0d0) /6.0d0 + ENDIF +* +* Calculate the binding energy of the convective envelope (multiplied by lambda). +* + EBINDI = EBINDI + M2*(M2-m2endstage1)/(LAMB2*R2) + ELSE +* +* Calculate the binding energy of the giant envelope (multiplied by lambda). +* + EBINDI = EBINDI + M2*(M2-MC2)/(LAMB2*R2) + ENDIF * * Calculate the initial orbital energy * @@ -167,7 +245,7 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, ENDIF EORBF = EORBI + EBINDI/ALPHA_CE * -* If the secondary is on the main sequence see if it fills its Roche lobe. +* If the secondary is on the main sequence. * IF(KW2.LE.1.OR.KW2.EQ.7)THEN SEPF = MC1*M2/(2.D0*EORBF) @@ -176,6 +254,40 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, RL1 = RL(Q1) RL2 = RL(Q2) * +* If degenerate or giant secondary. +* + ELSE + SEPF = MC1*MC2/(2.D0*EORBF) + Q1 = MC1/MC2 + Q2 = 1.D0/Q1 + RL1 = RL(Q1) + RL2 = RL(Q2) + ENDIF +* +* Two-stage energy formalism +* + IF(CE2STAGEFLAG.EQ.1)THEN + SEPF = m1endstage1*m2endstage1/(2.D0*EORBF) +* +* Second stage: stable mass transfer of the radiative intershell +* Binary hardening formula from Picker, Hirai & Mandel 2024 +* + IF(m1endstage1.GT.MC1)THEN + SEPF = SEPF*((m1endstage1+m2endstage1)/(MC1+m2endstage1)) + & *(m1endstage1/MC1)**2 + & *EXP(-2*(m1endstage1-MC1)/m2endstage1) + EORBF = MC1*m2endstage1/(2.D0*SEPF) + ENDIF + Q1 = MC1/m2endstage1 + Q2 = 1.D0/Q1 + RL1 = RL(Q1) + RL2 = RL(Q2) + ENDIF +* +* If the secondary is on the main sequence see if it fills its Roche lobe. +* + IF(KW2.LE.1.OR.KW2.EQ.7)THEN +* * If cemergeflag is set, cause kstars without clear core-envelope * structure to merge automatically if they enter a CE * diff --git a/src/cosmic/src/const_bse.h b/src/cosmic/src/const_bse.h index 3b3314f9d..022a8faed 100644 --- a/src/cosmic/src/const_bse.h +++ b/src/cosmic/src/const_bse.h @@ -18,7 +18,9 @@ REAL*8 don_lim,acc_lim(2),Mbh_initial,smt_periastron_check COMMON /MTVARS/ don_lim,acc_lim,Mbh_initial,smt_periastron_check + INTEGER ce2stageflag INTEGER ceflag,cekickflag,cemergeflag,cehestarflag,ussn + COMMON /CEFLAGS/ ce2stageflag COMMON /CEFLAGS/ ceflag,cekickflag,cemergeflag,cehestarflag,ussn INTEGER pisn_track(2) COMMON /TRACKERS/ pisn_track diff --git a/src/cosmic/src/test_bse.f b/src/cosmic/src/test_bse.f index 983d70cb0..47e2a6561 100644 --- a/src/cosmic/src/test_bse.f +++ b/src/cosmic/src/test_bse.f @@ -123,6 +123,7 @@ PROGRAM bse WRITE(*,*)mass0(1),mass0(2),tphysf,tb,kstar(1),kstar(2),z,ecc WRITE(*,*)neta,bwind,hewind,alpha1(1),alpha1(2),lambdaf,windflag + WRITE(*,*)ce2stageflag WRITE(*,*)rtmsflag,ceflag,tflag,ifflag,wdflag,bhflag,remnantflag WRITE(*,*)mxns,idum,pts1,pts2,pts3 WRITE(*,*)sigma,beta,xi,acc2,epsnov,eddfac,gamma diff --git a/src/cosmic/tests/data/Params.ini b/src/cosmic/tests/data/Params.ini index 9aba74431..8c0078ce0 100644 --- a/src/cosmic/tests/data/Params.ini +++ b/src/cosmic/tests/data/Params.ini @@ -194,6 +194,11 @@ alpha1 = [1.0, 1.0] ; default=0.5 lambdaf = 0.5 +; ce2stageflag=0 uses the standard energy formalism for CEE +; ce2stageflag=1 uses the two-stage formalism CEE from Hirai & Mandel 2022 +; default=0 +ce2stageflag=1 + ; ceflag=1 used the method from de Kool 1990 for setting the initial orbital energy ; ceflag=0 does not use this method (uses the core mass to calculate initial orbital energy) ; default=0 diff --git a/src/cosmic/tests/data/initial_conditions_for_testing.hdf5 b/src/cosmic/tests/data/initial_conditions_for_testing.hdf5 index eb582d7ab9eb034ee12dfd7c1ff78dd97758a275..9a5e289ff728c24858e29d9abda0c275d6c1abe6 100644 GIT binary patch delta 974 zcmYjPT}V?=9KUxT?(S}G_lnK!PPc0a4TQRGMi$Z(4TMPsQ(I`{5_1n_3xxy`CO()C zxn2G*_S7fXhmhApejqW4z$72)&CsC0&<8;!DD1_~xgX#;{P^=b|KI=Pod4PQ^MTG! zG~`0&=n5xS4x2iwl-CKfjiqRcq8>o$d=SvA7j`2e{GJiug2MxTwy^($cJzWvPdejLcj*jv-kTnR!4&yd5>qZ8v!AhUaVq=Ml_*lHiT-4 zpi522=vg_6Q)QMaO`U*Ds>D9ls7!^bc;>wHw`S}oIr zS5yN2X0jq^S7imPVRl@YdFdAS(e@YlotI0m6(c zO!OpC;QU2JVHQ!-439fk;>^R_Bz~vlvaj z2x|um6g*9xhAy*&Dym0I4Ili8htyOTpCL* k(HJHSGlm7jicx@J!(cG%7%YYZ!--Ld;hIk_x&Iyg2f=_1XaE2J delta 1852 zcmcIkYiJZ#6ux(6XP=XG*1WUb&7>7we848VrBT5~TwfU4!X`FFakGZVmc(Q+YafJC zC5lR+*~$4~u$64EU@Ai0PKznQ7=@BjDhdh;5k(aAhm`&jgo2(s_s;r2e+coyH)plUuAed{Cl|I4k7h?&^_X2oCsqVNqlA`zncA8z}MNXKwOEx_57PD))7^ zn}1iXp&RRK&IfXP2!V6%^=}$TvcMh8cOYG`GQ1G}R=m^?({&;5!9ISDpTPf04NKQ7 zgexg5cWrpnkcTGTM#z_(>mcjUq2SC`XaviP>z8(ZQU?n7(HS@0s!>*Egt+z-Cg4`> z6di|Sb!)wY`?>yx-%fEPcvN-Z?eVp_&M1bx+%q9}LX1U(?RpX!1m1}lGX?nYyrxbI z1Cx04wS{&}7wWl1qc9h`C<>XT(%s)c*4Ic1fyZq@*rE5p5xo*B+E>6;ecBw{ugd1} zMH^Mv-WZYoi4}`s3d0#UoNjVECfOPP#o~~co8aH3KKQwI6?6p`LoU|TQe*WlZ(M?P z;qeb=coZ#MkWK5?TIJ|ue*rG$lCAe9VU=SsoXkYbczmE5K5dOk+or&XR+(e14_q*+ z=r1+Om#=M<1e|fzz!hf;+;7t`2bL&}LrPz2rnXjTRDylc+S<>8{A4@gA56|ii^coY z9L_vYU~bs2XVbY%D%YRN?MZj_4A?A#|Cps9(4f{sHsChL`ZQ|Z>>mrjz2vOz~W5FxVdUf=Ddzia%fh9}wVie2}n9{4mA?-^4}fS(+-6 zWEI=qANBDqc@)nYF#RqmH6XVqKl(U9QD^(jyI3PO3;%y9jEUPN<{8F1ge~J??66>~ ze`o9~!M^*Iu?fM3{$OlWu=2CWShsmyTe+&k>QWkGmp9_kVN%)t;z%ef6{LQ}65&R) z7O}aV_>4&0O2`^EKZ#N|)8`5iu=lr$!bYI5<=f3Wd;7b3Vkv9D4&X_Gf1c~!%_!&b zm+bLW$!5hp+W-In diff --git a/src/cosmic/tests/data/kick_initial_conditions.h5 b/src/cosmic/tests/data/kick_initial_conditions.h5 index 0323919ebb57edf9a6d644ec4afeab56804d521b..10aeee830cb82d047efa768fb1fd854b2c29ecf5 100644 GIT binary patch delta 827 zcmY+BUr19?9LLY^ZoBuM>)h_D>(-LZOk#6kLro=h>!qbAp-))~B4H~ST58}!L2uPV z>FpPs2x6Cr@FB=?g7pw(6i5_3)luG}BYk>$~sJczH)_S{C#$72Trnrg)!viXA z7{{d{rRaRZn(EC zmQ!cT_hVDnB-SJ|lgT_d}VoFsy;d1JUug(9zT^COA|DoKg5+~J>Ns&YXC>vy!ka3^t!i3@%8(PUfRc| z&<->eKk9)xzM9t+bIGc{E9SJC?B~Me!R{iGz#(u7Tmo632vh+TXo4z%TTm_V2)y~h KZhdV-X#O8lx8NrL delta 1541 zcmbtUU1$_n6ux(UGP8F!$*k$-XGnI}6k)ScQM1OVF>AE0ilBi$xYpL)O%vMC)k=*e zqT)-zwjgFbfdoOjPnt)UdDtLIE+|AosrulP!IvV6E$D-YZK2yUGk2Cn-*lF5hwq*< z-~G;m__72kF>#AcDr{i zMhWH^M91F``ND~@?8M1qW9~9(F$s~~xW|uAd`}3R(e_Hi43@N`FxTmN6na)e$%?_~qc{v}lR;fL z&EU@Pg_^YIj7tHwtU9likW8?y)}P^JxQQxI4Ih-AriSme9598)U)&I4m;ji4YuSw<78g!rfjD0Z(tSZ zMnAM0aeJYs^bs_t)=C3g?V;3v8r+-KrMoIHrPgiQ1h=d%_DuA<>Fs{!oR`;tvoe14 zYj_;CV0lruoA*RaDitqubi4_BTek}br#E<(X3OQL)Yjo$FWxV6#p=W&jC>Hb^9*|=eC7hTP-BzTtUug~`F z|CHj|Zm}4H2dz=o@24ovth{X1<>~dbSngsO5qP6v!1zh9&W%pcm=v^-bm^R+qgz}$ zE9k3!hq}$n1+aM^(_TCzS8H>Jr?OVqriN%DA)XV{L458PoHM+H-Af>hiTpX z=Rf+#c&nZdJrMhf43^=)s~{kR;*=L8?07(3!nk0;-XC$qJvO<{kG_NDM@u+obOJqaRwte5^kKHE2kYu2pORuR73+pL+FSAf*}k<5YdQeLNr&( zY4b(3$Q4aULiBFeIgWxH1vv_G6yzw#QIMk`M?sE)90fTFaunn!$Wf4^AWxgi)y7M8 K1|3r$Zt#E3FM^i< From bf11ddd4f094843494b1850d2cc6aedc789c314d Mon Sep 17 00:00:00 2001 From: "github-actions[bot]" Date: Wed, 17 Jun 2026 15:14:07 +0000 Subject: [PATCH 2/4] update Params.ini --- examples/Params.ini | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/Params.ini b/examples/Params.ini index 1bfb135da..51081b9c0 100644 --- a/examples/Params.ini +++ b/examples/Params.ini @@ -420,7 +420,7 @@ alpha1 = [1.0, 1.0] lambdaf = 0.0 ; ce2stageflag -; Selects whether to use the two-stage CE formalism from Hirai & Mandel 2022 +; Sets whether to use the standard energy formalism or the two-stage CE formalism from Hirai & Mandel 2022 ; ; Options: ; 0 - Use the standard energy formalism From fe7fd22352b178eaa9ad1a2d8afc9a5f25f6719c Mon Sep 17 00:00:00 2001 From: TristanBruel Date: Wed, 17 Jun 2026 17:22:06 +0200 Subject: [PATCH 3/4] Add fitting functions for the convective envelope to SSE_zfuncs.f --- src/cosmic/src/SSE/SSE_zfuncs.f | 88 +++++++++++++++++++++++++++++++++ 1 file changed, 88 insertions(+) diff --git a/src/cosmic/src/SSE/SSE_zfuncs.f b/src/cosmic/src/SSE/SSE_zfuncs.f index 25b78578a..7bce39384 100644 --- a/src/cosmic/src/SSE/SSE_zfuncs.f +++ b/src/cosmic/src/SSE/SSE_zfuncs.f @@ -845,6 +845,94 @@ real*8 FUNCTION lmcgbf(mc,GB) * return end +*** + real*8 FUNCTION tonset(tmin,z) + implicit none + real*8 tmin,z +* +* A function to evaluate Tonset (Teff when Mconv reaches 1% of its max) +* from Mandel, Hirai & Picker 2024 +* + tonset = tmin /min(0.695d0-0.057d0*log10(z),0.95d0) +* + return + end +*** + real*8 FUNCTION menvmax(kw,m,z) + implicit none + integer kw + real*8 m,z + real*8 mcagbf,mifmcf + external mcagbf +* +* A function to evaluate the maximum mass of the convective envelope +* from Picker, Hirai & Mandel 2024 +* +* I'm not sure how HG stars should be treated ... +* + if(kw.ge.10)then + menvmax = 0.0d0 + elseif(kw.lt.2)then + if(m.gt.1.25d0)then + menvmax = 0.0d0 + else + menvmax = m + if(m.gt.0.35d0)then + menvmax = 0.35d0*(1.25d0-m)*(1.25d0-m)/0.81d0 + endif + endif + else + mifmcf = -0.023*log10(z)-0.0023 + menvmax = max(m-mcagbf(m)*(1.0d0+mifmcf), 0.0d0) + endif +* + return + end +*** + real*8 FUNCTION mconvenv(kw,m,z,teff,tmin,tset,aj,tm) + implicit none + integer kw + real*8 m,z,teff,tmin,tset,aj,tm + real*8 menvmax + external menvmax +* +* A function to evaluate the mass of the convective envelope +* from Picker, Hirai & Mandel 2024 +* + if(kw.ge.10)then + mconvenv = 0.0d0 + elseif(kw.lt.2)then + if(m.gt.1.25d0)then + mconvenv = 0.0d0 + else + mconvenv = menvmax(kw,m,z)*sqrt(sqrt(1.0d0-aj/tm)) + endif + else + mconvenv = menvmax(kw,m,z)/(1.0d0+exp(4.6d0* + & (tmin+tset-2.0d0*teff)/(tmin-tset))) + endif +* + return + end +*** + real*8 FUNCTION celamhe(m,z,mconv,mconvmax) + implicit none + real*8 m,z,mconv,mconvmax + real*8 m2,b1 +* +* A function to evaluate the binding energy factor of the convective envelope +* from Picker, Hirai & Mandel 2024 +* + m2 = 0.0023d0*log10(z)**2+0.0088d0*log10(z)+0.013d0 + b1 = m2*m-0.23d0 + if((mconv/mconvmax).gt.0.3d0)then + celamhe = exp(0.42*(mconv/mconvmax)+b1) + else + celamhe = exp(0.3*0.42+b1) + endif +* + return + end *** real*8 FUNCTION lHeIf(m,mhefl) implicit none From 7c90e6b334eaa3501e46fa16ddb1f7ac6483bd18 Mon Sep 17 00:00:00 2001 From: TristanBruel Date: Wed, 17 Jun 2026 23:51:09 +0200 Subject: [PATCH 4/4] log in bpp some physical values at the end of the first CE stage when ce2stageflag=1 --- src/cosmic/src/comenv.f | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/src/cosmic/src/comenv.f b/src/cosmic/src/comenv.f index d0a57c52b..748efc802 100644 --- a/src/cosmic/src/comenv.f +++ b/src/cosmic/src/comenv.f @@ -175,6 +175,7 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, * If the secondary star is also giant-like add its envelopes energy. * Determine EORBI based on CEFLAG (CEFLAG=1 for de Kool prescription) * + m2endstage1 = M2 IF(KW2.GE.2.AND.KW2.LE.9.AND.KW2.NE.7)THEN MENVD = MENV/(M2-MC2) IF (using_METISSE.eq.1) THEN @@ -268,6 +269,29 @@ SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1, * IF(CE2STAGEFLAG.EQ.1)THEN SEPF = m1endstage1*m2endstage1/(2.D0*EORBF) +* log end of the first stage with evolve_type=3 + q1_bpp = m1endstage1/m2endstage1 + q2_bpp = 1.d0/q1_bpp + rrl1_bpp = RC1/(RL(q1_bpp)*SEP_postCE) + rrl2_bpp = R2/(RL(q2_bpp)*SEP_postCE) + TB = (SEPF/AURSUN)* + & SQRT(SEPF/(AURSUN*(m1endstage1+m2endstage1))) + CALL writetab(jp,tphys,3.d0, + & m1endstage1,m2endstage1,KW1,KW2, + & SEPF,TB,ECC, + & rrl1_bpp,rrl2_bpp, + & AJ1,AJ2,tms1_bpp,tms2_bpp, + & mc_he(1),mc_he(2),mc_co(1),mc_co(2), + & rad1_bpp,rad2_bpp, + & M02,M01,lumin(1),lumin(2), + & teff1,teff2, + & RC2,RC1,menv_bpp(1),menv_bpp(2), + & renv_bpp(1), + & renv_bpp(2),OSPIN2,OSPIN1,B_0(1),B_0(2), + & bacc(1),bacc(2),tacc(1),tacc(2),epoch(1), + & epoch(2),bhspin2,bhspin1, + & deltam_2,deltam_1,formation2,formation1, + & binstate,mergertype,zpars(14)**2.d5,'bpp') * * Second stage: stable mass transfer of the radiative intershell * Binary hardening formula from Picker, Hirai & Mandel 2024