https://github.com/florentrenaud/nbody6tt
Tip revision: 8a4716382ead3ece116c48a4ae5c65f8c9534437 authored by Florent on 29 January 2015, 12:19:28 UTC
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Nbody6 - 29 January 2015 (added GPU2/Build/.gitkeep)
Tip revision: 8a47163
comenv.f
***
SUBROUTINE COMENV(M01,M1,MC1,AJ1,JSPIN1,KW1,
& M02,M2,MC2,AJ2,JSPIN2,KW2,ECC,SEP,COEL)
*
* Common Envelope Evolution.
*
INCLUDE 'common6.h'
*
INTEGER KW1,KW2,KW
* INTEGER IDUM
*
REAL*8 M01,M1,MC1,AJ1,JSPIN1,R1,L1
REAL*8 M02,M2,MC2,AJ2,JSPIN2,R2,L2,MC22
REAL*8 TSCLS1(20),TSCLS2(20),LUMS(10),GB(10),TM1,TM2,TN
REAL*8 EBINDI,EBINDF,EORBI,EORBF,ECIRC
REAL*8 AA,BB,CC,EFAC
REAL*8 ECC,SEP,SEPF,SEPL,MF,XX
REAL*8 CONST,DELY,DERI,DELMF,MC3,FAGE1,FAGE2
REAL*8 TB,OORB,OSPIN1,OSPIN2
REAL*8 RC1,RC2,Q1,Q2,RL1,RL2
REAL*8 MENV,RENV,MENVD,RZAMS
REAL*8 LAMB1,LAMB2,K21,K22
REAL*8 AURSUN,K3,LAMBDA,ALPHA1,MCH
PARAMETER (AURSUN = 214.95D0,K3 = 0.21D0)
PARAMETER (LAMBDA = 0.D0,ALPHA1 = 3.D0)
PARAMETER (MCH = 1.44D0)
LOGICAL COEL,ECCFLG
REAL*8 CELAMF,RL,RZAMSF
EXTERNAL CELAMF,RL,RZAMSF
*
* Common envelope evolution - entered only when KW1 = 2, 3, 4, 5, 6, 8 or 9.
*
* For simplicity energies are divided by -G.
*
COEL = .FALSE.
ECCFLG = .TRUE.
*
* Obtain the core masses and radii.
*
KW = KW1
CALL star(KW1,M01,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS)
IF(AJ1.GT.TN) AJ1 = TN
CALL hrdiag(M01,AJ1,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS,
& R1,L1,KW1,MC1,RC1,MENV,RENV,K21)
OSPIN1 = JSPIN1/(K21*R1*R1*(M1-MC1)+K3*RC1*RC1*MC1)
IF(M1.GT.MC1)THEN
MENVD = MENV/(M1-MC1)
RZAMS = RZAMSF(M01)
LAMB1 = CELAMF(KW,M01,L1,R1,RZAMS,MENVD,LAMBDA)
ELSE
LAMB1 = 0.5D0
ENDIF
IF(KW.NE.KW1) WRITE(38,*)' COMENV TYPE CHANGE *1'
*
KW = KW2
CALL star(KW2,M02,M2,TM2,TN,TSCLS2,LUMS,GB,ZPARS)
IF(AJ2.GT.TN) AJ2 = TN
CALL hrdiag(M02,AJ2,M2,TM2,TN,TSCLS2,LUMS,GB,ZPARS,
& R2,L2,KW2,MC2,RC2,MENV,RENV,K22)
OSPIN2 = JSPIN2/(K22*R2*R2*(M2-MC2)+K3*RC2*RC2*MC2)
IF(KW.NE.KW2) WRITE(38,*)' COMENV TYPE CHANGE *2'
*
IF(SEP.GT.0.D0)THEN
TB = (SEP/AURSUN)*SQRT(SEP/(AURSUN*(M1+M2)))
OORB = TWOPI/TB
ENDIF
*
WRITE(38,66)KW1,M01,M1,MC1,AJ1,SEP
66 FORMAT(' COMENV OLD *1:',I4,3F7.2,F9.2,F12.4)
WRITE(38,67)KW2,M02,M2,MC2,AJ2
67 FORMAT(' COMENV OLD *2:',I4,3F7.2,F9.2)
*
* Calculate the binding energy of the giant envelope (multiplied by lambda).
*
EBINDI = M1*(M1-MC1)/(LAMB1*R1)
*
* Calculate the initial orbital energy
*
EORBI = MC1*M2/(2.D0*SEP)
*
* Allow for an eccentric orbit.
*
* ECIRC = EORBI/(1.D0 - ECC*ECC)
*
IF(ECCFLG.AND.ECC.GT.0.01D0)THEN
*
* Assume 10% of the binding energy is lost on the first orbit and determine
* whether the star just passes through. For the moment put all the
* secondary's envelope with the primary's.
*
* DE = MIN(EBINDI/10.D0,(ECIRC-EORBI)*ALPHA1)
*
* Adopt new procedure based on a cubic fitting (May 2003).
* (Oct 2005: YP is maximum for a head-on collision and would reduce
* to zero for a grazing encounter, i.e. QP=1. However,
* in order to be consistent with Kochanek collision detection
* which allows QP values up to 1.7 on entry to expel/comenv
* we set a non-zero YP for QP > 1. This is a temporary fix
* until a better collision detection scheme comes along sometime
* in the near future, we hope.)
*
QP = SEP*(1.D0 - ECC)/(R1 + R2)
YP = 0.5D0*(QP - 0.5D0)**3 - 0.375D0*(QP - 0.5D0) + 0.125D0
YP = MAX(YP,0.05D0)
YP = MIN(YP,0.5D0)
DE = YP*EBINDI
* Remove all envelope mass if below 0.05 m_sun.
IF (M1-MC1.LT.0.05) DE = EBINDI
EBINDF = EBINDI - DE
MF = 0.5D0*(MC1 + SQRT(MC1*MC1 + 4.D0*LAMB1*R1*EBINDF))
EORBI = MF*M2/(2.D0*SEP)
ECIRC = EORBI/(1.D0 - ECC*ECC)
EORBF = EORBI + DE/ALPHA1
IF(EORBF.LT.ECIRC)THEN
M1 = MF
ECC = SQRT(1.D0 - EORBF/ECIRC)
ELSE
* EFAC = 1.d0 - 1.d0/(ECC*ECC)
EFAC = ECC*ECC/(1.d0 - ECC*ECC)
AA = 1.d0/(LAMB1*R1)
* BB = (2.d0*M1-MC1)/(LAMB1*R1) + ALPHA1*M2/(SEP*EFAC)
BB = (2.d0*M1-MC1)/(LAMB1*R1) + ALPHA1*EFAC*M2/(2.d0*SEP)
* CC = M1*MC1/(LAMB1*R1) + ALPHA1*M1*M2/(SEP*EFAC)
CC = ALPHA1*EFAC*M1*M2/(2.d0*SEP)
DM = (BB - SQRT(BB*BB - 4.d0*AA*CC))/(2.d0*AA)
IF(DM.LE.0.0.OR.DM.GT.(M1-MC1))THEN
WRITE(*,*)' WARNING: CE DM ',DM,M1,MC1
ENDIF
M1 = MAX(M1-DM,MC1)
EBINDF = M1*(M1-MC1)/(LAMB1*R1)
DE = EBINDI - EBINDF
EORBF = EORBI + DE/ALPHA1
ECC = 0.001D0
ENDIF
SEPF = M1*M2/(2.D0*EORBF)
IF((M1-MC1).LT.1.0D-07)THEN
M1 = MC1
CALL star(KW1,M01,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS)
CALL hrdiag(M01,AJ1,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS,
& R1,L1,KW1,MC1,RC1,MENV,RENV,K21)
ELSEIF(KW1.EQ.4)THEN
AJ1 = (AJ1 - TSCLS1(2))/TSCLS1(3)
CALL gntage(MC1,M1,KW1,ZPARS,M01,AJ1)
ENDIF
GOTO 30
ENDIF
*
* If the secondary star is also giant-like add its envelopes's energy.
*
IF(KW2.GE.2.AND.KW2.LE.9.AND.KW2.NE.7)THEN
IF(M2.GT.MC2)THEN
MENVD = MENV/(M2-MC2)
RZAMS = RZAMSF(M02)
LAMB2 = CELAMF(KW,M02,L2,R2,RZAMS,MENVD,LAMBDA)
ELSE
LAMB2 = 0.5D0
ENDIF
EBINDI = EBINDI + M2*(M2-MC2)/(LAMB2*R2)
*
* Re-calculate the initial orbital energy
*
EORBI = MC1*MC2/(2.D0*SEP)
ENDIF
*
*
* Calculate the final orbital energy without coalescence.
*
EORBF = EORBI + EBINDI/ALPHA1
*
* If the secondary is on the main sequence see if it fills its Roche lobe.
*
IF(KW2.LE.1.OR.KW2.EQ.7)THEN
SEPF = MC1*M2/(2.D0*EORBF)
***
* Assume the energy generated by forcing the secondary to
* co-rotate goes into the envelope (experimental).
* TB = (SEPF/AURSUN)*SQRT(SEPF/(AURSUN*(MC1+M2)))
* OORB = TWOPI/TB
* DELY = 0.5D0*M2*R2*R2*(OSPIN2*OSPIN2 - OORB*OORB)/3.91D+08
* DELY = K22*DELY
* EBINDI = MAX(0.D0,EBINDI - DELY)
* EORBF = EORBI + EBINDI/ALPHA1
* SEPF = MC1*M2/(2.D0*EORBF)
***
Q1 = MC1/M2
Q2 = 1.D0/Q1
RL1 = RL(Q1)
RL2 = RL(Q2)
IF(RC1/RL1.GE.R2/RL2)THEN
*
* The helium core of a very massive star of type 4 may actually fill
* its Roche lobe in a wider orbit with a very low-mass secondary.
*
IF(RC1.GT.RL1*SEPF)THEN
COEL = .TRUE.
SEPL = RC1/RL1
ENDIF
ELSE
IF(R2.GT.RL2*SEPF)THEN
COEL = .TRUE.
SEPL = R2/RL2
ENDIF
ENDIF
IF(COEL)THEN
*
KW = KTYPE(KW1,KW2) - 100
MC3 = MC1
IF(KW2.EQ.7.AND.KW.EQ.4) MC3 = MC3 + M2
*
* Coalescence - calculate final binding energy.
*
EORBF = MAX(MC1*M2/(2.D0*SEPL),EORBI)
EBINDF = EBINDI - ALPHA1*(EORBF - EORBI)
ELSE
*
* Primary becomes a black hole, neutron star, white dwarf or helium star.
*
MF = M1
M1 = MC1
CALL star(KW1,M01,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS)
CALL hrdiag(M01,AJ1,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS,
& R1,L1,KW1,MC1,RC1,MENV,RENV,K21)
ENDIF
ELSE
*
* Degenerate or giant secondary. Check if the least massive core fills its
* Roche lobe.
*
SEPF = MC1*MC2/(2.D0*EORBF)
***
* Assume the energy generated by forcing the secondary to
* co-rotate goes into the envelope (experimental).
* IF(KW2.GE.10.AND.KW2.LE.14)THEN
* TB = (SEPF/AURSUN)*SQRT(SEPF/(AURSUN*(MC1+MC2)))
* OORB = TWOPI/TB
* DELY = 0.5D0*M2*R2*R2*(OSPIN2*OSPIN2 - OORB*OORB)/3.91D+08
* DELY = K3*DELY
* EBINDI = MAX(0.D0,EBINDI - DELY)
* EORBF = EORBI + EBINDI/ALPHA1
* SEPF = MC1*MC2/(2.D0*EORBF)
* ENDIF
***
Q1 = MC1/MC2
Q2 = 1.D0/Q1
RL1 = RL(Q1)
RL2 = RL(Q2)
IF(RC1/RL1.GE.RC2/RL2)THEN
IF(RC1.GT.RL1*SEPF)THEN
COEL = .TRUE.
SEPL = RC1/RL1
ENDIF
ELSE
IF(RC2.GT.RL2*SEPF)THEN
COEL = .TRUE.
SEPL = RC2/RL2
ENDIF
ENDIF
*
IF(COEL)THEN
*
* If the secondary was a neutron star or black hole the outcome
* is an unstable Thorne-Zytkow object that leaves only the core.
*
SEPF = 0.D0
IF(KW2.GE.13)THEN
MC1 = MC2
M1 = MC1
MC2 = 0.D0
M2 = 0.D0
KW1 = KW2
KW2 = 15
AJ1 = 0.D0
*
* The envelope mass is not required in this case.
*
GOTO 30
ENDIF
*
KW = KTYPE(KW1,KW2) - 100
MC3 = MC1 + MC2
*
* Calculate the final envelope binding energy.
*
EORBF = MAX(MC1*MC2/(2.D0*SEPL),EORBI)
EBINDF = EBINDI - ALPHA1*(EORBF - EORBI)
*
* Check if we have the merging of two degenerate cores and if so
* then see if the resulting core will survive or change form.
*
IF(KW1.EQ.6.AND.(KW2.EQ.6.OR.KW2.GE.11))THEN
CALL dgcore(KW1,KW2,KW,MC1,MC2,MC3,EBINDF)
ENDIF
IF(KW1.LE.3.AND.M01.LE.ZPARS(2))THEN
IF((KW2.GE.2.AND.KW2.LE.3.AND.M02.LE.ZPARS(2)).OR.
& KW2.EQ.10)THEN
CALL dgcore(KW1,KW2,KW,MC1,MC2,MC3,EBINDF)
IF(KW.GE.10)THEN
KW1 = KW
M1 = MC3
MC1 = MC3
IF(KW.LT.15) M01 = MC3
AJ1 = 0.D0
MC2 = 0.D0
M2 = 0.D0
KW2 = 15
GOTO 30
ENDIF
ENDIF
ENDIF
*
ELSE
*
* The cores do not coalesce - assign the correct masses and ages.
*
MF = M1
M1 = MC1
CALL star(KW1,M01,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS)
CALL hrdiag(M01,AJ1,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS,
& R1,L1,KW1,MC1,RC1,MENV,RENV,K21)
MF = M2
KW = KW2
M2 = MC2
CALL star(KW2,M02,M2,TM2,TN,TSCLS2,LUMS,GB,ZPARS)
CALL hrdiag(M02,AJ2,M2,TM2,TN,TSCLS2,LUMS,GB,ZPARS,
& R2,L2,KW2,MC2,RC2,MENV,RENV,K22)
ENDIF
ENDIF
*
IF(COEL)THEN
MC22 = MC2
IF(KW.EQ.4.OR.KW.EQ.7)THEN
* If making a helium burning star calculate the fractional age
* depending on the amount of helium that has burnt.
IF(KW1.LE.3)THEN
FAGE1 = 0.D0
ELSEIF(KW1.GE.6)THEN
FAGE1 = 1.D0
ELSE
FAGE1 = (AJ1 - TSCLS1(2))/(TSCLS1(13) - TSCLS1(2))
ENDIF
IF(KW2.LE.3.OR.KW2.EQ.10)THEN
FAGE2 = 0.D0
ELSEIF(KW2.EQ.7)THEN
FAGE2 = AJ2/TM2
MC22 = M2
ELSEIF(KW2.GE.6)THEN
FAGE2 = 1.D0
ELSE
FAGE2 = (AJ2 - TSCLS2(2))/(TSCLS2(13) - TSCLS2(2))
ENDIF
ENDIF
ENDIF
*
* Now calculate the final mass following coelescence. This requires a
* Newton-Raphson iteration.
*
IF(COEL)THEN
*
* Calculate the orbital spin just before coalescence.
*
TB = (SEPL/AURSUN)*SQRT(SEPL/(AURSUN*(MC1+MC2)))
OORB = TWOPI/TB
*
XX = 1.D0 + ZPARS(7)
IF(EBINDF.LE.0.D0)THEN
MF = MC3
GOTO 20
ELSE
CONST = ((M1+M2)**XX)*(M1-MC1+M2-MC22)*EBINDF/EBINDI
ENDIF
*
* Initial Guess.
*
MF = MAX(MC1 + MC22,(M1 + M2)*(EBINDF/EBINDI)**(1.D0/XX))
10 DELY = (MF**XX)*(MF - MC1 - MC22) - CONST
* IF(ABS(DELY/MF**(1.D0+XX)).LE.1.0D-02) GOTO 20
IF(ABS(DELY/MF).LE.1.0D-03) GOTO 20
DERI = MF**ZPARS(7)*((1.D0+XX)*MF - XX*(MC1 + MC22))
DELMF = DELY/DERI
MF = MF - DELMF
GOTO 10
*
* Set the masses and separation.
*
20 IF(MC22.EQ.0.D0) MF = MAX(MF,MC1+M2)
M2 = 0.D0
M1 = MF
KW2 = 15
*
* Combine the core masses.
*
IF(KW.EQ.2)THEN
CALL star(KW,M1,M1,TM2,TN,TSCLS2,LUMS,GB,ZPARS)
IF(GB(9).GE.MC1)THEN
M01 = M1
AJ1 = TM2 + (TSCLS2(1) - TM2)*(AJ1-TM1)/(TSCLS1(1) - TM1)
CALL star(KW,M01,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS)
ENDIF
ELSEIF(KW.EQ.7)THEN
M01 = M1
CALL star(KW,M01,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS)
AJ1 = TM1*(FAGE1*MC1 + FAGE2*MC22)/(MC1 + MC22)
ELSEIF(KW.EQ.4.OR.MC2.GT.0.D0.OR.KW.NE.KW1)THEN
IF(KW.EQ.4) AJ1 = (FAGE1*MC1 + FAGE2*MC22)/(MC1 + MC22)
MC1 = MC1 + MC2
MC2 = 0.D0
*
* Obtain a new age for the giant.
*
CALL gntage(MC1,M1,KW,ZPARS,M01,AJ1)
CALL star(KW,M01,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS)
ENDIF
CALL hrdiag(M01,AJ1,M1,TM1,TN,TSCLS1,LUMS,GB,ZPARS,
& R1,L1,KW,MC1,RC1,MENV,RENV,K21)
JSPIN1 = OORB*(K21*R1*R1*(M1-MC1)+K3*RC1*RC1*MC1)
JSPIN2 = 0.D0
KW1 = KW
ECC = 0.001D0
ELSE
*
* Determine if any eccentricity remains in the orbit.
*
* IF(ECCFLG.AND.EORBF.LT.ECIRC)THEN
* ECC = SQRT(1.D0 - EORBF/ECIRC)
* ELSE
* ECC = 0.001D0
* ENDIF
*
* Set both cores in co-rotation with the orbit on exit of CE,
*
TB = (SEPF/AURSUN)*SQRT(SEPF/(AURSUN*(M1+M2)))
OORB = TWOPI/TB
* JSPIN1 = OORB*(K21*R1*R1*(M1-MC1)+K3*RC1*RC1*MC1)
* JSPIN2 = OORB*(K22*R2*R2*(M2-MC2)+K3*RC2*RC2*MC2)
*
* or, leave the spins of the cores as they were on entry.
* Tides will deal with any synchronization later.
*
JSPIN1 = OSPIN1*(K21*R1*R1*(M1-MC1)+K3*RC1*RC1*MC1)
JSPIN2 = OSPIN2*(K22*R2*R2*(M2-MC2)+K3*RC2*RC2*MC2)
ENDIF
30 SEP = SEPF
AJ1 = MAX(AJ1,1.0D-10)
AJ2 = MAX(AJ2,1.0D-10)
* AJ1 = MAX(AJ1,1.0D-10*TM1)
* AJ2 = MAX(AJ2,1.0D-10*TM2)
*
WRITE(38,68)KW1,M01,M1,MC1,AJ1
68 FORMAT(' COMENV NEW *1:',I4,3F7.2,F9.2)
WRITE(38,69)KW2,M02,M2,MC2,AJ2
69 FORMAT(' COMENV NEW *2:',I4,3F7.2,F9.2)
*
RETURN
END
***