Nbody6
xtpert.f
Go to the documentation of this file.
1  SUBROUTINE xtpert(ACC,CHTIME)
2 *
3 *
4 * External perturbations on chain members.
5 * ----------------------------------------
6 *
7  include 'common6.h'
8  REAL*8 m,mass,mc,mij,mkk
9  parameter(nmx=10,nmx3=3*nmx,nmxm=nmx*(nmx-1)/2)
10  common/chain1/ xch(nmx3),vch(nmx3),m(nmx),
11  & zz(nmx3),wc(nmx3),mc(nmx),
12  & xj(nmx3),pi(nmx3),mass,rinv(nmxm),rsum,mkk(nmx),
13  & mij(nmx,nmx),tkk(nmx),tk1(nmx),iname(nmx),nn
14  common/chainc/ xc(3,ncmax),uc(3,ncmax),bodyc(ncmax),ich,
15  & listc(lmax)
16  common/clump/ bodys(ncmax,5),t0s(5),ts(5),steps(5),rmaxs(5),
17  & names(ncmax,5),isys(5)
18  REAL*8 acc(nmx3)
19 *
20 *
21 * Predict current coordinates of perturbers & chain members.
22  time0 = time
23  isub = isys(5)
24  time = chtime + t0s(isub)
25  CALL xcpred(1)
26 *
27 * Initialize the external perturbations.
28  nk = 3*nn
29  DO 1 k = 1,nk
30  acc(k) = 0.0d0
31  1 CONTINUE
32 *
33 * Consider each perturber in turn.
34  npc = listc(1) + 1
35  DO 20 l = 2,npc
36  j = listc(l)
37 * Sum perturber contributions over each chain component.
38  ik = -3
39  ks = 0
40  DO 10 i = 1,nn
41  ik = ik + 3
42  a1 = x(1,j) - xc(1,i)
43  a2 = x(2,j) - xc(2,i)
44  a3 = x(3,j) - xc(3,i)
45  rij2 = a1*a1 + a2*a2 + a3*a3
46  IF (j.LE.n) go to 5
47  jpair = j - n
48 * Check c.m. approximation (only resolve KS components once).
49  IF (rij2.GT.cmsep2*r(jpair)**2) go to 5
50  IF (ks.EQ.0) THEN
51  CALL ksres(jpair,j1,j2,rij2)
52  ks = 1
53  END IF
54  kdum = j1
55  k = kdum
56  4 a1 = x(1,k) - xc(1,i)
57  a2 = x(2,k) - xc(2,i)
58  a3 = x(3,k) - xc(3,i)
59  rij2 = a1*a1 + a2*a2 + a3*a3
60  a6 = body(k)/(rij2*sqrt(rij2))
61  acc(ik+1) = acc(ik+1) + a1*a6
62  acc(ik+2) = acc(ik+2) + a2*a6
63  acc(ik+3) = acc(ik+3) + a3*a6
64 * See whether the second component has been included.
65  IF (k.EQ.kdum) THEN
66  k = k + 1
67  go to 4
68  END IF
69  go to 10
70 *
71  5 a6 = body(j)/(rij2*sqrt(rij2))
72  acc(ik+1) = acc(ik+1) + a1*a6
73  acc(ik+2) = acc(ik+2) + a2*a6
74  acc(ik+3) = acc(ik+3) + a3*a6
75  10 CONTINUE
76  20 CONTINUE
77 *
78 * Restore the global time.
79  time = time0
80 *
81  RETURN
82  END