Nbody6
orbit.f
Go to the documentation of this file.
1  SUBROUTINE orbit(I,J,SEMI,ECC,GI)
2 *
3 *
4 * Two-body elements.
5 * -----------------
6 *
7  include 'common6.h'
8 *
9 *
10 * Find the dominant neighbour if body #J not specified.
11  gi = 0.0
12  IF (j.LE.0) THEN
13  nnb = list(1,i)
14  fmax = 1.0d-10
15  jm = list(2,i)
16  DO 2 l = 1,nnb
17  jj = list(l+1,i)
18  jlist(l) = jj
19  rij2 = 0.0
20  DO 1 k = 1,3
21  rij2 = rij2 + (x(k,i) - x(k,jj))**2
22  1 CONTINUE
23  fij = (body(i) + body(jj))/rij2
24 * Exclude any c.m. bodies from dominant motion.
25  IF (fij.GT.fmax.AND.jj.LE.n) THEN
26  fmax = fij
27  jm = jj
28  END IF
29  2 CONTINUE
30 * Avoid rare case of halo orbit with zero neighbour number.
31  IF (jm.EQ.i) jm = i - 1
32 *
33 * Obtain the relative perturbation for decision-making.
34  CALL fpert(i,jm,nnb,pert)
35  gi = pert/fmax
36  j = jm
37  END IF
38 *
39 * Determine the semi-major axis and eccentricity.
40  rij2 = 0.0
41  vij2 = 0.0
42  rdot = 0.0
43  DO 5 k = 1,3
44  rij2 = rij2 + (x(k,i) - x(k,j))**2
45  vij2 = vij2 + (xdot(k,i) - xdot(k,j))**2
46  rdot = rdot + (x(k,i) - x(k,j))*(xdot(k,i) - xdot(k,j))
47  5 CONTINUE
48  rij = sqrt(rij2)
49  a1 = 2.0/rij - vij2/(body(i) + body(j))
50  semi = 1.0/a1
51  ecc2 = (1.0 - rij/semi)**2 + rdot**2/(semi*(body(i) + body(j)))
52  ecc = sqrt(ecc2)
53 *
54  RETURN
55 *
56  END