Nbody6
 All Files Functions Variables
ksrect.f
Go to the documentation of this file.
1  SUBROUTINE ksrect(IPAIR)
2 *
3 *
4 * Rectification of KS orbit.
5 * --------------------------
6 *
7  include 'common6.h'
8 *
9 *
10 * Include some diagnostic output.
11  upr2 = 0.0
12  DO 1 k = 1,4
13  upr2 = upr2 + udot(k,ipair)**2
14  1 CONTINUE
15 *
16 * Skip rectification for small eccentricity or large perturbation.
17  i = n + ipair
18  semi = -0.5*body(i)/h(ipair)
19  ecc2 = (1.0 - r(ipair)/semi)**2 + tdot2(ipair)**2/(semi*body(i))
20  ecc = sqrt(ecc2)
21  IF (ecc.LE.0.01) go to 50
22  IF (gamma(ipair).GT.0.1) go to 50
23 *
24 * Include Roche rectification for large ECC (additional to roche.f).
25  IF (kstar(i).GT.10.AND.ecc.GT.0.2) THEN
26  qp = semi*(1.0 - ecc)
27  icirc = -1
28  CALL tcirc(qp,ecc,2*ipair-1,2*ipair,icirc,tc)
29  WRITE (6,2) name(2*ipair-1), kstar(i), ecc, qp, tc,
30  & gamma(ipair)
31  2 FORMAT (' RECTIFY K* NM K* E QP TC G',i8,i4,f8.4,1p,3e10.2)
32  kstar(i) = 0
33  END IF
34 *
35  hi = (2.0*upr2 - body(n+ipair))/r(ipair)
36  err = (hi - h(ipair))/hi
37  zmu = body(2*ipair)*body(2*ipair-1)/body(n+ipair)
38  db = zmu*(hi - h(ipair))
39  IF (abs(db).GT.1.0d-08) THEN
40  semi = -0.5*body(n+ipair)/h(ipair)
41  ra = r(ipair)/semi
42  IF (semi.LT.0.0) ra = r(ipair)
43  WRITE (16,3) time+toff, ipair, ra, h(ipair), gamma(ipair),
44  & db, err
45  3 FORMAT (' KSRECT: T # R/A H G DB DH/H ',
46  & f8.2,i4,f8.4,f8.1,f7.3,1p,2e10.1)
47  CALL flush(16)
48  END IF
49 *
50 * Initialize iteration counter for difficult case (SJA 10/97).
51  iter = 0
52 *
53 * Form square regularized velocity for the explicit binding energy.
54  10 upr2 = 0.0
55  DO 15 k = 1,4
56  upr2 = upr2 + udot(k,ipair)**2
57  15 CONTINUE
58 *
59 * Form KS scaling factors from energy and angular momentum relation.
60  a1 = 0.25d0*body(n+ipair)/upr2
61 * Solve for C1 from H = (2*U'*U'*C1**2 - M)/(U*U*C2**2) with C2 = 1/C1.
62  a2 = a1**2 + 0.5d0*h(ipair)*r(ipair)/upr2
63 *
64 * Avoid negative round-off value on second try (NB! no change in CK).
65  IF (iter.EQ.2.AND.a2.LT.0.0) a2 = 0.0d0
66 *
67 * Check for undefined case (circular orbit or eccentric anomaly = 90).
68  IF (a2.GE.0.0d0) THEN
69  IF (a1.LT.1.0) THEN
70 * Choose square root sign from eccentric anomaly (e*cos(E) = 1 - R/a).
71  c1 = sqrt(a1 + sqrt(a2))
72  ELSE
73  c1 = sqrt(a1 - sqrt(a2))
74  END IF
75  ck = 1.0
76  ELSE
77 * Adopt C1*C2 = CK for difficult case (Seppo's suggestion of 1991).
78  c1 = 1.0
79  ck = body(n+ipair)/sqrt(-8.0d0*h(ipair)*r(ipair)*upr2)
80  WRITE (6,20) ipair, kstar(n+ipair), ecc, r(ipair), h(ipair),
81  & gamma(ipair), upr2, a2, ck-1.0
82  20 FORMAT (' WARNING! KSRECT KS K* E R H G UPR2 A2 CK-1 ',
83  & 2i4,f8.4,1p,6e10.2)
84  iter = iter + 1
85  END IF
86 *
87 * Specify KS coordinate scaling from angular momentum conservation.
88  c2 = ck/c1
89 *
90 * Transform KS variables to yield the prescribed elements.
91  r(ipair) = 0.0d0
92 * UPR2 = 0.0D0
93  td2 = 0.0d0
94  DO 25 k = 1,4
95  u(k,ipair) = c2*u(k,ipair)
96  udot(k,ipair) = c1*udot(k,ipair)
97  u0(k,ipair) = u(k,ipair)
98  r(ipair) = r(ipair) + u(k,ipair)**2
99  td2 = td2 + u(k,ipair)*udot(k,ipair)
100 * UPR2 = UPR2 + UDOT(K,IPAIR)**2
101  25 CONTINUE
102  tdot2(ipair) = 2.0*td2
103 *
104 * Include diagnostic output.
105 * HI = (2.0*UPR2 - BODY(N+IPAIR))/R(IPAIR)
106 * WRITE (16,30) IPAIR, IPHASE, KSTAR(N+IPAIR), R(IPAIR),
107 * & GAMMA(IPAIR), (HI-H(IPAIR))/HI
108 * 30 FORMAT (' KSRECT: KS IPH K* R G DH/H ',3I4,1P,3E10.2)
109 * CALL FLUSH(16)
110 *
111 * Improve solution by second iteration in case of CK procedure.
112  iter = iter + 1
113  IF (iter.EQ.2) go to 10
114 *
115  50 RETURN
116 *
117  END