Nbody6
 All Files Functions Variables
ksin2.f
Go to the documentation of this file.
1  SUBROUTINE ksin2(ICASE)
2 *
3 *
4 * Initialization of hierarchical KS.
5 * ----------------------------------
6 *
7  include 'common6.h'
8  REAL*8 q(3),rdot(3),ui(4),vi(4),a1(3,4)
9  CHARACTER*8 which1
10 *
11 *
12 * Copy KS pair index and set c.m. index.
13  ipair = kspair
14  icm = n + ipair
15 *
16 * Distinguish first or second stage (termination: ICASE > 1 only).
17  IF (icase.GT.1) go to 50
18 *
19 * Specify mass & name for new c.m. and initialize radius & type.
20  bodyi = body(icomp)
21  body(icm) = body(icomp) + body(jcomp)
22  name(icm) = nzero + name(icomp)
23  radius(icm) = 0.0
24  tev(icm) = 1.0e+10
25  tev0(icm) = 1.0e+10
26  body0(icm) = bodym
27  epoch(icm) = time*tstar
28  kstar(icm) = 0
29 *
30 * Define relative coordinates and velocities in physical units.
31  DO 20 k = 1,3
32  q(k) = x(k,icomp) - x(k,jcomp)
33  rdot(k) = x0dot(k,icomp) - x0dot(k,jcomp)
34  20 CONTINUE
35 *
36 * Introduce regularized variables using definition of 1985 paper.
37  r(ipair) = sqrt(q(1)**2 + q(2)**2 + q(3)**2)
38 *
39 * Initialize the regularized coordinates according to sign of Q(1).
40  IF (q(1).LE.0.0d0) THEN
41  ui(3) = 0.0d0
42  ui(2) = sqrt(0.5d0*(r(ipair) - q(1)))
43  ui(1) = 0.5d0*q(2)/ui(2)
44  ui(4) = 0.5d0*q(3)/ui(2)
45  ELSE
46  ui(4) = 0.0d0
47  ui(1) = sqrt(0.5d0*(r(ipair) + q(1)))
48  ui(2) = 0.5d0*q(2)/ui(1)
49  ui(3) = 0.5d0*q(3)/ui(1)
50  END IF
51 *
52 * Set current transformation matrix.
53  CALL matrix(ui,a1)
54 *
55 * Form regularized velocity and set initial KS coordinates & TDOT2.
56  tdot2(ipair) = 0.0d0
57  DO 30 k = 1,4
58  udot(k,ipair) = 0.50d0*(a1(1,k)*rdot(1) + a1(2,k)*rdot(2) +
59  & a1(3,k)*rdot(3))
60 * Note that A1(J,K) is the transpose of A1(K,J).
61  u(k,ipair) = ui(k)
62  u0(k,ipair) = u(k,ipair)
63  tdot2(ipair) = tdot2(ipair) + 2.0d0*ui(k)*udot(k,ipair)
64  30 CONTINUE
65 *
66 * Evaluate initial binding energy per unit mass and initialize T0.
67  h(ipair) = (2.0d0*(udot(1,ipair)**2 + udot(2,ipair)**2 +
68  & udot(3,ipair)**2 + udot(4,ipair)**2) -
69  & body(icm))/r(ipair)
70  t0(2*ipair-1) = time
71 *
72 * Copy old c.m. mass to first KS component.
73  body(2*ipair-1) = bodyi
74 *
75 * Define c.m. coordinates & velocities.
76  DO 40 k = 1,3
77  x(k,icm) = (bodyi*x(k,icomp) + body(jcomp)*x(k,jcomp))/
78  & body(icm)
79  x0dot(k,icm) = (bodyi*x0dot(k,icomp) + body(jcomp)*
80  & x0dot(k,jcomp))/body(icm)
81  xdot(k,icm) = x0dot(k,icm)
82  x0(k,icm) = x(k,icm)
83  40 CONTINUE
84 *
85 * Exit after initialization.
86  go to 100
87 *
88 * Check for sufficient perturbers (RMAX = apocentre set in IMPACT).
89  50 nnb = list(1,icm)
90  IF (sqrt(cmsep2)*rmax.GT.2.0*rs(icm)) THEN
91  fac = float(nnbmax)/float(nnb)
92  rsi = fac**0.3333*rs(icm)
93  IF (rsi.GT.rs(icm)) THEN
94  CALL nblist(icm,rsi)
95  END IF
96  END IF
97 *
98 * Initialize force polynomial for c.m. using new neighbour list.
99  CALL fpoly1(icm,icm,0)
100  CALL fpoly2(icm,icm,0)
101 *
102 * Form perturber list.
103  CALL kslist(ipair)
104 *
105 * Transform any unperturbed hard binary to apocentre and set time-step.
106  imod = 1
107  eb = h(ipair)*body(icomp)*body(jcomp)/body(icm)
108 * Suppress the following for now.
109  IF (list(1,icomp).LT.0.AND.eb.LT.ebh) THEN
110  semi = -0.5*body(icm)/h(ipair)
111  tk = twopi*abs(semi)*sqrt(abs(semi)/body(icm))
112  IF (iphase.NE.7) THEN
113  DO 55 k = 1,4
114  ui(k) = u(k,ipair)
115  vi(k) = udot(k,ipair)
116  55 CONTINUE
117 * Determine pericentre time (TP < 0 if TDOT2 < 0) and add TK/2.
118  CALL tperi(semi,ui,vi,body(icm),tp)
119  step(icomp) = 0.5*min(tk,step(icm)) - tp
120 * Transform KS variables to peri and by pi/2 to apocentre (skip apo).
121  IF (abs(tdot2(ipair)).GT.1.0e-12.OR.r(ipair).LT.semi) THEN
122  CALL ksperi(ipair)
123  CALL ksapo(ipair)
124  ELSE IF (tdot2(ipair).GT.0.0) THEN
125  tdot2(ipair) = -1.0e-20
126  END IF
127  END IF
128 * Estimate an appropriate KS slow-down index for G < GMIN.
129  IF (kz(26).GT.0.AND.step(icm).GT.tk) THEN
130  imod = 1 + log(step(icm)/tk)/0.69
131  imod = min(imod,5)
132  END IF
133  END IF
134 *
135 * Specify zero membership and large step for second component.
136  list(1,jcomp) = 0
137  step(jcomp) = 1.0e+06
138  stepr(jcomp) = 1.0e+06
139 *
140 * Obtain polynomials for perturbed KS motion (standard case).
141  CALL kspoly(ipair,imod)
142 *
143 * Increase regularization counters (NKSHYP for hyperbolic orbits).
144  nksreg = nksreg + 1
145 *
146  IF (kz(10).GT.0.AND.icase.GE.2) THEN
147  ri = sqrt((x(1,icm) - rdens(1))**2 +
148  & (x(2,icm) - rdens(2))**2 +
149  & (x(3,icm) - rdens(3))**2)
150  which1 = ' MERGE2 '
151  IF (icase.EQ.3) which1 = ' RESET2 '
152  WRITE (6,60) which1, time+toff, name(icomp), name(jcomp),
153  & dtau(ipair), r(ipair), ri, h(ipair), ipair,
154  & gamma(ipair), step(icm), list(1,icomp)
155  60 FORMAT (/,' NEW',a8,' T =',f8.2,2i6,f12.3,1p,e10.1,0p,
156  & f8.2,f9.2,i5,f8.3,1p,e10.1,0p,i5)
157  END IF
158 *
159 * See whether either component has been regularized recently.
160  nnb = listd(1) + 1
161  k = 0
162 * Check case of initial binary and loop over disrupted pairs.
163  IF (iabs(name(icomp) - name(jcomp)).EQ.1) k = -1
164  DO 80 l = 2,nnb
165  IF (name(icomp).EQ.listd(l).OR.name(jcomp).EQ.listd(l)) k = -1
166  80 CONTINUE
167 *
168 * Ensure that mergers are treated as new binaries.
169  IF (iphase.EQ.6) k = 0
170 * Set flags to distinguish primordial binaries & standard KS motion.
171  list(2,jcomp) = k
172  kslow(ipair) = 1
173 *
174 * Check diagnostic output of new hard binary.
175  IF (kz(8).GT.0.AND.k.EQ.0) THEN
176  IF (eb.GT.ebh) go to 100
177  semi = -0.5*body(icm)/h(ipair)
178  ri = sqrt((x(1,icm) - rdens(1))**2 +
179  & (x(2,icm) - rdens(2))**2 +
180  & (x(3,icm) - rdens(3))**2)
181  WRITE (8,90) time+toff, name(icomp), name(jcomp), k,
182  & body(icomp), body(jcomp), eb, semi, r(ipair),
183  & gamma(ipair), ri
184  90 FORMAT (' NEW BINARY T =',f7.1,' NAME = ',2i6,i3,
185  & ' M =',2f8.4,' EB =',f9.4,' A =',f8.5,
186  & ' R =',f8.5,' G =',f6.3,' RI =',f5.2)
187  END IF
188 *
189  100 RETURN
190 *
191  END