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