GNU Octave  3.8.0 A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
cwrsk.f
Go to the documentation of this file.
1  SUBROUTINE cwrsk(ZR, FNU, KODE, N, Y, NZ, CW, TOL, ELIM, ALIM)
2 C***BEGIN PROLOGUE CWRSK
3 C***REFER TO CBESI,CBESK
4 C
5 C CWRSK COMPUTES THE I BESSEL FUNCTION FOR RE(Z).GE.0.0 BY
6 C NORMALIZING THE I FUNCTION RATIOS FROM CRATI BY THE WRONSKIAN
7 C
8 C***ROUTINES CALLED CBKNU,CRATI,R1MACH
9 C***END PROLOGUE CWRSK
10  COMPLEX cinu, cscl, ct, cw, c1, c2, rct, st, y, zr
11  REAL act, acw, alim, ascle, elim, fnu, s1, s2, tol, yy
12  INTEGER i, kode, n, nw, nz
13  dimension y(n), cw(2)
14 C-----------------------------------------------------------------------
15 C I(FNU+I-1,Z) BY BACKWARD RECURRENCE FOR RATIOS
16 C Y(I)=I(FNU+I,Z)/I(FNU+I-1,Z) FROM CRATI NORMALIZED BY THE
17 C WRONSKIAN WITH K(FNU,Z) AND K(FNU+1,Z) FROM CBKNU.
18 C-----------------------------------------------------------------------
19  nz = 0
20  CALL cbknu(zr, fnu, kode, 2, cw, nw, tol, elim, alim)
21  IF (nw.NE.0) go to 50
22  CALL crati(zr, fnu, n, y, tol)
23 C-----------------------------------------------------------------------
24 C RECUR FORWARD ON I(FNU+1,Z) = R(FNU,Z)*I(FNU,Z),
25 C R(FNU+J-1,Z)=Y(J), J=1,...,N
26 C-----------------------------------------------------------------------
27  cinu = cmplx(1.0e0,0.0e0)
28  IF (kode.EQ.1) go to 10
29  yy = aimag(zr)
30  s1 = cos(yy)
31  s2 = sin(yy)
32  cinu = cmplx(s1,s2)
33  10 CONTINUE
34 C-----------------------------------------------------------------------
35 C ON LOW EXPONENT MACHINES THE K FUNCTIONS CAN BE CLOSE TO BOTH
36 C THE UNDER AND OVERFLOW LIMITS AND THE NORMALIZATION MUST BE
37 C SCALED TO PREVENT OVER OR UNDERFLOW. CUOIK HAS DETERMINED THAT
38 C THE RESULT IS ON SCALE.
39 C-----------------------------------------------------------------------
40  acw = cabs(cw(2))
41  ascle = 1.0e+3*r1mach(1)/tol
42  cscl = cmplx(1.0e0,0.0e0)
43  IF (acw.GT.ascle) go to 20
44  cscl = cmplx(1.0e0/tol,0.0e0)
45  go to 30
46  20 CONTINUE
47  ascle = 1.0e0/ascle
48  IF (acw.LT.ascle) go to 30
49  cscl = cmplx(tol,0.0e0)
50  30 CONTINUE
51  c1 = cw(1)*cscl
52  c2 = cw(2)*cscl
53  st = y(1)
54 C-----------------------------------------------------------------------
55 C CINU=CINU*(CONJG(CT)/CABS(CT))*(1.0E0/CABS(CT) PREVENTS
56 C UNDER- OR OVERFLOW PREMATURELY BY SQUARING CABS(CT)
57 C-----------------------------------------------------------------------
58  ct = zr*(c2+st*c1)
59  act = cabs(ct)
60  rct = cmplx(1.0e0/act,0.0e0)
61  ct = conjg(ct)*rct
62  cinu = cinu*rct*ct
63  y(1) = cinu*cscl
64  IF (n.EQ.1) RETURN
65  DO 40 i=2,n
66  cinu = st*cinu
67  st = y(i)
68  y(i) = cinu*cscl
69  40 CONTINUE
70  RETURN
71  50 CONTINUE
72  nz = -1
73  IF(nw.EQ.(-2)) nz=-2
74  RETURN
75  END