GNU Octave  9.1.0
A high-level interpreted language, primarily intended for numerical computations, mostly compatible with Matlab
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
subroutine xzexp(AR, AI, BR, BI)
Definition: xzexp.f:2
subroutine xzsqrt(AR, AI, BR, BI)
Definition: xzsqrt.f:2
subroutine zacai(ZR, ZI, FNU, KODE, MR, N, YR, YI, NZ, RL, TOL, ELIM, ALIM)
Definition: zacai.f:3
subroutine zairy(ZR, ZI, ID, KODE, AIR, AII, NZ, IERR)
Definition: zairy.f:2
subroutine zbknu(ZR, ZI, FNU, KODE, N, YR, YI, NZ, TOL, ELIM, ALIM)
Definition: zbknu.f:3