← Back to branch summary

~ubuntu-branches/ubuntu/karmic/python-scipy/karmic

« back to all changes in this revision

Viewing changes to Lib/special/cdflib/basym.f

  • Committer: Bazaar Package Importer
  • Author(s): Daniel T. Chen (new)
  • Date: 2005-03-16 02:15:29 UTC
  • Revision ID: james.westby@ubuntu.com-20050316021529-xrjlowsejs0cijig
Tags: upstream-0.3.2
Import upstream version 0.3.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
Lib/special/cdflib/basym.f
 
1
      DOUBLE PRECISION FUNCTION basym(a,b,lambda,eps)
 
2
C-----------------------------------------------------------------------
 
3
C     ASYMPTOTIC EXPANSION FOR IX(A,B) FOR LARGE A AND B.
 
4
C     LAMBDA = (A + B)*Y - B  AND EPS IS THE TOLERANCE USED.
 
5
C     IT IS ASSUMED THAT LAMBDA IS NONNEGATIVE AND THAT
 
6
C     A AND B ARE GREATER THAN OR EQUAL TO 15.
 
7
C-----------------------------------------------------------------------
 
8
C     .. Scalar Arguments ..
 
9
      DOUBLE PRECISION a,b,eps,lambda
 
10
C     ..
 
11
C     .. Local Scalars ..
 
12
      DOUBLE PRECISION bsum,dsum,e0,e1,f,h,h2,hn,j0,j1,r,r0,r1,s,sum,t,
 
13
     +                 t0,t1,u,w,w0,z,z0,z2,zn,znm1
 
14
      INTEGER i,im1,imj,j,m,mm1,mmj,n,np1,num
 
15
C     ..
 
16
C     .. Local Arrays ..
 
17
      DOUBLE PRECISION a0(21),b0(21),c(21),d(21)
 
18
C     ..
 
19
C     .. External Functions ..
 
20
      DOUBLE PRECISION bcorr,erfc1,rlog1
 
21
      EXTERNAL bcorr,erfc1,rlog1
 
22
C     ..
 
23
C     .. Intrinsic Functions ..
 
24
      INTRINSIC abs,exp,sqrt
 
25
C     ..
 
26
C     .. Data statements ..
 
27
C------------------------
 
28
C     ****** NUM IS THE MAXIMUM VALUE THAT N CAN TAKE IN THE DO LOOP
 
29
C            ENDING AT STATEMENT 50. IT IS REQUIRED THAT NUM BE EVEN.
 
30
C            THE ARRAYS A0, B0, C, D HAVE DIMENSION NUM + 1.
 
31
C
 
32
C------------------------
 
33
C     E0 = 2/SQRT(PI)
 
34
C     E1 = 2**(-3/2)
 
35
C------------------------
 
36
      DATA num/20/
 
37
      DATA e0/1.12837916709551D0/,e1/.353553390593274D0/
 
38
C     ..
 
39
C     .. Executable Statements ..
 
40
C------------------------
 
41
      basym = 0.0D0
 
42
      IF (a.GE.b) GO TO 10
 
43
      h = a/b
 
44
      r0 = 1.0D0/ (1.0D0+h)
 
45
      r1 = (b-a)/b
 
46
      w0 = 1.0D0/sqrt(a* (1.0D0+h))
 
47
      GO TO 20
 
48
 
 
49
   10 h = b/a
 
50
      r0 = 1.0D0/ (1.0D0+h)
 
51
      r1 = (b-a)/a
 
52
      w0 = 1.0D0/sqrt(b* (1.0D0+h))
 
53
C
 
54
   20 f = a*rlog1(-lambda/a) + b*rlog1(lambda/b)
 
55
      t = exp(-f)
 
56
      IF (t.EQ.0.0D0) RETURN
 
57
      z0 = sqrt(f)
 
58
      z = 0.5D0* (z0/e1)
 
59
      z2 = f + f
 
60
C
 
61
      a0(1) = (2.0D0/3.0D0)*r1
 
62
      c(1) = -0.5D0*a0(1)
 
63
      d(1) = -c(1)
 
64
      j0 = (0.5D0/e0)*erfc1(1,z0)
 
65
      j1 = e1
 
66
      sum = j0 + d(1)*w0*j1
 
67
C
 
68
      s = 1.0D0
 
69
      h2 = h*h
 
70
      hn = 1.0D0
 
71
      w = w0
 
72
      znm1 = z
 
73
      zn = z2
 
74
      DO 70 n = 2,num,2
 
75
          hn = h2*hn
 
76
          a0(n) = 2.0D0*r0* (1.0D0+h*hn)/ (n+2.0D0)
 
77
          np1 = n + 1
 
78
          s = s + hn
 
79
          a0(np1) = 2.0D0*r1*s/ (n+3.0D0)
 
80
C
 
81
          DO 60 i = n,np1
 
82
              r = -0.5D0* (i+1.0D0)
 
83
              b0(1) = r*a0(1)
 
84
              DO 40 m = 2,i
 
85
                  bsum = 0.0D0
 
86
                  mm1 = m - 1
 
87
                  DO 30 j = 1,mm1
 
88
                      mmj = m - j
 
89
                      bsum = bsum + (j*r-mmj)*a0(j)*b0(mmj)
 
90
   30             CONTINUE
 
91
                  b0(m) = r*a0(m) + bsum/m
 
92
   40         CONTINUE
 
93
              c(i) = b0(i)/ (i+1.0D0)
 
94
C
 
95
              dsum = 0.0D0
 
96
              im1 = i - 1
 
97
              DO 50 j = 1,im1
 
98
                  imj = i - j
 
99
                  dsum = dsum + d(imj)*c(j)
 
100
   50         CONTINUE
 
101
              d(i) = - (dsum+c(i))
 
102
   60     CONTINUE
 
103
C
 
104
          j0 = e1*znm1 + (n-1.0D0)*j0
 
105
          j1 = e1*zn + n*j1
 
106
          znm1 = z2*znm1
 
107
          zn = z2*zn
 
108
          w = w0*w
 
109
          t0 = d(n)*w*j0
 
110
          w = w0*w
 
111
          t1 = d(np1)*w*j1
 
112
          sum = sum + (t0+t1)
 
113
          IF ((abs(t0)+abs(t1)).LE.eps*sum) GO TO 80
 
114
   70 CONTINUE
 
115
C
 
116
   80 u = exp(-bcorr(a,b))
 
117
      basym = e0*t*u*sum
 
118
      RETURN
 
119
 
 
120
      END