spinup.f

Go to the documentation of this file.

      SUBROUTINE spinup(IPAIR,ITERM)
*
*
*       Spin synchronization of hierarchical binary.
*       -------------------------------------------
*
      include 'common6.h'
      common/binary/  cm(4,mmax),xrel(3,mmax),vrel(3,mmax),
     &                hm(mmax),um(4,mmax),umdot(4,mmax),tmdis(mmax),
     &                namem(mmax),nameg(mmax),kstarm(mmax),iflagm(mmax)
      REAL*8  m1,m2,mx,mb
      DATA  rg2 /0.1/
      EXTERNAL rl
*
*
*       Set basic pointers.
      i = n + ipair
      i1 = 2*ipair - 1
      iterm = 0
      im = 0
      DO 1 k = 1,nmerge
         IF (namem(k).EQ.name(i)) im = k
    1 CONTINUE
*
*       Find location of the secondary (i.e. ghost).
      CALL findj(i,i2,im)
      IF (i2.LT.0) go to 50
      m1 = cm(1,im)
      m2 = cm(2,im)
      mb = m1 + m2
      semi0 = -0.5*mb/hm(im)
      semi1 = semi0
*
*       Obtain current separation and t'' = 2*U*U'.
      rb = 0.0
      td2 = 0.0
      v20 = 0.0
      DO 5 k = 1,4
          rb = rb + um(k,im)**2
          td2 = td2 + 2.0*um(k,im)*umdot(k,im)
          v20 = v20 + umdot(k,im)**2
    5 CONTINUE
*
*       Evaluate inner eccentricity and exit non-circular orbit.
      ecc2 = (1.0 - rb/semi0)**2 + td2**2/(mb*semi0)
      ecc = sqrt(ecc2)
      IF (ecc.GT.0.1) THEN
          go to 50
      END IF
*
*       Determine largest angular momentum and Roche radius.
      IF (spin(i2).GT.spin(i1)) THEN
          rx = radius(i2)
          mx = m2
          j1 = i2
          q = m2/m1
          rl1 = rl(q)*semi0
      ELSE
          rx = radius(i1)
          mx = m1
          j1 = i1
          q = m1/m2
          rl1 = rl(q)*semi0
      END IF
*
      IF (rx.GT.rl1) THEN
          iterm = 1
          go to 50
      END IF
*
*       Solve for joint angular frequency for star and binary by iteration.
      zmu = m1*m2/mb
      amb0 = zmu*sqrt(mb*semi0)
      spin0 = spin(j1)
      DO 10 l = 1,10
          w = spin(j1)/(rg2*mx*rx**2)
          wb = sqrt(mb/semi0**3)
          spin(j1) = spin(j1)*(wb/w)
          fac = (amb0 + spin0 - spin(j1))**2
*       Obtain new semi-major axis from angular momentum conservation.
          semi = fac/(zmu**2*mb)
*       Copy new value for continued iteration.
          dw = (wb - w)/wb
          iter = l
          IF (abs(dw).LT.1.0d-03) THEN
*             WRITE (6,8)  L, SEMI, WB, W, DW
              go to 15
          END IF
          semi0 = semi
*         WRITE (6,8)  L, SEMI, WB, W, DW
*   8     FORMAT (' ITERATE    L A WB W DW/W ',I4,1P,E12.4,3E10.2)
   10 CONTINUE
*
*       Check minimum size.
   15 IF (semi.LT.rx) THEN
          semi = rx
      END IF
*
*       Update binding energy.
      hi = hm(im)
      hm(im) = -0.5*mb/semi
*
*       Correct EMERGE & ECOLL (consistency; no net effect).
      decorr = zmu*(hi - hm(im))
      emerge = emerge - decorr
      ecoll = ecoll + decorr
      egrav = egrav + decorr
*
*       Specify KS coordinate & velocity scaling factors at general point.
      c2 = sqrt(semi/semi1)
      v2 = 0.5*(mb + hm(im)*rb*(semi/semi1))
      c1 = sqrt(v2/v20)
*
*       Re-scale KS variables to new energy with constant eccentricity.
      rb = 0.0d0
      DO 20 k = 1,4
          um(k,im) = c2*um(k,im)
          umdot(k,im) = c1*umdot(k,im)
*         RB = RB + UM(K,IM)**2
   20 CONTINUE
*
*       Rectify the hierarchical KS variables.
      CALL hirect(im)
*
*       Set new look-up time slightly ahead of components.
      tev0(i) = tev(i)
      tev(i) = tev(j1) + 0.01
*
      WRITE (6,25)  time+toff, name(i1), kstar(j1), iter, semi, w, dw
   25 FORMAT (' SPINUP    T NM K* IT A W DW/W ',f9.3,i6,2i4,1p,3e10.2)
*
*     WRITE (6,30)  KSTAR(J1), KSTAR(I), MX*SMU, RX, RL1,
*    &              TEV(I1), TEV(I2)
*  30 FORMAT (' SPINCHECK   K*1 K*I MX RX RL TEV ',2I4,F7.2,1P,4E10.2)
*
   50 RETURN
*
      END