← Back to branch summary

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

« back to all changes in this revision

Viewing changes to scipy/special/cdflib/brcomp.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/brcomp.f
 
1
      DOUBLE PRECISION FUNCTION brcomp(a,b,x,y)
 
2
C-----------------------------------------------------------------------
 
3
C               EVALUATION OF X**A*Y**B/BETA(A,B)
 
4
C-----------------------------------------------------------------------
 
5
C     .. Scalar Arguments ..
 
6
      DOUBLE PRECISION a,b,x,y
 
7
C     ..
 
8
C     .. Local Scalars ..
 
9
      DOUBLE PRECISION a0,apb,b0,c,const,e,h,lambda,lnx,lny,t,u,v,x0,y0,
 
10
     +                 z
 
11
      INTEGER i,n
 
12
C     ..
 
13
C     .. External Functions ..
 
14
      DOUBLE PRECISION algdiv,alnrel,bcorr,betaln,gam1,gamln1,rlog1
 
15
      EXTERNAL algdiv,alnrel,bcorr,betaln,gam1,gamln1,rlog1
 
16
C     ..
 
17
C     .. Intrinsic Functions ..
 
18
      INTRINSIC abs,dble,dlog,dmax1,dmin1,exp,sqrt
 
19
C     ..
 
20
C     .. Data statements ..
 
21
C-----------------
 
22
C     CONST = 1/SQRT(2*PI)
 
23
C-----------------
 
24
      DATA const/.398942280401433D0/
 
25
C     ..
 
26
C     .. Executable Statements ..
 
27
C
 
28
      brcomp = 0.0D0
 
29
      IF (x.EQ.0.0D0 .OR. y.EQ.0.0D0) RETURN
 
30
      a0 = dmin1(a,b)
 
31
      IF (a0.GE.8.0D0) GO TO 130
 
32
C
 
33
      IF (x.GT.0.375D0) GO TO 10
 
34
      lnx = dlog(x)
 
35
      lny = alnrel(-x)
 
36
      GO TO 30
 
37
 
 
38
   10 IF (y.GT.0.375D0) GO TO 20
 
39
      lnx = alnrel(-y)
 
40
      lny = dlog(y)
 
41
      GO TO 30
 
42
 
 
43
   20 lnx = dlog(x)
 
44
      lny = dlog(y)
 
45
C
 
46
   30 z = a*lnx + b*lny
 
47
      IF (a0.LT.1.0D0) GO TO 40
 
48
      z = z - betaln(a,b)
 
49
      brcomp = exp(z)
 
50
      RETURN
 
51
C-----------------------------------------------------------------------
 
52
C              PROCEDURE FOR A .LT. 1 OR B .LT. 1
 
53
C-----------------------------------------------------------------------
 
54
   40 b0 = dmax1(a,b)
 
55
      IF (b0.GE.8.0D0) GO TO 120
 
56
      IF (b0.GT.1.0D0) GO TO 70
 
57
C
 
58
C                   ALGORITHM FOR B0 .LE. 1
 
59
C
 
60
      brcomp = exp(z)
 
61
      IF (brcomp.EQ.0.0D0) RETURN
 
62
C
 
63
      apb = a + b
 
64
      IF (apb.GT.1.0D0) GO TO 50
 
65
      z = 1.0D0 + gam1(apb)
 
66
      GO TO 60
 
67
 
 
68
   50 u = dble(a) + dble(b) - 1.D0
 
69
      z = (1.0D0+gam1(u))/apb
 
70
C
 
71
   60 c = (1.0D0+gam1(a))* (1.0D0+gam1(b))/z
 
72
      brcomp = brcomp* (a0*c)/ (1.0D0+a0/b0)
 
73
      RETURN
 
74
C
 
75
C                ALGORITHM FOR 1 .LT. B0 .LT. 8
 
76
C
 
77
   70 u = gamln1(a0)
 
78
      n = b0 - 1.0D0
 
79
      IF (n.LT.1) GO TO 90
 
80
      c = 1.0D0
 
81
      DO 80 i = 1,n
 
82
          b0 = b0 - 1.0D0
 
83
          c = c* (b0/ (a0+b0))
 
84
   80 CONTINUE
 
85
      u = dlog(c) + u
 
86
C
 
87
   90 z = z - u
 
88
      b0 = b0 - 1.0D0
 
89
      apb = a0 + b0
 
90
      IF (apb.GT.1.0D0) GO TO 100
 
91
      t = 1.0D0 + gam1(apb)
 
92
      GO TO 110
 
93
 
 
94
  100 u = dble(a0) + dble(b0) - 1.D0
 
95
      t = (1.0D0+gam1(u))/apb
 
96
  110 brcomp = a0*exp(z)* (1.0D0+gam1(b0))/t
 
97
      RETURN
 
98
C
 
99
C                   ALGORITHM FOR B0 .GE. 8
 
100
C
 
101
  120 u = gamln1(a0) + algdiv(a0,b0)
 
102
      brcomp = a0*exp(z-u)
 
103
      RETURN
 
104
C-----------------------------------------------------------------------
 
105
C              PROCEDURE FOR A .GE. 8 AND B .GE. 8
 
106
C-----------------------------------------------------------------------
 
107
  130 IF (a.GT.b) GO TO 140
 
108
      h = a/b
 
109
      x0 = h/ (1.0D0+h)
 
110
      y0 = 1.0D0/ (1.0D0+h)
 
111
      lambda = a - (a+b)*x
 
112
      GO TO 150
 
113
 
 
114
  140 h = b/a
 
115
      x0 = 1.0D0/ (1.0D0+h)
 
116
      y0 = h/ (1.0D0+h)
 
117
      lambda = (a+b)*y - b
 
118
C
 
119
  150 e = -lambda/a
 
120
      IF (abs(e).GT.0.6D0) GO TO 160
 
121
      u = rlog1(e)
 
122
      GO TO 170
 
123
 
 
124
  160 u = e - dlog(x/x0)
 
125
C
 
126
  170 e = lambda/b
 
127
      IF (abs(e).GT.0.6D0) GO TO 180
 
128
      v = rlog1(e)
 
129
      GO TO 190
 
130
 
 
131
  180 v = e - dlog(y/y0)
 
132
C
 
133
  190 z = exp(- (a*u+b*v))
 
134
      brcomp = const*sqrt(b*x0)*z*exp(-bcorr(a,b))
 
135
      RETURN
 
136
 
 
137
      END