← Back to branch summary

~ubuntu-branches/ubuntu/saucy/python-scipy/saucy

« back to all changes in this revision

Viewing changes to scipy/special/cdflib/bfrac.f

  • Committer: Bazaar Package Importer
  • Author(s): Ondrej Certik
  • Date: 2008-06-16 22:58:01 UTC
  • mfrom: (2.1.24 intrepid)
  • Revision ID: james.westby@ubuntu.com-20080616225801-irdhrpcwiocfbcmt
Tags: 0.6.0-12
http://bugs.debian.org/489149
* The description updated to match the current SciPy (Closes: #489149).
* Standards-Version bumped to 3.8.0 (no action needed)
* Build-Depends: netcdf-dev changed to libnetcdf-dev

Show diffs side-by-side

added added

removed removed

Lines of Context:
scipy/special/cdflib/bfrac.f
 
1
      DOUBLE PRECISION FUNCTION bfrac(a,b,x,y,lambda,eps)
 
2
C-----------------------------------------------------------------------
 
3
C     CONTINUED FRACTION EXPANSION FOR IX(A,B) WHEN A,B .GT. 1.
 
4
C     IT IS ASSUMED THAT  LAMBDA = (A + B)*Y - B.
 
5
C-----------------------------------------------------------------------
 
6
C     .. Scalar Arguments ..
 
7
      DOUBLE PRECISION a,b,eps,lambda,x,y
 
8
C     ..
 
9
C     .. Local Scalars ..
 
10
      DOUBLE PRECISION alpha,an,anp1,beta,bn,bnp1,c,c0,c1,e,n,p,r,r0,s,
 
11
     +                 t,w,yp1
 
12
C     ..
 
13
C     .. External Functions ..
 
14
      DOUBLE PRECISION brcomp
 
15
      EXTERNAL brcomp
 
16
C     ..
 
17
C     .. Intrinsic Functions ..
 
18
      INTRINSIC abs
 
19
C     ..
 
20
C     .. Executable Statements ..
 
21
C--------------------
 
22
      bfrac = brcomp(a,b,x,y)
 
23
      IF (bfrac.EQ.0.0D0) RETURN
 
24
C
 
25
      c = 1.0D0 + lambda
 
26
      c0 = b/a
 
27
      c1 = 1.0D0 + 1.0D0/a
 
28
      yp1 = y + 1.0D0
 
29
C
 
30
      n = 0.0D0
 
31
      p = 1.0D0
 
32
      s = a + 1.0D0
 
33
      an = 0.0D0
 
34
      bn = 1.0D0
 
35
      anp1 = 1.0D0
 
36
      bnp1 = c/c1
 
37
      r = c1/c
 
38
C
 
39
C        CONTINUED FRACTION CALCULATION
 
40
C
 
41
   10 n = n + 1.0D0
 
42
      t = n/a
 
43
      w = n* (b-n)*x
 
44
      e = a/s
 
45
      alpha = (p* (p+c0)*e*e)* (w*x)
 
46
      e = (1.0D0+t)/ (c1+t+t)
 
47
      beta = n + w/s + e* (c+n*yp1)
 
48
      p = 1.0D0 + t
 
49
      s = s + 2.0D0
 
50
C
 
51
C        UPDATE AN, BN, ANP1, AND BNP1
 
52
C
 
53
      t = alpha*an + beta*anp1
 
54
      an = anp1
 
55
      anp1 = t
 
56
      t = alpha*bn + beta*bnp1
 
57
      bn = bnp1
 
58
      bnp1 = t
 
59
C
 
60
      r0 = r
 
61
      r = anp1/bnp1
 
62
      IF (abs(r-r0).LE.eps*r) GO TO 20
 
63
C
 
64
C        RESCALE AN, BN, ANP1, AND BNP1
 
65
C
 
66
      an = an/bnp1
 
67
      bn = bn/bnp1
 
68
      anp1 = r
 
69
      bnp1 = 1.0D0
 
70
      GO TO 10
 
71
C
 
72
C                 TERMINATION
 
73
C
 
74
   20 bfrac = bfrac*r
 
75
      RETURN
 
76
 
 
77
      END