GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
zbesi.f
Go to the documentation of this file.
1  SUBROUTINE zbesi(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR)
2 C***BEGIN PROLOGUE ZBESI
3 C***DATE WRITTEN 830501 (YYMMDD)
4 C***REVISION DATE 890801 (YYMMDD)
5 C***CATEGORY NO. B5K
6 C***KEYWORDS I-BESSEL FUNCTION,COMPLEX BESSEL FUNCTION,
7 C MODIFIED BESSEL FUNCTION OF THE FIRST KIND
8 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
9 C***PURPOSE TO COMPUTE I-BESSEL FUNCTIONS OF COMPLEX ARGUMENT
10 C***DESCRIPTION
11 C
12 C ***A DOUBLE PRECISION ROUTINE***
13 C ON KODE=1, ZBESI COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
14 C BESSEL FUNCTIONS CY(J)=I(FNU+J-1,Z) FOR REAL, NONNEGATIVE
15 C ORDERS FNU+J-1, J=1,...,N AND COMPLEX Z IN THE CUT PLANE
16 C -PI.LT.ARG(Z).LE.PI. ON KODE=2, ZBESI RETURNS THE SCALED
17 C FUNCTIONS
18 C
19 C CY(J)=EXP(-ABS(X))*I(FNU+J-1,Z) J = 1,...,N , X=REAL(Z)
20 C
21 C WITH THE EXPONENTIAL GROWTH REMOVED IN BOTH THE LEFT AND
22 C RIGHT HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION
23 C ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS
24 C (REF. 1).
25 C
26 C INPUT ZR,ZI,FNU ARE DOUBLE PRECISION
27 C ZR,ZI - Z=CMPLX(ZR,ZI), -PI.LT.ARG(Z).LE.PI
28 C FNU - ORDER OF INITIAL I FUNCTION, FNU.GE.0.0D0
29 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION
30 C KODE= 1 RETURNS
31 C CY(J)=I(FNU+J-1,Z), J=1,...,N
32 C = 2 RETURNS
33 C CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X)), J=1,...,N
34 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
35 C
36 C OUTPUT CYR,CYI ARE DOUBLE PRECISION
37 C CYR,CYI- DOUBLE PRECISION VECTORS WHOSE FIRST N COMPONENTS
38 C CONTAIN REAL AND IMAGINARY PARTS FOR THE SEQUENCE
39 C CY(J)=I(FNU+J-1,Z) OR
40 C CY(J)=I(FNU+J-1,Z)*EXP(-ABS(X)) J=1,...,N
41 C DEPENDING ON KODE, X=REAL(Z)
42 C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW,
43 C NZ= 0 , NORMAL RETURN
44 C NZ.GT.0 , LAST NZ COMPONENTS OF CY SET TO ZERO
45 C TO UNDERFLOW, CY(J)=CMPLX(0.0D0,0.0D0)
46 C J = N-NZ+1,...,N
47 C IERR - ERROR FLAG
48 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
49 C IERR=1, INPUT ERROR - NO COMPUTATION
50 C IERR=2, OVERFLOW - NO COMPUTATION, REAL(Z) TOO
51 C LARGE ON KODE=1
52 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
53 C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
54 C REDUCTION PRODUCE LESS THAN HALF OF MACHINE
55 C ACCURACY
56 C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
57 C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
58 C CANCE BY ARGUMENT REDUCTION
59 C IERR=5, ERROR - NO COMPUTATION,
60 C ALGORITHM TERMINATION CONDITION NOT MET
61 C
62 C***LONG DESCRIPTION
63 C
64 C THE COMPUTATION IS CARRIED OUT BY THE POWER SERIES FOR
65 C SMALL CABS(Z), THE ASYMPTOTIC EXPANSION FOR LARGE CABS(Z),
66 C THE MILLER ALGORITHM NORMALIZED BY THE WRONSKIAN AND A
67 C NEUMANN SERIES FOR IMTERMEDIATE MAGNITUDES, AND THE
68 C UNIFORM ASYMPTOTIC EXPANSIONS FOR I(FNU,Z) AND J(FNU,Z)
69 C FOR LARGE ORDERS. BACKWARD RECURRENCE IS USED TO GENERATE
70 C SEQUENCES OR REDUCE ORDERS WHEN NECESSARY.
71 C
72 C THE CALCULATIONS ABOVE ARE DONE IN THE RIGHT HALF PLANE AND
73 C CONTINUED INTO THE LEFT HALF PLANE BY THE FORMULA
74 C
75 C I(FNU,Z*EXP(M*PI)) = EXP(M*PI*FNU)*I(FNU,Z) REAL(Z).GT.0.0
76 C M = +I OR -I, I**2=-1
77 C
78 C FOR NEGATIVE ORDERS,THE FORMULA
79 C
80 C I(-FNU,Z) = I(FNU,Z) + (2/PI)*SIN(PI*FNU)*K(FNU,Z)
81 C
82 C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO INTEGERS, THE
83 C THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE POSITIVE
84 C INTEGER,THE MAGNITUDE OF I(-FNU,Z)=I(FNU,Z) IS A LARGE
85 C NEGATIVE POWER OF TEN. BUT WHEN FNU IS NOT AN INTEGER,
86 C K(FNU,Z) DOMINATES IN MAGNITUDE WITH A LARGE POSITIVE POWER OF
87 C TEN AND THE MOST THAT THE SECOND TERM CAN BE REDUCED IS BY
88 C UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, WIDE CHANGES CAN
89 C OCCUR WITHIN UNIT ROUNDOFF OF A LARGE INTEGER FOR FNU. HERE,
90 C LARGE MEANS FNU.GT.CABS(Z).
91 C
92 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
93 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
94 C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
95 C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
96 C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
97 C IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
98 C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
99 C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
100 C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
101 C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
102 C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
103 C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
104 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
105 C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
106 C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
107 C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
108 C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
109 C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
110 C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
111 C
112 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
113 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
114 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
115 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
116 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
117 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
118 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
119 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
120 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
121 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
122 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
123 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
124 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
125 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
126 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
127 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
128 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
129 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
130 C OR -PI/2+P.
131 C
132 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
133 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
134 C COMMERCE, 1955.
135 C
136 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
137 C BY D. E. AMOS, SAND83-0083, MAY, 1983.
138 C
139 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
140 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
141 C
142 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
143 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
144 C 1018, MAY, 1985
145 C
146 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
147 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
148 C MATH. SOFTWARE, 1986
149 C
150 C***ROUTINES CALLED ZBINU,I1MACH,D1MACH
151 C***END PROLOGUE ZBESI
152 C COMPLEX CONE,CSGN,CW,CY,CZERO,Z,ZN
153  DOUBLE PRECISION AA, ALIM, ARG, CONEI, CONER, CSGNI, CSGNR, CYI,
154  * CYR, DIG, ELIM, FNU, FNUL, PI, RL, R1M5, STR, TOL, ZI, ZNI, ZNR,
155  * ZR, D1MACH, AZ, BB, FN, XZABS, ASCLE, RTOL, ATOL, STI
156  INTEGER I, IERR, INU, K, KODE, K1,K2,N,NZ,NN, I1MACH
157  dimension cyr(n), cyi(n)
158  DATA pi /3.14159265358979324d0/
159  DATA coner, conei /1.0d0,0.0d0/
160 C
161 C***FIRST EXECUTABLE STATEMENT ZBESI
162  ierr = 0
163  nz=0
164  IF (fnu.LT.0.0d0) ierr=1
165  IF (kode.LT.1 .OR. kode.GT.2) ierr=1
166  IF (n.LT.1) ierr=1
167  IF (ierr.NE.0) RETURN
168 C-----------------------------------------------------------------------
169 C SET PARAMETERS RELATED TO MACHINE CONSTANTS.
170 C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
171 C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
172 C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND
173 C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR
174 C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
175 C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
176 C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
177 C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU.
178 C-----------------------------------------------------------------------
179  tol = dmax1(d1mach(4),1.0d-18)
180  k1 = i1mach(15)
181  k2 = i1mach(16)
182  r1m5 = d1mach(5)
183  k = min0(iabs(k1),iabs(k2))
184  elim = 2.303d0*(dble(float(k))*r1m5-3.0d0)
185  k1 = i1mach(14) - 1
186  aa = r1m5*dble(float(k1))
187  dig = dmin1(aa,18.0d0)
188  aa = aa*2.303d0
189  alim = elim + dmax1(-aa,-41.45d0)
190  rl = 1.2d0*dig + 3.0d0
191  fnul = 10.0d0 + 6.0d0*(dig-3.0d0)
192 C-----------------------------------------------------------------------------
193 C TEST FOR PROPER RANGE
194 C-----------------------------------------------------------------------
195  az = xzabs(zr,zi)
196  fn = fnu+dble(float(n-1))
197  aa = 0.5d0/tol
198  bb=dble(float(i1mach(9)))*0.5d0
199  aa = dmin1(aa,bb)
200  IF (az.GT.aa) GO TO 260
201  IF (fn.GT.aa) GO TO 260
202  aa = dsqrt(aa)
203  IF (az.GT.aa) ierr=3
204  IF (fn.GT.aa) ierr=3
205  35 CONTINUE
206  znr = zr
207  zni = zi
208  csgnr = coner
209  csgni = conei
210  IF (zr.GE.0.0d0) GO TO 40
211  znr = -zr
212  zni = -zi
213 C-----------------------------------------------------------------------
214 C CALCULATE CSGN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
215 C WHEN FNU IS LARGE
216 C-----------------------------------------------------------------------
217  inu = int(sngl(fnu))
218  arg = (fnu-dble(float(inu)))*pi
219  IF (zi.LT.0.0d0) arg = -arg
220  csgnr = dcos(arg)
221  csgni = dsin(arg)
222  IF (mod(inu,2).EQ.0) GO TO 40
223  csgnr = -csgnr
224  csgni = -csgni
225  40 CONTINUE
226  CALL zbinu(znr, zni, fnu, kode, n, cyr, cyi, nz, rl, fnul, tol,
227  * elim, alim)
228  IF (nz.LT.0) GO TO 120
229  IF (zr.GE.0.0d0) RETURN
230 C-----------------------------------------------------------------------
231 C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE
232 C-----------------------------------------------------------------------
233  nn = n - nz
234  IF (nn.EQ.0) RETURN
235  rtol = 1.0d0/tol
236  ascle = d1mach(1)*rtol*1.0d+3
237  DO 50 i=1,nn
238 C STR = CYR(I)*CSGNR - CYI(I)*CSGNI
239 C CYI(I) = CYR(I)*CSGNI + CYI(I)*CSGNR
240 C CYR(I) = STR
241  aa = cyr(i)
242  bb = cyi(i)
243  atol = 1.0d0
244  IF (dmax1(dabs(aa),dabs(bb)).GT.ascle) GO TO 55
245  aa = aa*rtol
246  bb = bb*rtol
247  atol = tol
248  55 CONTINUE
249  str = aa*csgnr - bb*csgni
250  sti = aa*csgni + bb*csgnr
251  cyr(i) = str*atol
252  cyi(i) = sti*atol
253  csgnr = -csgnr
254  csgni = -csgni
255  50 CONTINUE
256  RETURN
257  120 CONTINUE
258  IF(nz.EQ.(-2)) GO TO 130
259  nz = 0
260  ierr=2
261  RETURN
262  130 CONTINUE
263  nz=0
264  ierr=5
265  RETURN
266  260 CONTINUE
267  ierr=4
268  GO TO 35
269  END
T mod(T x, T y)
Definition: lo-mappers.h:294
subroutine zbesi(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR)
Definition: zbesi.f:2
subroutine zbinu(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL, TOL, ELIM, ALIM)
Definition: zbinu.f:3