Nbody6
ksperi.f
Go to the documentation of this file.
1  SUBROUTINE ksperi(IPAIR)
2 *
3 *
4 * Pericentre KS variables.
5 * ------------------------
6 *
7  include 'common6.h'
8  common/slow0/ range,islow(10)
9 *
10 *
11 * Save current time and initialize integration indicator.
12  time0 = time
13  itime = 0
14  i1 = 2*ipair - 1
15  icm = n + ipair
16 *
17 * Re-calculate TDOT2 which may have misleading value.
18  1 td2 = 0.0d0
19  DO 5 k = 1,4
20  td2 = td2 + u0(k,ipair)*udot(k,ipair)
21  5 CONTINUE
22  td2 = 2.0*td2
23 *
24 * Define auxiliary quantities and obtain the eccentricity.
25  semi = -0.5d0*body(icm)/h(ipair)
26  zeta = 1.0 - r(ipair)/semi
27  psi = td2/sqrt(body(icm))
28  ecc = sqrt(zeta**2 + psi**2/semi)
29 *
30 * Avoid nearly circular orbits (undefined pericentre).
31  IF (ecc.LT.0.0001) go to 30
32 *
33 * Distinguish between near-collision, elliptic & hyperbolic case.
34  q = (psi/ecc)**2
35  IF (zeta.GT.0.0d0.AND.abs(q/semi).LT.0.1) THEN
36 * Expansion of THETA/(SIN(THETA)*COS(THETA)) - 1 for Q = SIN(THETA)**2.
37  sn1 = 1.0
38  s = 0.0d0
39  q = q/semi
40  DO 8 it = 1,10
41  sn1 = q*sn1*float(2*it)/float(2*it + 1)
42  s = s + sn1
43  8 CONTINUE
44  s = s + sn1*q*0.9d0/(1.0 - q)
45  dt = (r(ipair)*zeta - psi**2 + zeta*s*semi)*psi/
46  & (ecc**2*sqrt(body(icm)))
47  ELSE IF (semi.GT.0.0d0) THEN
48 * Determine the eccentric anomaly with respect to pericentre (0,PI).
49  theta = atan2(abs(psi)/sqrt(semi),zeta)
50 * Obtain pericentre time interval from Kepler's equation.
51  dt = semi*sqrt(semi/body(icm))*(theta - abs(psi)/sqrt(semi))
52  ELSE IF (semi.LT.0.0d0) THEN
53  a1 = psi/(ecc*sqrt(abs(semi)))
54 * Use EXP(F) = SINH(F) + COSH(F) to obtain the eccentric anomaly THETA.
55  a2 = abs(a1) + sqrt(a1**2 + 1.0d0)
56  theta = log(a2)
57  IF (a1.LT.0.0d0) theta = -theta
58  a0 = abs(semi)
59  dt = a0*sqrt(a0/body(icm))*(abs(psi)/sqrt(a0) - theta)
60  END IF
61 *
62 * Re-define current time (NB! not quantized; T0(I1) may be << TIME).
63  time = time - dt
64 *
65 * Integrate backwards for perturbed motion (reflection gives errors).
66  IF (list(1,i1).GT.0.AND.itime.EQ.0.AND.dt.GT.step(i1)) THEN
67  time = time0
68  imod = kslow(ipair)
69  zmod = float(islow(imod))
70  iph = iphase
71  iphase = -1
72 *
73 * Integrate step by step if interval is too large (note IPHASE < 0).
74  10 IF (dt.GT.step(i1)) THEN
75  time = time - step(i1)
76  dt = dt - step(i1)
77  h0(ipair) = h(ipair)
78  z = -0.5d0*h(ipair)*dtau(ipair)**2
79  CALL stumpf(ipair,z)
80  dtau(ipair) = -abs(dtau(ipair))
81  CALL ksint(i1)
82  dtu = dtau(ipair)
83 * Use negative DTU and treat STEP as positive (not used elsewhere).
84  step(i1) = ((one6*tdot3(ipair)*dtu + 0.5*tdot2(ipair))*dtu
85  & + r(ipair))*dtu
86  step(i1) = -zmod*step(i1)
87  itime = itime + 1
88  IF (itime.LT.200) go to 10
89  END IF
90 *
91  itime = itime + 1
92  dtu = dt/(r(ipair)*zmod)
93  dtu0 = dtau(ipair)
94  iter = 0
95 * Determine the regularized step by Newton-Raphson iteration (DT > 0).
96  20 y0 = dt - zmod*((one6*tdot3(ipair)*dtu +
97  & 0.5*tdot2(ipair))*dtu + r(ipair))*dtu
98  ypr = -((0.5*tdot3(ipair)*dtu + tdot2(ipair))*dtu + r(ipair))
99  ypr = zmod*ypr
100  dtu = dtu - y0/ypr
101  dt1 = ((one6*tdot3(ipair)*dtu + 0.5*tdot2(ipair))*dtu +
102  & r(ipair))*dtu
103  dt1 = zmod*dt1
104  iter = iter + 1
105  IF (abs(dt - dt1).GT.1.0d-10*step(i1).AND.iter.LT.5) go to 20
106 *
107 * Integrate back to pericentre using temporary indicator < 0 for exit.
108  time = time - dt
109  dtau(ipair) = -dtu
110 * Re-initialize Stumpff functions.
111  h0(ipair) = h(ipair)
112  z = -0.5d0*h(ipair)*dtau(ipair)**2
113  CALL stumpf(ipair,z)
114  CALL ksint(i1)
115 * Update energy and set positive step in case no KS initialization.
116  h0(ipair) = h(ipair)
117  dtau(ipair) = abs(dtu0)
118  step(i1) = zmod*dt
119  iphase = iph
120 * Use reflection procedure to improve provisional pericentre.
121  go to 1
122 * Note: typically 2 iterations and final TDOT2 > 0.
123  END IF
124 *
125 * Specify transformation coefficients (Seppo Mikkola's procedure).
126  IF (zeta.GE.0.0) THEN
127  xc = sqrt(0.5d0 + 0.5d0*zeta/ecc)
128  ys = psi/(ecc*xc*sqrt(body(icm)))
129  ELSE
130 * Employ well behaved expressions for R > A (SM 29/5/97).
131  xc = 0.5*abs(psi)/(sqrt(semi)*ecc)/sqrt(0.5d0-0.5d0*zeta/ecc)
132 * Avoid division by small XC near apocentre (ZETA < 0 only).
133  ys = 2.0*sqrt(semi/body(icm)*(0.5d0 - 0.5d0*zeta/ecc))
134  IF (psi.LT.0.0) ys = -ys
135  END IF
136 *
137  zz = body(icm)/(4.0*semi)
138  r(ipair) = 0.0d0
139  tdot2(ipair) = 0.0d0
140 *
141 * Generate analytical solutions for U & UDOT using old U0 & UDOT.
142  t0(i1) = time
143  DO 25 k = 1,4
144  u(k,ipair) = u0(k,ipair)*xc - udot(k,ipair)*ys
145  udot(k,ipair) = u0(k,ipair)*ys*zz + udot(k,ipair)*xc
146  u0(k,ipair) = u(k,ipair)
147  r(ipair) = r(ipair) + u(k,ipair)**2
148  tdot2(ipair) = tdot2(ipair) + 2.0d0*u(k,ipair)*udot(k,ipair)
149  25 CONTINUE
150 *
151 * Do not allow R' < 0 for repeated pericentre in routine KSINT.
152  IF (tdot2(ipair).LT.0.0d0) THEN
153  tdot2(ipair) = 0.0d0
154  END IF
155 *
156 * Predict c.m. coordinates & velocities.
157  IF (abs(dt).LT.step(icm).AND.time.GT.0.0) THEN
158  CALL xvpred(icm,0)
159  END IF
160 *
161  30 RETURN
162 *
163  END