Nbody6
kick.f
Go to the documentation of this file.
1  SUBROUTINE kick(I,ICASE,KW,DM)
2 *
3 *
4 * Velocity kick for WD, neutron stars or black holes.
5 * ---------------------------------------------------
6 *
7  include 'common6.h'
8  parameter(vfac=2.0d0)
9  REAL*8 ran2,vk(4)
10  SAVE ipair, kc, vdis, ri
11  DATA ipair,kc /0,0/
12 * WD velocity kicks: Fellhauer et al. 2003, ApJ Letters 595, L53.
13 *
14 *
15 * Save orbital parameters in case of KS binary (called from KSAPO).
16  IF (icase.EQ.0) THEN
17  ipair = i
18  kc = kstar(n+ipair)
19 * Identify the correct component (KSTAR reversed in MDOT).
20  i1 = 2*ipair - 1
21  IF (kstar(i1).LT.0) THEN
22  in = i1
23  ELSE IF (kstar(i1+1).LT.0) THEN
24  in = i1 + 1
25  END IF
26  kstar(in) = -kstar(in)
27 *
28 * Determine disruption velocity after mass loss.
29  vd2 = 2.0*body(n+ipair)/r(ipair)
30  vdis = sqrt(vd2)*vstar
31 * Set cluster escape velocity (add twice central potential).
32  vp2 = 2.0*body(n+ipair)/r(ipair)
33  vesc = sqrt(vp2 + 4.0)*vstar
34  semi = -0.5*body(n+ipair)/h(ipair)
35  zm1 = body(i1)*smu
36  zm2 = body(i1+1)*smu
37 * Skip case of massless binary.
38  IF (body(n+ipair).GT.0.0d0) THEN
39  eb = body(i1)*body(i1+1)/body(n+ipair)*h(ipair)
40  ELSE
41  eb = 0.0
42  END IF
43  ri = r(ipair)
44 * Skip on #25 = 0/1 for consistent net WD modification of EKICK.
45  IF ((kw.LT.13.AND.kz(25).EQ.0).OR.
46  & (kw.EQ.12.AND.kz(25).NE.2)) go to 30
47 * Sum whole binding energy (used by BINOUT for net change).
48  ekick = ekick + eb
49  egrav = egrav + eb
50  i2 = i1 + 1
51  WRITE (6,1) name(i1), name(i2), kstar(i1), kstar(i2), zm1,
52  & zm2, vesc, vdis, r(ipair)/semi, eb, r(ipair)
53  1 FORMAT (' BINARY KICK: NAM K* M1 M2 VESC VDIS R/A EB R ',
54  & 2i6,2i4,4f7.1,f6.2,f9.4,1p,e10.2)
55  nbkick = nbkick + 1
56 * Remove any circularized KS binary from the CHAOS/SYNCH table.
57  IF (kstar(n+ipair).GE.10.AND.nchaos.GT.0) THEN
58  ii = -(n + ipair)
59  CALL spiral(ii)
60  kstar(n+ipair) = 0
61  END IF
62  go to 30
63  END IF
64 *
65 * Generate velocity kick for neutron star (Gordon Drukier Tokyo paper).
66 * IT = 0
67 * V0 = 330.0
68 * VCUT = 1000.0
69 * 2 VT = VCUT/V0*RAN2(IDUM1)
70 * VP = VT*(2.0*RAN2(IDUM1) - 1.0)
71 * VN = SQRT(VT**2 - VP**2)
72 * FAC = 1.0/0.847*VN**0.3/(1.0 + VN**3.3)
73 * IF (FAC.LT.RAN2(IDUM1).AND.IT.LT.10) GO TO 2
74 * VKICK = V0*VT
75 *
76 * Adopt the Maxwellian of Hansen & Phinney (MN 291, 569, 1997).
77  disp = 190.0
78 *
79 * Include velocity dispersion in terms of VSTAR (Parameter statement!).
80  IF (vfac.GT.0.0d0) THEN
81  disp = vfac*vstar
82  END IF
83 *
84 * Allow for optional type-dependent WD kick.
85  IF (kw.LT.13) THEN
86  IF (kz(25).EQ.1.AND.(kw.EQ.10.OR.kw.EQ.11)) THEN
87  disp = 5.0
88  ELSE IF (kz(25).GT.1.AND.kw.EQ.12) THEN
89  disp = 5.0
90  ELSE
91  disp = 0.0
92  END IF
93  END IF
94 * IF (KW.EQ.14) DISP = 0.0
95 *
96 * Use Henon's method for pairwise components (Douglas Heggie 22/5/97).
97  DO 2 k = 1,2
98  x1 = ran2(idum1)
99  x2 = ran2(idum1)
100 * Generate two velocities from polar coordinates S & THETA.
101  s = disp*sqrt(-2.0*log(1.0 - x1))
102  theta = twopi*x2
103  vk(2*k-1) = s*cos(theta)
104  vk(2*k) = s*sin(theta)
105  2 CONTINUE
106  vkick = sqrt(vk(1)**2 + vk(2)**2 + vk(3)**2)
107  vk(4) = vkick
108 *
109 * Limit kick velocity to VDIS+10*VSTAR/10*VST for binary/single stars.
110  IF (ipair.GT.0) THEN
111  vbf = sqrt(vdis**2 + 100.0*vstar**2)
112  vkick = min(vkick,vbf)
113  ELSE
114  vkick = min(vkick,10.0d0*vstar)
115 * Ensure escape of massless star.
116  IF (body(i).EQ.0.0d0) vkick = 10.0*vstar
117  END IF
118  vkick = vkick/vstar
119 * Skip WD case with KZ(25) = 0.
120  IF (vkick.EQ.0.0d0) go to 30
121 *
122 * Randomize the velocity components.
123 * A(4) = 0.0
124 * DO 5 K = 1,3
125 * A(K) = 2.0*RAN2(IDUM1) - 1.0
126 * A(4) = A(4) + A(K)**2
127 * 5 CONTINUE
128 *
129 * Add truncated/full kick velocity and initialize X0DOT.
130  vi2 = 0.0
131  vf2 = 0.0
132  DO 10 k = 1,3
133  vi2 = vi2 + xdot(k,i)**2
134 * XDOT(K,I) = XDOT(K,I) + VKICK*A(K)/SQRT(A(4))
135  xdot(k,i) = xdot(k,i) + vkick*vk(k)/vk(4)
136  x0dot(k,i) = xdot(k,i)
137  vf2 = vf2 + xdot(k,i)**2
138  10 CONTINUE
139 *
140 * Modify energy loss due to increased velocity of single particle.
141  ecdot = ecdot - 0.5*body(i)*(vf2 - vi2)
142  nkick = nkick + 1
143 *
144 * Replace final velocity by relative velocity for binary kick.
145  IF (ipair.GT.0) THEN
146  jp = kvec(i)
147  j = i + 1
148  IF (i.EQ.2*jp) j = i - 1
149  vf2 = 0.0
150  DO 15 k = 1,3
151  vf2 = vf2 + (xdot(k,i) - xdot(k,j))**2
152  15 CONTINUE
153  hnew = 0.5*vf2 - (body(i) + body(j))/ri
154  eb1 = body(i)*body(j)/(body(i) + body(j))*hnew
155  IF (eb1.LT.0.0) THEN
156  ekick = ekick - eb1
157  egrav = egrav - eb1
158  END IF
159  ipair = 0
160  END IF
161 *
162  IF (nkick.LT.50.OR.name(i).LE.2*nbin0.OR.
163  & (kw.GE.13.AND.ttot*tstar.GT.100.0)) THEN
164  zm = body(i)*zmbar
165  WRITE (6,20) i, name(i), kstar(i), kc, body0(i)*zmbar, zm,
166  & sqrt(vi2)*vstar, vkick*vstar, sqrt(vf2)*vstar
167  20 FORMAT (' VELOCITY KICK: I NAM K* KC* M0 M VI VK VF ',
168  & 2i6,2i4,2f7.2,3f7.1)
169  kc = 0
170  END IF
171 *
172 * Highlight BH/NS velocities below 4 times rms velocity.
173  IF (vkick.LT.4.0*sqrt(0.5).AND.kw.GE.13) THEN
174  WRITE (6,25) i, name(i), kw, vkick*vstar, sqrt(vf2)*vstar
175  25 FORMAT (' LOW KICK: I NAM K* VK VF ',2i7,i4,2f7.2)
176  END IF
177 *
178 * Include optional list of high-velocity particles (double counting!).
179 * IF (KZ(37).GT.0) THEN
180 * CALL HIVEL(I)
181 * END IF
182 *
183  30 RETURN
184 *
185  END