kick.f

Go to the documentation of this file.

      SUBROUTINE kick(I,ICASE,KW,DM)
*
*
*       Velocity kick for WD, neutron stars or black holes.
*       ---------------------------------------------------
*
      include 'common6.h'
      parameter(vfac=2.0d0)
      REAL*8  ran2,vk(4)
      SAVE  ipair, kc, vdis, ri
      DATA  ipair,kc /0,0/
*       WD velocity kicks: Fellhauer et al. 2003, ApJ Letters 595, L53.
*
*
*       Save orbital parameters in case of KS binary (called from KSAPO).
      IF (icase.EQ.0) THEN
          ipair = i
          kc = kstar(n+ipair)
*       Identify the correct component (KSTAR reversed in MDOT).
          i1 = 2*ipair - 1
          IF (kstar(i1).LT.0) THEN
              in = i1
          ELSE IF (kstar(i1+1).LT.0) THEN
              in = i1 + 1
          END IF
          kstar(in) = -kstar(in)
*
*       Determine disruption velocity after mass loss.
          vd2 = 2.0*body(n+ipair)/r(ipair)
          vdis = sqrt(vd2)*vstar
*       Set cluster escape velocity (add twice central potential).
          vp2 = 2.0*body(n+ipair)/r(ipair)
          vesc = sqrt(vp2 + 4.0)*vstar
          semi = -0.5*body(n+ipair)/h(ipair)
          zm1 = body(i1)*smu
          zm2 = body(i1+1)*smu
*       Skip case of massless binary.
          IF (body(n+ipair).GT.0.0d0) THEN
              eb = body(i1)*body(i1+1)/body(n+ipair)*h(ipair)
          ELSE
              eb = 0.0
          END IF
          ri = r(ipair)
*       Skip on #25 = 0/1 for consistent net WD modification of EKICK.
          IF ((kw.LT.13.AND.kz(25).EQ.0).OR.
     &        (kw.EQ.12.AND.kz(25).NE.2)) go to 30
*       Sum whole binding energy (used by BINOUT for net change).
          ekick = ekick + eb
          egrav = egrav + eb
          i2 = i1 + 1
          WRITE (6,1)  name(i1), name(i2), kstar(i1), kstar(i2), zm1,
     &                 zm2, vesc, vdis, r(ipair)/semi, eb, r(ipair)
    1     FORMAT (' BINARY KICK:    NAM K* M1 M2 VESC VDIS R/A EB R ',
     &                              2i6,2i4,4f7.1,f6.2,f9.4,1p,e10.2)
          nbkick = nbkick + 1
*       Remove any circularized KS binary from the CHAOS/SYNCH table.
          IF (kstar(n+ipair).GE.10.AND.nchaos.GT.0) THEN
              ii = -(n + ipair)
              CALL spiral(ii)
              kstar(n+ipair) = 0
          END IF
          go to 30
      END IF
*
*       Generate velocity kick for neutron star (Gordon Drukier Tokyo paper).
*     IT = 0
*     V0 = 330.0
*     VCUT = 1000.0
*   2 VT = VCUT/V0*RAN2(IDUM1)
*     VP = VT*(2.0*RAN2(IDUM1) - 1.0)
*     VN = SQRT(VT**2 - VP**2)
*     FAC = 1.0/0.847*VN**0.3/(1.0 + VN**3.3)
*     IF (FAC.LT.RAN2(IDUM1).AND.IT.LT.10) GO TO 2
*     VKICK = V0*VT
*
*       Adopt the Maxwellian of Hansen & Phinney (MN 291, 569, 1997).
      disp = 190.0
*
*       Include velocity dispersion in terms of VSTAR (Parameter statement!).
      IF (vfac.GT.0.0d0) THEN
          disp = vfac*vstar
      END IF
*
*       Allow for optional type-dependent WD kick.
      IF (kw.LT.13) THEN
          IF (kz(25).EQ.1.AND.(kw.EQ.10.OR.kw.EQ.11)) THEN
              disp = 5.0
          ELSE IF (kz(25).GT.1.AND.kw.EQ.12) THEN
              disp = 5.0
          ELSE
              disp = 0.0
          END IF
      END IF
*     IF (KW.EQ.14) DISP = 0.0
*
*       Use Henon's method for pairwise components (Douglas Heggie 22/5/97).
      DO 2 k = 1,2
          x1 = ran2(idum1)
          x2 = ran2(idum1)
*       Generate two velocities from polar coordinates S & THETA.
          s = disp*sqrt(-2.0*log(1.0 - x1))
          theta = twopi*x2
          vk(2*k-1) = s*cos(theta)
          vk(2*k) = s*sin(theta)
    2 CONTINUE
      vkick = sqrt(vk(1)**2 + vk(2)**2 + vk(3)**2)
      vk(4) = vkick
*
*       Limit kick velocity to VDIS+10*VSTAR/10*VST for binary/single stars.
      IF (ipair.GT.0) THEN
          vbf = sqrt(vdis**2 + 100.0*vstar**2)
          vkick = min(vkick,vbf)
      ELSE
          vkick = min(vkick,10.0d0*vstar)
*       Ensure escape of massless star.
          IF (body(i).EQ.0.0d0) vkick = 10.0*vstar
      END IF
      vkick = vkick/vstar
*       Skip WD case with KZ(25) = 0.
      IF (vkick.EQ.0.0d0) go to 30
*
*       Randomize the velocity components.
*     A(4) = 0.0
*     DO 5 K = 1,3
*         A(K) = 2.0*RAN2(IDUM1) - 1.0
*         A(4) = A(4) + A(K)**2
*   5 CONTINUE
*
*       Add truncated/full kick velocity and initialize X0DOT.
      vi2 = 0.0
      vf2 = 0.0
      DO 10 k = 1,3
          vi2 = vi2 + xdot(k,i)**2
*         XDOT(K,I) = XDOT(K,I) + VKICK*A(K)/SQRT(A(4))
          xdot(k,i) = xdot(k,i) + vkick*vk(k)/vk(4)
          x0dot(k,i) = xdot(k,i)
          vf2 = vf2 + xdot(k,i)**2
   10 CONTINUE
*
*       Modify energy loss due to increased velocity of single particle.
      ecdot = ecdot - 0.5*body(i)*(vf2 - vi2)
      nkick = nkick + 1
*
*       Replace final velocity by relative velocity for binary kick.
      IF (ipair.GT.0) THEN
          jp = kvec(i)
          j = i + 1
          IF (i.EQ.2*jp) j = i - 1
          vf2 = 0.0
          DO 15 k = 1,3
              vf2 = vf2 + (xdot(k,i) - xdot(k,j))**2
   15     CONTINUE
          hnew = 0.5*vf2 - (body(i) + body(j))/ri
          eb1 = body(i)*body(j)/(body(i) + body(j))*hnew
          IF (eb1.LT.0.0) THEN
              ekick = ekick - eb1
              egrav = egrav - eb1
          END IF
          ipair = 0
      END IF
*
      IF (nkick.LT.50.OR.name(i).LE.2*nbin0.OR.
     &    (kw.GE.13.AND.ttot*tstar.GT.100.0)) THEN
          zm = body(i)*zmbar
          WRITE (6,20)  i, name(i), kstar(i), kc, body0(i)*zmbar, zm,
     &                  sqrt(vi2)*vstar, vkick*vstar, sqrt(vf2)*vstar
   20     FORMAT (' VELOCITY KICK:    I NAM K* KC* M0 M VI VK VF ',
     &                                2i6,2i4,2f7.2,3f7.1)
          kc = 0
      END IF
*
*       Highlight BH/NS velocities below 4 times rms velocity.
      IF (vkick.LT.4.0*sqrt(0.5).AND.kw.GE.13) THEN
          WRITE (6,25)  i, name(i), kw, vkick*vstar, sqrt(vf2)*vstar
   25     FORMAT (' LOW KICK:    I NAM K* VK VF ',2i7,i4,2f7.2)
      END IF
*
*       Include optional list of high-velocity particles (double counting!).
*     IF (KZ(37).GT.0) THEN
*         CALL HIVEL(I)
*     END IF
*
   30 RETURN
*
      END