GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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  35 CONTINUE
203  inu = int(fnu)
204  inuh = inu/2
205  ir = inu - 2*inuh
206  arg = (fnu-float(inu-ir))*hpi
207  r1 = cos(arg)
208  r2 = sin(arg)
209  csgn = cmplx(r1,r2)
210  IF (mod(inuh,2).EQ.1) csgn = -csgn
211 C-----------------------------------------------------------------------
212 C ZN IS IN THE RIGHT HALF PLANE
213 C-----------------------------------------------------------------------
214  zn = -z*ci
215  IF (yy.GE.0.0e0) GO TO 40
216  zn = -zn
217  csgn = conjg(csgn)
218  ci = conjg(ci)
219  40 CONTINUE
220  CALL cbinu(zn, fnu, kode, n, cy, nz, rl, fnul, tol, elim, alim)
221  IF (nz.LT.0) GO TO 120
222  nl = n - nz
223  IF (nl.EQ.0) RETURN
224  rtol = 1.0e0/tol
225  ascle = r1mach(1)*rtol*1.0e+3
226  DO 50 i=1,nl
227 C CY(I)=CY(I)*CSGN
228  zn=cy(i)
229  aa=real(zn)
230  bb=aimag(zn)
231  atol=1.0e0
232  IF (amax1(abs(aa),abs(bb)).GT.ascle) GO TO 55
233  zn = zn*cmplx(rtol,0.0e0)
234  atol = tol
235  55 CONTINUE
236  zn = zn*csgn
237  cy(i) = zn*cmplx(atol,0.0e0)
238  csgn = csgn*ci
239  50 CONTINUE
240  RETURN
241  120 CONTINUE
242  IF(nz.EQ.(-2)) GO TO 130
243  nz = 0
244  ierr = 2
245  RETURN
246  130 CONTINUE
247  nz=0
248  ierr=5
249  RETURN
250  140 CONTINUE
251  ierr=4
252  GO TO 35
253  END
double complex cmplx
Definition: Faddeeva.cc:217
subroutine cbesj(Z, FNU, KODE, N, CY, NZ, IERR)
Definition: cbesj.f:2
subroutine cbinu(Z, FNU, KODE, N, CY, NZ, RL, FNUL, TOL, ELIM, ALIM)
Definition: cbinu.f:3
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:137
octave_int< T > mod(const octave_int< T > &x, const octave_int< T > &y)
Definition: oct-inttypes.h:932
static T abs(T x)
Definition: pr-output.cc:1678