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
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  cii = 1.0d0
204  inu = int(sngl(fnu))
205  inuh = inu/2
206  ir = inu - 2*inuh
207  arg = (fnu-dble(float(inu-ir)))*hpi
208  csgnr = dcos(arg)
209  csgni = dsin(arg)
210  IF (mod(inuh,2).EQ.0) go to 40
211  csgnr = -csgnr
212  csgni = -csgni
213  40 CONTINUE
214 C-----------------------------------------------------------------------
215 C ZN IS IN THE RIGHT HALF PLANE
216 C-----------------------------------------------------------------------
217  znr = zi
218  zni = -zr
219  IF (zi.GE.0.0d0) go to 50
220  znr = -znr
221  zni = -zni
222  csgni = -csgni
223  cii = -cii
224  50 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 130
228  nl = n - nz
229  IF (nl.EQ.0) RETURN
230  rtol = 1.0d0/tol
231  ascle = d1mach(1)*rtol*1.0d+3
232  DO 60 i=1,nl
233 C STR = CYR(I)*CSGNR - CYI(I)*CSGNI
234 C CYI(I) = CYR(I)*CSGNI + CYI(I)*CSGNR
235 C CYR(I) = STR
236  aa = cyr(i)
237  bb = cyi(i)
238  atol = 1.0d0
239  IF (dmax1(dabs(aa),dabs(bb)).GT.ascle) go to 55
240  aa = aa*rtol
241  bb = bb*rtol
242  atol = tol
243  55 CONTINUE
244  str = aa*csgnr - bb*csgni
245  sti = aa*csgni + bb*csgnr
246  cyr(i) = str*atol
247  cyi(i) = sti*atol
248  str = -csgni*cii
249  csgni = csgnr*cii
250  csgnr = str
251  60 CONTINUE
252  RETURN
253  130 CONTINUE
254  IF(nz.EQ.(-2)) go to 140
255  nz = 0
256  ierr = 2
257  RETURN
258  140 CONTINUE
259  nz=0
260  ierr=5
261  RETURN
262  260 CONTINUE
263  nz=0
264  ierr=4
265  RETURN
266  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
subroutine zbesj(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, IERR)
Definition: zbesj.f:1
F77_RET_T const double const double double * d