|
3217
|
1 SUBROUTINE ZWRSK(ZRR, ZRI, FNU, KODE, N, YR, YI, NZ, CWR, CWI, |
|
|
2 * TOL, ELIM, ALIM) |
|
|
3 C***BEGIN PROLOGUE ZWRSK |
|
|
4 C***REFER TO ZBESI,ZBESK |
|
|
5 C |
|
|
6 C ZWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY |
|
|
7 C NORMALIZING THE I FUNCTION RATIOS FROM ZRATI BY THE WRONSKIAN |
|
|
8 C |
|
|
9 C***ROUTINES CALLED D1MACH,ZBKNU,ZRATI,XZABS |
|
|
10 C***END PROLOGUE ZWRSK |
|
|
11 C COMPLEX CINU,CSCL,CT,CW,C1,C2,RCT,ST,Y,ZR |
|
|
12 DOUBLE PRECISION ACT, ACW, ALIM, ASCLE, CINUI, CINUR, CSCLR, CTI, |
|
|
13 * CTR, CWI, CWR, C1I, C1R, C2I, C2R, ELIM, FNU, PTI, PTR, RACT, |
|
|
14 * STI, STR, TOL, YI, YR, ZRI, ZRR, XZABS, D1MACH |
|
|
15 INTEGER I, KODE, N, NW, NZ |
|
|
16 DIMENSION YR(N), YI(N), CWR(2), CWI(2) |
|
|
17 C----------------------------------------------------------------------- |
|
|
18 C I(FNU+I-1,Z) BY BACKWARD RECURRENCE FOR RATIOS |
|
|
19 C Y(I)=I(FNU+I,Z)/I(FNU+I-1,Z) FROM CRATI NORMALIZED BY THE |
|
|
20 C WRONSKIAN WITH K(FNU,Z) AND K(FNU+1,Z) FROM CBKNU. |
|
|
21 C----------------------------------------------------------------------- |
|
|
22 NZ = 0 |
|
|
23 CALL ZBKNU(ZRR, ZRI, FNU, KODE, 2, CWR, CWI, NW, TOL, ELIM, ALIM) |
|
|
24 IF (NW.NE.0) GO TO 50 |
|
|
25 CALL ZRATI(ZRR, ZRI, FNU, N, YR, YI, TOL) |
|
|
26 C----------------------------------------------------------------------- |
|
|
27 C RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z), |
|
|
28 C R(FNU+J-1,Z)=Y(J), J=1,...,N |
|
|
29 C----------------------------------------------------------------------- |
|
|
30 CINUR = 1.0D0 |
|
|
31 CINUI = 0.0D0 |
|
|
32 IF (KODE.EQ.1) GO TO 10 |
|
|
33 CINUR = DCOS(ZRI) |
|
|
34 CINUI = DSIN(ZRI) |
|
|
35 10 CONTINUE |
|
|
36 C----------------------------------------------------------------------- |
|
|
37 C ON LOW EXPONENT MACHINES THE K FUNCTIONS CAN BE CLOSE TO BOTH |
|
|
38 C THE UNDER AND OVERFLOW LIMITS AND THE NORMALIZATION MUST BE |
|
|
39 C SCALED TO PREVENT OVER OR UNDERFLOW. CUOIK HAS DETERMINED THAT |
|
|
40 C THE RESULT IS ON SCALE. |
|
|
41 C----------------------------------------------------------------------- |
|
|
42 ACW = XZABS(CWR(2),CWI(2)) |
|
|
43 ASCLE = 1.0D+3*D1MACH(1)/TOL |
|
|
44 CSCLR = 1.0D0 |
|
|
45 IF (ACW.GT.ASCLE) GO TO 20 |
|
|
46 CSCLR = 1.0D0/TOL |
|
|
47 GO TO 30 |
|
|
48 20 CONTINUE |
|
|
49 ASCLE = 1.0D0/ASCLE |
|
|
50 IF (ACW.LT.ASCLE) GO TO 30 |
|
|
51 CSCLR = TOL |
|
|
52 30 CONTINUE |
|
|
53 C1R = CWR(1)*CSCLR |
|
|
54 C1I = CWI(1)*CSCLR |
|
|
55 C2R = CWR(2)*CSCLR |
|
|
56 C2I = CWI(2)*CSCLR |
|
|
57 STR = YR(1) |
|
|
58 STI = YI(1) |
|
|
59 C----------------------------------------------------------------------- |
|
|
60 C CINU=CINU*(CONJG(CT)/CABS(CT))*(1.0D0/CABS(CT) PREVENTS |
|
|
61 C UNDER- OR OVERFLOW PREMATURELY BY SQUARING CABS(CT) |
|
|
62 C----------------------------------------------------------------------- |
|
|
63 PTR = STR*C1R - STI*C1I |
|
|
64 PTI = STR*C1I + STI*C1R |
|
|
65 PTR = PTR + C2R |
|
|
66 PTI = PTI + C2I |
|
|
67 CTR = ZRR*PTR - ZRI*PTI |
|
|
68 CTI = ZRR*PTI + ZRI*PTR |
|
|
69 ACT = XZABS(CTR,CTI) |
|
|
70 RACT = 1.0D0/ACT |
|
|
71 CTR = CTR*RACT |
|
|
72 CTI = -CTI*RACT |
|
|
73 PTR = CINUR*RACT |
|
|
74 PTI = CINUI*RACT |
|
|
75 CINUR = PTR*CTR - PTI*CTI |
|
|
76 CINUI = PTR*CTI + PTI*CTR |
|
|
77 YR(1) = CINUR*CSCLR |
|
|
78 YI(1) = CINUI*CSCLR |
|
|
79 IF (N.EQ.1) RETURN |
|
|
80 DO 40 I=2,N |
|
|
81 PTR = STR*CINUR - STI*CINUI |
|
|
82 CINUI = STR*CINUI + STI*CINUR |
|
|
83 CINUR = PTR |
|
|
84 STR = YR(I) |
|
|
85 STI = YI(I) |
|
|
86 YR(I) = CINUR*CSCLR |
|
|
87 YI(I) = CINUI*CSCLR |
|
|
88 40 CONTINUE |
|
|
89 RETURN |
|
|
90 50 CONTINUE |
|
|
91 NZ = -1 |
|
|
92 IF(NW.EQ.(-2)) NZ=-2 |
|
|
93 RETURN |
|
|
94 END |
