Nbody6
freeze.f
Go to the documentation of this file.
1  SUBROUTINE freeze(IPAIR)
2 *
3 *
4 * Partial reflection of KS binary.
5 * --------------------------------
6 *
7  include 'common6.h'
8  SAVE nksref
9  DATA nksref /0/
10 *
11 *
12 * Set c.m. & component indices and form semi-major axis.
13  icm = n + ipair
14  i2 = 2*ipair
15  i1 = i2 - 1
16  semi = -0.5d0*body(icm)/h(ipair)
17 *
18 * Define auxiliary quantities and obtain the eccentricity.
19  zeta = 1.0 - r(ipair)/semi
20  psi = tdot2(ipair)/sqrt(body(icm))
21  ecc = sqrt(zeta**2 + psi**2/semi)
22 *
23 * Skip reflection of hyperbolic orbit (just in case).
24  IF (ecc.GE.1.0) go to 100
25 *
26 * Determine the eccentric anomaly with respect to pericentre (0,PI).
27  theta = atan2(abs(psi)/sqrt(semi),zeta)
28 *
29 * Obtain total reflection time from Kepler's equation.
30  dt = 2.0d0*semi*sqrt(semi/body(icm))*(theta - abs(psi)/sqrt(semi))
31 *
32 * Specify regularized time (based on Baumgarte & Stielel, 1974).
33  dtu = -2.0d0*(h(ipair)*dt + tdot2(ipair))/body(icm)
34 *
35 * Skip reflection near pericentre.
36  IF (dtu.LT.4.0*dtau(ipair)) go to 100
37 *
38 * Predict c.m. coordinates and resolve unreflected KS components.
39  CALL xvpred(icm,0)
40 *
41 * Specify transformation coefficients (Mikkola's proceure).
42  xc = zeta/ecc
43  ys = 2.0d0*psi/(ecc*sqrt(body(icm)))
44  zz = body(icm)/(4.0*semi)
45  r(ipair) = 0.0d0
46 *
47 * Generate analytical solutions for U & UDOT using old U & UDOT.
48  DO 10 k = 1,4
49  u(k,ipair) = u0(k,ipair)*xc - udot(k,ipair)*ys
50  udot(k,ipair) = u0(k,ipair)*ys*zz + udot(k,ipair)*xc
51  u0(k,ipair) = u(k,ipair)
52  r(ipair) = r(ipair) + u(k,ipair)**2
53 * Consistent R = U*U for stabilization procedure.
54  10 CONTINUE
55 *
56 * Copy significant perturbers for tidal correction (R-dependent).
57  corr = 2.0 + 2.0*(semi - r(ipair))/(semi*(1.0 + ecc))
58  rcrit2 = cmsep2*(corr*r(ipair))**2
59  nnb1 = list(1,i1) + 1
60  np = 0
61  DO 20 l = 2,nnb1
62  j = list(l,i1)
63  rij2 = (x(1,icm) - x(1,j))**2 + (x(2,icm) - x(2,j))**2 +
64  & (x(3,icm) - x(3,j))**2
65  IF (rij2.LT.rcrit2) THEN
66  np = np + 1
67  jpert(np) = j
68  END IF
69  20 CONTINUE
70 *
71 * Ensure that at least one perturber is retained.
72  IF (np.EQ.0) THEN
73  np = 1
74  jpert(np) = list(2,i1)
75  END IF
76 *
77 * Reverse second time derivative and specify unperturbed motion.
78  tdot2(ipair) = -tdot2(ipair)
79  list(1,i1) = 0
80 *
81 * Set new step in physical units and increase counter.
82  step(i1) = dt
83  nksref = nksref + 1
84 *
85 * Check minimum two-body distance.
86  dmin2 = min(dmin2,semi*(1.0d0 - ecc))
87 *
88 * Obtain potential energy with respect to the components.
89  jlist(1) = i1
90  jlist(2) = i2
91  CALL nbpot(2,np,pot1)
92 *
93 * Transform to reflected coordinates and repeat the summation.
94  CALL ksres(ipair,j1,j2,rij2)
95  CALL nbpot(2,np,pot2)
96 *
97 * Form correction factor due to the tidal potential.
98  vi2 = x0dot(1,icm)**2 + x0dot(2,icm)**2 + x0dot(3,icm)**2
99  corr = 1.0 + 2.0*(pot2 - pot1)/(body(icm)*vi2)
100  etcorr = etcorr + (pot2 - pot1)
101  IF (corr.GT.0.0d0) THEN
102  corr = sqrt(corr)
103  ELSE
104  corr = 0.0d0
105  END IF
106 *
107 * Modify the c.m. velocity.
108  DO 30 k = 1,3
109  x0dot(k,icm) = corr*x0dot(k,icm)
110  30 CONTINUE
111 *
112  100 RETURN
113 *
114  END