ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
airy.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file airy.f
1 
2  SUBROUTINE aiz(IFUN,IFAC,X0,Y0,GAIR,GAII,IERRO)
3 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
4 C COMPUTATION OF THE AIRY FUNCTION AI(Z) OR ITS DERIVATIVE AI'(Z)
5 C THE CODE USES:
6 C 1. MACLAURIN SERIES FOR |Y|<3 AND -2.5<X<1.3 (Z=X+I*Y)
7 C 2. GAUSS-LAGUERRE QUADRATURE FOR |Z|<15 AND WHEN
8 C MACLAURIN SERIES ARE NOT USED.
9 C 3. ASYMPTOTIC EXPANSION FOR |Z|>15.
10 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
11 C INPUTS:
12 C IFUN:
13 C * IFUN=1, THE CODE COMPUTES AI(Z)
14 C * IFUN=2, THE CODE COMPUTES AI'(Z)
15 C IFAC:
16 C * IFAC=1, THE CODE COMPUTES AI(Z) OR AI'(Z)
17 C * IFAC=2, THE CODE COMPUTES NORMALIZED AI(Z) OR AI'(Z)
18 C X0: REAL PART OF THE ARGUMENT Z
19 C Y0: IMAGINARY PART OF THE ARGUMENT Z
20 C
21 C OUTPUTS:
22 C GAIR: REAL PART OF AI(Z) OR AI'(Z)
23 C GAII: IMAGINARY PART OF AI(Z) OR AI'(Z)
24 C
25 C IERRO: ERROR FLAG
26 C * IERRO=0, SUCCESSFUL COMPUTATION
27 C * IERRO=1, COMPUTATION OUT OF RANGE
28 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
29 C MACHINE DEPENDENT CONSTANTS: FUNCTION D1MACH
30 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
31 C ACCURACY:
32 C
33 C 1) SCALED AIRY FUNCTIONS:
34 C RELATIVE ACCURACY BETTER THAN 10**(-13) EXCEPT CLOSE TO
35 C THE ZEROS, WHERE 10**(-13) IS THE ABSOLUTE PRECISION.
36 C GRADUAL LOSS OF PRECISION TAKES PLACE FOR |Z|>1000
37 C (REACHING 10**(-8) ABSOLUTE ACCURACY FOR |Z| CLOSE
38 C TO 10**(6)) IN THE CASE OF PHASE(Z) CLOSE TO PI.
39 C 2) UNSCALED AIRY FUNCTIONS:
40 C THE FUNCTION OVERFLOWS/UNDERFLOWS FOR
41 C 3/2*|Z|**(3/2)>LOG(OVER).
42 C FOR |Z|<30:
43 C A) RELATIVE ACCURACY FOR THE MODULUS (EXCEPT AT THE
44 C ZEROS) BETTER THAN 10**(-13).
45 C B) ABSOLUTE ACCURACY FOR MIN(R(Z),1/R(Z)) BETTER
46 C THAN 10**(-13), WHERE R(Z)=REAL(AI)/IMAG(AI)
47 C OR R(Z)=REAL(AI')/IMAG(AI').
48 C FOR |Z|>30, GRADUAL LOSS OF PRECISION TAKES PLACE
49 C AS |Z| INCREASES.
50 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
51 C AUTHORS:
52 C AMPARO GIL (U. AUTONOMA DE MADRID, MADRID, SPAIN).
53 C E-MAIL: AMPARO.GIL@UAM.ES
54 C JAVIER SEGURA (U. CARLOS III DE MADRID, MADRID, SPAIN).
55 C E-MAIL: JSEGURA@MATH.UC3M.ES
56 C NICO M. TEMME (CWI, AMSTERDAM, THE NETHERLANDS).
57 C E-MAIL: NICO.TEMME@CWI.NL
58 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
59 C REFERENCES:
60 C COMPUTING AIRY FUNCTIONS BY NUMERICAL QUADRATURE.
61 C NUMERICAL ALGORITHMS (2001).
62 C A. GIL, J. SEGURA, N.M. TEMME
63 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
64  DOUBLE PRECISION x0,y0,gair,gaii,x,w,xd,wd
65  DOUBLE PRECISION over,under,dl1,dl2,cover,d1mach
66  DOUBLE PRECISION pi,pihal,pih3,pisr,a,alf,thet,r,th15,s1,c1,
67  * r32,facto,th025,s3,c3,f23,pi23,sqrt3,xa,ya,f23r,df1,df2,
68  * s11,c11,dex,dre,dima,gar,gai,c,s,u,v,v0,ar,ai,ar1,ai1,
69  * ro,coe1,coe2,rex,dfr,dfi,ar11,ai11,phase
70  INTEGER ifun,ifac,ierro,iexpf,iexpf2,n
71  dimension x(25),w(25)
72  dimension xd(25),wd(25)
73  common/param1/pi,pihal
74  common/param2/pih3,pisr,a,alf
75  common/param3/thet,r,th15,s1,c1,r32
76  common/param4/facto,th025,s3,c3
77  SAVE x,w
78  SAVE xd,wd
79  DATA x,w/.283891417994567679d-1,.170985378860034935d0,
80  *.435871678341770460d0,.823518257913030858d0,1.33452543254227372d0,
81  *1.96968293206435071d0,2.72998134002859938d0,3.61662161916100897d0,
82  *4.63102611052654146d0,5.77485171830547694d0,7.05000568630218682d0,
83  *8.45866437513237792d0,10.0032955242749393d0,11.6866845947722423d0,
84  *13.5119659344693551d0,15.4826596959377140d0,17.6027156808069112d0,
85  *19.8765656022785451d0,22.3091856773962780d0,24.9061720212974207d0,
86  *27.6738320739497190d0,30.6192963295084111d0,33.7506560850239946d0,
87  *37.0771349708391198d0,40.6093049694341322d0,.143720408803313866d0,
88  *.230407559241880881d0,.242253045521327626d0,.203636639103440807d0,
89  *.143760630622921410d0,.869128834706078120d-1,.4541750018329
90  * 15883d-1,.206118031206069497d-1,.814278821268606972d-2,.280266
91  *075663377634d-2,.840337441621719716d-3,.219303732907765020d-3,
92  *.497401659009257760d-4,.978508095920717661d-5,.166542824603725
93  *563d-5,.244502736801316287d-6,.308537034236207072d-7,.3332960
94  *72940112245d-8,.306781892316295828d-9,.239331309885375719d-10,
95  *.157294707710054952d-11,.864936011664392267d-13,.394819815
96  *638647111d-14,.148271173082850884d-15,.453390377327054458d-17/
97  DATA xd,wd/.435079659953445d-1,.205779160144678d0,
98  *.489916161318751d0,.896390483211727d0,1.42582496737580d0,
99  *2.07903190767599d0,2.85702335104978d0,3.76102058198275d0,
100  *4.79246521225895d0,5.95303247470003d0,7.24464710774066d0,
101  *8.66950223642504d0,10.2300817341775d0,11.9291866622602d0,
102  *13.7699665302828d0,15.7559563095946d0,17.8911203751898d0,
103  *20.1799048700978d0,22.6273004064466d0,25.2389175786164d0,
104  *28.0210785229929d0,30.9809287996116d0,34.1265753192057d0,
105  *37.4672580871163d0,41.0135664833476d0,.576354557898966d-1,
106  *.139560003272262d0,.187792315011311d0,.187446935256946d0,
107  *.150716717316301d0,.101069904453380d0,.575274105486025d-1,
108  *.280625783448681d-1,.117972164134041d-1,.428701743297432d-2,
109  *.134857915232883d-2,.367337337105948d-3,.865882267841931d-4,
110  *.176391622890609d-4,.309929190938078d-5,.468479653648208d-6,
111  *.607273267228907d-7,.672514812555074d-8,.633469931761606d-9,
112  *.504938861248542d-10,.338602527895834d-11,.189738532450555d-12,
113  *.881618802142698d-14,.336676636121976d-15,.104594827170761d-16/
114 CC CONSTANTS CCCCCCCCCCCCCCCCCCCCCCC
115  pi=3.1415926535897932385d0
116  pihal=1.5707963267948966192d0
117  pih3=4.71238898038469d0
118  f23=.6666666666666666d0
119  pi23=2.09439510239320d0
120  pisr=1.77245385090552d0
121  sqrt3=1.7320508075688772935d0
122 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
123  ya=y0
124  xa=x0
125  ierro=0
126  iexpf=0
127  iexpf2=0
128  IF (ya.LT.0.d0) ya=-ya
129  r=sqrt(xa*xa+ya*ya)
130  r32=r*sqrt(r)
131  thet=phase(xa,ya)
132  cover=2.d0/3.d0*r32*abs(cos(1.5d0*thet))
133 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
134  over=d1mach(2)*1.d-3
135 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
136  under=d1mach(1)*1.d+3
137  dl1=log(over)
138  dl2=-log(under)
139  IF (dl1.GT.dl2) over=1/under
140  IF (ifac.EQ.1) THEN
141  IF (cover.GE.log(over)) THEN
142 CCC OVERFLOW/UNDERFLOW PROBLEMS.
143 CCC CALCULATION ABORTED
144  ierro=1
145  gair=0
146  gaii=0
147  ENDIF
148  IF (cover.GE.(log(over)*0.2)) iexpf2=1
149  ELSE
150  IF (cover.GE.(log(over)*0.2)) iexpf=1
151  ENDIF
152  IF (ierro.EQ.0) THEN
153  IF (ifun.EQ.1) THEN
154 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
155 CCC CALCULATION OF AI(Z) CCCCCCCCCCCCCCC
156 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
157 CCCCCCCCCCCC SERIES, INTEGRALS OR EXPANSIONS CCCCCCCCCCCCCCCCCCCCCCCC
158  IF ((ya.LT.3.d0).AND.(xa.LT.1.3d0).AND.(xa.GT.-2.5d0)) THEN
159 CCC SERIES CCC
160  CALL serai(xa,ya,gar,gai)
161  IF (ifac.EQ.2) THEN
162  thet=phase(xa,ya)
163  th15=1.5d0*thet
164  s1=sin(th15)
165  c1=cos(th15)
166  f23r=f23*r32
167  df1=f23r*c1
168  df2=f23r*s1
169  s11=sin(df2)
170  c11=cos(df2)
171  dex=exp(df1)
172  dre=dex*c11
173  dima=dex*s11
174  gair=dre*gar-dima*gai
175  gaii=dre*gai+dima*gar
176  ELSE
177  gair=gar
178  gaii=gai
179  IF (y0.EQ.0.) gaii=0.d0
180  ENDIF
181  ELSE
182  IF (r.GT.15.d0) THEN
183 CCC ASYMPTOTIC EXPANSIONS CCC
184  thet=phase(xa,ya)
185  facto=0.5d0/pisr*r**(-0.25d0)
186  IF (thet.GT.pi23) THEN
187 CCCCCCCCCCC CONNECTION FORMULAE CCCCCCCCCCCCCCCCCCCCCCCCCC
188 CCC N= 1: TRANSFORM Z TO W= U+IV=Z EXP( 2 PI I/3)
189  n=1
190  c=-0.5d0
191  s=n*0.5*sqrt3
192  u=xa*c-ya*s
193  v=xa*s+ya*c
194  v0=v
195  IF (v.LT.0.d0) v=-v
196  thet=phase(u,v)
197  th15=1.5d0*thet
198  s1=sin(th15)
199  c1=cos(th15)
200  th025=thet*0.25d0
201  s3=sin(th025)
202  c3=cos(th025)
203  CALL expai(ar1,ai1)
204  IF (v0.LT.0.d0) ai1=-ai1
205  ar=-(c*ar1-s*ai1)
206  ai=-(s*ar1+c*ai1)
207  IF (iexpf.EQ.0) THEN
208  IF (iexpf2.EQ.0) THEN
209 CCC N=-1: TRANSFORM Z TO W= U+IV=Z EXP(-2 PI I/3)
210  n=-1
211  c=-0.5d0
212  s=n*0.5*sqrt3
213  u=xa*c-ya*s
214  v=xa*s+ya*c
215  v0=v
216  IF (v.LT.0.d0) v=-v
217  thet=phase(u,v)
218  th15=1.5d0*thet
219  s1=sin(th15)
220  c1=cos(th15)
221  th025=thet*0.25d0
222  s3=sin(th025)
223  c3=cos(th025)
224  CALL expai(ar1,ai1)
225  IF (v0.LT.0.d0) ai1=-ai1
226  thet=phase(xa,ya)
227  th15=1.5d0*thet
228  s1=sin(th15)
229  c1=cos(th15)
230  ro=1.333333333333333d0*r32
231  coe1=ro*c1
232  coe2=ro*s1
233  rex=exp(coe1)
234  dfr=rex*cos(coe2)
235  dfi=rex*sin(coe2)
236  ar11=dfr*ar1-dfi*ai1
237  ai11=dfr*ai1+dfi*ar1
238  gair=ar-(c*ar11-s*ai11)
239  gaii=ai-(s*ar11+c*ai11)
240  ELSE
241  thet=phase(xa,ya)
242  th15=1.5d0*thet
243  s1=sin(th15)
244  c1=cos(th15)
245  gair=ar
246  gaii=ai
247  ENDIF
248  ELSE
249  gair=ar
250  gaii=ai
251  ENDIF
252  ELSE
253 CCCCCCC ASYMPTOTIC EXPANSION CCCCCCCCCCCCCCC
254  thet=phase(xa,ya)
255  th15=1.5d0*thet
256  s1=sin(th15)
257  c1=cos(th15)
258  th025=thet*0.25d0
259  s3=sin(th025)
260  c3=cos(th025)
261  CALL expai(gair,gaii)
262  ENDIF
263  ELSE
264 CCC INTEGRALS
265  a=0.1666666666666666d0
266  alf=-a
267  facto=0.280514117723058d0*r**(-0.25d0)
268  thet=phase(xa,ya)
269  IF (thet.LE.pihal) THEN
270  th15=1.5d0*thet
271  s1=sin(th15)
272  c1=cos(th15)
273  th025=thet*0.25d0
274  s3=sin(th025)
275  c3=cos(th025)
276  CALL airy1(x,w,gair,gaii)
277  ENDIF
278  IF ((thet.GT.pihal).AND.(thet.LE.pi23)) THEN
279  th15=1.5d0*thet
280  s1=sin(th15)
281  c1=cos(th15)
282  th025=thet*0.25d0
283  s3=sin(th025)
284  c3=cos(th025)
285  CALL airy2(x,w,gair,gaii)
286  ENDIF
287  IF (thet.GT.pi23) THEN
288  n=1
289  c=-0.5d0
290  s=n*0.5*sqrt3
291  u=xa*c-ya*s
292  v=xa*s+ya*c
293  v0=v
294  IF (v.LT.0.d0) v=-v
295  thet=phase(u,v)
296  IF (thet.LE.pihal) THEN
297  th15=1.5d0*thet
298  s1=sin(th15)
299  c1=cos(th15)
300  th025=thet*0.25d0
301  s3=sin(th025)
302  c3=cos(th025)
303  CALL airy1(x,w,ar1,ai1)
304  ENDIF
305  IF ((thet.GT.pihal).AND.(thet.LE.pi23)) THEN
306  th15=1.5d0*thet
307  s1=sin(th15)
308  c1=cos(th15)
309  th025=thet*0.25d0
310  s3=sin(th025)
311  c3=cos(th025)
312  CALL airy2(x,w,ar1,ai1)
313  ENDIF
314  IF (v0.LT.0.d0) ai1=-ai1
315  ar=-(c*ar1-s*ai1)
316  ai=-(s*ar1+c*ai1)
317  n=-1
318  c=-0.5d0
319  s=n*0.5*sqrt3
320  u=xa*c-ya*s
321  v=xa*s+ya*c
322  v0=v
323  IF (v.LT.0.d0) v=-v
324  thet=phase(u,v)
325  IF (thet.LE.pihal) THEN
326  th15=1.5d0*thet
327  s1=sin(th15)
328  c1=cos(th15)
329  th025=thet*0.25d0
330  s3=sin(th025)
331  c3=cos(th025)
332  CALL airy1(x,w,ar1,ai1)
333  ENDIF
334  IF ((thet.GT.pihal).AND.(thet.LE.pi23)) THEN
335  th15=1.5d0*thet
336  s1=sin(th15)
337  c1=cos(th15)
338  th025=thet*0.25d0
339  s3=sin(th025)
340  c3=cos(th025)
341  CALL airy2(x,w,ar1,ai1)
342  ENDIF
343  IF (v0.LT.0.d0) ai1=-ai1
344  thet=phase(xa,ya)
345  th15=1.5d0*thet
346  s1=sin(th15)
347  c1=cos(th15)
348  ro=1.333333333333333d0*r32
349  coe1=ro*c1
350  coe2=ro*s1
351  rex=exp(coe1)
352  dfr=rex*cos(coe2)
353  dfi=rex*sin(coe2)
354  ar11=dfr*ar1-dfi*ai1
355  ai11=dfr*ai1+dfi*ar1
356  gair=ar-(c*ar11-s*ai11)
357  gaii=ai-(s*ar11+c*ai11)
358  ENDIF
359  ENDIF
360  IF (ifac.EQ.1) THEN
361 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
362 CCC CALCULATION OF THE UNSCALED AI(Z) CCCCCCCCC
363 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
364  f23r=f23*r32
365  df1=f23r*c1
366  df2=f23r*s1
367  s11=sin(df2)
368  c11=cos(df2)
369  dex=exp(-df1)
370  dre=dex*c11
371  dima=-dex*s11
372  gar=dre*gair-dima*gaii
373  gai=dre*gaii+dima*gair
374  gair=gar
375  gaii=gai
376  IF (y0.EQ.0.) gaii=0.d0
377  ENDIF
378  ENDIF
379  ELSE
380 CCCC CALCULATION OF AI´(Z) CCCCCCCCCCC
381  alf=0.1666666666666666d0
382  facto=-0.270898621247918d0*r**0.25d0
383 CCCCCCCCCCCCCCC SERIES OR INTEGRALS CCCCCCCCCCCCCCCCCCCCCCCCCC
384  IF ((ya.LT.3.d0).AND.(xa.LT.1.3d0).AND.(xa.GT.-2.5d0)) THEN
385 CCC SERIES
386  CALL seraid(xa,ya,gar,gai)
387  IF (ifac.EQ.2) THEN
388  thet=phase(xa,ya)
389  th15=1.5d0*thet
390  s1=sin(th15)
391  c1=cos(th15)
392  f23r=f23*r32
393  df1=f23r*c1
394  df2=f23r*s1
395  s11=sin(df2)
396  c11=cos(df2)
397  dex=exp(df1)
398  dre=dex*c11
399  dima=dex*s11
400  gair=dre*gar-dima*gai
401  gaii=dre*gai+dima*gar
402  ELSE
403  gair=gar
404  gaii=gai
405  IF (y0.EQ.0.) gaii=0.d0
406  ENDIF
407  ELSE
408  IF (r.GT.15.d0) THEN
409 CCC ASYMPTOTIC EXPANSIONS CCCCCCCCCCCCC
410  thet=phase(xa,ya)
411  facto=0.5d0/pisr*r**0.25d0
412  IF (thet.GT.pi23) THEN
413 CCCCCCC CONNECTION FORMULAE CCCCCCCCCCCCCCCCCCCCCCCCCCCCC
414 CCC N= 1: TRANSFORM Z TO W= U+IV=Z EXP( 2 PI I/3)
415  n=1
416  c=-0.5d0
417  s=n*0.5*sqrt3
418  u=xa*c-ya*s
419  v=xa*s+ya*c
420  v0=v
421  IF (v.LT.0.d0) v=-v
422  thet=phase(u,v)
423  th15=1.5d0*thet
424  s1=sin(th15)
425  c1=cos(th15)
426  th025=thet*0.25d0
427  s3=sin(th025)
428  c3=cos(th025)
429  CALL expaid(ar1,ai1)
430  IF (v0.LT.0.d0) ai1=-ai1
431  ar=-(c*ar1+s*ai1)
432  ai=-(-s*ar1+c*ai1)
433  IF (iexpf.EQ.0) THEN
434  IF (iexpf2.EQ.0) THEN
435 CCC N=-1: TRANSFORM Z TO W= U+IV=Z EXP(-2 PI I/3)
436  n=-1
437  c=-0.5d0
438  s=n*0.5*sqrt3
439  u=xa*c-ya*s
440  v=xa*s+ya*c
441  v0=v
442  IF (v.LT.0.d0) v=-v
443  thet=phase(u,v)
444  th15=1.5d0*thet
445  s1=sin(th15)
446  c1=cos(th15)
447  th025=thet*0.25d0
448  s3=sin(th025)
449  c3=cos(th025)
450  CALL expaid(ar1,ai1)
451  IF (v0.LT.0.d0) ai1=-ai1
452  thet=phase(xa,ya)
453  th15=1.5d0*thet
454  s1=sin(th15)
455  c1=cos(th15)
456  ro=1.333333333333333d0*r32
457  coe1=ro*c1
458  coe2=ro*s1
459  rex=exp(coe1)
460  dfr=rex*cos(coe2)
461  dfi=rex*sin(coe2)
462  ar11=dfr*ar1-dfi*ai1
463  ai11=dfr*ai1+dfi*ar1
464  gair=ar-(c*ar11+s*ai11)
465  gaii=ai-(-s*ar11+c*ai11)
466  ELSE
467  thet=phase(xa,ya)
468  th15=1.5d0*thet
469  s1=sin(th15)
470  c1=cos(th15)
471  gair=ar
472  gaii=ai
473  ENDIF
474  ELSE
475  gair=ar
476  gaii=ai
477  ENDIF
478  ELSE
479  th15=1.5d0*thet
480  s1=sin(th15)
481  c1=cos(th15)
482  th025=thet*0.25d0
483  s3=sin(th025)
484  c3=cos(th025)
485  CALL expaid(gair,gaii)
486  ENDIF
487  ELSE
488 CCC INTEGRALS CCCCCCCCCCCCCCCC
489  thet=phase(xa,ya)
490  IF (thet.LE.pihal) THEN
491  th15=1.5d0*thet
492  s1=sin(th15)
493  c1=cos(th15)
494  th025=thet*0.25d0
495  s3=sin(th025)
496  c3=cos(th025)
497  CALL airy1d(xd,wd,gair,gaii)
498  ENDIF
499  IF ((thet.GT.pihal).AND.(thet.LE.pi23)) THEN
500  th15=1.5d0*thet
501  s1=sin(th15)
502  c1=cos(th15)
503  th025=thet*0.25d0
504  s3=sin(th025)
505  c3=cos(th025)
506  CALL airy2d(xd,wd,gair,gaii)
507  ENDIF
508  IF (thet.GT.pi23) THEN
509  n=1
510  c=-0.5d0
511  s=n*0.5*sqrt3
512  u=xa*c-ya*s
513  v=xa*s+ya*c
514  v0=v
515  IF (v.LT.0.d0) v=-v
516  thet=phase(u,v)
517  IF (thet.LE.pihal) THEN
518  th15=1.5d0*thet
519  s1=sin(th15)
520  c1=cos(th15)
521  th025=thet*0.25d0
522  s3=sin(th025)
523  c3=cos(th025)
524  CALL airy1d(xd,wd,ar1,ai1)
525  ENDIF
526  IF ((thet.GT.pihal).AND.(thet.LE.pi23)) THEN
527  th15=1.5d0*thet
528  s1=sin(th15)
529  c1=cos(th15)
530  th025=thet*0.25d0
531  s3=sin(th025)
532  c3=cos(th025)
533  CALL airy2d(xd,wd,ar1,ai1)
534  ENDIF
535  IF (v0.LT.0.d0) ai1=-ai1
536  ar=-(c*ar1+s*ai1)
537  ai=-(-s*ar1+c*ai1)
538  n=-1
539  c=-0.5d0
540  s=n*0.5*sqrt3
541  u=xa*c-ya*s
542  v=xa*s+ya*c
543  v0=v
544  IF (v.LT.0.d0) v=-v
545  thet=phase(u,v)
546  IF (thet.LE.pihal) THEN
547  th15=1.5d0*thet
548  s1=sin(th15)
549  c1=cos(th15)
550  th025=thet*0.25d0
551  s3=sin(th025)
552  c3=cos(th025)
553  CALL airy1d(xd,wd,ar1,ai1)
554  ENDIF
555  IF ((thet.GT.pihal).AND.(thet.LE.pi23)) THEN
556  th15=1.5d0*thet
557  s1=sin(th15)
558  c1=cos(th15)
559  th025=thet*0.25d0
560  s3=sin(th025)
561  c3=cos(th025)
562  CALL airy2d(xd,wd,ar1,ai1)
563  ENDIF
564  IF (v0.LT.0.d0) ai1=-ai1
565  thet=phase(xa,ya)
566  th15=1.5d0*thet
567  s1=sin(th15)
568  c1=cos(th15)
569  ro=1.333333333333333d0*r32
570  coe1=ro*c1
571  coe2=ro*s1
572  rex=exp(coe1)
573  dfr=rex*cos(coe2)
574  dfi=rex*sin(coe2)
575  ar11=dfr*ar1-dfi*ai1
576  ai11=dfr*ai1+dfi*ar1
577  gair=ar-(c*ar11+s*ai11)
578  gaii=ai-(-s*ar11+c*ai11)
579  ENDIF
580  ENDIF
581  IF (ifac.EQ.1) THEN
582 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
583 CCC CALCULATION OF THE UNSCALED AI'(z) CCCCCCC
584 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
585  f23r=f23*r32
586  df1=f23r*c1
587  df2=f23r*s1
588  s11=sin(df2)
589  c11=cos(df2)
590  dex=exp(-df1)
591  dre=dex*c11
592  dima=-dex*s11
593  gar=dre*gair-dima*gaii
594  gai=dre*gaii+dima*gair
595  gair=gar
596  gaii=gai
597  IF (y0.EQ.0) gaii=0.d0
598  ENDIF
599  ENDIF
600  ENDIF
601  ENDIF
602  IF (y0.LT.0.d0) gaii=-gaii
603  RETURN
604  END
605  SUBROUTINE airy1(X,W,GAIR,GAII)
606 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
607 CCC COMPUTES AI(Z) BY GAUSS-LAGUERRE QUADRATURE IN THE SECTOR C
608 CCC 0 <= PHASE(Z) <= PI/2 C
609 CCC C
610 CCC INPUTS: C
611 CCC X,W, NODES AND WEIGHTS FOR THE GAUSSIAN C
612 CCC QUADRATURE C
613 CCC OUTPUTS: C
614 CCC GAIR, GAII, REAL AND IMAGINARY PARTS OF AI(Z) C
615 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
616  DOUBLE PRECISION x,w,gair,gaii
617  DOUBLE PRECISION pih3,pisr,a,alf,thet,r,th15,s1,c1,
618  * r32,facto,th025,s3,c3,sumar,sumai,df1,df1c1,phi,phi6,
619  * s2,c2,dmodu,dmodu2,funr,funi,fac1,fac2,phase
620  INTEGER i
621  dimension x(25),w(25)
622  common/param2/pih3,pisr,a,alf
623  common/param3/thet,r,th15,s1,c1,r32
624  common/param4/facto,th025,s3,c3
625  sumar=0.d0
626  sumai=0.d0
627  DO 1 i=1,25
628  df1=1.5d0*x(i)/r32
629  df1c1=df1*c1
630  phi=phase(2.d0+df1c1,df1*s1)
631  phi6=phi/6.d0
632  s2=sin(phi6)
633  c2=cos(phi6)
634  dmodu=sqrt(4.d0+df1*df1+4.d0*df1c1)
635  dmodu2=dmodu**alf
636  funr=dmodu2*c2
637  funi=dmodu2*s2
638  sumar=sumar+w(i)*funr
639  sumai=sumai+w(i)*funi
640  1 CONTINUE
641  fac1=facto*c3
642  fac2=facto*s3
643  gair=fac1*sumar+fac2*sumai
644  gaii=fac1*sumai-fac2*sumar
645  RETURN
646  END
647  SUBROUTINE airy2(X,W,GAIR,GAII)
648 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
649 CCC COMPUTES AI(Z) BY GAUSS-LAGUERRE QUADRATURE IN THE SECTOR C
650 CCC PI/2 < PHASE(Z) <= 2PI/3 C
651 CCC C
652 CCC INPUTS: C
653 CCC X,W, NODES AND WEIGHTS FOR THE GAUSSIAN C
654 CCC QUADRATURE C
655 CCC OUTPUTS: C
656 CCC GAIR, GAII, REAL AND IMAGINARY PARTS OF AI(Z) C
657 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
658  DOUBLE PRECISION x,w,gair,gaii
659  DOUBLE PRECISION pih3,pisr,a,alf,thet,r,th15,s1,c1,
660  * r32,facto,th025,s3,c3,sumar,sumai,df1,df1c1,phi,phi6,
661  * s2,c2,dmodu,dmodu2,funr,funi,fac1,fac2,phase
662  DOUBLE PRECISION sqr2,sqr2i,tau,tgtau,b,ang,ctau,cfac,ct,st,
663  * sumr,sumi,ttau,beta
664  INTEGER i
665  dimension x(25),w(25)
666  common/param2/pih3,pisr,a,alf
667  common/param3/thet,r,th15,s1,c1,r32
668  common/param4/facto,th025,s3,c3
669  sqr2=1.41421356237310d0
670  sqr2i=0.707106781186548d0
671  tau=th15-pih3*0.5d0
672  tgtau=dtan(tau)
673  b=5.d0*a
674  ang=tau*b
675  ctau=cos(tau)
676  cfac=ctau**(-b)
677  ct=cos(ang)
678  st=sin(ang)
679  sumr=0.d0
680  sumi=0.d0
681  DO 2 i=1,25
682  df1=3.d0*x(i)/(ctau*r32)
683  df1c1=df1*sqr2i*0.5d0
684  phi=phase(2.d0-df1c1,df1c1)
685  phi6=phi/6.d0
686  ttau=x(i)*tgtau
687  beta=phi6-ttau
688  s2=sin(beta)
689  c2=cos(beta)
690  dmodu=sqrt(4.d0+df1*df1*0.25d0-sqr2*df1)
691  dmodu2=dmodu**alf
692  funr=dmodu2*c2
693  funi=dmodu2*s2
694  sumr=sumr+w(i)*funr
695  sumi=sumi+w(i)*funi
696  2 CONTINUE
697  sumar=cfac*(ct*sumr-st*sumi)
698  sumai=cfac*(ct*sumi+st*sumr)
699  fac1=facto*c3
700  fac2=facto*s3
701  gair=fac1*sumar+fac2*sumai
702  gaii=fac1*sumai-fac2*sumar
703  RETURN
704  END
705  SUBROUTINE airy1d(X,W,GAIR,GAII)
706 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
707 CCC COMPUTES AI'(Z) BY GAUSS-LAGUERRE QUADRATURE IN THE SECTOR C
708 CCC 0 <= PHASE(Z) <= PI/2 C
709 CCC C
710 CCC INPUTS: C
711 CCC X,W, NODES AND WEIGHTS FOR THE GAUSSIAN C
712 CCC QUADRATURE C
713 CCC OUTPUTS: C
714 CCC GAIR,GAII, REAL AND IMAGINARY PARTS OF AI'(Z) C
715 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
716  DOUBLE PRECISION x,w,gair,gaii
717  DOUBLE PRECISION pih3,pisr,a,alf,thet,r,th15,s1,c1,
718  * r32,facto,th025,s3,c3,sumar,sumai,df1,df1c1,phi,phi6,
719  * s2,c2,dmodu,dmodu2,funr,funi,fac1,fac2,phase
720  INTEGER i
721  dimension x(25),w(25)
722  common/param2/pih3,pisr,a,alf
723  common/param3/thet,r,th15,s1,c1,r32
724  common/param4/facto,th025,s3,c3
725  sumar=0.d0
726  sumai=0.d0
727  DO 3 i=1,25
728  df1=1.5d0*x(i)/r32
729  df1c1=df1*c1
730  phi=phase(2.d0+df1c1,df1*s1)
731  phi6=-phi*alf
732  s2=sin(phi6)
733  c2=cos(phi6)
734  dmodu=sqrt(4.d0+df1*df1+4.d0*df1c1)
735  dmodu2=dmodu**alf
736  funr=dmodu2*c2
737  funi=dmodu2*s2
738  sumar=sumar+w(i)*funr
739  sumai=sumai+w(i)*funi
740  3 CONTINUE
741  fac1=facto*c3
742  fac2=facto*s3
743  gair=fac1*sumar-fac2*sumai
744  gaii=fac1*sumai+fac2*sumar
745  RETURN
746  END
747  SUBROUTINE airy2d(X,W,GAIR,GAII)
748 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
749 CCC COMPUTES AI'(Z) BY GAUSS-LAGUERRE QUADRATURE IN THE SECTOR C
750 CCC PI/2 < PHASE(Z) <= 3PI/2 C
751 CCC C
752 CCC INPUTS: C
753 CCC X,W, NODES AND WEIGHTS FOR THE GAUSSIAN C
754 CCC QUADRATURE C
755 CCC OUTPUTS: C
756 CCC GAIR,GAII, REAL AND IMAGINARY PARTS OF AI'(Z) C
757 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
758  DOUBLE PRECISION x,w,gair,gaii
759  DOUBLE PRECISION pih3,pisr,a,alf,thet,r,th15,s1,c1,
760  * r32,facto,th025,s3,c3,sumar,sumai,df1,df1c1,phi,phi6,
761  * s2,c2,dmodu,dmodu2,funr,funi,fac1,fac2,phase
762  DOUBLE PRECISION sqr2,sqr2i,tau,tgtau,b,ang,ctau,cfac,ct,st,
763  * sumr,sumi,ttau,beta
764  INTEGER i
765  dimension x(25),w(25)
766  common/param2/pih3,pisr,a,alf
767  common/param3/thet,r,th15,s1,c1,r32
768  common/param4/facto,th025,s3,c3
769  sqr2=1.41421356237310d0
770  sqr2i=0.707106781186548d0
771  tau=th15-pih3*0.5d0
772  tgtau=dtan(tau)
773  b=7.d0*alf
774  ang=tau*b
775  ctau=cos(tau)
776  cfac=ctau**(-b)
777  ct=cos(ang)
778  st=sin(ang)
779  sumr=0.d0
780  sumi=0.d0
781  DO 4 i=1,25
782  df1=3.d0*x(i)/(ctau*r32)
783  df1c1=df1*sqr2i*0.5d0
784  phi=phase(2.d0-df1c1,df1c1)
785  phi6=-phi/6.d0
786  ttau=x(i)*tgtau
787  beta=phi6-ttau
788  s2=sin(beta)
789  c2=cos(beta)
790  dmodu=sqrt(4.d0+df1*df1*0.25d0-sqr2*df1)
791  dmodu2=dmodu**alf
792  funr=dmodu2*c2
793  funi=dmodu2*s2
794  sumr=sumr+w(i)*funr
795  sumi=sumi+w(i)*funi
796  4 CONTINUE
797  sumar=cfac*(ct*sumr-st*sumi)
798  sumai=cfac*(ct*sumi+st*sumr)
799  fac1=facto*c3
800  fac2=facto*s3
801  gair=fac1*sumar-fac2*sumai
802  gaii=fac1*sumai+fac2*sumar
803  RETURN
804  END
805  DOUBLE PRECISION FUNCTION phase(X,Y)
806  DOUBLE PRECISION pi,pihal,x,y,ay,p
807 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
808 CCC COMPUTES THE PHASE OF Z = X + IY, IN (-PI,PI]
809 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
810  common/param1/pi,pihal
811  IF ((x.EQ.0).AND.(y.EQ.0)) THEN
812  p=0.d0
813  ELSE
814  ay=abs(y)
815  IF (x.GE.ay) THEN
816  p=atan(ay/x)
817  ELSEIF ((x+ay).GE.0.d0) THEN
818  p=pihal-atan(x/ay)
819  ELSE
820  p=pi+atan(ay/x)
821  ENDIF
822  IF (y.LT.0.d0) p=-p
823  ENDIF
824  phase=p
825  END
826  SUBROUTINE fgp(X,Y,EPS,FR,FI,GR,GI)
827 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
828 C COMPUTES THE FUNCTIONS F AND G FOR THE SERIES C
829 C OF AI'(Z). C
830 C THIS ROUTINE IS CALLED BY SERAID. C
831 C C
832 C INPUTS: C
833 C X,Y, REAL AND IMAGINARY PARTS OF Z C
834 C EPS, PRECISION FOR THE COMPUTATION OF C
835 C THE SERIES C
836 C OUTPUTS: C
837 C FR,FI, REAL AND IMAGINARY PARTS OF F C
838 C GR,GI, REAL AND IMAGINARY PARTS OF G C
839 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
840  DOUBLE PRECISION x,y,eps,fr,fi,gr,gi
841  INTEGER a,b,k3
842  DOUBLE PRECISION x2,y2,u,v,p,q,cr,ci,dr,di
843  x2=x*x
844  y2=y*y
845  k3=0
846  u=x*(x2-3*y2)
847  v=y*(3*x2-y2)
848  cr=0.5d0
849  ci=0.d0
850  dr=1.d0
851  di=0.d0
852  fr=0.5d0
853  fi=0.d0
854  gr=1.d0
855  gi=0.d0
856  70 a=(k3+5)*(k3+3)
857  b=(k3+1)*(k3+3)
858  p=(u*cr-v*ci)/a
859  q=(v*cr+u*ci)/a
860  cr=p
861  ci=q
862  p=(u*dr-v*di)/b
863  q=(v*dr+u*di)/b
864  dr=p
865  di=q
866  fr=fr+cr
867  fi=fi+ci
868  gr=gr+dr
869  gi=gi+di
870  k3=k3+3
871  IF ((abs(cr)+abs(dr)+abs(ci)+abs(di)).GE.eps) goto 70
872  u=x2-y2
873  v=2.d0*x*y
874  p=u*fr-v*fi
875  q=u*fi+v*fr
876  fr=p
877  fi=q
878  RETURN
879  END
880  SUBROUTINE fg(X,Y,EPS,FR,FI,GR,GI)
881 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
882 C COMPUTES THE FUNCTIONS F AND G IN EXPRESSION C
883 C 10.4.2 OF ABRAMOWITZ & STEGUN FOR THE SERIES C
884 C OF AI(Z). C
885 C THIS ROUTINE IS CALLED BY SERAI. C
886 C C
887 C INPUTS: C
888 C X,Y, REAL AND IMAGINARY PARTS OF Z C
889 C EPS, PRECISION FOR THE COMPUTATION C
890 C OF THE SERIES. C
891 C OUTPUTS: C
892 C FR,FI, REAL AND IMAGINARY PARTS OF F C
893 C GR,GI, REAL AND IMAGINARY PARTS OF G C
894 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
895  INTEGER a,b,k3
896  DOUBLE PRECISION x2,y2,u,v,p,q,cr,ci,dr,di
897  DOUBLE PRECISION x,y,eps,fr,fi,gr,gi
898  x2=x*x
899  y2=y*y
900  k3=0
901  u=x*(x2-3.d0*y2)
902  v=y*(3.d0*x2-y2)
903  cr=1.d0
904  ci=0.d0
905  dr=1.d0
906  di=0.d0
907  fr=1.d0
908  fi=0.d0
909  gr=1.d0
910  gi=0.d0
911  71 a=(k3+2)*(k3+3)
912  b=(k3+4)*(k3+3)
913  p=(u*cr-v*ci)/a
914  q=(v*cr+u*ci)/a
915  cr=p
916  ci=q
917  p=(u*dr-v*di)/b
918  q=(v*dr+u*di)/b
919  dr=p
920  di=q
921  fr=fr+cr
922  fi=fi+ci
923  gr=gr+dr
924  gi=gi+di
925  k3=k3+3
926  IF ((abs(cr)+abs(dr)+abs(ci)+abs(di)).GE.eps) goto 71
927  p=x*gr-y*gi
928  q=x*gi+y*gr
929  gr=p
930  gi=q
931  RETURN
932  END
933  SUBROUTINE serai(X,Y,AIR,AII)
934 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
935 CCC AIRY AI(Z), TAYLOR, COMPLEX Z CCC
936 CCC CCC
937 CCC INPUTS: CCC
938 CCC X,Y, REAL AND IMAGINARY CCC
939 CCC PARTS OF Z CCC
940 CCC OUTPUTS: CCC
941 CCC AIR,AII, REAL AND IMAGINARY CCC
942 CCC PARTS OF AI(Z) CCC
943 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
944  DOUBLE PRECISION x,y,eps,air,aii
945  DOUBLE PRECISION fzr,fzi,gzr,gzi,cons1,cons2
946  DOUBLE PRECISION d1mach
947  eps=d1mach(3)
948  IF (eps.LT.1.d-15) eps=1.d-15
949  cons1=0.355028053887817239260d0
950  cons2=0.258819403792806798405d0
951  CALL fg(x,y,eps,fzr,fzi,gzr,gzi)
952  air=cons1*fzr-cons2*gzr
953  aii=cons1*fzi-cons2*gzi
954  RETURN
955  END
956  SUBROUTINE seraid(X,Y,AIR,AII)
957 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
958 CCC AIRY AI'(Z), TAYLOR, COMPLEX Z CCC
959 CCC CCC
960 CCC INPUTS: CCC
961 CCC X,Y, REAL AND IMAGINARY CCC
962 CCC PARTS OF Z CCC
963 CCC OUTPUTS: CCC
964 CCC AIR,AII, REAL AND IMAGINARY CCC
965 CCC PARTS OF AI'(Z) CCC
966 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
967  DOUBLE PRECISION x,y,eps,air,aii
968  DOUBLE PRECISION fzr,fzi,gzr,gzi,cons1,cons2
969  DOUBLE PRECISION d1mach
970  eps=d1mach(3)
971  IF (eps.LT.1.d-15) eps=1.d-15
972  cons1=0.355028053887817239260d0
973  cons2=0.258819403792806798405d0
974  CALL fgp(x,y,eps,fzr,fzi,gzr,gzi)
975  air=cons1*fzr-cons2*gzr
976  aii=cons1*fzi-cons2*gzi
977  RETURN
978  END
979  SUBROUTINE expai(GAIR,GAII)
980 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
981 CCC AIRY AI(Z), ASYMPTOTIC EXPANSION, COMPLEX Z CCC
982 CCC CCC
983 CCC OUTPUTS: CCC
984 CCC GAIR, GAII, REAL AND IMAGINARY CCC
985 CCC PARTS OF AI(Z) CCC
986 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
987  DOUBLE PRECISION eps,gair,gaii
988  DOUBLE PRECISION thet,r,th15,s1,
989  * c1,r32,facto,th025,s3,c3
990  DOUBLE PRECISION df1,psiir,psiii,ck,dfrr,dfii,sumar,sumai,
991  * dfr,dfi,deltar,deltai,fac1,fac2
992  DOUBLE PRECISION co,df
993  DOUBLE PRECISION d1mach
994  INTEGER k
995  dimension co(20)
996  common/param3/thet,r,th15,s1,c1,r32
997  common/param4/facto,th025,s3,c3
998  SAVE co
999  DATA co/-6.944444444444445d-2,3.713348765432099d-2,
1000  * -3.799305912780064d-2,5.764919041266972d-2,-0.116099064025515d0,
1001  * 0.291591399230751d0,-0.877666969510017d0,3.07945303017317d0,
1002  * -12.3415733323452d0,55.6227853659171d0,-278.465080777603d0,
1003  * 1533.16943201280d0,-9207.20659972641d0,59892.5135658791d0,
1004  * -419524.875116551d0,3148257.41786683d0,-25198919.8716024d0,
1005  * 214288036.963680d0,-1929375549.18249d0,18335766937.8906d0/
1006  eps=d1mach(3)
1007  IF (eps.LT.1.d-15) eps=1.d-15
1008  df1=1.5d0/r32
1009  psiir=c1
1010  psiii=-s1
1011  k=0
1012  ck=1.d0
1013  df=1.d0
1014  dfrr=1.d0
1015  dfii=0.d0
1016  sumar=1.d0
1017  sumai=0.d0
1018 80 df=df*df1
1019  ck=co(k+1)*df
1020  dfr=dfrr
1021  dfi=dfii
1022  dfrr=dfr*psiir-dfi*psiii
1023  dfii=dfr*psiii+dfi*psiir
1024  deltar=dfrr*ck
1025  deltai=dfii*ck
1026  sumar=sumar+deltar
1027  sumai=sumai+deltai
1028  k=k+1
1029  IF (sumar.NE.0) THEN
1030  IF (abs(deltar/sumar).GT.eps) goto 80
1031  ENDIF
1032  IF (sumai.NE.0) THEN
1033  IF (abs(deltai/sumai).GT.eps) goto 80
1034  ENDIF
1035  fac1=facto*c3
1036  fac2=facto*s3
1037  gair=fac1*sumar+fac2*sumai
1038  gaii=fac1*sumai-fac2*sumar
1039  RETURN
1040  END
1041  SUBROUTINE expaid(GAIR,GAII)
1042 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1043 CCC AIRY AI'(Z), ASYMPTOTIC EXPANSION, COMPLEX Z CCC
1044 CCC CCC
1045 CCC OUTPUTS: CCC
1046 CCC GAIR, GAII, REAL AND IMAGINARY CCC
1047 CCC PARTS OF AI'(Z) CCC
1048 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1049  DOUBLE PRECISION eps,gair,gaii
1050  DOUBLE PRECISION thet,r,th15,s1,
1051  * c1,r32,facto,th025,s3,c3
1052  DOUBLE PRECISION df1,psiir,psiii,vk,dfrr,dfii,sumar,sumai,
1053  * dfr,dfi,deltar,deltai,fac1,fac2
1054  DOUBLE PRECISION co,df
1055  DOUBLE PRECISION d1mach
1056  INTEGER k
1057  dimension co(20)
1058  common/param3/thet,r,th15,s1,c1,r32
1059  common/param4/facto,th025,s3,c3
1060  SAVE co
1061  DATA co/9.722222222222222d-2,-4.388503086419753d-2,
1062  * 4.246283078989484d-2,-6.266216349203230d-2,
1063  * 0.124105896027275d0,-0.308253764901079d0,
1064  * 0.920479992412945d0,-3.21049358464862d0,
1065  * 12.8072930807356d0,-57.5083035139143d0,
1066  * 287.033237109221d0,-1576.35730333710d0,
1067  * 9446.35482309593d0,-61335.7066638521d0,
1068  * 428952.400400069d0,-3214536.52140086d0,
1069  * 25697908.3839113d0,-218293420.832160d0,
1070  * 1963523788.99103d0,-18643931088.1072d0/
1071  eps=d1mach(3)
1072  IF (eps.LT.1.d-15) eps=1.d-15
1073  df1=1.5d0/r32
1074  psiir=c1
1075  psiii=-s1
1076  k=0
1077  df=1.d0
1078  dfrr=1.d0
1079  dfii=0.d0
1080  sumar=1.d0
1081  sumai=0.d0
1082  81 df=df*df1
1083  vk=co(k+1)*df
1084  dfr=dfrr
1085  dfi=dfii
1086  dfrr=dfr*psiir-dfi*psiii
1087  dfii=dfr*psiii+dfi*psiir
1088  deltar=dfrr*vk
1089  deltai=dfii*vk
1090  sumar=sumar+deltar
1091  sumai=sumai+deltai
1092  k=k+1
1093  IF (sumar.NE.0) THEN
1094  IF (abs(deltar/sumar).GT.eps) goto 81
1095  ENDIF
1096  IF (sumai.NE.0) THEN
1097  IF (abs(deltai/sumai).GT.eps) goto 81
1098  ENDIF
1099  fac1=facto*c3
1100  fac2=facto*s3
1101  gair=-(fac1*sumar-fac2*sumai)
1102  gaii=-(fac1*sumai+fac2*sumar)
1103  RETURN
1104  END
1105 
1106  SUBROUTINE biz(IFUN,IFAC,X0,Y0,GBIR,GBII,IERRO)
1107 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1108 C COMPUTATION OF THE AIRY FUNCTION BI(Z) OR ITS DERIVATIVE BI'(Z)
1109 C THE CODE USES THE CONNECTION OF BI(Z) WITH AI(Z).
1110 C BIZ CALLS THE ROUTINE AIZ
1111 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1112 C INPUTS:
1113 C IFUN:
1114 C * IFUN=1, THE CODE COMPUTES BI(Z)
1115 C * IFUN=2, THE CODE COMPUTES BI'(Z)
1116 C IFAC:
1117 C * IFAC=1, THE CODE COMPUTES BI(Z) OR BI'(Z)
1118 C * IFAC=2, THE CODE COMPUTES NORMALIZED BI(Z) OR BI'(Z)
1119 C X0: REAL PART OF THE ARGUMENT Z
1120 C Y0: IMAGINARY PART OF THE ARGUMENT Z
1121 C
1122 C OUTPUTS:
1123 C GBIR: REAL PART OF BI(Z) OR BI'(Z)
1124 C GBII: IMAGINARY PART OF BI(Z) OR BI'(Z)
1125 C
1126 C IERRO: ERROR FLAG
1127 C * IERRO=0, SUCCESSFUL COMPUTATION
1128 C * IERRO=1, COMPUTATION OUT OF RANGE
1129 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1130 C MACHINE DEPENDENT CONSTANTS: FUNCTION D1MACH
1131 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1132 C ACCURACY:
1133 C
1134 C 1) SCALED AIRY FUNCTIONS:
1135 C RELATIVE ACCURACY BETTER THAN 10**(-13) EXCEPT CLOSE TO
1136 C THE ZEROS, WHERE 10**(-13) IS THE ABSOLUTE PRECISION.
1137 C GRADUAL LOSS OF PRECISION TAKES PLACE FOR |Z|>1000
1138 C IN THE CASE OF PHASE(Z) CLOSE TO +3*PI/2 OR -3*PI/2.
1139 C 2) UNSCALED AIRY FUNCTIONS:
1140 C THE FUNCTION OVERFLOWS/UNDERFLOWS FOR
1141 C 3/2*|Z|**(3/2)>LOG(OVER).
1142 C FOR |Z|<30:
1143 C A) RELATIVE ACCURACY FOR THE MODULUS (EXCEPT AT THE
1144 C ZEROS) BETTER THAN 10**(-13).
1145 C B) ABSOLUTE ACCURACY FOR MIN(R(Z),1/R(Z)) BETTER
1146 C THAN 10**(-13), WHERE R(Z)=REAL(BI)/IMAG(BI)
1147 C OR R(Z)=REAL(BI')/IMAG(BI').
1148 C FOR |Z|>30, GRADUAL LOSS OF PRECISION TAKES PLACE
1149 C AS |Z| INCREASES.
1150 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1151 C AUTHORS:
1152 C AMPARO GIL (U. AUTONOMA DE MADRID, MADRID, SPAIN).
1153 C E-MAIL: AMPARO.GIL@UAM.ES
1154 C JAVIER SEGURA (U. CARLOS III DE MADRID, MADRID, SPAIN).
1155 C E-MAIL: JSEGURA@MATH.UC3M.ES
1156 C NICO M. TEMME (CWI, AMSTERDAM, THE NETHERLANDS).
1157 C E-MAIL: NICO.TEMME@CWI.NL
1158 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1159 C REFERENCES:
1160 C COMPUTING AIRY FUNCTIONS BY NUMERICAL QUADRATURE.
1161 C NUMERICAL ALGORITHMS (2001).
1162 C A. GIL, J. SEGURA, N.M. TEMME
1163 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1164  DOUBLE PRECISION x0,y0,gbir,gbii
1165  DOUBLE PRECISION over,under,dl1,dl2,cover,d1mach
1166  DOUBLE PRECISION pi,pi3,pi23,sqrt3,c,s,c1,s1,u,v,x,y,ar,ai,
1167  * apr,api,br,bi,bpr,bpi,bbr,bbi,bbpr,bbpi,phase
1168  DOUBLE PRECISION thet,r,r32,thet32,a1,b1,df1,expo,expoi,
1169  * etar,etai,etagr,etagi,pihal
1170  INTEGER ifun,ifac,iexpf,ierr,ierro
1171  common/param1/pi,pihal
1172  sqrt3=1.7320508075688772935d0
1173  pi=3.1415926535897932385d0
1174  pihal=1.5707963267948966192d0
1175  pi3=pi/3.d0
1176  pi23=2.d0*pi3
1177  x=x0
1178  c=0.5d0*sqrt3
1179  s=0.5d0
1180  ierro=0
1181  iexpf=0
1182  IF (y0.LT.0.d0) THEN
1183  y=-y0
1184  ELSE
1185  y=y0
1186  ENDIF
1187  r=sqrt(x*x+y*y)
1188  r32=r*sqrt(r)
1189  thet=phase(x,y)
1190  cover=2.d0/3.d0*r32*abs(cos(1.5d0*thet))
1191 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
1192  over=d1mach(2)*1.d-3
1193 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
1194  under=d1mach(1)*1.d+3
1195  dl1=log(over)
1196  dl2=-log(under)
1197  IF (dl1.GT.dl2) over=1/under
1198  IF (ifac.EQ.1) THEN
1199  IF (cover.GE.log(over)) THEN
1200 CCC OVERFLOW/UNDERFLOW PROBLEMS.
1201 CCC CALCULATION ABORTED
1202  ierro=1
1203  gbir=0
1204  gbii=0
1205  ENDIF
1206  ELSE
1207  IF (cover.GE.(log(over)*0.2)) iexpf=1
1208  ENDIF
1209  IF (ierro.EQ.0) THEN
1210  IF (ifac.EQ.1) THEN
1211  IF (y.EQ.0.d0) THEN
1212 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1213 CCC TAKE TWICE THE REAL PART OF EXP(-PI I/6) AI_(1)(Z)
1214 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1215  c1=-0.5d0
1216  s1=-0.5d0*sqrt3
1217  u=x*c1-y*s1
1218  v=x*s1+y*c1
1219  CALL aiz(ifun,ifac,u,v,ar,ai,ierr)
1220  IF (ifun.EQ.1) THEN
1221  br=sqrt3*ar+ai
1222  bi=0.d0
1223  ELSE
1224  u=ar*c1-ai*s1
1225  v=ar*s1+ai*c1
1226  apr=u
1227  api=v
1228  bpr=sqrt3*apr+api
1229  bpi=0.d0
1230  ENDIF
1231  ELSE
1232  IF ((x.LT.0.d0).AND.(y.LT.-x*sqrt3)) THEN
1233 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1234 CCC 2 PI/3 < PHASE(Z) < PI
1235 CCC BI(Z)=EXP(I PI/6) AI_(-1)(Z) + EXP(-I PI/6) AI_(1)(Z)
1236 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1237  c1=-0.5d0
1238  s1=0.5d0*sqrt3
1239  u=x*c1-y*s1
1240  v=x*s1+y*c1
1241  CALL aiz(ifun,ifac,u,v,ar,ai,ierr)
1242  IF (ifun.EQ.1) THEN
1243  br=c*ar-s*ai
1244  bi=c*ai+s*ar
1245  ELSE
1246  u=ar*c1-ai*s1
1247  v=ar*s1+ai*c1
1248  apr=u
1249  api=v
1250  bpr=c*apr-s*api
1251  bpi=c*api+s*apr
1252  ENDIF
1253 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1254 CCC WE NEED ALSO AI_(1)(Z)
1255 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1256  c1=-0.5d0
1257  s1=-0.5d0*sqrt3
1258  u=x*c1-y*s1
1259  v=x*s1+y*c1
1260  CALL aiz(ifun,ifac,u,v,ar,ai,ierr)
1261  IF (ifun.EQ.1) THEN
1262  s=-s
1263  br=br+c*ar-s*ai
1264  bi=bi+c*ai+s*ar
1265  ELSE
1266  u=ar*c1-ai*s1
1267  v=ar*s1+ai*c1
1268  apr=u
1269  api=v
1270  s=-s
1271  bpr=bpr+c*apr-s*api
1272  bpi=bpi+c*api+s*apr
1273  ENDIF
1274  ELSE
1275 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1276 CCC BI(Z) = I AI(Z) + 2 EXP(-I PI/6) AI_(1)(Z)
1277 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1278  c1=-0.5d0
1279  s1=-0.5d0*sqrt3
1280  u=x*c1-y*s1
1281  v=x*s1+y*c1
1282  CALL aiz(ifun,ifac,u,v,ar,ai,ierr)
1283  IF (ifun.EQ.1) THEN
1284  br=sqrt3*ar+ai
1285  bi=-ar+sqrt3*ai
1286  ELSE
1287  u=ar*c1-ai*s1
1288  v=ar*s1+ai*c1
1289  apr=u
1290  api=v
1291  bpr=sqrt3*apr+api
1292  bpi=-apr+sqrt3*api
1293  ENDIF
1294  CALL aiz(ifun,ifac,x,y,ar,ai,ierr)
1295  IF (ifun.EQ.1) THEN
1296  br=br-ai
1297  bi=bi+ar
1298  ELSE
1299  bpr=bpr-ai
1300  bpi=bpi+ar
1301  ENDIF
1302  ENDIF
1303  ENDIF
1304  ELSE
1305 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1306 CCC SCALED BI AIRY FUNCTIONS C
1307 CCC WE USE THE FOLLOWING NORMALIZATION: C
1308 CCC LET ARGZ=ARG(Z), THEN: C
1309 CCC A) IF 0 <= ARGZ <= PI/3 C
1310 CCC BI=EXP(-2/3Z^3/2)BI C
1311 CCC B) IF PI/3 <= ARGZ <= PI C
1312 CCC BI=EXP(2/3Z^3/2)BI C
1313 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1314  thet=phase(x,y)
1315  IF (thet.LE.pi3) THEN
1316  c1=-0.5d0
1317  s1=-0.5d0*sqrt3
1318  u=x*c1-y*s1
1319  v=x*s1+y*c1
1320  CALL aiz(ifun,ifac,u,v,ar,ai,ierr)
1321  IF (ifun.EQ.1) THEN
1322  br=sqrt3*ar+ai
1323  bi=-ar+sqrt3*ai
1324  ELSE
1325  u=ar*c1-ai*s1
1326  v=ar*s1+ai*c1
1327  apr=u
1328  api=v
1329  bpr=sqrt3*apr+api
1330  bpi=-apr+sqrt3*api
1331  ENDIF
1332  IF (iexpf.EQ.0) THEN
1333  r=sqrt(x*x+y*y)
1334  r32=r*sqrt(r)
1335  thet32=thet*1.5d0
1336  a1=cos(thet32)
1337  b1=sin(thet32)
1338  df1=4.d0/3.d0*r32
1339  expo=exp(df1*a1)
1340  expoi=1.d0/expo
1341  etar=expo*cos(df1*b1)
1342  etai=expo*sin(df1*b1)
1343  etagr=expoi*cos(-df1*b1)
1344  etagi=expoi*sin(-df1*b1)
1345  CALL aiz(ifun,ifac,x,y,ar,ai,ierr)
1346  IF (ifun.EQ.1) THEN
1347  br=br-ar*etagi-etagr*ai
1348  bi=bi+ar*etagr-etagi*ai
1349  ELSE
1350  bpr=bpr-ar*etagi-etagr*ai
1351  bpi=bpi+ar*etagr-etagi*ai
1352  ENDIF
1353  ENDIF
1354  ENDIF
1355  IF ((thet.GT.pi3).AND.(thet.LE.pi23)) THEN
1356  IF (iexpf.EQ.0) THEN
1357  r=sqrt(x*x+y*y)
1358  r32=r*sqrt(r)
1359  thet32=thet*1.5d0
1360  a1=cos(thet32)
1361  b1=sin(thet32)
1362  df1=4.d0/3.d0*r32
1363  expo=exp(df1*a1)
1364  expoi=1.d0/expo
1365  etar=expo*cos(df1*b1)
1366  etai=expo*sin(df1*b1)
1367  etagr=expoi*cos(-df1*b1)
1368  etagi=expoi*sin(-df1*b1)
1369  c1=-0.5d0
1370  s1=-0.5d0*sqrt3
1371  u=x*c1-y*s1
1372  v=x*s1+y*c1
1373  CALL aiz(ifun,ifac,u,v,ar,ai,ierr)
1374  IF (ifun.EQ.1) THEN
1375  bbr=sqrt3*ar+ai
1376  bbi=-ar+sqrt3*ai
1377  br=bbr*etar-bbi*etai
1378  bi=bbi*etar+bbr*etai
1379  ELSE
1380  u=ar*c1-ai*s1
1381  v=ar*s1+ai*c1
1382  apr=u
1383  api=v
1384  bbpr=sqrt3*apr+api
1385  bbpi=-apr+sqrt3*api
1386  bpr=bbpr*etar-bbpi*etai
1387  bpi=bbpi*etar+bbpr*etai
1388  ENDIF
1389  ELSE
1390  IF (ifun.EQ.1) THEN
1391  br=0.d0
1392  bi=0.d0
1393  ELSE
1394  bpr=0.d0
1395  bpi=0.d0
1396  ENDIF
1397  ENDIF
1398  CALL aiz(ifun,ifac,x,y,ar,ai,ierr)
1399  IF (ifun.EQ.1) THEN
1400  br=br-ai
1401  bi=bi+ar
1402  ELSE
1403  bpr=bpr-ai
1404  bpi=bpi+ar
1405  ENDIF
1406  ENDIF
1407  IF (thet.GT.pi23) THEN
1408  c1=-0.5d0
1409  s1=0.5d0*sqrt3
1410  u=x*c1-y*s1
1411  v=x*s1+y*c1
1412  CALL aiz(ifun,ifac,u,v,ar,ai,ierr)
1413  IF (ifun.EQ.1) THEN
1414  br=c*ar-s*ai
1415  bi=c*ai+s*ar
1416  ELSE
1417  u=ar*c1-ai*s1
1418  v=ar*s1+ai*c1
1419  apr=u
1420  api=v
1421  bpr=c*apr-s*api
1422  bpi=c*api+s*apr
1423  ENDIF
1424  IF (iexpf.EQ.0) THEN
1425 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1426 CCC WE NEED ALSO AI_(1)(Z)
1427 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1428  r=sqrt(x*x+y*y)
1429  r32=r*sqrt(r)
1430  thet32=thet*1.5d0
1431  a1=cos(thet32)
1432  b1=sin(thet32)
1433  df1=4.d0/3.d0*r32
1434  expo=exp(df1*a1)
1435  expoi=1.d0/expo
1436  etar=expo*cos(df1*b1)
1437  etai=expo*sin(df1*b1)
1438  etagr=expoi*cos(-df1*b1)
1439  etagi=expoi*sin(-df1*b1)
1440  c1=-0.5d0
1441  s1=-0.5d0*sqrt3
1442  u=x*c1-y*s1
1443  v=x*s1+y*c1
1444  CALL aiz(ifun,ifac,u,v,ar,ai,ierr)
1445  IF (ifun.EQ.1) THEN
1446  s=-s
1447  bbr=c*ar-s*ai
1448  bbi=c*ai+s*ar
1449  br=br+etar*bbr-etai*bbi
1450  bi=bi+bbi*etar+etai*bbr
1451  ELSE
1452  u=ar*c1-ai*s1
1453  v=ar*s1+ai*c1
1454  apr=u
1455  api=v
1456  s=-s
1457  bbpr=c*apr-s*api
1458  bbpi=c*api+s*apr
1459  bpr=bpr+etar*bbpr-etai*bbpi
1460  bpi=bpi+bbpi*etar+etai*bbpr
1461  ENDIF
1462  ENDIF
1463  ENDIF
1464  ENDIF
1465  IF (y0.LT.0) THEN
1466  bi=-bi
1467  bpi=-bpi
1468  ENDIF
1469  IF (ifun.EQ.1) THEN
1470  gbir=br
1471  gbii=bi
1472  ELSE
1473  gbir=bpr
1474  gbii=bpi
1475  ENDIF
1476  ENDIF
1477  RETURN
1478  END
1479 
1480 
1481 
1482 
1483 
1484