Nbody6
pfac.f
Go to the documentation of this file.
1  REAL*8 FUNCTION pfac(A,Z)
2
3 * Precession factor for hierarchy.
4 * --------------------------------
5 *
6 * Maths by Douglas Heggie; function by Rosemary Mardling (11/96).
7 * @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
8 *
9  REAL*8 m,m1,a,z,c,c1,c2,c3,d1,d2,d3,z1
10
11 *The function is singular for m1=0, but the following holds:
12 *For 0 < m1 <=1 (m1 never > 1), \infty > PFAC(m1) > 1.6, with P(0.01)=3.6.
13 *So taking it as 1 isn't too bad for an order of magnitude estimate.
14
15
16  c1 = 0.5*(z + a**2)
17  z1 = 5.0 - z + 4*a**2
18  c2 = (z1 + sqrt(z1))/6.0
19  c3 = (z1 - sqrt(z1))/6.0
20  IF (c1.GE.c2) THEN
21  d1 = c1
22  d2 = c2
23  d3 = c3
24  ELSE IF (c2.GE.c1.AND.c1.GE.c3) THEN
25  d1 = c2
26  d2 = c1
27  d3 = c3
28  ELSE
29  d1 = c2
30  d2 = c3
31  d3 = c1
32  END IF
33 *
34  c = 1.0/sqrt(d1 - d3)
35  m = (d2 - d3)/(d1 - d3)
36 *
37 * Evaluate elliptic integral by approximate expression.
38  a0=1.3862944
39  a1=0.1119723
40  a2=0.0725296
41  b0=0.5
42  b1=0.1213478
43  b2=0.0288729
44  m1 = 1.0 - m
45  m1 = max(m1,1.0d-10)
46
47  pfac = (a0+a1*m1+a2*m1**2)+(b0+b1*m1+b2*m1**2)*log(1.0/m1)
48 *
49 * Include the square root of (D1 - D3).
50  pfac = c*pfac
51
52  END