GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
double precision function d1mach(i)
Definition: d1mach.f:23
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:932
double precision function xzabs(ZR, ZI)
Definition: xzabs.f:2
subroutine xzexp(AR, AI, BR, BI)
Definition: xzexp.f:2
subroutine xzsqrt(AR, AI, BR, BI)
Definition: xzsqrt.f:2
subroutine zasyi(ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM, ALIM)
Definition: zasyi.f:3
subroutine zdiv(AR, AI, BR, BI, CR, CI)
Definition: zdiv.f:2
subroutine zmlt(AR, AI, BR, BI, CR, CI)
Definition: zmlt.f:2