GNU Octave  3.8.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
cbesj.f
Go to the documentation of this file.
1  SUBROUTINE cbesj(Z, FNU, KODE, N, CY, NZ, IERR)
2 C***BEGIN PROLOGUE CBESJ
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 ON KODE=1, CBESJ COMPUTES AN N MEMBER SEQUENCE OF COMPLEX
13 C BESSEL FUNCTIONS CY(I)=J(FNU+I-1,Z) FOR REAL, NONNEGATIVE
14 C ORDERS FNU+I-1, I=1,...,N AND COMPLEX Z IN THE CUT PLANE
15 C -PI.LT.ARG(Z).LE.PI. ON KODE=2, CBESJ RETURNS THE SCALED
16 C FUNCTIONS
17 C
18 C CY(I)=EXP(-ABS(Y))*J(FNU+I-1,Z) I = 1,...,N , Y=AIMAG(Z)
19 C
20 C WHICH REMOVE THE EXPONENTIAL GROWTH IN BOTH THE UPPER AND
21 C LOWER HALF PLANES FOR Z TO INFINITY. DEFINITIONS AND NOTATION
22 C ARE FOUND IN THE NBS HANDBOOK OF MATHEMATICAL FUNCTIONS
23 C (REF. 1).
24 C
25 C INPUT
26 C Z - Z=CMPLX(X,Y), -PI.LT.ARG(Z).LE.PI
27 C FNU - ORDER OF INITIAL J FUNCTION, FNU.GE.0.0E0
28 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION
29 C KODE= 1 RETURNS
30 C CY(I)=J(FNU+I-1,Z), I=1,...,N
31 C = 2 RETURNS
32 C CY(I)=J(FNU+I-1,Z)*EXP(-ABS(Y)), I=1,...
33 C N - NUMBER OF MEMBERS OF THE SEQUENCE, N.GE.1
34 C
35 C OUTPUT
36 C CY - A COMPLEX VECTOR WHOSE FIRST N COMPONENTS CONTAIN
37 C VALUES FOR THE SEQUENCE
38 C CY(I)=J(FNU+I-1,Z) OR
39 C CY(I)=J(FNU+I-1,Z)*EXP(-ABS(Y)) I=1,...,N
40 C DEPENDING ON KODE, Y=AIMAG(Z).
41 C NZ - NUMBER OF COMPONENTS SET TO ZERO DUE TO UNDERFLOW,
42 C NZ= 0 , NORMAL RETURN
43 C NZ.GT.0 , LAST NZ COMPONENTS OF CY SET TO ZERO
44 C DUE TO UNDERFLOW, CY(I)=CMPLX(0.0,0.0),
45 C I = N-NZ+1,...,N
46 C IERR - ERROR FLAG
47 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
48 C IERR=1, INPUT ERROR - NO COMPUTATION
49 C IERR=2, OVERFLOW - NO COMPUTATION, AIMAG(Z)
50 C TOO LARGE ON KODE=1
51 C IERR=3, CABS(Z) OR FNU+N-1 LARGE - COMPUTATION DONE
52 C BUT LOSSES OF SIGNIFCANCE BY ARGUMENT
53 C REDUCTION PRODUCE LESS THAN HALF OF MACHINE
54 C ACCURACY
55 C IERR=4, CABS(Z) OR FNU+N-1 TOO LARGE - NO COMPUTA-
56 C TION BECAUSE OF COMPLETE LOSSES OF SIGNIFI-
57 C CANCE BY ARGUMENT REDUCTION
58 C IERR=5, ERROR - NO COMPUTATION,
59 C ALGORITHM TERMINATION CONDITION NOT MET
60 C
61 C***LONG DESCRIPTION
62 C
63 C THE COMPUTATION IS CARRIED OUT BY THE FORMULA
64 C
65 C J(FNU,Z)=EXP( FNU*PI*I/2)*I(FNU,-I*Z) AIMAG(Z).GE.0.0
66 C
67 C J(FNU,Z)=EXP(-FNU*PI*I/2)*I(FNU, I*Z) AIMAG(Z).LT.0.0
68 C
69 C WHERE I**2 = -1 AND I(FNU,Z) IS THE I BESSEL FUNCTION.
70 C
71 C FOR NEGATIVE ORDERS,THE FORMULA
72 C
73 C J(-FNU,Z) = J(FNU,Z)*COS(PI*FNU) - Y(FNU,Z)*SIN(PI*FNU)
74 C
75 C CAN BE USED. HOWEVER,FOR LARGE ORDERS CLOSE TO INTEGERS, THE
76 C THE FUNCTION CHANGES RADICALLY. WHEN FNU IS A LARGE POSITIVE
77 C INTEGER,THE MAGNITUDE OF J(-FNU,Z)=J(FNU,Z)*COS(PI*FNU) IS A
78 C LARGE NEGATIVE POWER OF TEN. BUT WHEN FNU IS NOT AN INTEGER,
79 C Y(FNU,Z) DOMINATES IN MAGNITUDE WITH A LARGE POSITIVE POWER OF
80 C TEN AND THE MOST THAT THE SECOND TERM CAN BE REDUCED IS BY
81 C UNIT ROUNDOFF FROM THE COEFFICIENT. THUS, WIDE CHANGES CAN
82 C OCCUR WITHIN UNIT ROUNDOFF OF A LARGE INTEGER FOR FNU. HERE,
83 C LARGE MEANS FNU.GT.CABS(Z).
84 C
85 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
86 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z OR FNU+N-1 IS
87 C LARGE, LOSSES OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR.
88 C CONSEQUENTLY, IF EITHER ONE EXCEEDS U1=SQRT(0.5/UR), THEN
89 C LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR FLAG
90 C IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF. ALSO
91 C IF EITHER IS LARGER THAN U2=0.5/UR, THEN ALL SIGNIFICANCE IS
92 C LOST AND IERR=4. IN ORDER TO USE THE INT FUNCTION, ARGUMENTS
93 C MUST BE FURTHER RESTRICTED NOT TO EXCEED THE LARGEST MACHINE
94 C INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF Z AND FNU+N-1 IS
95 C RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2, AND U3
96 C ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE PRECISION
97 C ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE PRECISION
98 C ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMITING IN
99 C THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT ONE CAN EXPECT
100 C TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES, NO DIGITS
101 C IN SINGLE AND ONLY 7 DIGITS IN DOUBLE PRECISION ARITHMETIC.
102 C SIMILAR CONSIDERATIONS HOLD FOR OTHER MACHINES.
103 C
104 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
105 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
106 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
107 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
108 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
109 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
110 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
111 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
112 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
113 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
114 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
115 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
116 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
117 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
118 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
119 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
120 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
121 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
122 C OR -PI/2+P.
123 C
124 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
125 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
126 C COMMERCE, 1955.
127 C
128 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
129 C BY D. E. AMOS, SAND83-0083, MAY, 1983.
130 C
131 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
132 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
133 C
134 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
135 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
136 C 1018, MAY, 1985
137 C
138 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
139 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
140 C MATH. SOFTWARE, 1986
141 C
142 C***ROUTINES CALLED CBINU,I1MACH,R1MACH
143 C***END PROLOGUE CBESJ
144 C
145  COMPLEX ci, csgn, cy, z, zn
146  REAL aa, alim, arg, dig, elim, fnu, fnul, hpi, rl, r1, r1m5, r2,
147  * tol, yy, r1mach, az, fn, bb, ascle, rtol, atol
148  INTEGER i, ierr, inu, inuh, ir, kode, k1, k2, n, nl, nz, i1mach, k
149  dimension cy(n)
150  DATA hpi /1.57079632679489662e0/
151 C
152 C***FIRST EXECUTABLE STATEMENT CBESJ
153  ierr = 0
154  nz=0
155  IF (fnu.LT.0.0e0) ierr=1
156  IF (kode.LT.1 .OR. kode.GT.2) ierr=1
157  IF (n.LT.1) ierr=1
158  IF (ierr.NE.0) RETURN
159 C-----------------------------------------------------------------------
160 C SET PARAMETERS RELATED TO MACHINE CONSTANTS.
161 C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
162 C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
163 C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND
164 C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR
165 C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
166 C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
167 C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
168 C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU.
169 C-----------------------------------------------------------------------
170  tol = amax1(r1mach(4),1.0e-18)
171  k1 = i1mach(12)
172  k2 = i1mach(13)
173  r1m5 = r1mach(5)
174  k = min0(iabs(k1),iabs(k2))
175  elim = 2.303e0*(float(k)*r1m5-3.0e0)
176  k1 = i1mach(11) - 1
177  aa = r1m5*float(k1)
178  dig = amin1(aa,18.0e0)
179  aa = aa*2.303e0
180  alim = elim + amax1(-aa,-41.45e0)
181  rl = 1.2e0*dig + 3.0e0
182  fnul = 10.0e0 + 6.0e0*(dig-3.0e0)
183  ci = cmplx(0.0e0,1.0e0)
184  yy = aimag(z)
185  az = cabs(z)
186 C-----------------------------------------------------------------------
187 C TEST FOR RANGE
188 C-----------------------------------------------------------------------
189  aa = 0.5e0/tol
190  bb=float(i1mach(9))*0.5e0
191  aa=amin1(aa,bb)
192  fn=fnu+float(n-1)
193  IF(az.GT.aa) go to 140
194  IF(fn.GT.aa) go to 140
195  aa=sqrt(aa)
196  IF(az.GT.aa) ierr=3
197  IF(fn.GT.aa) ierr=3
198 C-----------------------------------------------------------------------
199 C CALCULATE CSGN=EXP(FNU*HPI*I) TO MINIMIZE LOSSES OF SIGNIFICANCE
200 C WHEN FNU IS LARGE
201 C-----------------------------------------------------------------------
202  inu = int(fnu)
203  inuh = inu/2
204  ir = inu - 2*inuh
205  arg = (fnu-float(inu-ir))*hpi
206  r1 = cos(arg)
207  r2 = sin(arg)
208  csgn = cmplx(r1,r2)
209  IF (mod(inuh,2).EQ.1) csgn = -csgn
210 C-----------------------------------------------------------------------
211 C ZN IS IN THE RIGHT HALF PLANE
212 C-----------------------------------------------------------------------
213  zn = -z*ci
214  IF (yy.GE.0.0e0) go to 40
215  zn = -zn
216  csgn = conjg(csgn)
217  ci = conjg(ci)
218  40 CONTINUE
219  CALL cbinu(zn, fnu, kode, n, cy, nz, rl, fnul, tol, elim, alim)
220  IF (nz.LT.0) go to 120
221  nl = n - nz
222  IF (nl.EQ.0) RETURN
223  rtol = 1.0e0/tol
224  ascle = r1mach(1)*rtol*1.0e+3
225  DO 50 i=1,nl
226 C CY(I)=CY(I)*CSGN
227  zn=cy(i)
228  aa=REAL(zn)
229  bb=aimag(zn)
230  atol=1.0e0
231  IF (amax1(abs(aa),abs(bb)).GT.ascle) go to 55
232  zn = zn*cmplx(rtol,0.0e0)
233  atol = tol
234  55 CONTINUE
235  zn = zn*csgn
236  cy(i) = zn*cmplx(atol,0.0e0)
237  csgn = csgn*ci
238  50 CONTINUE
239  RETURN
240  120 CONTINUE
241  IF(nz.EQ.(-2)) go to 130
242  nz = 0
243  ierr = 2
244  RETURN
245  130 CONTINUE
246  nz=0
247  ierr=5
248  RETURN
249  140 CONTINUE
250  nz=0
251  ierr=4
252  RETURN
253  END