Nbody6
 All Files Functions Variables
ksmod.f
Go to the documentation of this file.
1  SUBROUTINE ksmod(IPAIR,KMOD)
2 *
3 *
4 * Modified KS motion.
5 * -------------------
6 *
7  include 'common6.h'
8  common/slow0/ range,islow(10)
9 *
10 *
11 * Set current modification level.
12  imod = kslow(ipair)
13 *
14 * Check transition type (KMOD <= 1 denotes standard restart).
15  IF (kmod.GT.1) THEN
16 * Determine provisional index for KS slow-down.
17  DO 5 k = 2,10
18  isbin = k - 1
19  IF (islow(k).GT.kmod) go to 10
20  5 CONTINUE
21 * Restrict increase to two levels.
22  10 isbin = min(isbin,imod+2)
23 * See whether standard solution is called for.
24  IF (isbin.EQ.1) go to 30
25  ELSE
26 * Include termination (IMOD > 1 & KMOD <= 1).
27  isbin = 1
28  go to 30
29  END IF
30 *
31 * Estimate time interval to reach largest permitted perturbation.
32  gx = range*gmin
33  CALL tpert(ipair,gx,dt)
34 *
35 * Evaluate the unmodified Kepler period.
36  icm = n + ipair
37  semi = -0.5*body(icm)/h(ipair)
38  tk = twopi*semi*sqrt(semi/body(icm))
39 *
40 * Reduce level if modification factor is too large.
41  DO 20 k = 2,10
42  IF (tk*float(islow(isbin)).LT.dt.OR.isbin.EQ.1) go to 30
43  isbin = isbin - 1
44  20 CONTINUE
45 *
46 * Exit if level is unchanged.
47  30 IF (isbin.EQ.imod) go to 100
48 *
49 * Ensure the Kepler period is defined for transition to standard case.
50  IF (isbin.EQ.1) THEN
51  icm = n + ipair
52  semi = -0.5*body(icm)/h(ipair)
53  tk = twopi*semi*sqrt(semi/body(icm))
54  END IF
55 *
56 * Define auxiliary quantities.
57  zeta = 1.0 - r(ipair)/semi
58  psi = tdot2(ipair)/sqrt(body(icm))
59 *
60 * Determine the eccentric anomaly with respect to pericentre (-PI,PI).
61  theta = atan2(psi/sqrt(semi),zeta)
62 *
63 * Obtain apocentre time interval from Kepler's equation and the period.
64  dt = semi*sqrt(semi/body(icm))*(theta - psi/sqrt(semi))
65  dt = 0.5d0*tk + dt
66  dt = -dt
67 *
68 * Evaluate regularized apocentre time (Baumgarte & Stiefel, 1974).
69 * DTU = -2.0D0*(H(IPAIR)*DT + TDOT2(IPAIR))/BODY(ICM)
70 *
71 * Determine the regularized step by Newton-Raphson iteration (DT < 0).
72  i1 = 2*ipair - 1
73  dtu = dt/r(ipair)
74  iter = 0
75 * Note: explicit relation agrees with iterated value (bug fix 9/99).
76  40 y0 = dt - ((one6*tdot3(ipair)*dtu +
77  & 0.5*tdot2(ipair))*dtu + r(ipair))*dtu
78  ypr = -((0.5*tdot3(ipair)*dtu + tdot2(ipair))*dtu + r(ipair))
79  dtu = dtu - y0/ypr
80  dt1 = ((one6*tdot3(ipair)*dtu + 0.5*tdot2(ipair))*dtu +
81  & r(ipair))*dtu
82  iter = iter + 1
83  IF (abs(dt - dt1).GT.1.0d-10*step(i1).AND.iter.LT.5) go to 40
84 *
85 * Reset reference energy and generate new Stumpff coefficients.
86  h0(ipair) = h(ipair)
87  z = -0.5*h(ipair)*dtu**2
88  CALL stumpf(ipair,z)
89 *
90 * Integrate small step back to apocentre using temporary indicator # 0.
91  dtau(ipair) = dtu
92  time = t0(i1) + dt
93  iphase = -1
94  CALL ksint(i1)
95 *
96 * Predict current coordinates & velocities of ICM.
97  CALL xvpred(icm,0)
98 *
99 * Set new KS level and increase restart counter (ISBIN > 1).
100  kslow(ipair) = isbin
101  IF (isbin.GT.1) nksmod = nksmod + 1
102 *
103 * Perform perturbed restart of KS motion.
104  CALL kspoly(ipair,isbin)
105 *
106 * Ensure that apocentre criterion will fail after next step.
107  tdot2(ipair) = 0.0d0
108 *
109 * Set indicator = -1 to skip perturber selection in routine KSINT.
110  kmod = -1
111  iphase = 0
112 *
113 * Enforce new KSLIST after small perturbation (CHAOS/SPIRAL case).
114  IF (isbin.GE.6) kmod = 0
115 *
116  100 RETURN
117 *
118  END