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
cairy.f
Go to the documentation of this file.
1  SUBROUTINE cairy(Z, ID, KODE, AI, NZ, IERR)
2 C***BEGIN PROLOGUE CAIRY
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 ON KODE=1, CAIRY COMPUTES THE COMPLEX AIRY FUNCTION AI(Z) OR
12 C ITS DERIVATIVE DAI(Z)/DZ ON ID=0 OR ID=1 RESPECTIVELY. ON
13 C KODE=2, A SCALING OPTION CEXP(ZTA)*AI(Z) OR CEXP(ZTA)*
14 C DAI(Z)/DZ IS PROVIDED TO REMOVE THE EXPONENTIAL DECAY IN
15 C -PI/3.LT.ARG(Z).LT.PI/3 AND THE EXPONENTIAL GROWTH IN
16 C PI/3.LT.ABS(ARG(Z)).LT.PI WHERE ZTA=(2/3)*Z*CSQRT(Z)
17 C
18 C WHILE THE AIRY FUNCTIONS AI(Z) AND DAI(Z)/DZ ARE ANALYTIC IN
19 C THE WHOLE Z PLANE, THE CORRESPONDING SCALED FUNCTIONS DEFINED
20 C FOR KODE=2 HAVE A CUT ALONG THE NEGATIVE REAL AXIS.
21 C DEFINITIONS AND NOTATION ARE FOUND IN THE NBS HANDBOOK OF
22 C MATHEMATICAL FUNCTIONS (REF. 1).
23 C
24 C INPUT
25 C Z - Z=CMPLX(X,Y)
26 C ID - ORDER OF DERIVATIVE, ID=0 OR ID=1
27 C KODE - A PARAMETER TO INDICATE THE SCALING OPTION
28 C KODE= 1 RETURNS
29 C AI=AI(Z) ON ID=0 OR
30 C AI=DAI(Z)/DZ ON ID=1
31 C = 2 RETURNS
32 C AI=CEXP(ZTA)*AI(Z) ON ID=0 OR
33 C AI=CEXP(ZTA)*DAI(Z)/DZ ON ID=1 WHERE
34 C ZTA=(2/3)*Z*CSQRT(Z)
35 C
36 C OUTPUT
37 C AI - COMPLEX ANSWER DEPENDING ON THE CHOICES FOR ID AND
38 C KODE
39 C NZ - UNDERFLOW INDICATOR
40 C NZ= 0 , NORMAL RETURN
41 C NZ= 1 , AI=CMPLX(0.0,0.0) DUE TO UNDERFLOW IN
42 C -PI/3.LT.ARG(Z).LT.PI/3 ON KODE=1
43 C IERR - ERROR FLAG
44 C IERR=0, NORMAL RETURN - COMPUTATION COMPLETED
45 C IERR=1, INPUT ERROR - NO COMPUTATION
46 C IERR=2, OVERFLOW - NO COMPUTATION, REAL(ZTA)
47 C TOO LARGE WITH KODE=1.
48 C IERR=3, CABS(Z) LARGE - COMPUTATION COMPLETED
49 C LOSSES OF SIGNIFCANCE BY ARGUMENT REDUCTION
50 C PRODUCE LESS THAN HALF OF MACHINE ACCURACY
51 C IERR=4, CABS(Z) TOO LARGE - NO COMPUTATION
52 C COMPLETE LOSS OF ACCURACY BY ARGUMENT
53 C REDUCTION
54 C IERR=5, ERROR - NO COMPUTATION,
55 C ALGORITHM TERMINATION CONDITION NOT MET
56 C
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=R1MACH(4)=UNIT ROUNDOFF.
75 C ALSO, IF THE MAGNITUDE OF ZETA IS LARGER THAN U2=0.5/UR, THEN
76 C ALL SIGNIFICANCE IS LOST AND IERR=4. IN ORDER TO USE THE INT
77 C FUNCTION, ZETA MUST BE FURTHER RESTRICTED NOT TO EXCEED THE
78 C LARGEST INTEGER, U3=I1MACH(9). THUS, THE MAGNITUDE OF ZETA
79 C MUST BE RESTRICTED BY MIN(U2,U3). ON 32 BIT MACHINES, U1,U2,
80 C AND U3 ARE APPROXIMATELY 2.0E+3, 4.2E+6, 2.1E+9 IN SINGLE
81 C PRECISION ARITHMETIC AND 1.3E+8, 1.8E+16, 2.1E+9 IN DOUBLE
82 C PRECISION ARITHMETIC RESPECTIVELY. THIS MAKES U2 AND U3 LIMIT-
83 C ING IN THEIR RESPECTIVE ARITHMETICS. THIS MEANS THAT THE MAG-
84 C NITUDE OF Z CANNOT EXCEED 3.1E+4 IN SINGLE AND 2.1E+6 IN
85 C DOUBLE PRECISION ARITHMETIC. THIS ALSO MEANS THAT ONE CAN
86 C EXPECT TO RETAIN, IN THE WORST CASES ON 32 BIT MACHINES,
87 C NO DIGITS IN SINGLE PRECISION AND ONLY 7 DIGITS IN DOUBLE
88 C PRECISION ARITHMETIC. SIMILAR CONSIDERATIONS HOLD FOR OTHER
89 C MACHINES.
90 C
91 C THE APPROXIMATE RELATIVE ERROR IN THE MAGNITUDE OF A COMPLEX
92 C BESSEL FUNCTION CAN BE EXPRESSED BY P*10**S WHERE P=MAX(UNIT
93 C ROUNDOFF,1.0E-18) IS THE NOMINAL PRECISION AND 10**S REPRE-
94 C SENTS THE INCREASE IN ERROR DUE TO ARGUMENT REDUCTION IN THE
95 C ELEMENTARY FUNCTIONS. HERE, S=MAX(1,ABS(LOG10(CABS(Z))),
96 C ABS(LOG10(FNU))) APPROXIMATELY (I.E. S=MAX(1,ABS(EXPONENT OF
97 C CABS(Z),ABS(EXPONENT OF FNU)) ). HOWEVER, THE PHASE ANGLE MAY
98 C HAVE ONLY ABSOLUTE ACCURACY. THIS IS MOST LIKELY TO OCCUR WHEN
99 C ONE COMPONENT (IN ABSOLUTE VALUE) IS LARGER THAN THE OTHER BY
100 C SEVERAL ORDERS OF MAGNITUDE. IF ONE COMPONENT IS 10**K LARGER
101 C THAN THE OTHER, THEN ONE CAN EXPECT ONLY MAX(ABS(LOG10(P))-K,
102 C 0) SIGNIFICANT DIGITS; OR, STATED ANOTHER WAY, WHEN K EXCEEDS
103 C THE EXPONENT OF P, NO SIGNIFICANT DIGITS REMAIN IN THE SMALLER
104 C COMPONENT. HOWEVER, THE PHASE ANGLE RETAINS ABSOLUTE ACCURACY
105 C BECAUSE, IN COMPLEX ARITHMETIC WITH PRECISION P, THE SMALLER
106 C COMPONENT WILL NOT (AS A RULE) DECREASE BELOW P TIMES THE
107 C MAGNITUDE OF THE LARGER COMPONENT. IN THESE EXTREME CASES,
108 C THE PRINCIPAL PHASE ANGLE IS ON THE ORDER OF +P, -P, PI/2-P,
109 C OR -PI/2+P.
110 C
111 C***REFERENCES HANDBOOK OF MATHEMATICAL FUNCTIONS BY M. ABRAMOWITZ
112 C AND I. A. STEGUN, NBS AMS SERIES 55, U.S. DEPT. OF
113 C COMMERCE, 1955.
114 C
115 C COMPUTATION OF BESSEL FUNCTIONS OF COMPLEX ARGUMENT
116 C AND LARGE ORDER BY D. E. AMOS, SAND83-0643, MAY, 1983
117 C
118 C A SUBROUTINE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
119 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, SAND85-
120 C 1018, MAY, 1985
121 C
122 C A PORTABLE PACKAGE FOR BESSEL FUNCTIONS OF A COMPLEX
123 C ARGUMENT AND NONNEGATIVE ORDER BY D. E. AMOS, TRANS.
124 C MATH. SOFTWARE, 1986
125 C
126 C***ROUTINES CALLED CACAI,CBKNU,I1MACH,R1MACH
127 C***END PROLOGUE CAIRY
128  COMPLEX ai, cone, csq, cy, s1, s2, trm1, trm2, z, zta, z3
129  REAL aa, ad, ak, alim, atrm, az, az3, bk, ck, coef, c1, c2, dig,
130  * dk, d1, d2, elim, fid, fnu, rl, r1m5, sfac, tol, tth, zi, zr,
131  * z3i, z3r, r1mach, bb, alaz
132  INTEGER id, ierr, iflag, k, kode, k1, k2, mr, nn, nz, i1mach
133  dimension cy(1)
134  DATA tth, c1, c2, coef /6.66666666666666667e-01,
135  * 3.55028053887817240e-01,2.58819403792806799e-01,
136  * 1.83776298473930683e-01/
137  DATA cone / (1.0e0,0.0e0) /
138 C***FIRST EXECUTABLE STATEMENT CAIRY
139  ierr = 0
140  nz=0
141  IF (id.LT.0 .OR. id.GT.1) ierr=1
142  IF (kode.LT.1 .OR. kode.GT.2) ierr=1
143  IF (ierr.NE.0) RETURN
144  az = cabs(z)
145  tol = amax1(r1mach(4),1.0e-18)
146  fid = float(id)
147  IF (az.GT.1.0e0) go to 60
148 C-----------------------------------------------------------------------
149 C POWER SERIES FOR CABS(Z).LE.1.
150 C-----------------------------------------------------------------------
151  s1 = cone
152  s2 = cone
153  IF (az.LT.tol) go to 160
154  aa = az*az
155  IF (aa.LT.tol/az) go to 40
156  trm1 = cone
157  trm2 = cone
158  atrm = 1.0e0
159  z3 = z*z*z
160  az3 = az*aa
161  ak = 2.0e0 + fid
162  bk = 3.0e0 - fid - fid
163  ck = 4.0e0 - fid
164  dk = 3.0e0 + fid + fid
165  d1 = ak*dk
166  d2 = bk*ck
167  ad = amin1(d1,d2)
168  ak = 24.0e0 + 9.0e0*fid
169  bk = 30.0e0 - 9.0e0*fid
170  z3r = REAL(z3)
171  z3i = aimag(z3)
172  DO 30 k=1,25
173  trm1 = trm1*cmplx(z3r/d1,z3i/d1)
174  s1 = s1 + trm1
175  trm2 = trm2*cmplx(z3r/d2,z3i/d2)
176  s2 = s2 + trm2
177  atrm = atrm*az3/ad
178  d1 = d1 + ak
179  d2 = d2 + bk
180  ad = amin1(d1,d2)
181  IF (atrm.LT.tol*ad) go to 40
182  ak = ak + 18.0e0
183  bk = bk + 18.0e0
184  30 CONTINUE
185  40 CONTINUE
186  IF (id.EQ.1) go to 50
187  ai = s1*cmplx(c1,0.0e0) - z*s2*cmplx(c2,0.0e0)
188  IF (kode.EQ.1) RETURN
189  zta = z*csqrt(z)*cmplx(tth,0.0e0)
190  ai = ai*cexp(zta)
191  RETURN
192  50 CONTINUE
193  ai = -s2*cmplx(c2,0.0e0)
194  IF (az.GT.tol) ai = ai + z*z*s1*cmplx(c1/(1.0e0+fid),0.0e0)
195  IF (kode.EQ.1) RETURN
196  zta = z*csqrt(z)*cmplx(tth,0.0e0)
197  ai = ai*cexp(zta)
198  RETURN
199 C-----------------------------------------------------------------------
200 C CASE FOR CABS(Z).GT.1.0
201 C-----------------------------------------------------------------------
202  60 CONTINUE
203  fnu = (1.0e0+fid)/3.0e0
204 C-----------------------------------------------------------------------
205 C SET PARAMETERS RELATED TO MACHINE CONSTANTS.
206 C TOL IS THE APPROXIMATE UNIT ROUNDOFF LIMITED TO 1.0E-18.
207 C ELIM IS THE APPROXIMATE EXPONENTIAL OVER- AND UNDERFLOW LIMIT.
208 C EXP(-ELIM).LT.EXP(-ALIM)=EXP(-ELIM)/TOL AND
209 C EXP(ELIM).GT.EXP(ALIM)=EXP(ELIM)*TOL ARE INTERVALS NEAR
210 C UNDERFLOW AND OVERFLOW LIMITS WHERE SCALED ARITHMETIC IS DONE.
211 C RL IS THE LOWER BOUNDARY OF THE ASYMPTOTIC EXPANSION FOR LARGE Z.
212 C DIG = NUMBER OF BASE 10 DIGITS IN TOL = 10**(-DIG).
213 C-----------------------------------------------------------------------
214  k1 = i1mach(12)
215  k2 = i1mach(13)
216  r1m5 = r1mach(5)
217  k = min0(iabs(k1),iabs(k2))
218  elim = 2.303e0*(float(k)*r1m5-3.0e0)
219  k1 = i1mach(11) - 1
220  aa = r1m5*float(k1)
221  dig = amin1(aa,18.0e0)
222  aa = aa*2.303e0
223  alim = elim + amax1(-aa,-41.45e0)
224  rl = 1.2e0*dig + 3.0e0
225  alaz=alog(az)
226 C-----------------------------------------------------------------------
227 C TEST FOR RANGE
228 C-----------------------------------------------------------------------
229  aa=0.5e0/tol
230  bb=float(i1mach(9))*0.5e0
231  aa=amin1(aa,bb)
232  aa=aa**tth
233  IF (az.GT.aa) go to 260
234  aa=sqrt(aa)
235  IF (az.GT.aa) ierr=3
236  csq=csqrt(z)
237  zta=z*csq*cmplx(tth,0.0e0)
238 C-----------------------------------------------------------------------
239 C RE(ZTA).LE.0 WHEN RE(Z).LT.0, ESPECIALLY WHEN IM(Z) IS SMALL
240 C-----------------------------------------------------------------------
241  iflag = 0
242  sfac = 1.0e0
243  zi = aimag(z)
244  zr = REAL(z)
245  ak = aimag(zta)
246  IF (zr.GE.0.0e0) go to 70
247  bk = REAL(zta)
248  ck = -abs(bk)
249  zta = cmplx(ck,ak)
250  70 CONTINUE
251  IF (zi.NE.0.0e0) go to 80
252  IF (zr.GT.0.0e0) go to 80
253  zta = cmplx(0.0e0,ak)
254  80 CONTINUE
255  aa = REAL(zta)
256  IF (aa.GE.0.0e0 .AND. zr.GT.0.0e0) go to 100
257  IF (kode.EQ.2) go to 90
258 C-----------------------------------------------------------------------
259 C OVERFLOW TEST
260 C-----------------------------------------------------------------------
261  IF (aa.GT.(-alim)) go to 90
262  aa = -aa + 0.25e0*alaz
263  iflag = 1
264  sfac = tol
265  IF (aa.GT.elim) go to 240
266  90 CONTINUE
267 C-----------------------------------------------------------------------
268 C CBKNU AND CACAI RETURN EXP(ZTA)*K(FNU,ZTA) ON KODE=2
269 C-----------------------------------------------------------------------
270  mr = 1
271  IF (zi.LT.0.0e0) mr = -1
272  CALL cacai(zta, fnu, kode, mr, 1, cy, nn, rl, tol, elim, alim)
273  IF (nn.LT.0) go to 250
274  nz = nz + nn
275  go to 120
276  100 CONTINUE
277  IF (kode.EQ.2) go to 110
278 C-----------------------------------------------------------------------
279 C UNDERFLOW TEST
280 C-----------------------------------------------------------------------
281  IF (aa.LT.alim) go to 110
282  aa = -aa - 0.25e0*alaz
283  iflag = 2
284  sfac = 1.0e0/tol
285  IF (aa.LT.(-elim)) go to 180
286  110 CONTINUE
287  CALL cbknu(zta, fnu, kode, 1, cy, nz, tol, elim, alim)
288  120 CONTINUE
289  s1 = cy(1)*cmplx(coef,0.0e0)
290  IF (iflag.NE.0) go to 140
291  IF (id.EQ.1) go to 130
292  ai = csq*s1
293  RETURN
294  130 ai = -z*s1
295  RETURN
296  140 CONTINUE
297  s1 = s1*cmplx(sfac,0.0e0)
298  IF (id.EQ.1) go to 150
299  s1 = s1*csq
300  ai = s1*cmplx(1.0e0/sfac,0.0e0)
301  RETURN
302  150 CONTINUE
303  s1 = -s1*z
304  ai = s1*cmplx(1.0e0/sfac,0.0e0)
305  RETURN
306  160 CONTINUE
307  aa = 1.0e+3*r1mach(1)
308  s1 = cmplx(0.0e0,0.0e0)
309  IF (id.EQ.1) go to 170
310  IF (az.GT.aa) s1 = cmplx(c2,0.0e0)*z
311  ai = cmplx(c1,0.0e0) - s1
312  RETURN
313  170 CONTINUE
314  ai = -cmplx(c2,0.0e0)
315  aa = sqrt(aa)
316  IF (az.GT.aa) s1 = z*z*cmplx(0.5e0,0.0e0)
317  ai = ai + s1*cmplx(c1,0.0e0)
318  RETURN
319  180 CONTINUE
320  nz = 1
321  ai = cmplx(0.0e0,0.0e0)
322  RETURN
323  240 CONTINUE
324  nz = 0
325  ierr=2
326  RETURN
327  250 CONTINUE
328  IF(nn.EQ.(-1)) go to 240
329  nz=0
330  ierr=5
331  RETURN
332  260 CONTINUE
333  ierr=4
334  nz=0
335  RETURN
336  END
Generated on Mon Dec 30 2013 03:04:45 for GNU Octave by   doxygen 1.8.1.2