GNU Octave  8.1.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
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