Nbody6
flyby.f
Go to the documentation of this file.
1  SUBROUTINE flyby(I,ITERM)
2 *
3 *
4 * Flyby check of KS orbit.
5 * -----------------------
6 *
7  include 'common6.h'
8 *
9 *
10 * Set index of KS pair & first component of c.m. body #I.
11  ipair = i - n
12  i1 = 2*ipair - 1
13  rjmin2 = 1000.0
14  iterm = 0
15  jcomp = list(2,i1)
16 *
17 * Find the closest body (JCOMP) and copy perturbers to JLIST.
18  nnb = list(1,i1)
19  DO 10 l = 2,nnb+1
20  j = list(l,i1)
21  jlist(l-1) = j
22  rij2 = (x(1,i) - x(1,j))**2 + (x(2,i) - x(2,j))**2 +
23  & (x(3,i) - x(3,j))**2
24  IF (rij2.LT.rjmin2) THEN
25  rjmin2 = rij2
26  jcomp = j
27  END IF
28  10 CONTINUE
29 *
30 * Skip rare case of merged binary or chain c.m. (denoted by NAME <= 0).
31  IF (name(i).LE.0.OR.name(jcomp).LE.0) go to 20
32 *
33  rdot = (x(1,i) - x(1,jcomp))*(xdot(1,i) - xdot(1,jcomp)) +
34  & (x(2,i) - x(2,jcomp))*(xdot(2,i) - xdot(2,jcomp)) +
35  & (x(3,i) - x(3,jcomp))*(xdot(3,i) - xdot(3,jcomp))
36 *
37  vrel2 = (xdot(1,i) - xdot(1,jcomp))**2 +
38  & (xdot(2,i) - xdot(2,jcomp))**2 +
39  & (xdot(3,i) - xdot(3,jcomp))**2
40
41 * Evaluate outer semi-major axis, eccentricity & impact parameter.
42  rij = sqrt(rjmin2)
43  semi1 = 2.0/rij - vrel2/(body(i) + body(jcomp))
44  semi1 = 1.0/semi1
45  ecc1 = sqrt((1.0d0 - rij/semi1)**2 +
46  & rdot**2/(semi1*(body(i) + body(jcomp))))
47  pmin = semi1*(1.0d0 - ecc1)
48 *
49 * Set semi-major axis & eccentricity of inner binary.
50  semi = -0.5d0*body(i)/h(ipair)
51  ecc2 = (1.0d0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(body(i)*semi)
52  ecc = sqrt(ecc2)
53 *
54  i2 = i1 + 1
55  rji1 = (x(1,jcomp) - x(1,i1))**2 + (x(2,jcomp) - x(2,i1))**2 +
56  & (x(3,jcomp) - x(3,i1))**2
57  rji2 = (x(1,jcomp) - x(1,i2))**2 + (x(2,jcomp) - x(2,i2))**2 +
58  & (x(3,jcomp) - x(3,i2))**2
59 *
60 * Identify the dominant interaction (#I1 or #I2).
61  fji1 = (body(i1) + body(jcomp))/rji1
62  fji2 = (body(i2) + body(jcomp))/rji2
63  nnb = nnb + 1
64  IF (fji1.GT.fji2) THEN
65  j = i1
66  rj2 = rji1
67 * Add the other component to the perturber list.
68  jlist(nnb) = i2
69  ELSE
70  j = i2
71  rj2 = rji2
72  jlist(nnb) = i1
73  END IF
74 *
75 * Evaluate radial velocity (routine SEARCH requires RD < 0).
76  rd = 0.0
77  rjj = 0.0
78  vjj = 0.0
79  DO 15 k = 1,3
80  rd = rd + (x(k,jcomp) - x(k,j))*(xdot(k,jcomp) - xdot(k,j))
81  rjj = rjj + (x(k,jcomp) - x(k,j))**2
82  vjj = vjj + (xdot(k,jcomp) - xdot(k,j))**2
83  15 CONTINUE
84 *
85 * Determine vectorial perturbation on intruder and closest component.
86  CALL fpert(jcomp,j,nnb,pert)
87  gj = pert*rj2/(body(jcomp) + body(j))
88  rj = sqrt(rjj)
89 * WRITE (6,18) NAME(J),NAME(JCOMP),GAMMA(IPAIR),GJ,R(IPAIR),RJ,RD,
90 * & RDOT
91 * 18 FORMAT (' NMJ NMJC GI GJ R RJ RD RDI ',2I5,2F6.2,1P,4E9.1)
92 *
93 * Terminate if perturber & I1/I2 provides dominant force.
94  apo = abs(semi)*(1.0 + ecc)
95  rsum = r(ipair) + sqrt(rj2)
96  xf = 2.0
97 * Compare current radial velocity with #JCOMP & J.
98  IF (tdot2(ipair)/r(ipair).LT.rd/rj) xf = 1.0
99 *
100 * Increase the cross section for quadruples.
101  IF (jcomp.GT.n) THEN
102  semi2 = -0.5*body(jcomp)/h(jcomp-n)
103  apo = apo + abs(semi2)
104  END IF
105 *
106  semij = 2.0/rj - vjj/(body(jcomp) + body(j))
107  semij = 1.0/semij
108 * Check for chain regularization test or standard termination.
109  IF (pmin.LT.apo.AND.rsum.LT.xf*rmin.AND.rdot.LT.0.0) THEN
110  IF (semi.LT.0.0.AND.gamma(ipair).LT.0.5) THEN
111  iterm = 0
112  go to 20
113  END IF
114  iterm = 1
115  ELSE IF (gj.LT.0.7*gamma(ipair).AND.rd.LE.0.0) THEN
116  iterm = 2
117  END IF
118 *
119  20 RETURN
120 *
121  END