GNU Octave  3.8.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
zwrsk.f
Go to the documentation of this file.
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