GNU Octave  3.8.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
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