Nbody6
 All Files Functions Variables
subsys.f
Go to the documentation of this file.
1  SUBROUTINE subsys(NSYS,CM)
2 *
3 *
4 * Initialization of subsystem.
5 * ----------------------------
6 *
7  include 'common6.h'
8  common/clump/ bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
9  & names(ncmax,5),isys(5)
10  REAL*8 cm(7)
11 *
12 *
13 * Increase subsystem counter and set current index.
14  nsub = nsub + 1
15  isub = nsub
16 *
17 * Set zero name & mass in last component to distinguish triple case.
18  IF (nsys.EQ.3) THEN
19  names(4,isub) = 0
20  bodys(4,isub) = 0.0d0
21  END IF
22 *
23 * Save global masses & names of new subsystem components.
24  bx = 0.0
25  DO 10 l = 1,nsys
26  j = jlist(l)
27  bodys(l,isub) = body(j)
28  names(l,isub) = name(j)
29  IF (body(j).GT.bx) THEN
30  bx = body(j)
31  lx = l
32  END IF
33  10 CONTINUE
34 *
35 * Form ghosts and initialize integration variables (NSYS = 3 or 4).
36  DO 20 l = 1,nsys
37  j = jlist(l)
38  body(j) = 0.0d0
39  t0(j) = 1.0e+06
40  list(1,j) = 0
41  DO 15 k = 1,3
42  x0dot(k,j) = 0.0d0
43  xdot(k,j) = 0.0d0
44  f(k,j) = 0.0d0
45  fdot(k,j) = 0.0d0
46  d0(k,j) = 0.0d0
47  d1(k,j) = 0.0d0
48  d2(k,j) = 0.0d0
49  d3(k,j) = 0.0d0
50  d0r(k,j) = 0.0d0
51  d1r(k,j) = 0.0d0
52  d2r(k,j) = 0.0d0
53  d3r(k,j) = 0.0d0
54  15 CONTINUE
55 * Set large X0 & X to avoid perturber selection (no escape).
56  x0(1,j) = 1.0e+06
57  x(1,j) = 1.0e+06
58  20 CONTINUE
59 *
60 * Place c.m. of subsystem in first location (ICOMP may switch!).
61  icomp = jlist(1)
62 * Ensure heaviest body is selected in case of ARchain.
63  IF (isys(isub).EQ.4) THEN
64  icomp = jlist(lx)
65  END IF
66 * Define zero name for identification (only chain or ARchain c.m.).
67  IF (isys(isub).GE.3) name(icomp) = 0
68 *
69 * Initialize c.m in ICOMP (redefined above).
70  t0(icomp) = time
71  body(icomp) = cm(7)
72  DO 30 k = 1,3
73  x(k,icomp) = cm(k)
74  x0(k,icomp) = cm(k)
75  xdot(k,icomp) = cm(k+3)
76  x0dot(k,icomp) = cm(k+3)
77  30 CONTINUE
78 *
79 * Predict coordinates & velocities for all neighbours (order FDOT).
80  nnb = list(1,icomp)
81  CALL xvpred(icomp,nnb)
82 *
83 * Obtain new neighbour list (to ensure membership > 0).
84  rs0 = rs(icomp)
85  CALL nblist(icomp,rs0)
86 *
87 * Construct force polynomials for c.m. motion.
88  CALL fpoly1(icomp,icomp,0)
89  CALL fpoly2(icomp,icomp,0)
90 *
91 * Initialize decision-making variables for multiple regularization.
92  t0s(isub) = time
93  ts(isub) = time
94  steps(isub) = step(icomp)
95 *
96 * Obtain maximum unperturbed separation based on dominant neighbour.
97  CALL extend(isub)
98 *
99  RETURN
100 *
101  END