mix.f

Go to the documentation of this file.

      SUBROUTINE mix(J1,J2,DM)
*
*     Author : J. R. Hurley
*     Date :   7th July 1998
*     
*     Updated by A. D. Railton 23th Feb 2011 to include preMS evolution.
*
*       Evolution parameters for mixed star.
*       ------------------------------------
*
      include 'common6.h'
*
      REAL*8 tscls(20),lums(10),gb(10),tms1,tms2,tms3,tn
      REAL*8 m01,m02,m03,m1,m2,m3,age1,age2,age3,mc3,lum,rm,rcc
      REAL*8 rr,mf,rzams,tprems,left,right,err,xx
      REAL*8 menv,renv,k2e
      REAL*8 mch,mxns
      parameter(mch=1.44d0,mxns=3.0d0,err=1d-5)
      EXTERNAL rzamsf,premsf
      LOGICAL  first
      SAVE  first
      DATA  first /.true./
      DATA  rg2 /0.1d0/
*
*
*       Define global indices with body #I1 being most evolved.
      IF(kstar(j1).GE.kstar(j2))THEN
          i1 = j1
          i2 = j2
          IF(kstar(j1).LE.1.AND.body(j1).LT.body(j2))THEN
             i1 = j2
             i2 = j1
          ENDIF
      ELSE
          i1 = j2
          i2 = j1
      END IF
*
*       Specify case index for collision treatment.
      k1 = kstar(i1)
      k2 = kstar(i2)
*       PreMS collision scheme: -1, -1 ==> -1; -1, 0 ==> -1; -1, x ==> x.
      IF (k1.LT.0.OR.k2.LT.0) THEN
          icase = -1
          IF (max(k1,k2).GT.0) THEN
             icase = max(k1,k2)
             IF (k1 .EQ. -1) epoch(i1) = time*tstar
             IF (k2 .EQ. -1) epoch(i2) = time*tstar
          END IF
      ELSE
          icase = ktype(k1,k2)
      END IF
      IF(icase.GT.100)THEN
         WRITE(38,*)' MIX ERROR ICASE>100 ',icase,k1,k2
      ENDIF
*
*       Record name and stellar type of the colliding stars.
      jc = ncoll + 1
      itype(1) = name(i1)
      itype(2) = name(i2)
      itype(3) = kstar(i1)
      itype(4) = kstar(i2)
*
*       Set physical time and initialize mass loss & time increment.
      ttot = time + toff
      tphys = ttot*tstar
      dm = 0.d0
      dt = 0.d0
      rs1 = radius(i1)*su
      rs2 = radius(i2)*su
*
*       Evolve the stars to the current time unless they have been
*       evolved further by recent Roche interaction.
      tev1 = max(time,tev0(i1))
*       Determine evolution time scales for first star.
      m01 = body0(i1)*zmbar
      m1 = body(i1)*zmbar
      age1 = tev1*tstar - epoch(i1)
      CALL star(k1,m01,m1,tms1,tn,tscls,lums,gb,zpars)
*
*       Obtain time scales for second star.
      m02 = body0(i2)*zmbar
      m2 = body(i2)*zmbar
      age2 = tev1*tstar - epoch(i2)
      CALL star(k2,m02,m2,tms2,tn,tscls,lums,gb,zpars)
*
*       Check for planetary systems - defined as HeWDs - and low-mass WDs!
      IF(k1.EQ.10.AND.m1.LT.0.05)THEN
         icase = k2
         IF(k2.LE.1)THEN
            icase = 1
            age1 = 0.d0
         ENDIF
      ELSEIF(k1.GE.11.AND.m1.LT.0.15.AND.icase.EQ.6)THEN
         icase = 1
      ENDIF
      IF(k2.EQ.10.AND.m2.LT.0.05)THEN
         icase = k1
         IF(k1.LE.1)THEN
            icase = 1
            age2 = 0.d0
         ENDIF
      ENDIF
*
      WRITE(38,67)k1,m01,m1
 67   FORMAT(' MIX OLD *1:',i4,2f10.6)
      WRITE(38,68)k2,m02,m2
 68   FORMAT(' MIX OLD *2:',i4,2f10.6)
*
*       Specify total mass.
      m3 = m1 + m2
      m03 = m01 + m02
      mc3 = 0.d0
      kw = icase
      age3 = 0.d0
      tms3 = 0.d0
*
*       Check energy budget for partial disruption.
      zmb = body(i1) + body(i2)
      rij2 = 0.0
      vij2 = 0.0
      rdot = 0.0
      DO 50 k = 1,3
          rij2 = rij2 + (x(k,i1) - x(k,i2))**2
          vij2 = vij2 + (xdot(k,i1) - xdot(k,i2))**2
          rdot = rdot + (x(k,i1) - x(k,i2))*(xdot(k,i1) - xdot(k,i2))
   50 CONTINUE
      rij = sqrt(rij2)
*       Form binding energy of secondary and orbital kinetic energy.
      eb = -0.5*6.68d-08*4.0d+66*m2**2/(radius(i2)*su*7.0d+10)
      zmu = smu*body(i1)*body(i2)/(body(i1) + body(i2))
      zke = 0.5*2.0d+33*zmu*1.0d+10*vij2*vstar**2
      IF (zke + eb.GT.0.0.AND.rij.LT.0.5*(radius(i1)+radius(i2)).AND.
     &    max(k1,k2).LT.2) THEN
          zm = abs(zke/eb)
          fac = 1.0/sqrt(zm)
          WRITE (6,55)  fac, rij*su, sqrt(vij2)*vstar, eb, zke
   55     FORMAT (' ENFORCED MASS LOSS    FAC RIJ VIJ EB KE ',
     &                                    f7.3,f8.2,f8.2,1p,2e10.2)
*         DM = FAC*M2
          dm = 0.3d0*m2
          m3 = m3 - dm
          m03 = m3
          m02 = m02 - dm
      END IF
*
*       Evaluate apparent age and other parameters.
*
      IF (icase.EQ.-1) THEN
*       Collisions between -1 and -1/0 stars.
         mf = 0.1 
         m3 = (1.0 - mf)*(m1 + m2)
         dm = m1 + m2 - m3
*       Treat case that collided star is greater than 8 solar masses.
*        -> start on ZAMS.   
         IF(m3.GT.8.d0)THEN
            age3 = 0.0d0
         ELSE
*       From internal energy arguments, treating both stars as
*       n = 3/2 polytropes, get new radius:
         rr = m03**2*(1.0-mf)*2.33333
         rr = rr/(m1*m1/rs1 + m2*m2/rs2)*(3.0 - 4.0*mf)
*       Define the preMS variables.
         tprems = 10d0**(43.628d0-(35.835d0*m3**1.5029d-2)
     &            *exp(m3*3.9608d-3))
         rzams = rzamsf(m3)
*
*       Do bisection on the function PREMSR to find the time at
*       which the preMS star is at this radius.
*
         ic = 0
         left = 0
         right = 1
*       Make sure new radius is in right range.
         IF(rr.LE.rzams)THEN
            xx=0
            IF(m3.LT.1.25) kw = 0
            IF(m3.GE.1.25) kw = 1
            goto 103
         ENDIF
         xx = 1d0
         IF(premsf(m3,xx,rr).LE.0d0)THEN
            xx=1d0
            goto 103
         ENDIF      
*
 101     xx=0.5*(left + right)
         ic = ic + 1
         IF (ic.GT.100) stop
         IF(abs(premsf(m3,xx,rr)).LE.err)goto 103
         IF(abs(premsf(m3,xx,rr)).GT.err)goto 102
*
 102     IF(sign(1.d0,premsf(m3,xx,rr)).EQ.
     &        sign(1.d0,premsf(m3,left,rr)))left=xx
         IF(sign(1.d0,premsf(m3,xx,rr)).EQ.
     &        sign(1.d0,premsf(m3,right,rr)))right=xx
         goto 101
*     
 103     age3 = -tprems*xx            
         END IF
      IF (m3.GT.2.0) WRITE (6,720)  m3, rzams, xx, tprems, age3
  720 FORMAT (' WATCH!   M3 RZAMS XX TP AGE3 ',f7.2,1p,4e10.2)
         CALL star(kw,m3,m3,tms3,tn,tscls,lums,gb,zpars)
*     
      ELSE IF(icase.EQ.1)THEN
*       Specify new age based on complete mixing.
         IF(k1.EQ.7) kw = 7
         CALL star(kw,m03,m3,tms3,tn,tscls,lums,gb,zpars)
         age3 = 0.1d0*tms3*(age1*m01/tms1 + age2*m02/tms2)/m03
      ELSEIF(icase.EQ.3.OR.icase.EQ.6.OR.icase.EQ.9)THEN
         mc3 = m1
         CALL gntage(mc3,m3,kw,zpars,m03,age3)
      ELSEIF(icase.EQ.4)THEN
         mc3 = m1
         age3 = age1/tms1
         IF(age3.GT.1.d0)THEN
            WRITE(6,*)' WARNING! MIX AGE WRONG FOR KW=4 ',age1,tms1
            WRITE(6,*)k1,m01,m1,tev1,epoch(i1)
            age3 = 0.99d0
         ENDIF
         CALL gntage(mc3,m3,kw,zpars,m03,age3)
      ELSEIF(icase.EQ.7)THEN
         CALL star(kw,m03,m3,tms3,tn,tscls,lums,gb,zpars)
         age3 = tms3*(age2*m2/tms2)/m3
      ELSEIF(icase.LE.12)THEN
*       Ensure that a new WD has the initial mass set correctly.
         m03 = m3
         dt = 1.0d+04
      ELSEIF(icase.EQ.13.OR.icase.EQ.14)THEN
*       Set unstable Thorne-Zytkow object with fast mass loss of envelope
*       unless the less evolved star is a WD, NS or BH.
         IF(k2.LT.10)THEN
            m03 = m1
            m3 = m1
            dm = m2
            ntz = ntz + 1
            WRITE (6,2)  k1, k2, name(i1), name(i2)
    2       FORMAT (' NEW TZ    K* NM ',2i4,2i6)
         ELSEIF(icase.EQ.13)THEN
            ngb = ngb + 1
            IF(m3.GT.mxns) kw = 14
         ENDIF
         dt = 1.0d+04
      ELSEIF(icase.EQ.15)THEN
         dm = m3
         m3 = 0.d0
      ELSEIF(icase.GT.100)THEN
*       Common envelope case which should only be used after COMENV.
         kw = k1
         age3 = age1
         m3 = m1
         m03 = m01
         dm = m2
      ELSE
*       This should not be reached.
        WRITE(6,*)' ERROR MIX: ICASE NOT CAUGHT!!!'
        WRITE(6,*)' K1 K2 -> K3 ',k2,k2,kw
        stop
      ENDIF
*
      WRITE(38,69)kw,m03,m3
 69   FORMAT(' MIX NEW *3:',i4,2f10.6)
*
*       Determine consistent stellar type and specify mass loss.
      IF(kw.LE.14)THEN
         kw0 = kw 
         CALL star(kw,m03,m3,tms3,tn,tscls,lums,gb,zpars)
         CALL hrdiag(m03,age3,m3,tms3,tn,tscls,lums,gb,zpars,
     &               rm,lum,kw,mc3,rcc,menv,renv,k2e)
         IF(kw.NE.kw0)then
            dm = m1 + m2 - m3
            WRITE (6,5)  name(i1), name(i2), kw0, kw, dm
    5       FORMAT (' MIX TYPE CHANGE:  NAM K* DM ',2i6,2i4,f6.2)
         ENDIF
      ENDIF
      kstar(i1) = kw
      radius(i1) = rm/su
      zlmsty(i1) = lum
*       Update initial mass, epoch & evolution times.
      body0(i1) = m03/zmbar
      epoch(i1) = tev1*tstar - age3
      tev(i1) = tev1 + dt
      tev0(i1) = tev1
      dm = dm/zmbar
*
*       Record final type of mixed star.
      itype(5) = kstar(i1)
*
*       Check for blue straggler formation (TM < TPHYS & KSTAR <= 1).
      IF(tms3.LT.tphys.AND.kstar(i1).LE.1)THEN
          WRITE (6,8)  name(i1), m3, tms3, tphys, age3
    8     FORMAT (' NEW BS (MIX):    NAM M3 TM TP AGE ',
     &                               i6,f6.2,3f8.1)
*       Prepare spin diagnostics for fort.91 (rg2=0.1 is assumed).
          semi = 2.0/rij - vij2/zmb
          semi = 1.0/semi
          pd = days*semi*sqrt(abs(semi)/zmb)
          pd = max(pd,0.0d0)
          ecc2 = (1.0 - rij/semi)**2 + rdot**2/(semi*zmb)
          ecc = sqrt(ecc2)
          spin1 = spin(i1)/(rg2*body(i1)*radius(i1)**2)
          spin2 = spin(i2)/(rg2*body(i2)*radius(i2)**2)
          ts = 1.0d+06*365.0*twopi*tstar
          nbs = nbs + 1
          WRITE (91,9)  time+ttot, name(i1), name(i2), m3, ecc, pd,
     &                  ts/spin1, ts/spin2
    9     FORMAT (' NEW BS    T NAM M3 ECC P ROT ',
     &                        f8.1,2i7,f6.2,f8.4,1p,3e9.1)
          CALL flush(91)
      ENDIF
*
      WRITE (6,10)  iqcoll, k1, k2, kstar(i1), m1, m2, m3, rs1, rs2,
     &              radius(i1)*su, dm*zmbar
   10 FORMAT (' NEW STAR:    IQ K1 K2 K* M1 M2 M3 R1 R2 R* DM ',
     &                       4i3,3f6.1,3f7.1,f5.1)
*
*       Open unit #13 the first time.
      IF(first.AND.iqcoll.NE.3)THEN
          OPEN (unit=13,status='NEW',form='FORMATTED',file='COLL')
          first = .false.
*
*       Print cluster scaling parameters at start of the run.
          IF(ncoll.EQ.0)THEN
              WRITE (13,20)  rbar, bodym*zmbar, body1*zmbar, tscale,
     &                       nbin0, nzero
   20         FORMAT (/,6x,'MODEL:    RBAR =',f5.1,'  <M> =',f6.2,
     &                     '  M1 =',f6.1,'  TSCALE =',f6.2,
     &                     '  NB =',i4,'  N0 =',i6,//)
              WRITE (13,25)
   25         FORMAT ('    TIME  NAME  NAME  K1  K2  KC  M1   M2   MC',
     &                '   DM    R1     R2    r/Rc   R     ECC      P',/)
          ENDIF
      ENDIF
*
*       Form central distance (scaled by RC) and period in days.
      ri2 = 0.d0
      rij2 = 0.d0
      vij2 = 0.d0
      DO 30 k = 1,3
          ri2 = ri2 + (x(k,i1) - rdens(k))**2
          rij2 = rij2 + (x(k,i1) - x(k,i2))**2
          vij2 = vij2 + (xdot(k,i1) - xdot(k,i2))**2
   30 CONTINUE
      ri = sqrt(ri2)
      rij = sqrt(rij2)
      semi = 2.d0/rij - vij2/zmb
      semi = 1.d0/semi
      tk = days*semi*sqrt(abs(semi)/zmb)
      ecc = 1.0 - rij/semi
*
*       Accumulate collision diagnostics on unit #13 (filename COLL).
      IF (iqcoll.NE.3) THEN
          WRITE (13,35)  ttot, (itype(k),k=1,5), m1, m2, m3,
     &                   dm*zmbar, rs1, rs2, ri/rc, rij*su, ecc, tk
   35     FORMAT (1x,f7.1,2i6,3i4,4f5.1,2f7.2,f6.2,f7.2,f9.5,1p,e9.1)
          CALL flush(13)
      END IF
*
*       Re-define indices of colliding bodies with J1 as new c.m.
      j1 = i1
      j2 = i2
*
      IF(kstar(i1).GT.12)THEN
          WRITE (15,40)  k1, k2, kstar(i1), body(i1)*zmbar,
     &                   body(i2)*zmbar, m3
   40     FORMAT (' MIX:    K1 K2 K* M1 M2 M3 ',3i4,3f7.2)
          CALL flush(15)
      ENDIF
*
      RETURN
      END
***