GNU Octave  4.0.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
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  znr = zr
206  zni = zi
207  csgnr = coner
208  csgni = conei
209  IF (zr.GE.0.0d0) go to 40
210  znr = -zr
211  zni = -zi
212 C-----------------------------------------------------------------------
213 C CALCULATE CSGN=EXP(FNU*PI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
214 C WHEN FNU IS LARGE
215 C-----------------------------------------------------------------------
216  inu = int(sngl(fnu))
217  arg = (fnu-dble(float(inu)))*pi
218  IF (zi.LT.0.0d0) arg = -arg
219  csgnr = dcos(arg)
220  csgni = dsin(arg)
221  IF (mod(inu,2).EQ.0) go to 40
222  csgnr = -csgnr
223  csgni = -csgni
224  40 CONTINUE
225  CALL zbinu(znr, zni, fnu, kode, n, cyr, cyi, nz, rl, fnul, tol,
226  * elim, alim)
227  IF (nz.LT.0) go to 120
228  IF (zr.GE.0.0d0) RETURN
229 C-----------------------------------------------------------------------
230 C ANALYTIC CONTINUATION TO THE LEFT HALF PLANE
231 C-----------------------------------------------------------------------
232  nn = n - nz
233  IF (nn.EQ.0) RETURN
234  rtol = 1.0d0/tol
235  ascle = d1mach(1)*rtol*1.0d+3
236  DO 50 i=1,nn
237 C STR = CYR(I)*CSGNR - CYI(I)*CSGNI
238 C CYI(I) = CYR(I)*CSGNI + CYI(I)*CSGNR
239 C CYR(I) = STR
240  aa = cyr(i)
241  bb = cyi(i)
242  atol = 1.0d0
243  IF (dmax1(dabs(aa),dabs(bb)).GT.ascle) go to 55
244  aa = aa*rtol
245  bb = bb*rtol
246  atol = tol
247  55 CONTINUE
248  str = aa*csgnr - bb*csgni
249  sti = aa*csgni + bb*csgnr
250  cyr(i) = str*atol
251  cyi(i) = sti*atol
252  csgnr = -csgnr
253  csgni = -csgni
254  50 CONTINUE
255  RETURN
256  120 CONTINUE
257  IF(nz.EQ.(-2)) go to 130
258  nz = 0
259  ierr=2
260  RETURN
261  130 CONTINUE
262  nz=0
263  ierr=5
264  RETURN
265  260 CONTINUE
266  nz=0
267  ierr=4
268  RETURN
269  END
subroutine zbinu(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL, TOL, ELIM, ALIM)
Definition: zbinu.f:1
std::string str(char sep= 'x') const
Definition: dim-vector.cc:63
std::string dimension(void) const
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:959
double precision function d1mach(i)
Definition: d1mach.f:1
F77_RET_T const double const double double * d
subroutine zbesi(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR)
Definition: zbesi.f:1