← Back to branch summary

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

« back to all changes in this revision

Viewing changes to scipy/special/amos/zuni1.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/amos/zuni1.f
 
1
      SUBROUTINE ZUNI1(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL,
 
2
     * TOL, ELIM, ALIM)
 
3
C***BEGIN PROLOGUE  ZUNI1
 
4
C***REFER TO  ZBESI,ZBESK
 
5
C
 
6
C     ZUNI1 COMPUTES I(FNU,Z)  BY MEANS OF THE UNIFORM ASYMPTOTIC
 
7
C     EXPANSION FOR I(FNU,Z) IN -PI/3.LE.ARG Z.LE.PI/3.
 
8
C
 
9
C     FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
 
10
C     EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
 
11
C     NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
 
12
C     FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
 
13
C     Y(I)=CZERO FOR I=NLAST+1,N
 
14
C
 
15
C***ROUTINES CALLED  ZUCHK,ZUNIK,ZUOIK,D1MACH,AZABS
 
16
C***END PROLOGUE  ZUNI1
 
17
C     COMPLEX CFN,CONE,CRSC,CSCL,CSR,CSS,CWRK,CZERO,C1,C2,PHI,RZ,SUM,S1,
 
18
C    *S2,Y,Z,ZETA1,ZETA2
 
19
      DOUBLE PRECISION ALIM, APHI, ASCLE, BRY, CONER, CRSC,
 
20
     * CSCL, CSRR, CSSR, CWRKI, CWRKR, C1R, C2I, C2M, C2R, ELIM, FN,
 
21
     * FNU, FNUL, PHII, PHIR, RAST, RS1, RZI, RZR, STI, STR, SUMI,
 
22
     * SUMR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZEROI, ZEROR, ZETA1I,
 
23
     * ZETA1R, ZETA2I, ZETA2R, ZI, ZR, CYR, CYI, D1MACH, AZABS
 
24
      INTEGER I, IFLAG, INIT, K, KODE, M, N, ND, NLAST, NN, NUF, NW, NZ
 
25
      DIMENSION BRY(3), YR(N), YI(N), CWRKR(16), CWRKI(16), CSSR(3),
 
26
     * CSRR(3), CYR(2), CYI(2)
 
27
      DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 /
 
28
C
 
29
      NZ = 0
 
30
      ND = N
 
31
      NLAST = 0
 
32
C-----------------------------------------------------------------------
 
33
C     COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
 
34
C     NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
 
35
C     EXP(ALIM)=EXP(ELIM)*TOL
 
36
C-----------------------------------------------------------------------
 
37
      CSCL = 1.0D0/TOL
 
38
      CRSC = TOL
 
39
      CSSR(1) = CSCL
 
40
      CSSR(2) = CONER
 
41
      CSSR(3) = CRSC
 
42
      CSRR(1) = CRSC
 
43
      CSRR(2) = CONER
 
44
      CSRR(3) = CSCL
 
45
      BRY(1) = 1.0D+3*D1MACH(1)/TOL
 
46
C-----------------------------------------------------------------------
 
47
C     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
 
48
C-----------------------------------------------------------------------
 
49
      FN = DMAX1(FNU,1.0D0)
 
50
      INIT = 0
 
51
      CALL ZUNIK(ZR, ZI, FN, 1, 1, TOL, INIT, PHIR, PHII, ZETA1R,
 
52
     * ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
 
53
      IF (KODE.EQ.1) GO TO 10
 
54
      STR = ZR + ZETA2R
 
55
      STI = ZI + ZETA2I
 
56
      RAST = FN/AZABS(STR,STI)
 
57
      STR = STR*RAST*RAST
 
58
      STI = -STI*RAST*RAST
 
59
      S1R = -ZETA1R + STR
 
60
      S1I = -ZETA1I + STI
 
61
      GO TO 20
 
62
   10 CONTINUE
 
63
      S1R = -ZETA1R + ZETA2R
 
64
      S1I = -ZETA1I + ZETA2I
 
65
   20 CONTINUE
 
66
      RS1 = S1R
 
67
      IF (DABS(RS1).GT.ELIM) GO TO 130
 
68
   30 CONTINUE
 
69
      NN = MIN0(2,ND)
 
70
      DO 80 I=1,NN
 
71
        FN = FNU + DBLE(FLOAT(ND-I))
 
72
        INIT = 0
 
73
        CALL ZUNIK(ZR, ZI, FN, 1, 0, TOL, INIT, PHIR, PHII, ZETA1R,
 
74
     *   ZETA1I, ZETA2R, ZETA2I, SUMR, SUMI, CWRKR, CWRKI)
 
75
        IF (KODE.EQ.1) GO TO 40
 
76
        STR = ZR + ZETA2R
 
77
        STI = ZI + ZETA2I
 
78
        RAST = FN/AZABS(STR,STI)
 
79
        STR = STR*RAST*RAST
 
80
        STI = -STI*RAST*RAST
 
81
        S1R = -ZETA1R + STR
 
82
        S1I = -ZETA1I + STI + ZI
 
83
        GO TO 50
 
84
   40   CONTINUE
 
85
        S1R = -ZETA1R + ZETA2R
 
86
        S1I = -ZETA1I + ZETA2I
 
87
   50   CONTINUE
 
88
C-----------------------------------------------------------------------
 
89
C     TEST FOR UNDERFLOW AND OVERFLOW
 
90
C-----------------------------------------------------------------------
 
91
        RS1 = S1R
 
92
        IF (DABS(RS1).GT.ELIM) GO TO 110
 
93
        IF (I.EQ.1) IFLAG = 2
 
94
        IF (DABS(RS1).LT.ALIM) GO TO 60
 
95
C-----------------------------------------------------------------------
 
96
C     REFINE  TEST AND SCALE
 
97
C-----------------------------------------------------------------------
 
98
        APHI = AZABS(PHIR,PHII)
 
99
        RS1 = RS1 + DLOG(APHI)
 
100
        IF (DABS(RS1).GT.ELIM) GO TO 110
 
101
        IF (I.EQ.1) IFLAG = 1
 
102
        IF (RS1.LT.0.0D0) GO TO 60
 
103
        IF (I.EQ.1) IFLAG = 3
 
104
   60   CONTINUE
 
105
C-----------------------------------------------------------------------
 
106
C     SCALE S1 IF CABS(S1).LT.ASCLE
 
107
C-----------------------------------------------------------------------
 
108
        S2R = PHIR*SUMR - PHII*SUMI
 
109
        S2I = PHIR*SUMI + PHII*SUMR
 
110
        STR = DEXP(S1R)*CSSR(IFLAG)
 
111
        S1R = STR*DCOS(S1I)
 
112
        S1I = STR*DSIN(S1I)
 
113
        STR = S2R*S1R - S2I*S1I
 
114
        S2I = S2R*S1I + S2I*S1R
 
115
        S2R = STR
 
116
        IF (IFLAG.NE.1) GO TO 70
 
117
        CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
 
118
        IF (NW.NE.0) GO TO 110
 
119
   70   CONTINUE
 
120
        CYR(I) = S2R
 
121
        CYI(I) = S2I
 
122
        M = ND - I + 1
 
123
        YR(M) = S2R*CSRR(IFLAG)
 
124
        YI(M) = S2I*CSRR(IFLAG)
 
125
   80 CONTINUE
 
126
      IF (ND.LE.2) GO TO 100
 
127
      RAST = 1.0D0/AZABS(ZR,ZI)
 
128
      STR = ZR*RAST
 
129
      STI = -ZI*RAST
 
130
      RZR = (STR+STR)*RAST
 
131
      RZI = (STI+STI)*RAST
 
132
      BRY(2) = 1.0D0/BRY(1)
 
133
      BRY(3) = D1MACH(2)
 
134
      S1R = CYR(1)
 
135
      S1I = CYI(1)
 
136
      S2R = CYR(2)
 
137
      S2I = CYI(2)
 
138
      C1R = CSRR(IFLAG)
 
139
      ASCLE = BRY(IFLAG)
 
140
      K = ND - 2
 
141
      FN = DBLE(FLOAT(K))
 
142
      DO 90 I=3,ND
 
143
        C2R = S2R
 
144
        C2I = S2I
 
145
        S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I)
 
146
        S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R)
 
147
        S1R = C2R
 
148
        S1I = C2I
 
149
        C2R = S2R*C1R
 
150
        C2I = S2I*C1R
 
151
        YR(K) = C2R
 
152
        YI(K) = C2I
 
153
        K = K - 1
 
154
        FN = FN - 1.0D0
 
155
        IF (IFLAG.GE.3) GO TO 90
 
156
        STR = DABS(C2R)
 
157
        STI = DABS(C2I)
 
158
        C2M = DMAX1(STR,STI)
 
159
        IF (C2M.LE.ASCLE) GO TO 90
 
160
        IFLAG = IFLAG + 1
 
161
        ASCLE = BRY(IFLAG)
 
162
        S1R = S1R*C1R
 
163
        S1I = S1I*C1R
 
164
        S2R = C2R
 
165
        S2I = C2I
 
166
        S1R = S1R*CSSR(IFLAG)
 
167
        S1I = S1I*CSSR(IFLAG)
 
168
        S2R = S2R*CSSR(IFLAG)
 
169
        S2I = S2I*CSSR(IFLAG)
 
170
        C1R = CSRR(IFLAG)
 
171
   90 CONTINUE
 
172
  100 CONTINUE
 
173
      RETURN
 
174
C-----------------------------------------------------------------------
 
175
C     SET UNDERFLOW AND UPDATE PARAMETERS
 
176
C-----------------------------------------------------------------------
 
177
  110 CONTINUE
 
178
      IF (RS1.GT.0.0D0) GO TO 120
 
179
      YR(ND) = ZEROR
 
180
      YI(ND) = ZEROI
 
181
      NZ = NZ + 1
 
182
      ND = ND - 1
 
183
      IF (ND.EQ.0) GO TO 100
 
184
      CALL ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, NUF, TOL, ELIM, ALIM)
 
185
      IF (NUF.LT.0) GO TO 120
 
186
      ND = ND - NUF
 
187
      NZ = NZ + NUF
 
188
      IF (ND.EQ.0) GO TO 100
 
189
      FN = FNU + DBLE(FLOAT(ND-1))
 
190
      IF (FN.GE.FNUL) GO TO 30
 
191
      NLAST = ND
 
192
      RETURN
 
193
  120 CONTINUE
 
194
      NZ = -1
 
195
      RETURN
 
196
  130 CONTINUE
 
197
      IF (RS1.GT.0.0D0) GO TO 120
 
198
      NZ = N
 
199
      DO 140 I=1,N
 
200
        YR(I) = ZEROR
 
201
        YI(I) = ZEROI
 
202
  140 CONTINUE
 
203
      RETURN
 
204
      END