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