GNU Octave  6.2.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
cbiry.f
Go to the documentation of this file.
1  SUBROUTINE cbiry(Z, ID, KODE, BI, IERR)
2 C***BEGIN PROLOGUE CBIRY
3 C***DATE WRITTEN 830501 (YYMMDD)
4 C***REVISION DATE 890801 (YYMMDD)
5 C***CATEGORY NO. B5K
6 C***KEYWORDS AIRY FUNCTION,BESSEL FUNCTIONS OF ORDER ONE THIRD
7 C***AUTHOR AMOS, DONALD E., SANDIA NATIONAL LABORATORIES
8 C***PURPOSE TO COMPUTE AIRY FUNCTIONS BI(Z) AND DBI(Z) FOR COMPLEX Z
9 C***DESCRIPTION
10 C
11 C ON KODE=1, CBIRY COMPUTES THE COMPLEX AIRY FUNCTION BI(Z) OR
12 C ITS DERIVATIVE DBI(Z)/DZ ON ID=0 OR ID=1 RESPECTIVELY. ON
13 C KODE=2, A SCALING OPTION CEXP(-AXZTA)*BI(Z) OR CEXP(-AXZTA)*
14 C DBI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL BEHAVIOR IN
15 C BOTH THE LEFT AND RIGHT HALF PLANES WHERE
16 C ZTA=(2/3)*Z*CSQRT(Z)=CMPLX(XZTA,YZTA) AND AXZTA=ABS(XZTA).
17 C DEFINITIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF
18 C MATHEMATICAL FUNCTIONS (REF. 1).
19 C
20 C INPUT
21 C Z - Z=CMPLX(X,Y)
22 C ID - ORDER OF DERIVATIVE, ID=0 OR ID=1
23 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION
24 C KODE= 1 RETURNS
25 C BI=BI(Z) ON ID=0 OR
26 C BI=DBI(Z)/DZ ON ID=1
27 C = 2 RETURNS
28 C BI=CEXP(-AXZTA)*BI(Z) ON ID=0 OR
29 C BI=CEXP(-AXZTA)*DBI(Z)/DZ ON ID=1 WHERE
30 C ZTA=(2/3)*Z*CSQRT(Z)=CMPLX(XZTA,YZTA)
31 C AND AXZTA=ABS(XZTA)
32 C
33 C OUTPUT
34 C BI - COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND
35 C KODE
36 C IERR - ERROR FLAG
37 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
38 C IERR=1, INPUT ERROR - NO COMPUTATION
39 C IERR=2, OVERFLOW - NO COMPUTATION, REAL(Z)
40 C TOO LARGE WITH KODE=1
41 C IERR=3, CABS(Z) LARGE - COMPUTATION COMPLETED
42 C LOSSES OF SIGNIFCANCE BY ARGUMENT REDUCTION
43 C PRODUCE LESS THAN HALF OF MACHINE ACCURACY
44 C IERR=4, CABS(Z) TOO LARGE - NO COMPUTATION
45 C COMPLETE LOSS OF ACCURACY BY ARGUMENT
46 C REDUCTION
47 C IERR=5, ERROR - NO COMPUTATION,
48 C ALGORITHM TERMINATION CONDITION NOT MET
49 C
50 C***LONG DESCRIPTION
51 C
52 C BI AND DBI ARE COMPUTED FOR CABS(Z).GT.1.0 FROM THE I BESSEL
53 C FUNCTIONS BY
54 C
55 C BI(Z)=C*SQRT(Z)*( I(-1/3,ZTA) + I(1/3,ZTA) )
56 C DBI(Z)=C * Z * ( I(-2/3,ZTA) + I(2/3,ZTA) )
57 C C=1.0/SQRT(3.0)
58 C ZTA=(2/3)*Z**(3/2)
59 C
60 C WITH THE POWER SERIES FOR CABS(Z).LE.1.0.
61 C
62 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
63 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z IS LARGE, LOSSES
64 C OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. CONSEQUENTLY, IF
65 C THE MAGNITUDE OF ZETA=(2/3)*Z**1.5 EXCEEDS U1=SQRT(0.5/UR),
66 C THEN LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR
67 C FLAG IERR=3 IS TRIGGERED WHERE UR=R1MACH(4)=UNIT ROUNDOFF.
68 C ALSO, IF THE MAGNITUDE OF ZETA IS LARGER THAN U2=0.5/UR, THEN
69 C ALL SIGNIFICANCE IS LOST AND IERR=4. IN ORDER TO USE THE INT
70 C FUNCTION, ZETA MUST BE FURTHER RESTRICTED NOT TO EXCEED THE
71 C LARGEST INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF ZETA
72 C MUST BE RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2,
73 C AND U3 ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE
74 C PRECISION ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE
75 C PRECISION ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMIT-
76 C ING IN THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT THE MAG-
77 C NITUDE OF Z CANNOT EXCEED 3.1E+4 IN SINGLE AND 2.1E+6 IN
78 C DOUBLE PRECISION ARITHMETIC. THIS ALSO MEANS THAT ONE CAN
79 C EXPECT TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES,
80 C NO DIGITS IN SINGLE PRECISION AND ONLY 7 DIGITS IN DOUBLE
81 C PRECISION ARITHMETIC.
82 C
83 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
84 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
85 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
86 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
87 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
88 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
89 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
90 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
91 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
92 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
93 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
94 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
95 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
96 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
97 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
98 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
99 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
100 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
101 C OR -PI/2+P.
102 C
103 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
104 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
105 C COMMERCE, 1955.
106 C
107 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
108 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
109 C
110 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
111 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
112 C 1018, MAY, 1985
113 C
114 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
115 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
116 C MATH. SOFTWARE, 1986
117 C
118 C***ROUTINES CALLED CBINU,I1MACH,R1MACH
119 C***END PROLOGUE CBIRY
120  COMPLEX BI, CONE, CSQ, CY, S1, S2, TRM1, TRM2, Z, ZTA, Z3
121  REAL AA, AD, AK, ALIM, ATRM, AZ, AZ3, BB, BK, CK, COEF, C1, C2,
122  * DIG, DK, D1, D2, ELIM, FID, FMR, FNU, FNUL, PI, RL, R1M5, SFAC,
123  * TOL, TTH, ZI, ZR, Z3I, Z3R, R1MACH
124  INTEGER ID, IERR, K, KODE, K1, K2, NZ, I1MACH
125  dimension cy(2)
126  DATA tth, c1, c2, coef, pi /6.66666666666666667e-01,
127  * 6.14926627446000736e-01,4.48288357353826359e-01,
128  * 5.77350269189625765e-01,3.14159265358979324e+00/
129  DATA cone / (1.0e0,0.0e0) /
130 C***FIRST EXECUTABLE STATEMENT CBIRY
131  ierr = 0
132  nz=0
133  IF (id.LT.0 .OR. id.GT.1) ierr=1
134  IF (kode.LT.1 .OR. kode.GT.2) ierr=1
135  IF (ierr.NE.0) RETURN
136  az = cabs(z)
137  tol = amax1(r1mach(4),1.0e-18)
138  fid = float(id)
139  IF (az.GT.1.0e0) GO TO 60
140 C-----------------------------------------------------------------------
141 C POWER SERIES FOR CABS(Z).LE.1.
142 C-----------------------------------------------------------------------
143  s1 = cone
144  s2 = cone
145  IF (az.LT.tol) GO TO 110
146  aa = az*az
147  IF (aa.LT.tol/az) GO TO 40
148  trm1 = cone
149  trm2 = cone
150  atrm = 1.0e0
151  z3 = z*z*z
152  az3 = az*aa
153  ak = 2.0e0 + fid
154  bk = 3.0e0 - fid - fid
155  ck = 4.0e0 - fid
156  dk = 3.0e0 + fid + fid
157  d1 = ak*dk
158  d2 = bk*ck
159  ad = amin1(d1,d2)
160  ak = 24.0e0 + 9.0e0*fid
161  bk = 30.0e0 - 9.0e0*fid
162  z3r = real(z3)
163  z3i = aimag(z3)
164  DO 30 k=1,25
165  trm1 = trm1*cmplx(z3r/d1,z3i/d1)
166  s1 = s1 + trm1
167  trm2 = trm2*cmplx(z3r/d2,z3i/d2)
168  s2 = s2 + trm2
169  atrm = atrm*az3/ad
170  d1 = d1 + ak
171  d2 = d2 + bk
172  ad = amin1(d1,d2)
173  IF (atrm.LT.tol*ad) GO TO 40
174  ak = ak + 18.0e0
175  bk = bk + 18.0e0
176  30 CONTINUE
177  40 CONTINUE
178  IF (id.EQ.1) GO TO 50
179  bi = s1*cmplx(c1,0.0e0) + z*s2*cmplx(c2,0.0e0)
180  IF (kode.EQ.1) RETURN
181  zta = z*csqrt(z)*cmplx(tth,0.0e0)
182  aa = real(zta)
183  aa = -abs(aa)
184  bi = bi*cmplx(exp(aa),0.0e0)
185  RETURN
186  50 CONTINUE
187  bi = s2*cmplx(c2,0.0e0)
188  IF (az.GT.tol) bi = bi + z*z*s1*cmplx(c1/(1.0e0+fid),0.0e0)
189  IF (kode.EQ.1) RETURN
190  zta = z*csqrt(z)*cmplx(tth,0.0e0)
191  aa = real(zta)
192  aa = -abs(aa)
193  bi = bi*cmplx(exp(aa),0.0e0)
194  RETURN
195 C-----------------------------------------------------------------------
196 C CASE FOR CABS(Z).GT.1.0
197 C-----------------------------------------------------------------------
198  60 CONTINUE
199  fnu = (1.0e0+fid)/3.0e0
200 C-----------------------------------------------------------------------
201 C SET PARAMETERS RELATED TO MACHINE CONSTANTS.
202 C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
203 C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
204 C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND
205 C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR
206 C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
207 C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
208 C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
209 C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU.
210 C-----------------------------------------------------------------------
211  k1 = i1mach(12)
212  k2 = i1mach(13)
213  r1m5 = r1mach(5)
214  k = min0(iabs(k1),iabs(k2))
215  elim = 2.303e0*(float(k)*r1m5-3.0e0)
216  k1 = i1mach(11) - 1
217  aa = r1m5*float(k1)
218  dig = amin1(aa,18.0e0)
219  aa = aa*2.303e0
220  alim = elim + amax1(-aa,-41.45e0)
221  rl = 1.2e0*dig + 3.0e0
222  fnul = 10.0e0 + 6.0e0*(dig-3.0e0)
223 C-----------------------------------------------------------------------
224 C TEST FOR RANGE
225 C-----------------------------------------------------------------------
226  aa=0.5e0/tol
227  bb=float(i1mach(9))*0.5e0
228  aa=amin1(aa,bb)
229  aa=aa**tth
230  IF (az.GT.aa) GO TO 190
231  aa=sqrt(aa)
232  IF (az.GT.aa) ierr=3
233  csq=csqrt(z)
234  zta=z*csq*cmplx(tth,0.0e0)
235 C-----------------------------------------------------------------------
236 C RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL
237 C-----------------------------------------------------------------------
238  sfac = 1.0e0
239  zi = aimag(z)
240  zr = real(z)
241  ak = aimag(zta)
242  IF (zr.GE.0.0e0) GO TO 70
243  bk = real(zta)
244  ck = -abs(bk)
245  zta = cmplx(ck,ak)
246  70 CONTINUE
247  IF (zi.EQ.0.0e0 .AND. zr.LE.0.0e0) zta = cmplx(0.0e0,ak)
248  aa = real(zta)
249  IF (kode.EQ.2) GO TO 80
250 C-----------------------------------------------------------------------
251 C OVERFLOW TEST
252 C-----------------------------------------------------------------------
253  bb = abs(aa)
254  IF (bb.LT.alim) GO TO 80
255  bb = bb + 0.25e0*alog(az)
256  sfac = tol
257  IF (bb.GT.elim) GO TO 170
258  80 CONTINUE
259  fmr = 0.0e0
260  IF (aa.GE.0.0e0 .AND. zr.GT.0.0e0) GO TO 90
261  fmr = pi
262  IF (zi.LT.0.0e0) fmr = -pi
263  zta = -zta
264  90 CONTINUE
265 C-----------------------------------------------------------------------
266 C AA=FACTOR FOR ANALYTIC CONTINUATION OF I(FNU,ZTA)
267 C KODE=2 RETURNS EXP(-ABS(XZTA))*I(FNU,ZTA) FROM CBINU
268 C-----------------------------------------------------------------------
269  CALL cbinu(zta, fnu, kode, 1, cy, nz, rl, fnul, tol, elim, alim)
270  IF (nz.LT.0) GO TO 180
271  aa = fmr*fnu
272  z3 = cmplx(sfac,0.0e0)
273  s1 = cy(1)*cmplx(cos(aa),sin(aa))*z3
274  fnu = (2.0e0-fid)/3.0e0
275  CALL cbinu(zta, fnu, kode, 2, cy, nz, rl, fnul, tol, elim, alim)
276  cy(1) = cy(1)*z3
277  cy(2) = cy(2)*z3
278 C-----------------------------------------------------------------------
279 C BACKWARD RECUR ONE STEP FOR ORDERS -1/3 OR -2/3
280 C-----------------------------------------------------------------------
281  s2 = cy(1)*cmplx(fnu+fnu,0.0e0)/zta + cy(2)
282  aa = fmr*(fnu-1.0e0)
283  s1 = (s1+s2*cmplx(cos(aa),sin(aa)))*cmplx(coef,0.0e0)
284  IF (id.EQ.1) GO TO 100
285  s1 = csq*s1
286  bi = s1*cmplx(1.0e0/sfac,0.0e0)
287  RETURN
288  100 CONTINUE
289  s1 = z*s1
290  bi = s1*cmplx(1.0e0/sfac,0.0e0)
291  RETURN
292  110 CONTINUE
293  aa = c1*(1.0e0-fid) + fid*c2
294  bi = cmplx(aa,0.0e0)
295  RETURN
296  170 CONTINUE
297  nz=0
298  ierr=2
299  RETURN
300  180 CONTINUE
301  IF(nz.EQ.(-1)) GO TO 170
302  nz=0
303  ierr=5
304  RETURN
305  190 CONTINUE
306  ierr=4
307  nz=0
308  RETURN
309  END
cmplx
double complex cmplx
Definition: Faddeeva.cc:217
cbinu
subroutine cbinu(Z, FNU, KODE, N, CY, NZ, RL, FNUL, TOL, ELIM, ALIM)
Definition: cbinu.f:3
cbiry
subroutine cbiry(Z, ID, KODE, BI, IERR)
Definition: cbiry.f:2
real
ColumnVector real(const ComplexColumnVector &a)
Definition: dColVector.cc:137
abs
static T abs(T x)
Definition: pr-output.cc:1678