← Back to branch summary

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

« back to all changes in this revision

Viewing changes to scipy/special/amos/zunk1.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/zunk1.f
 
1
      SUBROUTINE ZUNK1(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, TOL, ELIM,
 
2
     * ALIM)
 
3
C***BEGIN PROLOGUE  ZUNK1
 
4
C***REFER TO  ZBESK
 
5
C
 
6
C     ZUNK1 COMPUTES K(FNU,Z) AND ITS ANALYTIC CONTINUATION FROM THE
 
7
C     RIGHT HALF PLANE TO THE LEFT HALF PLANE BY MEANS OF THE
 
8
C     UNIFORM ASYMPTOTIC EXPANSION.
 
9
C     MR INDICATES THE DIRECTION OF ROTATION FOR ANALYTIC CONTINUATION.
 
10
C     NZ=-1 MEANS AN OVERFLOW WILL OCCUR
 
11
C
 
12
C***ROUTINES CALLED  ZKSCL,ZS1S2,ZUCHK,ZUNIK,D1MACH,AZABS
 
13
C***END PROLOGUE  ZUNK1
 
14
C     COMPLEX CFN,CK,CONE,CRSC,CS,CSCL,CSGN,CSPN,CSR,CSS,CWRK,CY,CZERO,
 
15
C    *C1,C2,PHI,PHID,RZ,SUM,SUMD,S1,S2,Y,Z,ZETA1,ZETA1D,ZETA2,ZETA2D,ZR
 
16
      DOUBLE PRECISION ALIM, ANG, APHI, ASC, ASCLE, BRY, CKI, CKR,
 
17
     * CONER, CRSC, CSCL, CSGNI, CSPNI, CSPNR, CSR, CSRR, CSSR,
 
18
     * CWRKI, CWRKR, CYI, CYR, C1I, C1R, C2I, C2M, C2R, ELIM, FMR, FN,
 
19
     * FNF, FNU, PHIDI, PHIDR, PHII, PHIR, PI, RAST, RAZR, RS1, RZI,
 
20
     * RZR, SGN, STI, STR, SUMDI, SUMDR, SUMI, SUMR, S1I, S1R, S2I,
 
21
     * S2R, TOL, YI, YR, ZEROI, ZEROR, ZETA1I, ZETA1R, ZETA2I, ZETA2R,
 
22
     * ZET1DI, ZET1DR, ZET2DI, ZET2DR, ZI, ZR, ZRI, ZRR, D1MACH, AZABS
 
23
      INTEGER I, IB, IFLAG, IFN, IL, INIT, INU, IUF, K, KDFLG, KFLAG,
 
24
     * KK, KODE, MR, N, NW, NZ, INITD, IC, IPARD, J
 
25
      DIMENSION BRY(3), INIT(2), YR(N), YI(N), SUMR(2), SUMI(2),
 
26
     * ZETA1R(2), ZETA1I(2), ZETA2R(2), ZETA2I(2), CYR(2), CYI(2),
 
27
     * CWRKR(16,3), CWRKI(16,3), CSSR(3), CSRR(3), PHIR(2), PHII(2)
 
28
      DATA ZEROR,ZEROI,CONER / 0.0D0, 0.0D0, 1.0D0 /
 
29
      DATA PI / 3.14159265358979324D0 /
 
30
C
 
31
      KDFLG = 1
 
32
      NZ = 0
 
33
C-----------------------------------------------------------------------
 
34
C     EXP(-ALIM)=EXP(-ELIM)/TOL=APPROX. ONE PRECISION GREATER THAN
 
35
C     THE UNDERFLOW LIMIT
 
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
      BRY(2) = 1.0D0/BRY(1)
 
47
      BRY(3) = D1MACH(2)
 
48
      ZRR = ZR
 
49
      ZRI = ZI
 
50
      IF (ZR.GE.0.0D0) GO TO 10
 
51
      ZRR = -ZR
 
52
      ZRI = -ZI
 
53
   10 CONTINUE
 
54
      J = 2
 
55
      DO 70 I=1,N
 
56
C-----------------------------------------------------------------------
 
57
C     J FLIP FLOPS BETWEEN 1 AND 2 IN J = 3 - J
 
58
C-----------------------------------------------------------------------
 
59
        J = 3 - J
 
60
        FN = FNU + DBLE(FLOAT(I-1))
 
61
        INIT(J) = 0
 
62
        CALL ZUNIK(ZRR, ZRI, FN, 2, 0, TOL, INIT(J), PHIR(J), PHII(J),
 
63
     *   ZETA1R(J), ZETA1I(J), ZETA2R(J), ZETA2I(J), SUMR(J), SUMI(J),
 
64
     *   CWRKR(1,J), CWRKI(1,J))
 
65
        IF (KODE.EQ.1) GO TO 20
 
66
        STR = ZRR + ZETA2R(J)
 
67
        STI = ZRI + ZETA2I(J)
 
68
        RAST = FN/AZABS(STR,STI)
 
69
        STR = STR*RAST*RAST
 
70
        STI = -STI*RAST*RAST
 
71
        S1R = ZETA1R(J) - STR
 
72
        S1I = ZETA1I(J) - STI
 
73
        GO TO 30
 
74
   20   CONTINUE
 
75
        S1R = ZETA1R(J) - ZETA2R(J)
 
76
        S1I = ZETA1I(J) - ZETA2I(J)
 
77
   30   CONTINUE
 
78
        RS1 = S1R
 
79
C-----------------------------------------------------------------------
 
80
C     TEST FOR UNDERFLOW AND OVERFLOW
 
81
C-----------------------------------------------------------------------
 
82
        IF (DABS(RS1).GT.ELIM) GO TO 60
 
83
        IF (KDFLG.EQ.1) KFLAG = 2
 
84
        IF (DABS(RS1).LT.ALIM) GO TO 40
 
85
C-----------------------------------------------------------------------
 
86
C     REFINE  TEST AND SCALE
 
87
C-----------------------------------------------------------------------
 
88
        APHI = AZABS(PHIR(J),PHII(J))
 
89
        RS1 = RS1 + DLOG(APHI)
 
90
        IF (DABS(RS1).GT.ELIM) GO TO 60
 
91
        IF (KDFLG.EQ.1) KFLAG = 1
 
92
        IF (RS1.LT.0.0D0) GO TO 40
 
93
        IF (KDFLG.EQ.1) KFLAG = 3
 
94
   40   CONTINUE
 
95
C-----------------------------------------------------------------------
 
96
C     SCALE S1 TO KEEP INTERMEDIATE ARITHMETIC ON SCALE NEAR
 
97
C     EXPONENT EXTREMES
 
98
C-----------------------------------------------------------------------
 
99
        S2R = PHIR(J)*SUMR(J) - PHII(J)*SUMI(J)
 
100
        S2I = PHIR(J)*SUMI(J) + PHII(J)*SUMR(J)
 
101
        STR = DEXP(S1R)*CSSR(KFLAG)
 
102
        S1R = STR*DCOS(S1I)
 
103
        S1I = STR*DSIN(S1I)
 
104
        STR = S2R*S1R - S2I*S1I
 
105
        S2I = S1R*S2I + S2R*S1I
 
106
        S2R = STR
 
107
        IF (KFLAG.NE.1) GO TO 50
 
108
        CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
 
109
        IF (NW.NE.0) GO TO 60
 
110
   50   CONTINUE
 
111
        CYR(KDFLG) = S2R
 
112
        CYI(KDFLG) = S2I
 
113
        YR(I) = S2R*CSRR(KFLAG)
 
114
        YI(I) = S2I*CSRR(KFLAG)
 
115
        IF (KDFLG.EQ.2) GO TO 75
 
116
        KDFLG = 2
 
117
        GO TO 70
 
118
   60   CONTINUE
 
119
        IF (RS1.GT.0.0D0) GO TO 300
 
120
C-----------------------------------------------------------------------
 
121
C     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
 
122
C-----------------------------------------------------------------------
 
123
        IF (ZR.LT.0.0D0) GO TO 300
 
124
        KDFLG = 1
 
125
        YR(I)=ZEROR
 
126
        YI(I)=ZEROI
 
127
        NZ=NZ+1
 
128
        IF (I.EQ.1) GO TO 70
 
129
        IF ((YR(I-1).EQ.ZEROR).AND.(YI(I-1).EQ.ZEROI)) GO TO 70
 
130
        YR(I-1)=ZEROR
 
131
        YI(I-1)=ZEROI
 
132
        NZ=NZ+1
 
133
   70 CONTINUE
 
134
      I = N
 
135
   75 CONTINUE
 
136
      RAZR = 1.0D0/AZABS(ZRR,ZRI)
 
137
      STR = ZRR*RAZR
 
138
      STI = -ZRI*RAZR
 
139
      RZR = (STR+STR)*RAZR
 
140
      RZI = (STI+STI)*RAZR
 
141
      CKR = FN*RZR
 
142
      CKI = FN*RZI
 
143
      IB = I + 1
 
144
      IF (N.LT.IB) GO TO 160
 
145
C-----------------------------------------------------------------------
 
146
C     TEST LAST MEMBER FOR UNDERFLOW AND OVERFLOW. SET SEQUENCE TO ZERO
 
147
C     ON UNDERFLOW.
 
148
C-----------------------------------------------------------------------
 
149
      FN = FNU + DBLE(FLOAT(N-1))
 
150
      IPARD = 1
 
151
      IF (MR.NE.0) IPARD = 0
 
152
      INITD = 0
 
153
      CALL ZUNIK(ZRR, ZRI, FN, 2, IPARD, TOL, INITD, PHIDR, PHIDI,
 
154
     * ZET1DR, ZET1DI, ZET2DR, ZET2DI, SUMDR, SUMDI, CWRKR(1,3),
 
155
     * CWRKI(1,3))
 
156
      IF (KODE.EQ.1) GO TO 80
 
157
      STR = ZRR + ZET2DR
 
158
      STI = ZRI + ZET2DI
 
159
      RAST = FN/AZABS(STR,STI)
 
160
      STR = STR*RAST*RAST
 
161
      STI = -STI*RAST*RAST
 
162
      S1R = ZET1DR - STR
 
163
      S1I = ZET1DI - STI
 
164
      GO TO 90
 
165
   80 CONTINUE
 
166
      S1R = ZET1DR - ZET2DR
 
167
      S1I = ZET1DI - ZET2DI
 
168
   90 CONTINUE
 
169
      RS1 = S1R
 
170
      IF (DABS(RS1).GT.ELIM) GO TO 95
 
171
      IF (DABS(RS1).LT.ALIM) GO TO 100
 
172
C----------------------------------------------------------------------------
 
173
C     REFINE ESTIMATE AND TEST
 
174
C-------------------------------------------------------------------------
 
175
      APHI = AZABS(PHIDR,PHIDI)
 
176
      RS1 = RS1+DLOG(APHI)
 
177
      IF (DABS(RS1).LT.ELIM) GO TO 100
 
178
   95 CONTINUE
 
179
      IF (DABS(RS1).GT.0.0D0) GO TO 300
 
180
C-----------------------------------------------------------------------
 
181
C     FOR ZR.LT.0.0, THE I FUNCTION TO BE ADDED WILL OVERFLOW
 
182
C-----------------------------------------------------------------------
 
183
      IF (ZR.LT.0.0D0) GO TO 300
 
184
      NZ = N
 
185
      DO 96 I=1,N
 
186
        YR(I) = ZEROR
 
187
        YI(I) = ZEROI
 
188
   96 CONTINUE
 
189
      RETURN
 
190
C---------------------------------------------------------------------------
 
191
C     FORWARD RECUR FOR REMAINDER OF THE SEQUENCE
 
192
C----------------------------------------------------------------------------
 
193
  100 CONTINUE
 
194
      S1R = CYR(1)
 
195
      S1I = CYI(1)
 
196
      S2R = CYR(2)
 
197
      S2I = CYI(2)
 
198
      C1R = CSRR(KFLAG)
 
199
      ASCLE = BRY(KFLAG)
 
200
      DO 120 I=IB,N
 
201
        C2R = S2R
 
202
        C2I = S2I
 
203
        S2R = CKR*C2R - CKI*C2I + S1R
 
204
        S2I = CKR*C2I + CKI*C2R + S1I
 
205
        S1R = C2R
 
206
        S1I = C2I
 
207
        CKR = CKR + RZR
 
208
        CKI = CKI + RZI
 
209
        C2R = S2R*C1R
 
210
        C2I = S2I*C1R
 
211
        YR(I) = C2R
 
212
        YI(I) = C2I
 
213
        IF (KFLAG.GE.3) GO TO 120
 
214
        STR = DABS(C2R)
 
215
        STI = DABS(C2I)
 
216
        C2M = DMAX1(STR,STI)
 
217
        IF (C2M.LE.ASCLE) GO TO 120
 
218
        KFLAG = KFLAG + 1
 
219
        ASCLE = BRY(KFLAG)
 
220
        S1R = S1R*C1R
 
221
        S1I = S1I*C1R
 
222
        S2R = C2R
 
223
        S2I = C2I
 
224
        S1R = S1R*CSSR(KFLAG)
 
225
        S1I = S1I*CSSR(KFLAG)
 
226
        S2R = S2R*CSSR(KFLAG)
 
227
        S2I = S2I*CSSR(KFLAG)
 
228
        C1R = CSRR(KFLAG)
 
229
  120 CONTINUE
 
230
  160 CONTINUE
 
231
      IF (MR.EQ.0) RETURN
 
232
C-----------------------------------------------------------------------
 
233
C     ANALYTIC CONTINUATION FOR RE(Z).LT.0.0D0
 
234
C-----------------------------------------------------------------------
 
235
      NZ = 0
 
236
      FMR = DBLE(FLOAT(MR))
 
237
      SGN = -DSIGN(PI,FMR)
 
238
C-----------------------------------------------------------------------
 
239
C     CSPN AND CSGN ARE COEFF OF K AND I FUNCTIONS RESP.
 
240
C-----------------------------------------------------------------------
 
241
      CSGNI = SGN
 
242
      INU = INT(SNGL(FNU))
 
243
      FNF = FNU - DBLE(FLOAT(INU))
 
244
      IFN = INU + N - 1
 
245
      ANG = FNF*SGN
 
246
      CSPNR = DCOS(ANG)
 
247
      CSPNI = DSIN(ANG)
 
248
      IF (MOD(IFN,2).EQ.0) GO TO 170
 
249
      CSPNR = -CSPNR
 
250
      CSPNI = -CSPNI
 
251
  170 CONTINUE
 
252
      ASC = BRY(1)
 
253
      IUF = 0
 
254
      KK = N
 
255
      KDFLG = 1
 
256
      IB = IB - 1
 
257
      IC = IB - 1
 
258
      DO 270 K=1,N
 
259
        FN = FNU + DBLE(FLOAT(KK-1))
 
260
C-----------------------------------------------------------------------
 
261
C     LOGIC TO SORT OUT CASES WHOSE PARAMETERS WERE SET FOR THE K
 
262
C     FUNCTION ABOVE
 
263
C-----------------------------------------------------------------------
 
264
        M=3
 
265
        IF (N.GT.2) GO TO 175
 
266
  172   CONTINUE
 
267
        INITD = INIT(J)
 
268
        PHIDR = PHIR(J)
 
269
        PHIDI = PHII(J)
 
270
        ZET1DR = ZETA1R(J)
 
271
        ZET1DI = ZETA1I(J)
 
272
        ZET2DR = ZETA2R(J)
 
273
        ZET2DI = ZETA2I(J)
 
274
        SUMDR = SUMR(J)
 
275
        SUMDI = SUMI(J)
 
276
        M = J
 
277
        J = 3 - J
 
278
        GO TO 180
 
279
  175   CONTINUE
 
280
        IF ((KK.EQ.N).AND.(IB.LT.N)) GO TO 180
 
281
        IF ((KK.EQ.IB).OR.(KK.EQ.IC)) GO TO 172
 
282
        INITD = 0
 
283
  180   CONTINUE
 
284
        CALL ZUNIK(ZRR, ZRI, FN, 1, 0, TOL, INITD, PHIDR, PHIDI,
 
285
     *   ZET1DR, ZET1DI, ZET2DR, ZET2DI, SUMDR, SUMDI,
 
286
     *   CWRKR(1,M), CWRKI(1,M))
 
287
        IF (KODE.EQ.1) GO TO 200
 
288
        STR = ZRR + ZET2DR
 
289
        STI = ZRI + ZET2DI
 
290
        RAST = FN/AZABS(STR,STI)
 
291
        STR = STR*RAST*RAST
 
292
        STI = -STI*RAST*RAST
 
293
        S1R = -ZET1DR + STR
 
294
        S1I = -ZET1DI + STI
 
295
        GO TO 210
 
296
  200   CONTINUE
 
297
        S1R = -ZET1DR + ZET2DR
 
298
        S1I = -ZET1DI + ZET2DI
 
299
  210   CONTINUE
 
300
C-----------------------------------------------------------------------
 
301
C     TEST FOR UNDERFLOW AND OVERFLOW
 
302
C-----------------------------------------------------------------------
 
303
        RS1 = S1R
 
304
        IF (DABS(RS1).GT.ELIM) GO TO 260
 
305
        IF (KDFLG.EQ.1) IFLAG = 2
 
306
        IF (DABS(RS1).LT.ALIM) GO TO 220
 
307
C-----------------------------------------------------------------------
 
308
C     REFINE  TEST AND SCALE
 
309
C-----------------------------------------------------------------------
 
310
        APHI = AZABS(PHIDR,PHIDI)
 
311
        RS1 = RS1 + DLOG(APHI)
 
312
        IF (DABS(RS1).GT.ELIM) GO TO 260
 
313
        IF (KDFLG.EQ.1) IFLAG = 1
 
314
        IF (RS1.LT.0.0D0) GO TO 220
 
315
        IF (KDFLG.EQ.1) IFLAG = 3
 
316
  220   CONTINUE
 
317
        STR = PHIDR*SUMDR - PHIDI*SUMDI
 
318
        STI = PHIDR*SUMDI + PHIDI*SUMDR
 
319
        S2R = -CSGNI*STI
 
320
        S2I = CSGNI*STR
 
321
        STR = DEXP(S1R)*CSSR(IFLAG)
 
322
        S1R = STR*DCOS(S1I)
 
323
        S1I = STR*DSIN(S1I)
 
324
        STR = S2R*S1R - S2I*S1I
 
325
        S2I = S2R*S1I + S2I*S1R
 
326
        S2R = STR
 
327
        IF (IFLAG.NE.1) GO TO 230
 
328
        CALL ZUCHK(S2R, S2I, NW, BRY(1), TOL)
 
329
        IF (NW.EQ.0) GO TO 230
 
330
        S2R = ZEROR
 
331
        S2I = ZEROI
 
332
  230   CONTINUE
 
333
        CYR(KDFLG) = S2R
 
334
        CYI(KDFLG) = S2I
 
335
        C2R = S2R
 
336
        C2I = S2I
 
337
        S2R = S2R*CSRR(IFLAG)
 
338
        S2I = S2I*CSRR(IFLAG)
 
339
C-----------------------------------------------------------------------
 
340
C     ADD I AND K FUNCTIONS, K SEQUENCE IN Y(I), I=1,N
 
341
C-----------------------------------------------------------------------
 
342
        S1R = YR(KK)
 
343
        S1I = YI(KK)
 
344
        IF (KODE.EQ.1) GO TO 250
 
345
        CALL ZS1S2(ZRR, ZRI, S1R, S1I, S2R, S2I, NW, ASC, ALIM, IUF)
 
346
        NZ = NZ + NW
 
347
  250   CONTINUE
 
348
        YR(KK) = S1R*CSPNR - S1I*CSPNI + S2R
 
349
        YI(KK) = CSPNR*S1I + CSPNI*S1R + S2I
 
350
        KK = KK - 1
 
351
        CSPNR = -CSPNR
 
352
        CSPNI = -CSPNI
 
353
        IF (C2R.NE.0.0D0 .OR. C2I.NE.0.0D0) GO TO 255
 
354
        KDFLG = 1
 
355
        GO TO 270
 
356
  255   CONTINUE
 
357
        IF (KDFLG.EQ.2) GO TO 275
 
358
        KDFLG = 2
 
359
        GO TO 270
 
360
  260   CONTINUE
 
361
        IF (RS1.GT.0.0D0) GO TO 300
 
362
        S2R = ZEROR
 
363
        S2I = ZEROI
 
364
        GO TO 230
 
365
  270 CONTINUE
 
366
      K = N
 
367
  275 CONTINUE
 
368
      IL = N - K
 
369
      IF (IL.EQ.0) RETURN
 
370
C-----------------------------------------------------------------------
 
371
C     RECUR BACKWARD FOR REMAINDER OF I SEQUENCE AND ADD IN THE
 
372
C     K FUNCTIONS, SCALING THE I SEQUENCE DURING RECURRENCE TO KEEP
 
373
C     INTERMEDIATE ARITHMETIC ON SCALE NEAR EXPONENT EXTREMES.
 
374
C-----------------------------------------------------------------------
 
375
      S1R = CYR(1)
 
376
      S1I = CYI(1)
 
377
      S2R = CYR(2)
 
378
      S2I = CYI(2)
 
379
      CSR = CSRR(IFLAG)
 
380
      ASCLE = BRY(IFLAG)
 
381
      FN = DBLE(FLOAT(INU+IL))
 
382
      DO 290 I=1,IL
 
383
        C2R = S2R
 
384
        C2I = S2I
 
385
        S2R = S1R + (FN+FNF)*(RZR*C2R-RZI*C2I)
 
386
        S2I = S1I + (FN+FNF)*(RZR*C2I+RZI*C2R)
 
387
        S1R = C2R
 
388
        S1I = C2I
 
389
        FN = FN - 1.0D0
 
390
        C2R = S2R*CSR
 
391
        C2I = S2I*CSR
 
392
        CKR = C2R
 
393
        CKI = C2I
 
394
        C1R = YR(KK)
 
395
        C1I = YI(KK)
 
396
        IF (KODE.EQ.1) GO TO 280
 
397
        CALL ZS1S2(ZRR, ZRI, C1R, C1I, C2R, C2I, NW, ASC, ALIM, IUF)
 
398
        NZ = NZ + NW
 
399
  280   CONTINUE
 
400
        YR(KK) = C1R*CSPNR - C1I*CSPNI + C2R
 
401
        YI(KK) = C1R*CSPNI + C1I*CSPNR + C2I
 
402
        KK = KK - 1
 
403
        CSPNR = -CSPNR
 
404
        CSPNI = -CSPNI
 
405
        IF (IFLAG.GE.3) GO TO 290
 
406
        C2R = DABS(CKR)
 
407
        C2I = DABS(CKI)
 
408
        C2M = DMAX1(C2R,C2I)
 
409
        IF (C2M.LE.ASCLE) GO TO 290
 
410
        IFLAG = IFLAG + 1
 
411
        ASCLE = BRY(IFLAG)
 
412
        S1R = S1R*CSR
 
413
        S1I = S1I*CSR
 
414
        S2R = CKR
 
415
        S2I = CKI
 
416
        S1R = S1R*CSSR(IFLAG)
 
417
        S1I = S1I*CSSR(IFLAG)
 
418
        S2R = S2R*CSSR(IFLAG)
 
419
        S2I = S2I*CSSR(IFLAG)
 
420
        CSR = CSRR(IFLAG)
 
421
  290 CONTINUE
 
422
      RETURN
 
423
  300 CONTINUE
 
424
      NZ = -1
 
425
      RETURN
 
426
      END