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
zasyi.f
Go to the documentation of this file.
1  SUBROUTINE zasyi(ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM,
2  * alim)
3 C***BEGIN PROLOGUE ZASYI
4 C***REFER TO ZBESI,ZBESK
5 C
6 C ZASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
7 C MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z) IN THE
8 C REGION CABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN.
9 C NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1.
10 C
11 C***ROUTINES CALLED D1MACH,XZABS,ZDIV,XZEXP,ZMLT,XZSQRT
12 C***END PROLOGUE ZASYI
13 C COMPLEX AK1,CK,CONE,CS1,CS2,CZ,CZERO,DK,EZ,P1,RZ,S2,Y,Z
14  DOUBLE PRECISION aa, aez, ak, ak1i, ak1r, alim, arg, arm, atol,
15  * az, bb, bk, cki, ckr, conei, coner, cs1i, cs1r, cs2i, cs2r, czi,
16  * czr, dfnu, dki, dkr, dnu2, elim, ezi, ezr, fdn, fnu, pi, p1i,
17  * p1r, raz, rl, rtpi, rtr1, rzi, rzr, s, sgn, sqk, sti, str, s2i,
18  * s2r, tol, tzi, tzr, yi, yr, zeroi, zeror, zi, zr, d1mach, xzabs
19  INTEGER i, ib, il, inu, j, jl, k, kode, koded, m, n, nn, nz
20  dimension yr(n), yi(n)
21  DATA pi, rtpi /3.14159265358979324d0 , 0.159154943091895336d0 /
22  DATA zeror,zeroi,coner,conei / 0.0d0, 0.0d0, 1.0d0, 0.0d0 /
23 C
24  nz = 0
25  az = xzabs(zr,zi)
26  arm = 1.0d+3*d1mach(1)
27  rtr1 = dsqrt(arm)
28  il = min0(2,n)
29  dfnu = fnu + dble(float(n-il))
30 C-----------------------------------------------------------------------
31 C OVERFLOW TEST
32 C-----------------------------------------------------------------------
33  raz = 1.0d0/az
34  str = zr*raz
35  sti = -zi*raz
36  ak1r = rtpi*str*raz
37  ak1i = rtpi*sti*raz
38  CALL xzsqrt(ak1r, ak1i, ak1r, ak1i)
39  czr = zr
40  czi = zi
41  IF (kode.NE.2) go to 10
42  czr = zeror
43  czi = zi
44  10 CONTINUE
45  IF (dabs(czr).GT.elim) go to 100
46  dnu2 = dfnu + dfnu
47  koded = 1
48  IF ((dabs(czr).GT.alim) .AND. (n.GT.2)) go to 20
49  koded = 0
50  CALL xzexp(czr, czi, str, sti)
51  CALL zmlt(ak1r, ak1i, str, sti, ak1r, ak1i)
52  20 CONTINUE
53  fdn = 0.0d0
54  IF (dnu2.GT.rtr1) fdn = dnu2*dnu2
55  ezr = zr*8.0d0
56  ezi = zi*8.0d0
57 C-----------------------------------------------------------------------
58 C WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE
59 C FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE
60 C EXPANSION FOR THE IMAGINARY PART.
61 C-----------------------------------------------------------------------
62  aez = 8.0d0*az
63  s = tol/aez
64  jl = int(sngl(rl+rl)) + 2
65  p1r = zeror
66  p1i = zeroi
67  IF (zi.EQ.0.0d0) go to 30
68 C-----------------------------------------------------------------------
69 C CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF
70 C SIGNIFICANCE WHEN FNU OR N IS LARGE
71 C-----------------------------------------------------------------------
72  inu = int(sngl(fnu))
73  arg = (fnu-dble(float(inu)))*pi
74  inu = inu + n - il
75  ak = -dsin(arg)
76  bk = dcos(arg)
77  IF (zi.LT.0.0d0) bk = -bk
78  p1r = ak
79  p1i = bk
80  IF (mod(inu,2).EQ.0) go to 30
81  p1r = -p1r
82  p1i = -p1i
83  30 CONTINUE
84  DO 70 k=1,il
85  sqk = fdn - 1.0d0
86  atol = s*dabs(sqk)
87  sgn = 1.0d0
88  cs1r = coner
89  cs1i = conei
90  cs2r = coner
91  cs2i = conei
92  ckr = coner
93  cki = conei
94  ak = 0.0d0
95  aa = 1.0d0
96  bb = aez
97  dkr = ezr
98  dki = ezi
99  DO 40 j=1,jl
100  CALL zdiv(ckr, cki, dkr, dki, str, sti)
101  ckr = str*sqk
102  cki = sti*sqk
103  cs2r = cs2r + ckr
104  cs2i = cs2i + cki
105  sgn = -sgn
106  cs1r = cs1r + ckr*sgn
107  cs1i = cs1i + cki*sgn
108  dkr = dkr + ezr
109  dki = dki + ezi
110  aa = aa*dabs(sqk)/bb
111  bb = bb + aez
112  ak = ak + 8.0d0
113  sqk = sqk - ak
114  IF (aa.LE.atol) go to 50
115  40 CONTINUE
116  go to 110
117  50 CONTINUE
118  s2r = cs1r
119  s2i = cs1i
120  IF (zr+zr.GE.elim) go to 60
121  tzr = zr + zr
122  tzi = zi + zi
123  CALL xzexp(-tzr, -tzi, str, sti)
124  CALL zmlt(str, sti, p1r, p1i, str, sti)
125  CALL zmlt(str, sti, cs2r, cs2i, str, sti)
126  s2r = s2r + str
127  s2i = s2i + sti
128  60 CONTINUE
129  fdn = fdn + 8.0d0*dfnu + 4.0d0
130  p1r = -p1r
131  p1i = -p1i
132  m = n - il + k
133  yr(m) = s2r*ak1r - s2i*ak1i
134  yi(m) = s2r*ak1i + s2i*ak1r
135  70 CONTINUE
136  IF (n.LE.2) RETURN
137  nn = n
138  k = nn - 2
139  ak = dble(float(k))
140  str = zr*raz
141  sti = -zi*raz
142  rzr = (str+str)*raz
143  rzi = (sti+sti)*raz
144  ib = 3
145  DO 80 i=ib,nn
146  yr(k) = (ak+fnu)*(rzr*yr(k+1)-rzi*yi(k+1)) + yr(k+2)
147  yi(k) = (ak+fnu)*(rzr*yi(k+1)+rzi*yr(k+1)) + yi(k+2)
148  ak = ak - 1.0d0
149  k = k - 1
150  80 CONTINUE
151  IF (koded.EQ.0) RETURN
152  CALL xzexp(czr, czi, ckr, cki)
153  DO 90 i=1,nn
154  str = yr(i)*ckr - yi(i)*cki
155  yi(i) = yr(i)*cki + yi(i)*ckr
156  yr(i) = str
157  90 CONTINUE
158  RETURN
159  100 CONTINUE
160  nz = -1
161  RETURN
162  110 CONTINUE
163  nz=-2
164  RETURN
165  END