pfac.f

Go to the documentation of this file.

      REAL*8 FUNCTION pfac(A,Z)

*       Precession factor for hierarchy.
*       --------------------------------
*
*       Maths by Douglas Heggie; function by Rosemary Mardling (11/96).
*       @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
*
      REAL*8  m,m1,a,z,c,c1,c2,c3,d1,d2,d3,z1

*The function is singular for m1=0, but the following holds:
*For 0 < m1 <=1 (m1 never > 1),  \infty > PFAC(m1) > 1.6, with P(0.01)=3.6.
*So taking it as 1 isn't too bad for an order of magnitude estimate. 

 
      c1 = 0.5*(z + a**2)
      z1 = 5.0 - z + 4*a**2
      c2 = (z1 + sqrt(z1))/6.0
      c3 = (z1 - sqrt(z1))/6.0
      IF (c1.GE.c2) THEN
          d1 = c1
          d2 = c2
          d3 = c3
      ELSE IF (c2.GE.c1.AND.c1.GE.c3) THEN
          d1 = c2
          d2 = c1
          d3 = c3
      ELSE
          d1 = c2
          d2 = c3
          d3 = c1
      END IF
*
      c = 1.0/sqrt(d1 - d3)
      m = (d2 - d3)/(d1 - d3)
*
*       Evaluate elliptic integral by approximate expression.
      a0=1.3862944
      a1=0.1119723
      a2=0.0725296
      b0=0.5
      b1=0.1213478
      b2=0.0288729
      m1 = 1.0 - m
      m1 = max(m1,1.0d-10)

      pfac = (a0+a1*m1+a2*m1**2)+(b0+b1*m1+b2*m1**2)*log(1.0/m1)
*
*       Include the square root of (D1 - D3).
      pfac = c*pfac

      END