comenv.f

Go to the documentation of this file.

***
      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
***