← Back to branch summary

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

« back to all changes in this revision

Viewing changes to Lib/special/amos/zuni2.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/amos/zuni2.f
 
1
      SUBROUTINE ZUNI2(ZR, ZI, FNU, KODE, N, YR, YI, NZ, NLAST, FNUL,
 
2
     * TOL, ELIM, ALIM)
 
3
C***BEGIN PROLOGUE  ZUNI2
 
4
C***REFER TO  ZBESI,ZBESK
 
5
C
 
6
C     ZUNI2 COMPUTES I(FNU,Z) IN THE RIGHT HALF PLANE BY MEANS OF
 
7
C     UNIFORM ASYMPTOTIC EXPANSION FOR J(FNU,ZN) WHERE ZN IS Z*I
 
8
C     OR -Z*I AND ZN IS IN THE RIGHT HALF PLANE ALSO.
 
9
C
 
10
C     FNUL IS THE SMALLEST ORDER PERMITTED FOR THE ASYMPTOTIC
 
11
C     EXPANSION. NLAST=0 MEANS ALL OF THE Y VALUES WERE SET.
 
12
C     NLAST.NE.0 IS THE NUMBER LEFT TO BE COMPUTED BY ANOTHER
 
13
C     FORMULA FOR ORDERS FNU TO FNU+NLAST-1 BECAUSE FNU+NLAST-1.LT.FNUL.
 
14
C     Y(I)=CZERO FOR I=NLAST+1,N
 
15
C
 
16
C***ROUTINES CALLED  ZAIRY,ZUCHK,ZUNHJ,ZUOIK,D1MACH,AZABS
 
17
C***END PROLOGUE  ZUNI2
 
18
C     COMPLEX AI,ARG,ASUM,BSUM,CFN,CI,CID,CIP,CONE,CRSC,CSCL,CSR,CSS,
 
19
C    *CZERO,C1,C2,DAI,PHI,RZ,S1,S2,Y,Z,ZB,ZETA1,ZETA2,ZN
 
20
      DOUBLE PRECISION AARG, AIC, AII, AIR, ALIM, ANG, APHI, ARGI,
 
21
     * ARGR, ASCLE, ASUMI, ASUMR, BRY, BSUMI, BSUMR, CIDI, CIPI, CIPR,
 
22
     * CONER, CRSC, CSCL, CSRR, CSSR, C1R, C2I, C2M, C2R, DAII,
 
23
     * DAIR, ELIM, FN, FNU, FNUL, HPI, PHII, PHIR, RAST, RAZ, RS1, RZI,
 
24
     * RZR, STI, STR, S1I, S1R, S2I, S2R, TOL, YI, YR, ZBI, ZBR, ZEROI,
 
25
     * ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R, ZI, ZNI, ZNR, ZR, CYR,
 
26
     * CYI, D1MACH, AZABS, CAR, SAR
 
27
      INTEGER I, IFLAG, IN, INU, J, K, KODE, N, NAI, ND, NDAI, NLAST,
 
28
     * NN, NUF, NW, NZ, IDUM
 
29
      DIMENSION BRY(3), YR(N), YI(N), CIPR(4), CIPI(4), CSSR(3),
 
30
     * CSRR(3), CYR(2), CYI(2)
 
31
      DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 /
 
32
      DATA CIPR(1),CIPI(1),CIPR(2),CIPI(2),CIPR(3),CIPI(3),CIPR(4),
 
33
     * CIPI(4)/ 1.0D0,0.0D0, 0.0D0,1.0D0, -1.0D0,0.0D0, 0.0D0,-1.0D0/
 
34
      DATA HPI, AIC  /
 
35
     1      1.57079632679489662D+00,     1.265512123484645396D+00/
 
36
C
 
37
      NZ = 0
 
38
      ND = N
 
39
      NLAST = 0
 
40
C-----------------------------------------------------------------------
 
41
C     COMPUTED VALUES WITH EXPONENTS BETWEEN ALIM AND ELIM IN MAG-
 
42
C     NITUDE ARE SCALED TO KEEP INTERMEDIATE ARITHMETIC ON SCALE,
 
43
C     EXP(ALIM)=EXP(ELIM)*TOL
 
44
C-----------------------------------------------------------------------
 
45
      CSCL = 1.0D0/TOL
 
46
      CRSC = TOL
 
47
      CSSR(1) = CSCL
 
48
      CSSR(2) = CONER
 
49
      CSSR(3) = CRSC
 
50
      CSRR(1) = CRSC
 
51
      CSRR(2) = CONER
 
52
      CSRR(3) = CSCL
 
53
      BRY(1) = 1.0D+3*D1MACH(1)/TOL
 
54
C-----------------------------------------------------------------------
 
55
C     ZN IS IN THE RIGHT HALF PLANE AFTER ROTATION BY CI OR -CI
 
56
C-----------------------------------------------------------------------
 
57
      ZNR = ZI
 
58
      ZNI = -ZR
 
59
      ZBR = ZR
 
60
      ZBI = ZI
 
61
      CIDI = -CONER
 
62
      INU = INT(SNGL(FNU))
 
63
      ANG = HPI*(FNU-DBLE(FLOAT(INU)))
 
64
      C2R = DCOS(ANG)
 
65
      C2I = DSIN(ANG)
 
66
      CAR = C2R
 
67
      SAR = C2I
 
68
      IN = INU + N - 1
 
69
      IN = MOD(IN,4) + 1
 
70
      STR = C2R*CIPR(IN) - C2I*CIPI(IN)
 
71
      C2I = C2R*CIPI(IN) + C2I*CIPR(IN)
 
72
      C2R = STR
 
73
      IF (ZI.GT.0.0D0) GO TO 10
 
74
      ZNR = -ZNR
 
75
      ZBI = -ZBI
 
76
      CIDI = -CIDI
 
77
      C2I = -C2I
 
78
   10 CONTINUE
 
79
C-----------------------------------------------------------------------
 
80
C     CHECK FOR UNDERFLOW AND OVERFLOW ON FIRST MEMBER
 
81
C-----------------------------------------------------------------------
 
82
      FN = DMAX1(FNU,1.0D0)
 
83
      CALL ZUNHJ(ZNR, ZNI, FN, 1, TOL, PHIR, PHII, ARGR, ARGI, ZETA1R,
 
84
     * ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
 
85
      IF (KODE.EQ.1) GO TO 20
 
86
      STR = ZBR + ZETA2R
 
87
      STI = ZBI + ZETA2I
 
88
      RAST = FN/AZABS(STR,STI)
 
89
      STR = STR*RAST*RAST
 
90
      STI = -STI*RAST*RAST
 
91
      S1R = -ZETA1R + STR
 
92
      S1I = -ZETA1I + STI
 
93
      GO TO 30
 
94
   20 CONTINUE
 
95
      S1R = -ZETA1R + ZETA2R
 
96
      S1I = -ZETA1I + ZETA2I
 
97
   30 CONTINUE
 
98
      RS1 = S1R
 
99
      IF (DABS(RS1).GT.ELIM) GO TO 150
 
100
   40 CONTINUE
 
101
      NN = MIN0(2,ND)
 
102
      DO 90 I=1,NN
 
103
        FN = FNU + DBLE(FLOAT(ND-I))
 
104
        CALL ZUNHJ(ZNR, ZNI, FN, 0, TOL, PHIR, PHII, ARGR, ARGI,
 
105
     *   ZETA1R, ZETA1I, ZETA2R, ZETA2I, ASUMR, ASUMI, BSUMR, BSUMI)
 
106
        IF (KODE.EQ.1) GO TO 50
 
107
        STR = ZBR + ZETA2R
 
108
        STI = ZBI + ZETA2I
 
109
        RAST = FN/AZABS(STR,STI)
 
110
        STR = STR*RAST*RAST
 
111
        STI = -STI*RAST*RAST
 
112
        S1R = -ZETA1R + STR
 
113
        S1I = -ZETA1I + STI + DABS(ZI)
 
114
        GO TO 60
 
115
   50   CONTINUE
 
116
        S1R = -ZETA1R + ZETA2R
 
117
        S1I = -ZETA1I + ZETA2I
 
118
   60   CONTINUE
 
119
C-----------------------------------------------------------------------
 
120
C     TEST FOR UNDERFLOW AND OVERFLOW
 
121
C-----------------------------------------------------------------------
 
122
        RS1 = S1R
 
123
        IF (DABS(RS1).GT.ELIM) GO TO 120
 
124
        IF (I.EQ.1) IFLAG = 2
 
125
        IF (DABS(RS1).LT.ALIM) GO TO 70
 
126
C-----------------------------------------------------------------------
 
127
C     REFINE  TEST AND SCALE
 
128
C-----------------------------------------------------------------------
 
129
C-----------------------------------------------------------------------
 
130
        APHI = AZABS(PHIR,PHII)
 
131
        AARG = AZABS(ARGR,ARGI)
 
132
        RS1 = RS1 + DLOG(APHI) - 0.25D0*DLOG(AARG) - AIC
 
133
        IF (DABS(RS1).GT.ELIM) GO TO 120
 
134
        IF (I.EQ.1) IFLAG = 1
 
135
        IF (RS1.LT.0.0D0) GO TO 70
 
136
        IF (I.EQ.1) IFLAG = 3
 
137
   70   CONTINUE
 
138
C-----------------------------------------------------------------------
 
139
C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
 
140
C     EXPONENT EXTREMES
 
141
C-----------------------------------------------------------------------
 
142
        CALL ZAIRY(ARGR, ARGI, 0, 2, AIR, AII, NAI, IDUM)
 
143
        CALL ZAIRY(ARGR, ARGI, 1, 2, DAIR, DAII, NDAI, IDUM)
 
144
        STR = DAIR*BSUMR - DAII*BSUMI
 
145
        STI = DAIR*BSUMI + DAII*BSUMR
 
146
        STR = STR + (AIR*ASUMR-AII*ASUMI)
 
147
        STI = STI + (AIR*ASUMI+AII*ASUMR)
 
148
        S2R = PHIR*STR - PHII*STI
 
149
        S2I = PHIR*STI + PHII*STR
 
150
        STR = DEXP(S1R)*CSSR(IFLAG)
 
151
        S1R = STR*DCOS(S1I)
 
152
        S1I = STR*DSIN(S1I)
 
153
        STR = S2R*S1R - S2I*S1I
 
154
        S2I = S2R*S1I + S2I*S1R
 
155
        S2R = STR
 
156
        IF (IFLAG.NE.1) GO TO 80
 
157
        CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
 
158
        IF (NW.NE.0) GO TO 120
 
159
   80   CONTINUE
 
160
        IF (ZI.LE.0.0D0) S2I = -S2I
 
161
        STR = S2R*C2R - S2I*C2I
 
162
        S2I = S2R*C2I + S2I*C2R
 
163
        S2R = STR
 
164
        CYR(I) = S2R
 
165
        CYI(I) = S2I
 
166
        J = ND - I + 1
 
167
        YR(J) = S2R*CSRR(IFLAG)
 
168
        YI(J) = S2I*CSRR(IFLAG)
 
169
        STR = -C2I*CIDI
 
170
        C2I = C2R*CIDI
 
171
        C2R = STR
 
172
   90 CONTINUE
 
173
      IF (ND.LE.2) GO TO 110
 
174
      RAZ = 1.0D0/AZABS(ZR,ZI)
 
175
      STR = ZR*RAZ
 
176
      STI = -ZI*RAZ
 
177
      RZR = (STR+STR)*RAZ
 
178
      RZI = (STI+STI)*RAZ
 
179
      BRY(2) = 1.0D0/BRY(1)
 
180
      BRY(3) = D1MACH(2)
 
181
      S1R = CYR(1)
 
182
      S1I = CYI(1)
 
183
      S2R = CYR(2)
 
184
      S2I = CYI(2)
 
185
      C1R = CSRR(IFLAG)
 
186
      ASCLE = BRY(IFLAG)
 
187
      K = ND - 2
 
188
      FN = DBLE(FLOAT(K))
 
189
      DO 100 I=3,ND
 
190
        C2R = S2R
 
191
        C2I = S2I
 
192
        S2R = S1R + (FNU+FN)*(RZR*C2R-RZI*C2I)
 
193
        S2I = S1I + (FNU+FN)*(RZR*C2I+RZI*C2R)
 
194
        S1R = C2R
 
195
        S1I = C2I
 
196
        C2R = S2R*C1R
 
197
        C2I = S2I*C1R
 
198
        YR(K) = C2R
 
199
        YI(K) = C2I
 
200
        K = K - 1
 
201
        FN = FN - 1.0D0
 
202
        IF (IFLAG.GE.3) GO TO 100
 
203
        STR = DABS(C2R)
 
204
        STI = DABS(C2I)
 
205
        C2M = DMAX1(STR,STI)
 
206
        IF (C2M.LE.ASCLE) GO TO 100
 
207
        IFLAG = IFLAG + 1
 
208
        ASCLE = BRY(IFLAG)
 
209
        S1R = S1R*C1R
 
210
        S1I = S1I*C1R
 
211
        S2R = C2R
 
212
        S2I = C2I
 
213
        S1R = S1R*CSSR(IFLAG)
 
214
        S1I = S1I*CSSR(IFLAG)
 
215
        S2R = S2R*CSSR(IFLAG)
 
216
        S2I = S2I*CSSR(IFLAG)
 
217
        C1R = CSRR(IFLAG)
 
218
  100 CONTINUE
 
219
  110 CONTINUE
 
220
      RETURN
 
221
  120 CONTINUE
 
222
      IF (RS1.GT.0.0D0) GO TO 140
 
223
C-----------------------------------------------------------------------
 
224
C     SET UNDERFLOW AND UPDATE PARAMETERS
 
225
C-----------------------------------------------------------------------
 
226
      YR(ND) = ZEROR
 
227
      YI(ND) = ZEROI
 
228
      NZ = NZ + 1
 
229
      ND = ND - 1
 
230
      IF (ND.EQ.0) GO TO 110
 
231
      CALL ZUOIK(ZR, ZI, FNU, KODE, 1, ND, YR, YI, NUF, TOL, ELIM, ALIM)
 
232
      IF (NUF.LT.0) GO TO 140
 
233
      ND = ND - NUF
 
234
      NZ = NZ + NUF
 
235
      IF (ND.EQ.0) GO TO 110
 
236
      FN = FNU + DBLE(FLOAT(ND-1))
 
237
      IF (FN.LT.FNUL) GO TO 130
 
238
C      FN = CIDI
 
239
C      J = NUF + 1
 
240
C      K = MOD(J,4) + 1
 
241
C      S1R = CIPR(K)
 
242
C      S1I = CIPI(K)
 
243
C      IF (FN.LT.0.0D0) S1I = -S1I
 
244
C      STR = C2R*S1R - C2I*S1I
 
245
C      C2I = C2R*S1I + C2I*S1R
 
246
C      C2R = STR
 
247
      IN = INU + ND - 1
 
248
      IN = MOD(IN,4) + 1
 
249
      C2R = CAR*CIPR(IN) - SAR*CIPI(IN)
 
250
      C2I = CAR*CIPI(IN) + SAR*CIPR(IN)
 
251
      IF (ZI.LE.0.0D0) C2I = -C2I
 
252
      GO TO 40
 
253
  130 CONTINUE
 
254
      NLAST = ND
 
255
      RETURN
 
256
  140 CONTINUE
 
257
      NZ = -1
 
258
      RETURN
 
259
  150 CONTINUE
 
260
      IF (RS1.GT.0.0D0) GO TO 140
 
261
      NZ = N
 
262
      DO 160 I=1,N
 
263
        YR(I) = ZEROR
 
264
        YI(I) = ZEROI
 
265
  160 CONTINUE
 
266
      RETURN
 
267
      END