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
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
Generated on Mon Dec 30 2013 03:04:45 for GNU Octave by   doxygen 1.8.1.2