GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
zbiry.f
Go to the documentation of this file.
1  SUBROUTINE zbiry(ZR, ZI, ID, KODE, BIR, BII, IERR)
2 C***BEGIN PROLOGUE ZBIRY
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 ***A DOUBLE PRECISION ROUTINE***
12 C ON KODE=1, CBIRY COMPUTES THE COMPLEX AIRY FUNCTION BI(Z) OR
13 C ITS DERIVATIVE DBI(Z)/DZ ON ID=0 OR ID=1 RESPECTIVELY. ON
14 C KODE=2, A SCALING OPTION CEXP(-AXZTA)*BI(Z) OR CEXP(-AXZTA)*
15 C DBI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL BEHAVIOR IN
16 C BOTH THE LEFT AND RIGHT HALF PLANES WHERE
17 C ZTA=(2/3)*Z*CSQRT(Z)=CMPLX(XZTA,YZTA) AND AXZTA=ABS(XZTA).
18 C DEFINTIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF
19 C MATHEMATICAL FUNCTIONS (REF. 1).
20 C
21 C INPUT ZR,ZI ARE DOUBLE PRECISION
22 C ZR,ZI - Z=CMPLX(ZR,ZI)
23 C ID - ORDER OF DERIVATIVE, ID=0 OR ID=1
24 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION
25 C KODE= 1 RETURNS
26 C BI=BI(Z) ON ID=0 OR
27 C BI=DBI(Z)/DZ ON ID=1
28 C = 2 RETURNS
29 C BI=CEXP(-AXZTA)*BI(Z) ON ID=0 OR
30 C BI=CEXP(-AXZTA)*DBI(Z)/DZ ON ID=1 WHERE
31 C ZTA=(2/3)*Z*CSQRT(Z)=CMPLX(XZTA,YZTA)
32 C AND AXZTA=ABS(XZTA)
33 C
34 C OUTPUT BIR,BII ARE DOUBLE PRECISION
35 C BIR,BII- COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND
36 C KODE
37 C IERR - ERROR FLAG
38 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
39 C IERR=1, INPUT ERROR - NO COMPUTATION
40 C IERR=2, OVERFLOW - NO COMPUTATION, REAL(Z)
41 C TOO LARGE ON KODE=1
42 C IERR=3, CABS(Z) LARGE - COMPUTATION COMPLETED
43 C LOSSES OF SIGNIFCANCE BY ARGUMENT REDUCTION
44 C PRODUCE LESS THAN HALF OF MACHINE ACCURACY
45 C IERR=4, CABS(Z) TOO LARGE - NO COMPUTATION
46 C COMPLETE LOSS OF ACCURACY BY ARGUMENT
47 C REDUCTION
48 C IERR=5, ERROR - NO COMPUTATION,
49 C ALGORITHM TERMINATION CONDITION NOT MET
50 C
51 C***LONG DESCRIPTION
52 C
53 C BI AND DBI ARE COMPUTED FOR CABS(Z).GT.1.0 FROM THE I BESSEL
54 C FUNCTIONS BY
55 C
56 C BI(Z)=C*SQRT(Z)*( I(-1/3,ZTA) + I(1/3,ZTA) )
57 C DBI(Z)=C * Z * ( I(-2/3,ZTA) + I(2/3,ZTA) )
58 C C=1.0/SQRT(3.0)
59 C ZTA=(2/3)*Z**(3/2)
60 C
61 C WITH THE POWER SERIES FOR CABS(Z).LE.1.0.
62 C
63 C IN MOST COMPLEX VARIABLE COMPUTATION, ONE MUST EVALUATE ELE-
64 C MENTARY FUNCTIONS. WHEN THE MAGNITUDE OF Z IS LARGE, LOSSES
65 C OF SIGNIFICANCE BY ARGUMENT REDUCTION OCCUR. CONSEQUENTLY, IF
66 C THE MAGNITUDE OF ZETA=(2/3)*Z**1.5 EXCEEDS U1=SQRT(0.5/UR),
67 C THEN LOSSES EXCEEDING HALF PRECISION ARE LIKELY AND AN ERROR
68 C FLAG IERR=3 IS TRIGGERED WHERE UR=DMAX1(D1MACH(4),1.0D-18) IS
69 C DOUBLE PRECISION UNIT ROUNDOFF LIMITED TO 18 DIGITS PRECISION.
70 C ALSO, IF THE MAGNITUDE OF ZETA IS LARGER THAN U2=0.5/UR, THEN
71 C ALL SIGNIFICANCE IS LOST AND IERR=4. IN ORDER TO USE THE INT
72 C FUNCTION, ZETA MUST BE FURTHER RESTRICTED NOT TO EXCEED THE
73 C LARGEST INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF ZETA
74 C MUST BE RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2,
75 C AND U3 ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE
76 C PRECISION ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE
77 C PRECISION ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMIT-
78 C ING IN THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT THE MAG-
79 C NITUDE OF Z CANNOT EXCEED 3.1E+4 IN SINGLE AND 2.1E+6 IN
80 C DOUBLE PRECISION ARITHMETIC. THIS ALSO MEANS THAT ONE CAN
81 C EXPECT TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES,
82 C NO DIGITS IN SINGLE PRECISION AND ONLY 7 DIGITS IN DOUBLE
83 C PRECISION ARITHMETIC. SIMILAR CONSIDERATIONS HOLD FOR OTHER
84 C MACHINES.
85 C
86 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
87 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
88 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
89 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
90 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
91 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
92 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
93 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
94 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
95 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
96 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
97 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
98 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
99 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
100 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
101 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
102 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
103 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
104 C OR -PI/2+P.
105 C
106 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
107 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
108 C COMMERCE, 1955.
109 C
110 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
111 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
112 C
113 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
114 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
115 C 1018, MAY, 1985
116 C
117 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
118 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
119 C MATH. SOFTWARE, 1986
120 C
121 C***ROUTINES CALLED ZBINU,XZABS,ZDIV,XZSQRT,D1MACH,I1MACH
122 C***END PROLOGUE ZBIRY
123 C COMPLEX BI,CONE,CSQ,CY,S1,S2,TRM1,TRM2,Z,ZTA,Z3
124  DOUBLE PRECISION AA, AD, AK, ALIM, ATRM, AZ, AZ3, BB, BII, BIR,
125  * BK, CC, CK, COEF, CONEI, CONER, CSQI, CSQR, CYI, CYR, C1, C2,
126  * DIG, DK, D1, D2, EAA, ELIM, FID, FMR, FNU, FNUL, PI, RL, R1M5,
127  * SFAC, STI, STR, S1I, S1R, S2I, S2R, TOL, TRM1I, TRM1R, TRM2I,
128  * TRM2R, TTH, ZI, ZR, ZTAI, ZTAR, Z3I, Z3R, D1MACH, XZABS
129  INTEGER ID, IERR, K, KODE, K1, K2, NZ, I1MACH
130  dimension cyr(2), cyi(2)
131  DATA tth, c1, c2, coef, pi /6.66666666666666667d-01,
132  * 6.14926627446000736d-01,4.48288357353826359d-01,
133  * 5.77350269189625765d-01,3.14159265358979324d+00/
134  DATA coner, conei /1.0d0,0.0d0/
135 C***FIRST EXECUTABLE STATEMENT ZBIRY
136  ierr = 0
137  nz=0
138  IF (id.LT.0 .OR. id.GT.1) ierr=1
139  IF (kode.LT.1 .OR. kode.GT.2) ierr=1
140  IF (ierr.NE.0) RETURN
141  az = xzabs(zr,zi)
142  tol = dmax1(d1mach(4),1.0d-18)
143  fid = dble(float(id))
144  IF (az.GT.1.0e0) GO TO 70
145 C-----------------------------------------------------------------------
146 C POWER SERIES FOR CABS(Z).LE.1.
147 C-----------------------------------------------------------------------
148  s1r = coner
149  s1i = conei
150  s2r = coner
151  s2i = conei
152  IF (az.LT.tol) GO TO 130
153  aa = az*az
154  IF (aa.LT.tol/az) GO TO 40
155  trm1r = coner
156  trm1i = conei
157  trm2r = coner
158  trm2i = conei
159  atrm = 1.0d0
160  str = zr*zr - zi*zi
161  sti = zr*zi + zi*zr
162  z3r = str*zr - sti*zi
163  z3i = str*zi + sti*zr
164  az3 = az*aa
165  ak = 2.0d0 + fid
166  bk = 3.0d0 - fid - fid
167  ck = 4.0d0 - fid
168  dk = 3.0d0 + fid + fid
169  d1 = ak*dk
170  d2 = bk*ck
171  ad = dmin1(d1,d2)
172  ak = 24.0d0 + 9.0d0*fid
173  bk = 30.0d0 - 9.0d0*fid
174  DO 30 k=1,25
175  str = (trm1r*z3r-trm1i*z3i)/d1
176  trm1i = (trm1r*z3i+trm1i*z3r)/d1
177  trm1r = str
178  s1r = s1r + trm1r
179  s1i = s1i + trm1i
180  str = (trm2r*z3r-trm2i*z3i)/d2
181  trm2i = (trm2r*z3i+trm2i*z3r)/d2
182  trm2r = str
183  s2r = s2r + trm2r
184  s2i = s2i + trm2i
185  atrm = atrm*az3/ad
186  d1 = d1 + ak
187  d2 = d2 + bk
188  ad = dmin1(d1,d2)
189  IF (atrm.LT.tol*ad) GO TO 40
190  ak = ak + 18.0d0
191  bk = bk + 18.0d0
192  30 CONTINUE
193  40 CONTINUE
194  IF (id.EQ.1) GO TO 50
195  bir = c1*s1r + c2*(zr*s2r-zi*s2i)
196  bii = c1*s1i + c2*(zr*s2i+zi*s2r)
197  IF (kode.EQ.1) RETURN
198  CALL xzsqrt(zr, zi, str, sti)
199  ztar = tth*(zr*str-zi*sti)
200  ztai = tth*(zr*sti+zi*str)
201  aa = ztar
202  aa = -dabs(aa)
203  eaa = dexp(aa)
204  bir = bir*eaa
205  bii = bii*eaa
206  RETURN
207  50 CONTINUE
208  bir = s2r*c2
209  bii = s2i*c2
210  IF (az.LE.tol) GO TO 60
211  cc = c1/(1.0d0+fid)
212  str = s1r*zr - s1i*zi
213  sti = s1r*zi + s1i*zr
214  bir = bir + cc*(str*zr-sti*zi)
215  bii = bii + cc*(str*zi+sti*zr)
216  60 CONTINUE
217  IF (kode.EQ.1) RETURN
218  CALL xzsqrt(zr, zi, str, sti)
219  ztar = tth*(zr*str-zi*sti)
220  ztai = tth*(zr*sti+zi*str)
221  aa = ztar
222  aa = -dabs(aa)
223  eaa = dexp(aa)
224  bir = bir*eaa
225  bii = bii*eaa
226  RETURN
227 C-----------------------------------------------------------------------
228 C CASE FOR CABS(Z).GT.1.0
229 C-----------------------------------------------------------------------
230  70 CONTINUE
231  fnu = (1.0d0+fid)/3.0d0
232 C-----------------------------------------------------------------------
233 C SET PARAMETERS RELATED TO MACHINE CONSTANTS.
234 C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
235 C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
236 C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND
237 C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR
238 C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
239 C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
240 C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
241 C FNUL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC SERIES FOR LARGE FNU.
242 C-----------------------------------------------------------------------
243  k1 = i1mach(15)
244  k2 = i1mach(16)
245  r1m5 = d1mach(5)
246  k = min0(iabs(k1),iabs(k2))
247  elim = 2.303d0*(dble(float(k))*r1m5-3.0d0)
248  k1 = i1mach(14) - 1
249  aa = r1m5*dble(float(k1))
250  dig = dmin1(aa,18.0d0)
251  aa = aa*2.303d0
252  alim = elim + dmax1(-aa,-41.45d0)
253  rl = 1.2d0*dig + 3.0d0
254  fnul = 10.0d0 + 6.0d0*(dig-3.0d0)
255 C-----------------------------------------------------------------------
256 C TEST FOR RANGE
257 C-----------------------------------------------------------------------
258  aa=0.5d0/tol
259  bb=dble(float(i1mach(9)))*0.5d0
260  aa=dmin1(aa,bb)
261  aa=aa**tth
262  IF (az.GT.aa) GO TO 260
263  aa=dsqrt(aa)
264  IF (az.GT.aa) ierr=3
265  CALL xzsqrt(zr, zi, csqr, csqi)
266  ztar = tth*(zr*csqr-zi*csqi)
267  ztai = tth*(zr*csqi+zi*csqr)
268 C-----------------------------------------------------------------------
269 C RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL
270 C-----------------------------------------------------------------------
271  sfac = 1.0d0
272  ak = ztai
273  IF (zr.GE.0.0d0) GO TO 80
274  bk = ztar
275  ck = -dabs(bk)
276  ztar = ck
277  ztai = ak
278  80 CONTINUE
279  IF (zi.NE.0.0d0 .OR. zr.GT.0.0d0) GO TO 90
280  ztar = 0.0d0
281  ztai = ak
282  90 CONTINUE
283  aa = ztar
284  IF (kode.EQ.2) GO TO 100
285 C-----------------------------------------------------------------------
286 C OVERFLOW TEST
287 C-----------------------------------------------------------------------
288  bb = dabs(aa)
289  IF (bb.LT.alim) GO TO 100
290  bb = bb + 0.25d0*dlog(az)
291  sfac = tol
292  IF (bb.GT.elim) GO TO 190
293  100 CONTINUE
294  fmr = 0.0d0
295  IF (aa.GE.0.0d0 .AND. zr.GT.0.0d0) GO TO 110
296  fmr = pi
297  IF (zi.LT.0.0d0) fmr = -pi
298  ztar = -ztar
299  ztai = -ztai
300  110 CONTINUE
301 C-----------------------------------------------------------------------
302 C AA=FACTOR FOR ANALYTIC CONTINUATION OF I(FNU,ZTA)
303 C KODE=2 RETURNS EXP(-ABS(XZTA))*I(FNU,ZTA) FROM CBESI
304 C-----------------------------------------------------------------------
305  CALL zbinu(ztar, ztai, fnu, kode, 1, cyr, cyi, nz, rl, fnul, tol,
306  * elim, alim)
307  IF (nz.LT.0) GO TO 200
308  aa = fmr*fnu
309  z3r = sfac
310  str = dcos(aa)
311  sti = dsin(aa)
312  s1r = (str*cyr(1)-sti*cyi(1))*z3r
313  s1i = (str*cyi(1)+sti*cyr(1))*z3r
314  fnu = (2.0d0-fid)/3.0d0
315  CALL zbinu(ztar, ztai, fnu, kode, 2, cyr, cyi, nz, rl, fnul, tol,
316  * elim, alim)
317  cyr(1) = cyr(1)*z3r
318  cyi(1) = cyi(1)*z3r
319  cyr(2) = cyr(2)*z3r
320  cyi(2) = cyi(2)*z3r
321 C-----------------------------------------------------------------------
322 C BACKWARD RECUR ONE STEP FOR ORDERS -1/3 OR -2/3
323 C-----------------------------------------------------------------------
324  CALL zdiv(cyr(1), cyi(1), ztar, ztai, str, sti)
325  s2r = (fnu+fnu)*str + cyr(2)
326  s2i = (fnu+fnu)*sti + cyi(2)
327  aa = fmr*(fnu-1.0d0)
328  str = dcos(aa)
329  sti = dsin(aa)
330  s1r = coef*(s1r+s2r*str-s2i*sti)
331  s1i = coef*(s1i+s2r*sti+s2i*str)
332  IF (id.EQ.1) GO TO 120
333  str = csqr*s1r - csqi*s1i
334  s1i = csqr*s1i + csqi*s1r
335  s1r = str
336  bir = s1r/sfac
337  bii = s1i/sfac
338  RETURN
339  120 CONTINUE
340  str = zr*s1r - zi*s1i
341  s1i = zr*s1i + zi*s1r
342  s1r = str
343  bir = s1r/sfac
344  bii = s1i/sfac
345  RETURN
346  130 CONTINUE
347  aa = c1*(1.0d0-fid) + fid*c2
348  bir = aa
349  bii = 0.0d0
350  RETURN
351  190 CONTINUE
352  ierr=2
353  nz=0
354  RETURN
355  200 CONTINUE
356  IF(nz.EQ.(-1)) GO TO 190
357  nz=0
358  ierr=5
359  RETURN
360  260 CONTINUE
361  ierr=4
362  nz=0
363  RETURN
364  END
xzsqrt
subroutine xzsqrt(AR, AI, BR, BI)
Definition: xzsqrt.f:2
zbinu
subroutine zbinu(ZR, ZI, FNU, KODE, N, CYR, CYI, NZ, RL, FNUL, TOL, ELIM, ALIM)
Definition: zbinu.f:3
zbiry
subroutine zbiry(ZR, ZI, ID, KODE, BIR, BII, IERR)
Definition: zbiry.f:2
zdiv
subroutine zdiv(AR, AI, BR, BI, CR, CI)
Definition: zdiv.f:2