Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
src.f
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file src.f
1  SUBROUTINE dkia(IFAC,X,A,DKI,DKID,IERRO)
2 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
3 cccc the purpouse of the routine is the calculation of the
4 cccc modified bessel functions k_ia(x) and k'_IA(X),
5 CCCC WHERE I IS THE IMAGINARY UNIT AND X IS THE ARGUMENT OF THE
6 CCCC FUNCTIONS. WE WILL REFER TO A AS THE ORDER OF THE FUNCTIONS.
7 CCCC
8 CCCC THE ROUTINE HAS THE OPTION OF COMPUTING SCALED FUNCTIONS.
9 CCCC THIS SCALING CAN BE USED TO ENLARGE THE RANGE OF
10 CCCC COMPUTATION.
11 CCCC THE SCALED FUNCTIONS ARE DEFINED AS FOLLOWS
12 CCCC (S STANDS FOR SCALED AND
13 CCCC L=SQRT{X**2-A**2} + A*ARCSIN(A/X)):
14 CCCC
15 CCCC EXP(L)*K_IA(X), IF X >=ABS(A)
16 CCCC SK_IA(X) =
17 CCCC EXP(ABS(A)*PI/2)*K_IA(X), IF X < ABS(A)
18 CCCC
19 CCCC
20 CCCC EXP(L)*K'_ia(x), IF x >=abs(a)
21 cccc sk'_IA(X) =
22 CCCC EXP(ABS(A)*PI/2)*K'_ia(x), IF x < abs(a)
23 cccc
24 cccc the range of the parameters(x,a) for the computation of
25 cccc scaled functions is:
26 cccc 0 < x <= 1500, -1500 <= a <= 1500.
27 cccc for functions without scaling, the range is usually larger
28 cccc than
29 cccc 0 < x <= 500, -400 <= a <= 400,
30 cccc depending on the machine overflow/underflow parameters, which
31 cccc are set up by the routine d1mach.
32 cccc
33 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
34 cccc methods of computation:
35 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
36 cccc 1) series, IF abs(a) > 0.044*abs(x-3.1)**1.9+(x-3.1)
37 cccc 2) continued fraction,
38 cccc IF x>3 and abs(a) < 380*(abs((x-3)/2300))**0.572
39 cccc 3) airy-TYPE asymptotic expansion,
40 cccc IF abs(a) > 0.4*x+7.5 and abs(a) < 3.7*x-10
41 cccc 4) quadratures,
42 cccc in the rest of cases
43 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
44 cccc description of input/output variables:
45 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
46 cccc inputs:
47 cccc x : argument of the functions
48 ccc a : order of the functions
49 cccc ifac : INTEGER flag for the scaling
50 cccc * IF ifac=1, the code computes kia(x), kia'(X)
51 CCCC * OTHERWISE, THE CODE COMPUTES SCALED KIA(X), KIA'(x)
52 cccc outputs
53 cccc dki : kia(x) function
54 cccc dkid : derivative of the kia(x) FUNCTION
55 cccc ierro: error flag
56 cccc * IF ierro=0, computation successful.
57 cccc * IF ierro=1, computation out of range.
58 cccc * IF ierro=2, variables x and/or a, out of range.
59 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
60 cccc accuracy:
61 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
62 cccc the relative accuracy is:
63 cccc better than 10**(-13) for(x,a) in(0,200]x[-200,200];
64 cccc better than 5.10**(-13) for(x,a) in(0,500]x[-500,500];
65 cccc CLOSE to 10**(-12) for(x,a) in(0,1500]x[-1500,1500].
66 cccc near the zeros of the functions(there are infinitely
67 cccc many of them in the oscillatory region) relative precision
68 cccc looses meaning and only absolute precision makes sense.
69 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
70 c authors:
71 c amparo gil(u. cantabria, santander, spain).
72 c e-mail: amparo.gil@unican.es
73 c javier segura(u. cantabria, santander, spain).
74 c e-mail: segurajj@unican.es
75 c nico m. temme(cwi, amsterdam, the netherlands).
76 c e-mail: nico.temme@cwi.nl
77 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
78 c references:
79 c this is the companion software of the articles
80 c 1)'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL DIFFERENTIAL
81 C EQUATION FOR IMAGINARY ORDERS AND POSITIVE ARGUMENTS',
82 c a. gil, j. segura, n.m. temme
83 c acm trans. math. soft. (2004)
84 c 2)'MODIFIED BESSEL FUNCTIONS OF IMAGINARY ORDER AND
85 C POSITIVE ARGUMENT',
86 c a. gil, j. segura, n.m. temme
87 c acm trans. math. soft. (2004)
88 c additional references:
89 c - 'COMPUTATION OF THE MODIFIED BESSEL FUNCTION OF THE
90 C THIRD KIND FOR IMAGINARY ORDERS'
91 c a. gil, j. segura, n.m. temme
92 c j. comput. phys. 175 (2002) 398-411
93 c - 'COMPUTATION OF THE MODIFIED BESSEL FUNCTIONS OF THE
94 C THIRD KIND OF IMAGINARY ORDERS:
95 C UNIFORM AIRY-TYPE ASYMPTOTIC EXPANSION'
96 c a. gil, j. segura, n.m. temme, proceedings opsfa 2001
97 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
98 ccc local variables:
99 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
100 c d: x-3.1
101 c df1: 0.044*abs(d)**1.9+(x-3.1)
102 c df2: 380*(abs((x-3)/2300))**0.572
103 c df3: 0.4*x+7.5
104 c df4: 3.7*x-10
105 c pnu: order of the FUNCTION
106 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
107  DOUBLE PRECISION a,d,df1,df2,df3,df4,dki,
108  + dkid,pnu,x
109  INTEGER ierro,ifac
110  ierro=0
111  pnu=a
112  IF ((x.GT.1500.d0).OR.(x.LE.0.d0)) THEN
113  ierro=2
114  dki=0.d0
115  dkid=0.d0
116  ENDIF
117  IF (abs(pnu).GT.1500.d0) THEN
118  ierro=2
119  dki=0.d0
120  dkid=0.d0
121  ENDIF
122  IF (ierro.EQ.0) THEN
123 cccccccccccccccccccccccccccccccccccccccccccccccc
124 ccc these functions are even functions in the c
125 ccc PARAMETER a(=pnu) c
126 cccccccccccccccccccccccccccccccccccccccccccccccc
127  IF (pnu.LT.0.d0) pnu=-pnu
128  d=x-3.1d0
129  df1=0.044d0*abs(d)**1.9d0+d
130  df2=380.d0*(abs((x-3.d0)/2300.d0))**0.572d0
131  df3=0.4d0*x+7.5d0
132  df4=3.7d0*x-10.d0
133  IF (pnu.GT.df1) THEN
134  CALL serkia(ifac,x,pnu,dki,dkid,ierro)
135  ELSEIF ((x.GT.3.d0).AND.(pnu.LT.df2)) THEN
136  CALL frakia(ifac,x,pnu,dki,dkid,ierro)
137  ELSEIF ((pnu.GT.df3).AND.(pnu.LT.df4)) THEN
138  CALL aiexki(ifac,x,pnu,dki,dkid,ierro)
139  ELSE
140  CALL dkint(ifac,x,pnu,dki,dkid,ierro)
141  ENDIF
142  ENDIF
143  RETURN
144  END
145 
146  SUBROUTINE dkint(IFAC,XX,PNUA,DKINF,DKIND,IERRO)
147 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
148 ccc calculation of the functions k,k' USING NON-OSCILLATING C
149 CCC INTEGRAL REPRESENTATIONS C
150 CCC C
151 CCC INPUTS: C
152 CCC XX: ARGUMENT OF THE FUNCTIONS C
153 CCC PNUA: ORDER OF THE FUNCTIONS C
154 CCC IFAC: C
155 CCC =1, NON SCALED FUNCTIONS C
156 CCC OTYHERWISE, SCALED FUNCTIONS C
157 CCC OUTPUTS: C
158 CCC DKINF: K FUNCTION C
159 CCC DKIND: DERIVATIVE OF THE K FUNCTION C
160 CCC IERRO: ERROR FLAG C
161 CCC * IF IERRO=0, COMPUTATION SUCCESSFUL. C
162 CCC * IF IERRO=1, COMPUTATION OUT OF RANGE. C
163 CCC * IF IERRO=2, ARGUMENT AND/OR ORDER, OUT OF RANGE. C
164 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
165 CCC METHOD OF COMPUTATION:
166 CCC * IF XX>=PNUA, THE NON-OSCILLATING INTEGRAL REPRESENTATIONS
167 CCC FOR THE MONOTONIC REGION ARE USED
168 CCC * IF XX<PNUA, THE NON-OSCILLATING INTEGRAL REPRESENTATIONS
169 CCC FOR THE OSCILLATORY REGION ARE USED
170 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
171 CCC LOCAL VARIABLES:
172 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
173 C CHI : X*SINH(MU)-PNU*MU
174 C CONTR1 : CONTRIBUTION OF THE SEMI-INFINITE INTEGRAL
175 C IN THE OSCILLATORY CASE (INCLUDING ADDITIONAL
176 C FACTORS: CONTR1=FOS1*FACTORS).
177 C COSCHI : COS(CHI)
178 C COSH2M : COSH(2*MU)
179 C COSHM : COSH(MU)
180 C COSTH : COS(THET), THET=ASIN(PNU/X)
181 C DF1 : FACTOR (DEPENDING ON IFAC)
182 C DMAIN : PI*PNU*0.5
183 C DMU : SOLUTION MU OF COSH(MU)=PNU/X
184 C DMU2 : 2*MU
185 C DMU3 : 3*MU
186 C DMU5 : 5*MU
187 C DMU7 : 7*MU
188 C DMU9 : 9*MU
189 C DMUFAC : MU*COSH(MU)-SINH(MU)
190 C DMUTAN : MU-TANH(MU)
191 C FDOMIN : X*(COS(THET)+THET*SIN(THET))
192 C FOS1 : CONTRIBUTION OF THE SEMI-INFINITE INTEGRAL
193 C IN THE OSCILLATORY CASE (KIA(X)).
194 C FOSD1 : CONTRIBUTION OF THE SEMI-INFINITE INTEGRAL
195 C IN THE OSCILLATORY CASE (KIA'(x)).
196 c hir : monotonic CASE, contribution of the integral
197 c(kia(x)).
198 c hird : monotonic CASE, contribution of the integral
199 c(kia'(X)).
200 C PI : 3.1415...
201 C PINU : PI*PNU
202 C PNU : ORDER OF THE FUNCTION
203 C SINCHI : SIN(CHI)
204 C SINH2M : SINH(2*MU)
205 C SINHM : SINH(MU)
206 C SINTH : SIN(THET)
207 C THET : ASIN(PNU/X)
208 C UNDER : UNDERFLOW NUMBER
209 C X : ARGUMENT OF THE FUNCTION
210 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
211  DOUBLE PRECISION CHI,CONTR1,COSCHI,COSH2M,COSHM,COSTH,
212  + D1MACH,DF1,DKIND,DKINF,DMAIN,DMU,DMU2,DMU3,DMU5,DMU7,
213  + DMU9,DMUFAC,DMUTAN,FDOMIN,FOS1,FOSD1,HIR,HIRD,PI,
214  + PINU,PNU,PNUA,SINCHI,SINH2M,SINHM,SINTH,THET,UNDER,
215  + X,XX
216  INTEGER IERRO,IFAC
217  COMMON/ARGU/X,PNU
218  COMMON/PARMON/THET,SINTH,COSTH,FDOMIN
219  COMMON/PAROS1/DMU,COSHM,SINHM,DMUFAC,COSCHI,SINCHI
220  COMMON/PAROS2/COSH2M,SINH2M,DMAIN
221  COMMON/PAROS3/DMUTAN
222 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
223  UNDER=D1MACH(1)*1.D+8
224  X=XX
225  PNU=PNUA
226  IERRO=0
227 .GE. IF (XPNU) THEN
228 CCCCCCCCCCCCCCCCCCCCCCCCCC
229 CCC MONOTONIC REGION CC
230 CCCCCCCCCCCCCCCCCCCCCCCCCC
231  SINTH=PNU/X
232  THET=ASIN(SINTH)
233  COSTH=COS(THET)
234  FDOMIN=X*(COSTH+THET*SINTH)
235 .EQ. IF (IFAC1) THEN
236 .LE. IF (-FDOMINLOG(UNDER)) IERRO=1
237  ENDIF
238 .EQ. IF (IERRO0) THEN
239 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
240 CCC CALCULATION OF KIA:
241 CCC CALL TO THE TRAPEZOIDAL ROUTINE C
242 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
243  CALL TRAPRE(1,HIR)
244 .EQ. IF (IFAC1) THEN
245  DKINF=0.5D0*HIR*EXP(-FDOMIN)
246  ELSE
247  DKINF=0.5D0*HIR
248  ENDIF
249 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
250 CCC CALCULATION OF KIA':
251 ccc CALL to the trapezoidal routine c
252 ccccccccccccccccccccccccccccccccccccccccc
253  CALL trapre(10,hird)
254  IF (ifac.EQ.1) THEN
255  dkind=-hird*exp(-fdomin)
256  ELSE
257  dkind=-hird
258  ENDIF
259  ELSE
260  dkinf=0.d0
261  dkind=0.d0
262  ENDIF
263  ELSE
264 cccccccccccccccccccccccccccc
265 ccc oscillatory region cc
266 cccccccccccccccccccccccccccc
267  pi=acos(-1.d0)
268  IF (ifac.EQ.1) THEN
269  IF ((-pi*pnu*0.5d0).LE.log(under)) ierro=1
270  ENDIF
271  IF (ierro.EQ.0) THEN
272  coshm=pnu/x
273  dmu=log(coshm+sqrt((coshm-1.d0)*(coshm+1.d0)))
274  dmu2=2.d0*dmu
275  cosh2m=cosh(dmu2)
276  sinhm=sinh(dmu)
277  sinh2m=sinh(dmu2)
278  IF (dmu.GT.0.1d0) THEN
279  chi=x*sinhm-pnu*dmu
280  dmufac=dmu*coshm-sinhm
281  ELSE
282  dmu2=dmu*dmu
283  dmu3=dmu2*dmu
284  dmu5=dmu3*dmu2
285  dmu7=dmu5*dmu2
286  dmu9=dmu7*dmu2
287  chi=-2.d0*x*(1.d0/6.d0*dmu3+1.d0/60.d0*dmu5+
288  + 3.d0/5040.d0*dmu7+4.d0/362880.d0*dmu9)
289  dmufac=dmu3/3.d0+dmu5/30.d0+dmu7/840.d0+dmu9/45360.d0
290  ENDIF
291  dmutan=dmu-tanh(dmu)
292  coschi=cos(chi)
293  sinchi=sin(chi)
294  pinu=pi*pnu
295  dmain=pinu*0.5d0
296 cccccccccccccccccc
297 ccc integrals cc
298 cccccccccccccccccc
299  CALL trapre(2,fos1)
300 ccccccccccccccccccccccccccccccccccccccccccccc
301 ccc THEN, the kia FUNCTION is given by ... cc
302 ccccccccccccccccccccccccccccccccccccccccccccc
303  IF (ifac.EQ.1) THEN
304  df1=exp(-dmain)
305  contr1=df1*fos1
306  ELSE
307  df1=1.d0
308  contr1=df1*fos1
309  ENDIF
310  dkinf=contr1
311 cccccccccccccccccccccccccc
312 ccc calculation of kia' CC
313 CCCCCCCCCCCCCCCCCCCCCCCCCC
314 CCCCCCCCCCCCCCCCCC
315 CCC INTEGRALS CC
316 CCCCCCCCCCCCCCCCCC
317  CALL TRAPRE(20,FOSD1)
318 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
319 CCC THEN, KIA' is given by ...
320 cccccccccccccccccccccccccccccccccc
321  IF (ifac.EQ.1) THEN
322  df1=exp(-dmain)
323  contr1=df1*fosd1
324  ELSE
325  df1=1.d0
326  contr1=df1*fosd1
327  ENDIF
328  dkind=contr1
329  ELSE
330  dkinf=0.d0
331  dkind=0.d0
332  ENDIF
333  ENDIF
334  RETURN
335  END
336 
337  DOUBLE PRECISION FUNCTION fa(U)
338 cccccccccccccccccccccccccccccccccccccccccccccccc
339 ccc computation of the integrand for the
340 ccc integral representation of the k FUNCTION
341 ccc in the monotonic region.
342 cccccccccccccccccccccccccccccccccccccccccccccccc
343 ccc input variable:
344 cccccccccccccccccccccccccccccccccccccccccccccccc
345 c u : argument of the FUNCTION
346 c(variable of integration).
347 cccccccccccccccccccccccccccccccccccccccccccccccc
348 ccc local variables:
349 cccccccccccccccccccccccccccccccccccccccccccccccc
350 c phir : x*cosh(u)*cos(v(u))+pnu*v(u)
351 c -fdomin, WHERE
352 c v(u)=asin(u/sinh(u)*sin(thet))
353 c and fdomin=x*(cos(thet)+
354 c thet*sin(thet))
355 c under : underflow number
356 cccccccccccccccccccccccccccccccccccccccccccccccc
357  DOUBLE PRECISION d1mach,phir,u,under
358 ccc machine dependent constant(underflow number)
359  under=d1mach(1)*1.d+8
360  IF ((-phir(u)).LE.log(under)) THEN
361  fa=0.d0
362  ELSE
363  fa=exp(-phir(u))
364  ENDIF
365  RETURN
366  END
367 
368  DOUBLE PRECISION FUNCTION fad(T)
369 cccccccccccccccccccccccccccccccccccccccccccccccccc
370 ccc computation of the integrand for the
371 ccc integral representation of the derivative of
372 ccc the k FUNCTION in the monotonic region.
373 cccccccccccccccccccccccccccccccccccccccccccccccccc
374 ccc input variable:
375 cccccccccccccccccccccccccccccccccccccccccccccccccc
376 c t : argument of the FUNCTION
377 c(variable of integration).
378 cccccccccccccccccccccccccccccccccccccccccccccccccc
379 ccc local variables:
380 cccccccccccccccccccccccccccccccccccccccccccccccccc
381 c ang : thet-v(u)
382 c angh : ang/2
383 c costh : cos(thet)
384 c cosvu : cos(v(u))
385 c dfac : cos(thet)+(fu1+2.d0*sin(ang/2)
386 c *sin(ang/2))/cos(vu)
387 c djaco : cosh(t)*exp(s)/(1.d0+exp(s))
388 c expo : exp(-phir(u))
389 c fdomin : x*(cos(thet)+thet*sin(thet))
390 c fu1 : 2*sinh(u/2)**2
391 c fuac : u**3/6+u**5/120+u**7/5040
392 c phir(u) : x*cosh(u)*cos(v(u))+pnu*v(u)-fdomin
393 c s : sinh(t)
394 c sinanh : sin(ang/2)
395 c sinhu : sinh(u)
396 c sinhuh : sinh(u/2)
397 c sinth : sin(thet)
398 c thet : asin(pnu/x)
399 c u : log(1+exp(s))
400 c u2 : u**2
401 c u3 : u**3
402 c u5 : u**5
403 c u7 : u**7
404 c uh : u/2
405 c under : underflow number
406 c v(u) : asin(u/sinh(u)*sin(thet))
407 c vu : v(u)
408 c y : exp(s)
409 cccccccccccccccccccccccccccccccccccccccccccccccccccc
410  DOUBLE PRECISION ang,angh,costh,cosvu,d1mach,
411  + dfac,djaco,expo,fdomin,fu1,fuac,phir,s,sinanh,
412  + sinhu,sinhuh,sinth,t,thet,u,u2,u3,u5,u7,uh,under,
413  + v,vu,y
414  common/parmon/thet,sinth,costh,fdomin
415 ccc machine dependent constant(underflow number)
416  under=d1mach(1)*1.d+8
417  s=sinh(t)
418  y=exp(s)
419  u=log(1.d0+y)
420  djaco=cosh(t)*y/(1.d0+y)
421  IF (abs(u).LT.1.d-1) THEN
422  IF (abs(u).LT.under) THEN
423  dfac=costh
424  ELSE
425  uh=u*0.5d0
426  sinhuh=sinh(uh)
427  fu1=2.d0*sinhuh*sinhuh
428  u2=u*u
429  u3=u2*u
430  u5=u2*u3
431  u7=u5*u2
432  fuac=u3/6.d0+u5/120.d0+u7/5040.d0
433  sinhu=sinh(u)
434  vu=v(u)
435  cosvu=cos(vu)
436  ang=asin(-sinth/(costh*u/sinhu+cosvu)*
437  + (sinhu+u)*fuac/(sinhu*sinhu))
438  angh=0.5d0*ang
439  sinanh=sin(angh)
440  dfac=costh+(fu1+2.d0*sinanh*sinanh)/cosvu
441  ENDIF
442  ELSE
443  vu=v(u)
444  ang=thet-vu
445  fu1=cosh(u)-1.d0
446  angh=0.5d0*ang
447  sinanh=sin(angh)
448  dfac=costh+(fu1+2.d0*sinanh*sinanh)/cos(vu)
449  ENDIF
450  IF ((-phir(u)).LE.log(under)) THEN
451  expo=0.d0
452  ELSE
453  expo=exp(-phir(u))
454  ENDIF
455  fad=expo*dfac*djaco
456  RETURN
457  END
458 
459  DOUBLE PRECISION FUNCTION fdtau2(T)
460 cccccccccccccccccccccccccccccccccccccccccccccccccccc
461 ccc computation of the integrand for the
462 ccc integral for the derivative of the k function
463 ccc in eq.(37) of the reference:
464 ccc 'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL
465 CCC DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS
466 CCC AND POSITIVE ARGUMENTS',
467 ccc a. gil, j. segura, n.m. temme
468 ccc acm trans. math. soft. (2004)
469 cccccccccccccccccccccccccccccccccccccccccccccccccccc
470 ccc input variable:
471 cccccccccccccccccccccccccccccccccccccccccccccccccccc
472 c t : argument of the FUNCTION
473 c(variable of integration).
474 cccccccccccccccccccccccccccccccccccccccccccccccccccc
475 ccc local variables:
476 cccccccccccccccccccccccccccccccccccccccccccccccccccc
477 c argu : (cosh(mu)*u-dmufac)/sinh(u)
478 c argu2 : argu**2
479 c coschi : cosh(chi), chi=x*sinh(mu)-pnu*mu
480 c cosh2m : cosh(2*mu)
481 c coshm : cosh(mu)
482 c coshu : cosh(u)
483 c coss : cos(sigma(u))
484 c d1 : cosh(mu)-dmufac
485 c delta : -sin(chi)*cos(sigma(u))*cosh(u)
486 c +cos(chi)*sin(sigma(u))*sinh(u)
487 c deri : (cosh(mu)/sinh(u)-d1*cosh(u)/sinh(u)**2)
488 c /sqrt(1-argu2)
489 c djaco : cosh(t)*exp(s)/sqrt(1+exp(s)**2)
490 c dmain : pi*pnu*0.5
491 c dmu : solution mu of cosh(mu)=pnu/x
492 c dmu2 : 2*mu
493 c dmufac : mu*cosh(mu)-sinh(mu)
494 c expon : exp(-(phib(u)-pi*pnu*0.5))
495 c f1 : sinh(u-mu)/(u-mu)
496 c g1 : z/6+z3/120+z5/5040+z7/362880, z=2*y
497 c gamma : sin(chi)*sin(sigma(u))*sinh(u)+
498 c cos(chi)*cos(sigma(u))*cosh(u)
499 c phib(u) : x*cosh(u)*cos(sigma(u))+pnu*sigma(u),
500 c WHERE x=argument of the function,
501 c sigma(u)=arcsin((cosh(mu)*u-dmufac)/sinh(u)
502 c resto : -sin(chi)*cos(sigma(u))*cosh(u)
503 c +cos(chi)*sin(sigma(u))*sinh(u)+
504 c(sin(chi)*sin(sigma(u))*sinh(u)+
505 c cos(chi)*cos(sigma(u))*cosh(u))*deri
506 c s : sinh(t)
507 c sigma(u): asin((cosh(mu)*u-dmufac)/sinh(u))
508 c sigmau : sigma(u)
509 c sinchi : sin(chi)
510 c sinh2m : sinh(2*mu)
511 c sinhm : sinh(mu)
512 c sinhu : sinh(u)
513 c sinhu2 : 2*sinh(u)
514 c sins : sin(sigma(u))
515 c u : mu+log(x+sqrt(x**2+1))
516 c under : underflow number
517 c x : exp(s)
518 c y : u-mu
519 c z : 2*y
520 c z2 : z**2
521 c z3 : z**3
522 c z5 : z**5
523 c z7 : z**7
524 cccccccccccccccccccccccccccccccccccccccccccccccccccc
525  DOUBLE PRECISION argu,argu2,coschi,cosh2m,coshm,
526  + coshu,coss,d1,d1mach,delta,deri,djaco,dmain,dmu,
527  + dmu2,dmufac,dmutan,expon,f1,g1,gamma,phib,resto,
528  + s,sigma,sigmau,sinchi,sinh2m,sinhm,sinhu,sinhu2,
529  + sins,t,u,under,x,y,z,z2,z3,z5,z7
530  common/paros1/dmu,coshm,sinhm,dmufac,coschi,sinchi
531  common/paros2/cosh2m,sinh2m,dmain
532  common/paros3/dmutan
533 ccc machine dependent constant(underflow number)
534  under=d1mach(1)*1.d+8
535  s=sinh(t)
536  x=exp(s)
537  u=dmutan+log(x+sqrt(x*x+1.d0))
538  y=u-dmu
539  coshu=cosh(u)
540  sinhu=sinh(u)
541  IF (abs(y).LE.1.d-1) THEN
542  IF (abs(y).GT.under) THEN
543  f1=sinh(y)/y
544  ELSE
545  f1=1.d0
546  ENDIF
547  z=y*2.d0
548  z2=z*z
549  z3=z2*z
550  z5=z3*z2
551  z7=z5*z2
552  g1=z/6.d0+z3/120.d0+z5/5040.d0+z7/362880.d0
553  dmu2=2.d0*dmu
554  deri=sinhu/(f1-coshm*coshu)*sqrt(cosh(dmu2)*f1*f1+
555  + 2.d0*sinh(dmu2)*g1-coshm*coshm)
556  deri=1.d0/deri
557  ELSE
558  d1=coshm*u-dmufac
559  argu=d1/sinhu
560  argu2=argu*argu
561  sinhu2=sinhu*sinhu
562  deri=1.d0/sqrt(1.d0-argu2)*(coshm/sinhu-
563  + d1*coshu/sinhu2)
564  IF (u.LT.dmu) deri=-deri
565  ENDIF
566  djaco=cosh(t)*x/sqrt(1.d0+x*x)
567  IF ((-(phib(u)-dmain)).LE.log(under)) THEN
568  expon=0.d0
569  fdtau2=0.d0
570  ELSE
571  expon=exp(-(phib(u)-dmain))
572  sigmau=sigma(u)
573  sins=sin(sigmau)
574  coss=cos(sigmau)
575  gamma=coschi*sins*sinhu-sinchi*coss*coshu
576  delta=-coschi*coss*coshu-sinchi*sins*sinhu
577  resto=delta+gamma*deri
578  fdtau2=expon*resto*djaco
579  ENDIF
580  RETURN
581  END
582 
583  DOUBLE PRECISION FUNCTION fstau2(T)
584 ccccccccccccccccccccccccccccccccccccccccccccccccc
585 ccc computation of the integrand for the
586 ccc integral for the k FUNCTION in eq.(37) OF
587 ccc the reference:
588 ccc 'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL
589 CCC DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS
590 CCC AND POSITIVE ARGUMENTS',
591 ccc a. gil, j. segura, n.m. temme
592 ccc acm trans. math. soft. (2004)
593 cccccccccccccccccccccccccccccccccccccccccccccccccccc
594 ccc input variable:
595 cccccccccccccccccccccccccccccccccccccccccccccccccccc
596 c t : argument of the FUNCTION
597 c(variable of integration).
598 ccccccccccccccccccccccccccccccccccccccccccccccccc
599 ccc local variables:
600 ccccccccccccccccccccccccccccccccccccccccccccccccc
601 c argu : (cosh(mu)*u-dmufac)/sinh(u)
602 c argu2 : argu**2
603 c coschi : cosh(chi), chi=x*sinh(mu)-pnu*mu
604 c cosh2m : cosh(2*mu)
605 c coshm : cosh(mu)
606 c d1 : cosh(mu)-dmufac
607 c deri : (cosh(mu)/sinh(u)-d1*cosh(u)/sinh(u)**2)
608 c /sqrt(1-argu2)
609 c djaco : cosh(t)*exp(s)/sqrt(1+exp(s)**2)
610 c dmain : pi*pnu*0.5
611 c dmu : solution mu of cosh(mu)=pnu/x
612 c dmu2 : 2*mu
613 c dmufac : mu*cosh(mu)-sinh(mu)
614 c expon : exp(-(phib(u)-pi*pnu*0.5))
615 c f1 : sinh(u-mu)/(u-mu)
616 c g1 : z/6+z3/120+z5/5040+z7/362880, z=2*y
617 c phib(u) : x*cosh(u)*cos(sigma(u))+pnu*sigma(u),
618 c WHERE x=argument of the function,
619 c sigma(u)=arcsin((cosh(mu)*u-dmufac)/sinh(u)
620 c resto : cos(chi)+sin(chi)*deri
621 c s : sinh(t)
622 c sinchi : sin(chi)
623 c sinh2m : sinh(2*mu)
624 c sinhm : sinh(mu)
625 c sinhu : sinh(u)
626 c sinhu2 : 2*sinh(u)
627 c u : mu+log(x+sqrt(x**2+1))
628 c under : underflow number
629 c x : exp(s)
630 c y : u-mu
631 c z : 2*y
632 c z2 : z**2
633 c z3 : z**3
634 c z5 : z**5
635 c z7 : z**7
636 ccccccccccccccccccccccccccccccccccccccccccccccccc
637  DOUBLE PRECISION argu,argu2,coschi,cosh2m,
638  + coshm,d1,d1mach,deri,djaco,dmain,dmu,dmu2,
639  + dmufac,dmutan,expon,f1,g1,phib,resto,s,sinchi,
640  + sinh2m,sinhm,sinhu,sinhu2,t,u,under,x,y,z,z2,
641  + z3,z5,z7
642  common/paros1/dmu,coshm,sinhm,dmufac,coschi,sinchi
643  common/paros2/cosh2m,sinh2m,dmain
644  common/paros3/dmutan
645 ccc machine dependent constant(underflow number)
646  under=d1mach(1)*1.d+8
647  s=sinh(t)
648  x=exp(s)
649  u=dmutan+log(x+sqrt(x*x+1.d0))
650  y=u-dmu
651  IF (abs(y).GT.under) THEN
652  f1=sinh(y)/y
653  ELSE
654  f1=1.d0
655  ENDIF
656  z=y*2.d0
657  z2=z*z
658  IF (abs(y).LE.0.1d0) THEN
659  z3=z2*z
660  z5=z3*z2
661  z7=z5*z2
662  g1=z/6.d0+z3/120.d0+z5/5040.d0+z7/362880.d0
663  dmu2=2.d0*dmu
664  deri=sinh(u)/(f1-coshm*cosh(u))*
665  + sqrt(cosh(dmu2)*f1*f1+2.d0*sinh(dmu2)*g1-
666  + coshm*coshm)
667  deri=1.d0/deri
668  ELSE
669  d1=coshm*u-dmufac
670  sinhu=sinh(u)
671  argu=d1/sinhu
672  argu2=argu*argu
673  sinhu2=sinhu*sinhu
674  deri=1.d0/sqrt(1.d0-argu2)*
675  + (coshm/sinhu-d1*cosh(u)/sinhu2)
676  IF (u.LT.dmu) deri=-deri
677  ENDIF
678  djaco=cosh(t)*x/sqrt(1.d0+x*x)
679  IF ((-(phib(u)-dmain)).LE.log(under)) THEN
680  expon=0.d0
681  fstau2=0.d0
682  ELSE
683  expon=exp(-(phib(u)-dmain))
684  resto=(coschi+sinchi*deri)
685  fstau2=expon*resto*djaco
686  ENDIF
687  RETURN
688  END
689 
690  SUBROUTINE trapre(IC,TI)
691 cccccccccccccccccccccccccccccccccccccccccccccccccc
692 ccc implementation of an adaptive trapezoidal rule
693 ccc for computing the integral representations of
694 ccc the functions, for the different regions
695 ccc(monotonic or oscillatory).
696 ccc
697 ccc input variable:
698 ccc ic: depending on the values of ic,
699 ccc different integrals are computed:
700 ccc *ic=1, k function, monotonic region
701 ccc *ic=2, k function, oscillatory region
702 ccc *ic=10, derivative of the k function,
703 ccc monotonic region.
704 ccc *ic=20, derivative of the k function,
705 ccc oscillatory region.
706 ccc output variable:
707 ccc ti, result of the integral
708 cccccccccccccccccccccccccccccccccccccccccccccccccccc
709 ccc local variables
710 cccccccccccccccccccccccccccccccccccccccccccccccccccc
711 c a : lower integration limit
712 c b : upper integration limit
713 c delta : calculates the relative precision
714 c eps : relative precision PARAMETER used in
715 c the calculation
716 c h : integration step
717 c pnu : order of the function
718 c sum : accumulates the elementary
719 c contributions
720 c tin : evaluated integral
721 c x : argument
722 c xac : integration abcissa
723 c z : x/pnu
724 cccccccccccccccccccccccccccccccccccccccccccccccccccc
725  DOUBLE PRECISION a,b,d1mach,delta,eps,finte,h,
726  + pnu,sum,ti,tin,x,xac,z
727  INTEGER i,ic,ifi,n
728  common/argu/x,pnu
729  eps=d1mach(3)
730  IF (eps.LT.1.d-14) eps=1.d-14
731  n=0
732 ccccc integration limits: a,b
733  z=x/pnu
734  IF ((z.GT.0.999d0).AND.(z.LT.1.001d0)) THEN
735  IF ((ic.EQ.1).OR.(ic.EQ.10)) THEN
736  a=-3.5d0
737  ELSE
738  a=-4.5d0
739  ENDIF
740  ELSE
741  IF ((ic.EQ.2).OR.(ic.EQ.20)) THEN
742  a=-2.5d0
743  ELSE
744  a=-4.5d0
745  ENDIF
746  ENDIF
747  b=-a
748  h=b-a
749  ti=0.5d0*h*(finte(ic,a)+finte(ic,b))
750  delta=1.d0+eps
751 11 n=n+1
752  h=0.5d0*h
753  IF (n.EQ.1) THEN
754  ifi=1
755  ELSE
756  ifi=2*ifi
757  ENDIF
758  sum=0.d0
759  DO 3 i=1,ifi
760  xac=a+dble(2*i-1)*h
761  sum=sum+finte(ic,xac)
762  3 CONTINUE
763  tin=0.5d0*ti+h*sum
764  IF ((tin.NE.0.d0).AND.(n.GT.4)) THEN
765  delta=abs(1.d0-ti/tin)
766  ENDIF
767  ti=tin
768  IF ((delta.GT.eps).AND.(n.LT.9)) goto 11
769  RETURN
770  END
771 
772  DOUBLE PRECISION FUNCTION finte(IC,T)
773 ccccccccccccccccccccccccccccccccccccccccccccccccccc
774 ccc auxiliary FUNCTION for the subroutine trapre.
775 ccc this FUNCTION calls the different functions
776 ccc containing the integrands which are integrated
777 ccc by trapre.
778 ccccccccccccccccccccccccccccccccccccccccccccccccccc
779 ccc input: c
780 ccc ic: INTEGER parameter. its admissible c
781 ccc values are the same as in the c
782 ccc SUBROUTINE trapre. c
783 ccc t: integration abscissa c
784 ccccccccccccccccccccccccccccccccccccccccccccccccccc
785 ccc local variables:
786 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
787 c fa : monotonic part contribution
788 c(kia(x) function).
789 c fad : monotonic part contribution
790 c(kia'(X) FUNCTION).
791 C FDTAU2 : OSCILLATORY PART: TAU CONTRIBUTION ROUTINE
792 C (KIA'(x) function). semi-infinite integral.
793 c used for large pnu
794 c fstau2 : oscillatory part: tau contribution routine
795 c(kia(x) function). semi-infinite integral.
796 c used for large pnu
797 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
798  DOUBLE PRECISION fa,fad,fdtau2,fstau2,t
799  INTEGER ic
800  IF (ic.EQ.1) THEN
801  finte=fa(t)
802  ENDIF
803  IF (ic.EQ.2) THEN
804  finte=fstau2(t)
805  ENDIF
806  IF (ic.EQ.10) THEN
807  finte=fad(t)
808  ENDIF
809  IF (ic.EQ.20) THEN
810  finte=fdtau2(t)
811  ENDIF
812  END
813 
814  SUBROUTINE frakia(IFAC,X,PNU,PSER,PSERD,IERRO)
815 cccccccccccccccccccccccccccccccccccccccccccccccccccc
816 ccc implementation of the continued fraction c
817 ccc method for the calculation of k and k' C
818 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
819 CCC INPUT VARIABLES: C
820 CCC X: ARGUMENT OF THE FUNCTIONS C
821 CCC PNU: ORDER OF THE FUNCTIONS C
822 CCC IFAC: C
823 CCC =1, NON SCALED FUNCTIONS C
824 CCC OTHERWISE, SCALED FUNCTIONS C
825 CCC OUTPUT VARIABLES: C
826 CCC PSER: K FUNCTION C
827 CCC PSERD: DERIVATIVE OF THE K FUNCTION C
828 CCC IERRO: ERROR FLAG C
829 CCC * IF IERRO=0, COMPUTATION SUCCESSFUL. C
830 CCC * IF IERRO=1, COMPUTATION OUT OF RANGE. C
831 CCC * IF IERRO=2, ARGUMENT AND/OR ORDER, C
832 CCC OUT OF RANGE. C
833 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
834 CCC LOCAL VARIABLES: C
835 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
836 C A : (PNU), ORDER OF THE FUNCTIONS
837 C A2 : A**2
838 C AA : AUXILIARY VARIABLE IN THE IMPLEMENTATION
839 C OF THE LENZ-THOMPSON ALGORITHM
840 C AB : AUXILIARY VARIABLE IN THE IMPLEMENTATION
841 C OF THE LENZ-THOMPSON ALGORITHM
842 C B : AUXILIARY VARIABLE IN THE IMPLEMENTATION
843 C OF THE LENZ-THOMPSON ALGORITHM
844 C CC : AUXILIARY VARIABLE IN THE IMPLEMENTATION
845 C OF THE LENZ-THOMPSON ALGORITHM
846 C COSTH : COS(THET)
847 C D : AUXILIARY VARIABLE IN THE IMPLEMENTATION
848 C OF THE LENZ-THOMPSON ALGORITHM
849 C DELS : AUXILIARY VARIABLE IN THE IMPLEMENTATION
850 C OF THE LENZ-THOMPSON ALGORITHM
851 C DELTA : AUXILIARY VARIABLE IN THE IMPLEMENTATION
852 C OF THE LENZ-THOMPSON ALGORITHM
853 C FC : CONTINUED FRACTION FOR K(PNU+1)/K(PNU)
854 C FDOMIN : X*(COS(THET)+THET*SIN(THET))
855 C K : KIA FUNCTION
856 C KP : KIA' function
857 c pi : 3.1415..
858 c pia : pi*pnu
859 c pisq : sqrt(pi)
860 c preci : relative precision used in the calculation
861 c q0b : auxiliary variable in the implementation
862 c of the lenz-thompson algorithm
863 c q1b : auxiliary variable in the implementation
864 c of the lenz-thompson algorithm
865 c qb : auxiliary variable in the implementation
866 c of the lenz-thompson algorithm
867 c qqb : auxiliary variable in the implementation
868 c of the lenz-thompson algorithm
869 c s : auxiliary variable in the implementation
870 c of the lenz-thompson algorithm
871 c sinth : sin(thet)
872 c thet : asin(pnu/x)
873 c under : underflow number
874 c z0 : 1/(1+s)
875 ccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
876  DOUBLE PRECISION a,a2,aa,ab,b,cc,costh,d,d1mach,
877  + dels,delta,fc,fdomin,k,kp,pi,pia,pisq,pnu,preci,
878  + pser,pserd,q0b,q1b,qb,qqb,s,sinth,thet,under,x,z0
879  INTEGER ierro,ifac,mm
880 ccc machine dependent constant(underflow number)
881  under=d1mach(1)*1.d+8
882  pi=acos(-1.d0)
883  ierro=0
884  IF (ifac.EQ.1) THEN
885  IF (x.GE.pnu) THEN
886  sinth=pnu/x
887  thet=asin(sinth)
888  costh=cos(thet)
889  fdomin=x*(costh+thet*sinth)
890  IF (-fdomin.LE.log(under)) ierro=1
891  ELSE
892  IF ((-pi*pnu*0.5d0).LE.log(under)) ierro=1
893  ENDIF
894  ENDIF
895  IF (ierro.EQ.0) THEN
896  preci=d1mach(3)*10
897  a=pnu
898  a2=a*a
899  pisq=pi**0.5d0
900 ccc cf for k(nu+1)/k(nu)
901  mm=2
902  ab=-(0.25d0+a2)
903  aa=ab
904  cc=-ab
905  q0b=0.d0
906  q1b=1.d0
907  qqb=cc
908  b=2.d0*(1.d0+x)
909  d=1.d0/b
910  delta=d
911  fc=delta
912  s=qqb*delta
913  ab=-2.d0+ab
914  b=b+2.d0
915  91 d=1.d0/(b+ab*d)
916  delta=(b*d-1.d0)*delta
917  fc=fc+delta
918  cc=-ab*cc/mm
919  qb=(q0b-(b-2.d0)*q1b)/ab
920  q0b=q1b
921  q1b=qb
922  qqb=qqb+cc*q1b
923  dels=qqb*delta
924  s=s+dels
925  b=b+2.d0
926  ab=-2.d0*mm+ab
927  mm=mm+1
928  IF (mm.LT.10000) THEN
929  IF (abs(dels/s).GT.preci) goto 91
930  ENDIF
931  z0=1.d0/(1.d0+s)
932  IF (ifac.EQ.1) THEN
933  k=pisq*(2.d0*x)**(-0.5d0)*exp(-x)*z0
934  kp=-k/x*(0.5d0+x-(a2+0.25d0)*fc)
935  ELSE
936  IF (x.LT.a) THEN
937  pia=pi*a
938  k=pisq*(2.d0*x)**(-0.5d0)*exp(-x+pia*0.5d0)*z0
939  kp=-k/x*(0.5d0+x-(a2+0.25d0)*fc)
940  ELSE
941  sinth=a/x
942  thet=asin(sinth)
943  costh=cos(thet)
944  fdomin=x*(costh+thet*sinth)
945  k=pisq*(2.d0*x)**(-0.5d0)*exp(-x+fdomin)*z0
946  kp=-k/x*(0.5d0+x-(a2+0.25d0)*fc)
947  ENDIF
948  ENDIF
949  pser=k
950  pserd=kp
951  ELSE
952  pser=0.d0
953  pserd=0.d0
954  ENDIF
955  RETURN
956  END
957 
958  SUBROUTINE serkia(IFAC,X,PNU,PSER,PSERD,IERRO)
959 ccccccccccccccccccccccccccccccccccccccccccccccccccccc
960 ccc calculation of power series for k, k'
961 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
962 CCC INPUT VARIABLES: C
963 CCC X: ARGUMENT OF THE FUNCTIONS C
964 CCC PNU: ORDER OF THE FUNCTIONS C
965 CCC IFAC: C
966 CCC =1, NON SCALED FUNCTIONS C
967 CCC OTHERWISE, SCALED FUNCTIONS C
968 CCC OUTPUT VARIABLES: C
969 CCC PSER: K FUNCTION C
970 CCC PSERD: DERIVATIVE OF THE K FUNCTION C
971 CCC IERRO: ERROR FLAG C
972 CCC * IF IERRO=0, COMPUTATION SUCCESSFUL. C
973 CCC * IF IERRO=1, COMPUTATION OUT OF RANGE. C
974 CCC * IF IERRO=2, ARGUMENT AND/OR ORDER, C
975 CCC OUT OF RANGE. C
976 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
977 CCC LOCAL VARIABLES: C
978 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
979 C A : (PNU) ORDER OF THE FUNCTIONS
980 C A2 : A**2
981 C A2H : A**2/2
982 C A2N : A**(2*N)
983 C ACCP : ACCUMULATES THE P COEFFICIENTS
984 C ACCQ : ACCUMULATES THE Q COEFFICIENTS
985 C ARGU : SIG(0)-A*LOG(X/2)
986 C C : X**(2*K)/K!
987 C COCI : 1/(K**2+A**2)
988 C COSTH : COS(THET)
989 C DELTAK : ACCUMULATES THE CONTRIBUTION FOR THE
990 C KIA(X) FUNCTION
991 C DELTKP : ACCUMULATES THE CONTRIBUTION FOR THE
992 C KIA'(x) function
993 c df1 : factor(depending on ifac)
994 c eta0 : PARAMETER for the calculation of the
995 c coulomb phase shift
996 c eta02 : eta0**2
997 c f(k) : sin(phi(a,k)-a*log(x/2))
998 c /(a**2*(1+a**2)...(k**2+a**2))**1/2,
999 c WHERE phi(a,k)=phase(gamma(1+k+ia))
1000 c fdomin : x*(cos(thet)+thet*sin(thet))
1001 c k : contribution to the kia(x) function
1002 c kp : contribution to the kia'(X) FUNCTION
1003 C OVER : OVERFLOW NUMBER
1004 C P0 : PARAMETERS FOR THE CALCULATION OF THE COULOMB
1005 C PHASE SHIFT
1006 C P1 : PARAMETERS FOR THE CALCULATION OF THE COULOMB
1007 C PHASE SHIFT
1008 C P2 : PARAMETERS FOR THE CALCULATION OF THE COULOMB
1009 C PHASE SHIFT
1010 C PI : 3.1415...
1011 C PIA : PI*A
1012 C PIA2 : 2*PI*A
1013 C PRECI : RELATIVE PRECISION USED IN THE CALCULATION
1014 C Q0 : PARAMETERS FOR THE CALCULATION OF THE COULOMB
1015 C PHASE SHIFT
1016 C Q1 : PARAMETERS FOR THE CALCULATION OF THE COULOMB
1017 C PHASE SHIFT
1018 C Q2 : PARAMETERS FOR THE CALCULATION OF THE COULOMB
1019 C PHASE SHIFT
1020 C R(K) : F(K)*A/TAN(PHI(A,K)-A*LOG(X/2))
1021 C SIG0 : COULOMB PHASE SHIFT
1022 C SINTH : SIN(THET)
1023 C THET : ASIN(A/X)
1024 C UNDER : UNDERFLOW NUMBER
1025 C X2 : X*X
1026 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1027  DOUBLE PRECISION A,A2,A2N,ACCP,ACCQ,ARGU,C,COCI,COSTH,
1028  + D1MACH,DDS,DEE,DELTAK,DELTKP,DF1,DSMALL,ETA0,ETA02,
1029  + EULER,F(0:500),FDOMIN,K,KP,OVER,P0(0:9),P1(0:9),P2(0:6),
1030  + PI,PIA,PIA2,PNU,PRECI,PSER,PSERD,Q0(0:9),Q1(0:9),Q2(0:6),
1031  + R(0:500),SIG0,SINTH,THET,UNDER,X,X2
1032  INTEGER IERRO,IFAC,L,M,N
1033  SAVE P0,Q0,P1,Q1,P2,Q2
1034  DATA P0/1.08871504904797411683D5,3.64707573081160914640D5,
1035  + 4.88801471582878013158D5,3.36275736298197324009D5,
1036  + 1.26899226277838479804D5,
1037  + 2.60795543527084582682D4,2.73352480554497990544D3,
1038  + 1.26447543569902963184D2,
1039  + 1.85446022125533909390D0,1.90716219990037648146D-3/
1040  DATA Q0/6.14884786346071135090D5,2.29801588515708014282D6,
1041  + 3.50310844128424021934D6,
1042  + 2.81194990286041080264D6,1.28236441994358406742D6,
1043  + 3.35209348711803753154D5,
1044  + 4.84319580247948701171D4,3.54877039006873206531D3,
1045  + 1.11207201299804390166D2,1.D0/
1046  DATA P1/-1.044100987526487618670D10,-1.508574107180079913696D10,
1047  + -5.582652833355901160542D9,4.052529174369477275446D8,
1048  + 5.461712273118594275192D8,
1049  + 9.510404403068169395714D7,6.281126609997342119416D6,
1050  + 1.651178048950518520416D5,
1051  + 1.498824421329341285521D3,2.974686506595477984776D0/
1052  DATA Q1/1.808868161493543887787D10,3.869142051704700267785D10,
1053  + 3.003264575147162634046D10,1.075554651494601843525D10,
1054  + 1.901298501823290694245D9,
1055  + 1.665999832151229472632D8,6.952188089169487375936D6,
1056  + 1.253235080625688652718D5,7.904420414560291396996D2,1.D0/
1057  DATA P2/7.08638611024520906826D-3,-6.54026368947801591128D-2,
1058  + 2.92684143106158043933D-1,4.66821392319665609167D0,
1059  + -3.43943790382690949054D0,
1060  + -7.72786486869252994370D0,-9.88841771200290647461D-01/
1061  DATA Q2/-7.08638611024520908189D-3,6.59931690706339630254D-2,
1062  + -2.98754421632058618922D-1,-4.63752355513412248006D0,
1063  + 3.79700454098863541593D0,7.06184065426336718524D0,1.D0/
1064  PI=ACOS(-1.D0)
1065  ETA0=1.8055470716051069198764D0
1066  EULER=0.577215664901532860606512D0
1067  PRECI=D1MACH(3)*10
1068 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
1069  OVER=D1MACH(2)*1.D-8
1070 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
1071  UNDER=D1MACH(1)*1.D+8
1072  IERRO=0
1073 .EQ. IF (IFAC1) THEN
1074 .GE. IF (XPNU) THEN
1075  SINTH=PNU/X
1076  THET=ASIN(SINTH)
1077  COSTH=COS(THET)
1078  FDOMIN=X*(COSTH+THET*SINTH)
1079 .LE. IF (-FDOMINLOG(UNDER)) IERRO=1
1080  ELSE
1081 .LE. IF ((-PI*PNU*0.5D0)LOG(UNDER)) IERRO=1
1082  ENDIF
1083  ENDIF
1084 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1085 CCC COEFFICIENTS FOR THE CALCULATION OF THE COULOMB PHASE SHIFT
1086 CCC FROM CODY & HILLSTROM, MATH. COMPUT. 24(111) 1970
1087 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1088 .EQ. IF (IERRO0) THEN
1089  A=PNU
1090  ETA02=ETA0*ETA0
1091  A2=A*A
1092  N=0
1093  ACCP=0.D0
1094  ACCQ=0.D0
1095 .LE. IF (A2.D0) THEN
1096  33 A2N=A2**N
1097  ACCP=ACCP+P0(N)*A2N
1098  ACCQ=ACCQ+Q0(N)*A2N
1099  N=N+1
1100 .LE. IF (N9) GOTO 33
1101  SIG0=A*(A2-ETA02)*ACCP/ACCQ
1102  ELSE
1103 .GT..AND..LE. IF ((A2.D0)(A4.D0)) THEN
1104  44 A2N=A2**N
1105  ACCP=ACCP+P1(N)*A2N
1106  ACCQ=ACCQ+Q1(N)*A2N
1107  N=N+1
1108 .LE. IF (N9) GOTO 44
1109  SIG0=A*ACCP/ACCQ
1110  ELSE
1111  55 A2N=A2**N
1112  ACCP=ACCP+P2(N)/A2N
1113  ACCQ=ACCQ+Q2(N)/A2N
1114  N=N+1
1115 .LE. IF (N6) GOTO 55
1116  SIG0=ATAN(A)*0.5D0+A*(LOG(1.D0+A2)*0.5D0+ACCP/ACCQ)
1117  ENDIF
1118  ENDIF
1119 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1120 CCC EVALUATION OF F(0), R(0), R(1)
1121 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1122  PIA=PI*A
1123  ARGU=SIG0-A*LOG(X*0.5D0)
1124  R(0)=COS(ARGU)
1125  R(1)=1.D0/(1.D0+A2)*(R(0)-A*SIN(ARGU))
1126 .LT. IF (AUNDER) THEN
1127  F(0)=-(EULER+LOG(X*0.5D0))
1128  ELSE
1129  F(0)=1.D0/A*SIN(ARGU)
1130  ENDIF
1131  C=1.D0
1132  K=F(0)
1133  KP=-0.5D0*R(0)
1134 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1135 CCC RECURSION FOR F(K), R(K), C(K)
1136 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1137  X2=0.25D0*X*X
1138  M=1
1139  COCI=1.D0/(M*M+A2)
1140  DELTAK=K*10
1141  DELTKP=KP*10
1142  66 F(M)=(M*F(M-1)+R(M-1))*COCI
1143  C=X2*C/M
1144  DELTAK=F(M)*C
1145  K=K+DELTAK
1146  DELTKP=(M*F(M)-0.5D0*R(M))*C
1147  KP=KP+DELTKP
1148  M=M+1
1149  COCI=1.D0/(M*M+A2)
1150  R(M)=((2.D0*M-1.D0)*R(M-1)-R(M-2))*COCI
1151 .GT. IF (KOVER) M=500
1152 .GT. IF (KPOVER) M=500
1153 .LT. IF (M500) THEN
1154 .GT..OR. IF ((ABS(DELTAK/K)PRECI)
1155 .GT. + (ABS(DELTKP/KP)PRECI)) GOTO 66
1156  ENDIF
1157  PIA2=2.D0*PIA
1158 .EQ. IF (IFAC1) THEN
1159 .LE. IF (-PIA2LOG(UNDER)) THEN
1160  DEE=0.D0
1161  ELSE
1162  DEE=EXP(-PIA2)
1163  ENDIF
1164 .LT. IF (A1.D-1) THEN
1165 .LT. IF (AUNDER) THEN
1166  DF1=1.D0
1167  ELSE
1168  L=0
1169  DDS=1.D0
1170  DSMALL=1.D0
1171  47 L=L+1
1172  DDS=-DDS*PIA2/(L+1.D0)
1173  DSMALL=DSMALL+DDS
1174 .GT. IF (ABS(DDS/DSMALL)PRECI) GOTO 47
1175  DF1=EXP(PIA*0.5D0)*SQRT(DSMALL)
1176  ENDIF
1177  ELSE
1178  DF1=EXP(PIA*0.5D0)*SQRT((1.D0-DEE)/PIA2)
1179  ENDIF
1180  PSER=K/DF1
1181  PSERD=KP*2.D0/X/DF1
1182  ELSE
1183 .LT. IF (XA) THEN
1184 .LE. IF (-PIA2LOG(UNDER)) THEN
1185  DEE=0.D0
1186  ELSE
1187  DEE=EXP(-PIA2)
1188  ENDIF
1189 .LT. IF (A1.D-1) THEN
1190 .LT. IF (AUNDER) THEN
1191  DF1=1.D0
1192  ELSE
1193  L=0
1194  DDS=1.D0
1195  DSMALL=1.D0
1196  48 L=L+1
1197  DDS=-DDS*PIA2/(L+1.D0)
1198  DSMALL=DSMALL+DDS
1199 .GT. IF (ABS(DDS/DSMALL)PRECI) GOTO 48
1200  DF1=SQRT(DSMALL)
1201  ENDIF
1202  ELSE
1203  DF1=SQRT((1.D0-DEE)/PIA2)
1204  ENDIF
1205  ELSE
1206  SINTH=A/X
1207  THET=ASIN(SINTH)
1208  COSTH=COS(THET)
1209  FDOMIN=X*(COSTH+THET*SINTH)
1210 .LE. IF (-PIA2LOG(UNDER)) THEN
1211  DEE=0.D0
1212  ELSE
1213  DEE=EXP(-PIA2)
1214  ENDIF
1215 .LT. IF (A1.D-1) THEN
1216 .LT. IF (AUNDER) THEN
1217  DF1=1.D0
1218  ELSE
1219  L=0
1220  DDS=1.D0
1221  DSMALL=1.D0
1222  49 L=L+1
1223  DDS=-DDS*PIA2/(L+1.D0)
1224  DSMALL=DSMALL+DDS
1225 .GT. IF (ABS(DDS/DSMALL)PRECI) GOTO 49
1226  DF1=EXP(PIA*0.5D0-FDOMIN)*SQRT(DSMALL)
1227  ENDIF
1228  ELSE
1229  DF1=EXP(PIA*0.5D0-FDOMIN)*SQRT((1.D0-DEE)/PIA2)
1230  ENDIF
1231  ENDIF
1232  PSER=K/DF1
1233  PSERD=KP*2.D0/X/DF1
1234  ENDIF
1235  ELSE
1236  PSER=0.D0
1237  PSERD=0.D0
1238  ENDIF
1239  RETURN
1240  END
1241 
1242  SUBROUTINE AIEXKI(IFAC,X,A,DKAI,DKAID,IERROK)
1243 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1244 CCC AIRY-TYPE ASYMPTOTIC EXPANSION FOR THE K AND K' c
1245 ccc functions c
1246 ccc c
1247 ccc input variables: c
1248 ccc x: argument of the functions c
1249 ccc a: order of the functions c
1250 ccc ifac: c
1251 ccc =1, non scaled functions c
1252 ccc otherwise, scaled functions c
1253 ccc output variables: c
1254 ccc dkai: k FUNCTION c
1255 ccc dkaid: derivative of the k FUNCTION c
1256 ccc ierrok: error flag c
1257 ccc * IF ierrok=0, computation successful. c
1258 ccc * IF ierrok=1, computation out of range. c
1259 ccc * IF ierrok=2, argument and/or order, c
1260 ccc out of range. c
1261 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1262 ccc local variables: c
1263 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1264 c a0ex : exact value of the coefficient a_0
1265 c a13 : a**(1/3)
1266 c a2 : a**2
1267 c a23 : a**(2/3)
1268 c a2k : a**(2*k)
1269 c ac : coefficients a_s(psi) used in the expansions
1270 c aii : imaginary part of the airy FUNCTION ai(Z)
1271 c air : REAL part of the airy function ai(z)
1272 c apihal : a*pi/2
1273 c arg : -a**(2/3)*psi
1274 c as : accumulates the contribution of the ac coefficients
1275 c(for the kia(x) function)
1276 c asp : accumulates the contribution of the ac coefficients
1277 c(for the kia'(X) FUNCTION)
1278 C B0EX : EXACT VALUE OF THE COEFFICIENT B_0
1279 C B0PEX : EXACT VALUE OF THE COEFFICIENT B'_0
1280 c bc : coefficients b_s(psi) used in the expansions
1281 c bs : accumulates the contribution of the bc coefficients
1282 c(for the kia(x) function)
1283 c bso : accumulates the old values of bs
1284 c bsp : accumulates the contribution of the bc coefficients
1285 c(for the kia'(X) FUNCTION)
1286 C BSPO : ACCUMULATES THE OLD VALUES OF BSP
1287 C C0EX : EXACT VALUE OF THE COEFFICIENT C_0
1288 C CHI : ACCUMULATES THE CONTRIBUTION OF THE CHIN COEFFICIENTS
1289 C CHIEX : EXACT VALUE OF THE FUNCTION CHI(PSI)
1290 C CHIN : COEFFICIENTS FOR THE EXPANSION OF THE FUNCTION
1291 C CHI(PSI)
1292 C COSTH : COS(THET)
1293 C D0EX : EXACT VALUE OF THE COEFFICIENT D_0
1294 C DAII : IMAGINARY PART OF THE AIRY FUNCTION AI'(z)
1295 c dair : REAL part of the airy function ai'(Z)
1296 C DF : 2**(1/3)
1297 C DZZ : (ABS(1-Z**2))**(1/2)
1298 C ETA : PSI/2**(1/3)
1299 C ETAJ : ETA**J
1300 C ETAK : ETA**K
1301 C ETAL : ETA**L
1302 C EXPAM : EXP(A*PI/2-FDOMIN)
1303 C EXPAPI : EXP(A*PI/2)
1304 C F2 : (-)**K/A**(2*K)
1305 C F4 : (A**(1/3))**4
1306 C FAC : FACTOR FOR THE KIA(X) FUNCTION
1307 C FACD : FACTOR FOR THE KIA'(x) function
1308 c fdomin : x*(cos(thet)+thet*sin(thet))
1309 c ierro : error flag for the airy functions
1310 c ifaca :
1311 c *IF ifaca=1, computation of unscaled airy ai(z),ai'(Z)
1312 C *IF IFACA=2, COMPUTATION OF SCALED AIRY AI(Z),AI'(z)
1313 c ifun :
1314 c *IF ifun=1, computation of airy ai(z)
1315 c *IF ifun=2, computation of airy ai'(Z)
1316 C PHI : COEFFICIENTS FOR THE EXPANSION OF THE FUNCTION
1317 C PHI(PSI)
1318 C PHIEX : EXACT VALUE OF THE FUNCTION PHI(PSI)
1319 C PHIEX2 : PHI(PSI)**2
1320 C PHIS : VALUE OF THE FUNCTION PHI(PSI)
1321 C PI : 3.1415...
1322 C PIHALF : PI/2
1323 C PSI :
1324 C *IF 0<=Z<=1, 2/3*PSI**3/2=LOG((1+SQRT(1-Z**2))/Z)
1325 C -SQRT(1-Z**2)
1326 C *IF Z>=1, 2/3*(-PSI)**3/2=SQRT(Z**2-1)-ARCCOS(1/Z)
1327 C PSI12 : SQRT(PSI)
1328 C PSI2 : PSI*PSI
1329 C PSI3 : PSI**3
1330 C SAS : ACCUMULATES THE CONTRIBUTION OF A_S(PSI)
1331 C SBS : ACCUMULATES THE CONTRIBUTION OF B_S(PSI)
1332 C SCS : ACCUMULATES THE CONTRIBUTION OF C_S(PSI)
1333 C SDS : ACCUMULATES THE CONTRIBUTION OF D_S(PSI)
1334 C SIG : (-)**K
1335 C SINTH : SIN(THET)
1336 C THET : ASIN(A/X)
1337 C UNDER : UNDERFLOW NUMBER
1338 C Y : Z-1
1339 C Z : X/A
1340 C Z2 : Z**2
1341 C Z21M : 1-Z**2
1342 C ZDS : SQRT(1-Z**2)
1343 C ZMASF : EXPANSION IN TERMS OF (Z-1) FOR THE
1344 C CALCULATION OF PSI(Z)
1345 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1346  DOUBLE PRECISION A,A0EX,A13,A2,A23,A2K,AC(0:5,0:20),
1347  + AII,AIR,APIHAL,ARG,AS,ASP,B0EX,B0PEX,BC(0:5,0:20),
1348  + BS,BSO,BSP,BSPO,C0EX,CHI,CHIEX,CHIN(0:26),COSTH,
1349  + COZMAS(0:12),D0EX,D1MACH,DAII,DAIR,DF,DKAI,DKAID
1350  DOUBLE PRECISION DZZ,ETA,ETAJ,ETAK,ETAL,EXPAM,
1351  + EXPAPI,F2,F4,FAC,FACD,FDOMIN,PHI(0:35),PHIEX,PHIEX2,
1352  + PHIS,PI,PIHALF,PSI,PSI12,PSI2,PSI3,SAS,SBS,SCS,SDS,
1353  + SIG,SINTH,THET,UNDER,X,Y,Z,Z2,Z21M,ZDS,ZMASF
1354  INTEGER IERRO,IERROK,IFAC,IFACA,IFUN,INDA(0:5),INDB(0:5),
1355  + J,K,L
1356  SAVE PHI,CHIN,AC,BC
1357 CCCC VALUES OF THE COEFFICIENTS
1358  DATA PHI/1.D0,0.2D0,.25714285714285714286D-1,
1359  + -.56507936507936507937D-2,-.39367965367965367965D-2,
1360  + -.5209362066504924D-3,.3708541264731741D-3,
1361  + .2123827840293627D-3,.2150629788753145D-4,
1362  + -.2636904062981799D-4,-.1405469826493129D-4,
1363  + -.1149328899029441D-5,.1972641193938624D-5,
1364  + .1014324305593329D-5,.7034331100192200D-7,
1365  + -.1525044777392676D-6,-.7677866256900572D-7,
1366  + -.4669842638693018D-8,.1206673645965329D-7,
1367  + .5990877668092344D-8,.3269102150077715D-9,
1368  + -.97138350244428335095D-9,-.47745934295232233834D-9,
1369  + -.23750318642839155779D-10,.79244598109106655567D-10,
1370  + .38653584230817865528D-10,.17734105846426807873D-11,
1371  + -.65332110030864751956D-11,-.31673512094772575686D-11,
1372  + -.13524195177121030660D-12,.54324103217903338951D-12,
1373  + .26204918647967626464D-12,.10488134973664584898D-13,
1374  + -.45490420121539001435D-13,-.21851238232690420878D-13,
1375  + -.82456080260379042800D-15/
1376  DATA CHIN/.2D0,.1142857142857142857142857D-1,
1377  + -.2438095238095238095238095D-1,-.1003471449185734900020614D-1,
1378  + .8811404468547325690182833D-3,.2318339655809043564145605D-2,
1379  + .7794494413564441575646057D-3,-.1373952504077228744949558D-3,
1380  + -.2162301322540308393022684D-3,-.6585004634375583559042795D-4,
1381  + .1502851367097217578058824D-4,.1991904617871647675455262D-4,
1382  + .5731972706719818525615427D-5,-.1496427612747891044606885D-5,
1383  + -.1821231830428939670133992D-5,-.5054254875930882534538217D-6,
1384  + .1433283859497625931203415D-6,.1657681319078678321113361D-6,
1385  + .4485592642702540575627044D-7,-.1346138242826094098161839D-7,
1386  + -.1504601862773535117708677D-7,-.3995660198654805921651406D-8,
1387  + .1250124931952495738882300D-8,.1363187174221864073749532D-8,
1388  + .3567608459777506132029204D-9,-.1152749290139859167119863D-9,
1389  + -.1233547289799408517691696D-9/
1390  DATA COZMAS/.9428090415820634D0,-.4242640687119285D0,
1391  + .2904188565587606D0,-.2234261009999160D0,
1392  + .1821754153944745D0,-.1539625154198624D0,
1393  + .1333737583085222D0,-.1176617834148007D0,
1394  + .1052687194772381D0,-.9524025714638822D-1,
1395  + .8695738197073783D-1,-.8000034897653656D-1,
1396  + .7407421797273086D-1/
1397  DATA AC(0,0)/1.D0/
1398  DATA AC(1,0)/-.44444444444444444445D-2/
1399  DATA AC(1,1)/-.18441558441558441558D-2/
1400  DATA AC(1,2)/.11213675213675213675D-2/
1401  DATA AC(1,3)/.13457752124418791086D-2/
1402  DATA AC(1,4)/.0003880626562979504D0/
1403  DATA AC(1,5)/-.0001830686723781799D0/
1404  DATA AC(1,6)/-.0001995460887806733D0/
1405  DATA AC(1,7)/-.00005256191234041587D0/
1406  DATA AC(1,8)/.00002460619652459158D0/
1407  DATA AC(1,9)/.00002519246780924541D0/
1408  DATA AC(1,10)/.6333157376533242D-5/
1409  DATA AC(1,11)/-.2957485733830202D-5/
1410  DATA AC(1,12)/-.2925255920564838D-5/
1411  DATA AC(1,13)/-.7159702610502009D-6/
1412  DATA AC(1,14)/.3331510720390949D-6/
1413  DATA AC(1,15)/.3227670475692310D-6/
1414  DATA AC(1,16)/.7767729381664199D-7/
1415  DATA AC(1,17)/-.3600954237921120D-7/
1416  DATA AC(1,18)/-.3441724449034226D-7/
1417  DATA AC(1,19)/-.8188194356398772D-8/
1418  DATA AC(1,20)/.3783148485152038D-8/
1419  DATA AC(2,0)/.69373554135458897365D-3/
1420  DATA AC(2,1)/.46448349036584330703D-3/
1421  DATA AC(2,2)/-.42838130171535112460D-3/
1422  DATA AC(2,3)/-.0007026702868771135D0/
1423  DATA AC(2,4)/-.0002632580046778811D0/
1424  DATA AC(2,5)/.0001663853666288703D0/
1425  DATA AC(2,6)/.0002212087687818584D0/
1426  DATA AC(2,7)/.00007020345615329662D0/
1427  DATA AC(2,8)/-.00004000421782540614D0/
1428  DATA AC(2,9)/-.00004786324966453962D0/
1429  DATA AC(2,10)/-.00001394600741473631D0/
1430  DATA AC(2,11)/.7536186591273727D-5/
1431  DATA AC(2,12)/.8478502161067667D-5/
1432  DATA AC(2,13)/.2345355228453912D-5/
1433  DATA AC(2,14)/-.1225943294710883D-5/
1434  DATA AC(2,15)/-.1325082343401027D-5/
1435  DATA AC(2,16)/-.3539954776569997D-6/
1436  DATA AC(2,17)/.1808291719376674D-6/
1437  DATA AC(2,18)/.1900383515233655D-6/
1438  DATA AC(3,0)/-.35421197145774384076D-3/
1439  DATA AC(3,1)/-.31232252789031883276D-3/
1440  DATA AC(3,2)/.3716442237502298D-3/
1441  DATA AC(3,3)/.0007539269155977733D0/
1442  DATA AC(3,4)/.0003408300059444739D0/
1443  DATA AC(3,5)/-.0002634968172069594D0/
1444  DATA AC(3,6)/-.0004089275726648432D0/
1445  DATA AC(3,7)/-.0001501108759563460D0/
1446  DATA AC(3,8)/.00009964015205538056D0/
1447  DATA AC(3,9)/.0001352492955751283D0/
1448  DATA AC(3,10)/.00004443117087272903D0/
1449  DATA AC(3,11)/-.00002713205071914117D0/
1450  DATA AC(3,12)/-.00003396796969771860D0/
1451  DATA AC(3,13)/-.00001040708865273043D0/
1452  DATA AC(3,14)/.6024639065414414D-5/
1453  DATA AC(3,15)/.7143919607846883D-5/
1454  DATA AC(4,0)/.378194199201773D-3/
1455  DATA AC(4,1)/.000404943905523634D0/
1456  DATA AC(4,2)/-.000579130526946453D0/
1457  DATA AC(4,3)/-.00138017901171011D0/
1458  DATA AC(4,4)/-.000722520056780091D0/
1459  DATA AC(4,5)/.000651265924036825D0/
1460  DATA AC(4,6)/.00114674563328389D0/
1461  DATA AC(4,7)/.000474423189340405D0/
1462  DATA AC(4,8)/-.000356495172735468D0/
1463  DATA AC(4,9)/-.000538157791035111D0/
1464  DATA AC(4,10)/-.000195687390661225D0/
1465  DATA AC(4,11)/.000132563525210293D0/
1466  DATA AC(4,12)/.000181949256267291D0/
1467  DATA AC(5,0)/-.69114139728829416760D-3/
1468  DATA AC(5,1)/-.00085995326611774383285D0/
1469  DATA AC(5,2)/.0014202335568143511489D0/
1470  DATA AC(5,3)/.0038535426995603052443D0/
1471  DATA AC(5,4)/.0022752811642901374595D0/
1472  DATA AC(5,5)/-.0023219572034556988366D0/
1473  DATA AC(5,6)/-.0045478643394434635622D0/
1474  DATA AC(5,7)/-.0020824431758272449919D0/
1475  DATA AC(5,8)/.0017370443573195808719D0/
1476  DATA BC(0,0)/.14285714285714285714D-1/
1477  DATA BC(0,1)/.88888888888888888889D-2/
1478  DATA BC(0,2)/.20482374768089053803D-2/
1479  DATA BC(0,3)/-.57826617826617826618D-3/
1480  DATA BC(0,4)/-.60412089799844901886D-3/
1481  DATA BC(0,5)/-.0001472685745626922D0/
1482  DATA BC(0,6)/.00005324102148009784D0/
1483  DATA BC(0,7)/.00005206561006583416D0/
1484  DATA BC(0,8)/.00001233115050894939D0/
1485  DATA BC(0,9)/-.4905932728531366D-5/
1486  DATA BC(0,10)/-.4632230987136350D-5/
1487  DATA BC(0,11)/-.1077174523455235D-5/
1488  DATA BC(0,12)/.4475963978932822D-6/
1489  DATA BC(0,13)/.4152586188464624D-6/
1490  DATA BC(0,14)/.9555819293589234D-7/
1491  DATA BC(0,15)/-.4060599208403059D-7/
1492  DATA BC(0,16)/-.3731367187988482D-7/
1493  DATA BC(0,17)/-.8532670645553778D-8/
1494  DATA BC(0,18)/.3673017245573624D-8/
1495  DATA BC(0,19)/.3355960460784536D-8/
1496  DATA BC(0,20)/.7643107095110475D-9/
1497  DATA BC(1,0)/-.11848595848595848596D-2/
1498  DATA BC(1,1)/-.13940630797773654917D-2/
1499  DATA BC(1,2)/-.48141005586383737645D-3/
1500  DATA BC(1,3)/.26841705366016142958D-3/
1501  DATA BC(1,4)/.0003419706982709903D0/
1502  DATA BC(1,5)/.0001034548234902078D0/
1503  DATA BC(1,6)/-.00005418191982095504D0/
1504  DATA BC(1,7)/-.00006202184830690167D0/
1505  DATA BC(1,8)/-.00001724885886056087D0/
1506  DATA BC(1,9)/.8744675992887053D-5/
1507  DATA BC(1,10)/.9420684216180929D-5/
1508  DATA BC(1,11)/.2494922112085850D-5/
1509  DATA BC(1,12)/-.1238458608836357D-5/
1510  DATA BC(1,13)/-.1285461713809769D-5/
1511  DATA BC(1,14)/-.3299710862537507D-6/
1512  DATA BC(1,15)/.1613441105788315D-6/
1513  DATA BC(1,16)/.1633623194402374D-6/
1514  DATA BC(1,17)/.4104252949605779D-7/
1515  DATA BC(1,18)/-.1984317042326989D-7/
1516  DATA BC(1,19)/-.1973948142769706D-7/
1517  DATA BC(1,20)/-.4882194808588752D-8/
1518  DATA BC(2,0)/.43829180944898810994D-3/
1519  DATA BC(2,1)/.71104865116708668943D-3/
1520  DATA BC(2,2)/.31858383945387580576D-3/
1521  DATA BC(2,3)/-.0002404809426804458D0/
1522  DATA BC(2,4)/-.0003722966038621536D0/
1523  DATA BC(2,5)/-.0001352752059595618D0/
1524  DATA BC(2,6)/.00008691694372704142D0/
1525  DATA BC(2,7)/.0001158750753591377D0/
1526  DATA BC(2,8)/.00003724965927846447D0/
1527  DATA BC(2,9)/-.00002198334949606935D0/
1528  DATA BC(2,10)/-.00002686449633870452D0/
1529  DATA BC(2,11)/-.8023061612032524D-5/
1530  DATA BC(2,12)/.4494756592180126D-5/
1531  DATA BC(2,13)/.5193504763856015D-5/
1532  DATA BC(2,14)/.1477156191529617D-5/
1533  DATA BC(2,15)/-.7988793826096784D-6/
1534  DATA BC(3,0)/-.37670439477105454219D-3/
1535  DATA BC(3,1)/-.75856271658798642365D-3/
1536  DATA BC(3,2)/-.0004103253968775048D0/
1537  DATA BC(3,3)/.0003791263310429010D0/
1538  DATA BC(3,4)/.0006850981673903450D0/
1539  DATA BC(3,5)/.0002878310571932216D0/
1540  DATA BC(3,6)/-.0002157010636115705D0/
1541  DATA BC(3,7)/-.0003260863991373500D0/
1542  DATA BC(3,8)/-.0001181317008748678D0/
1543  DATA BC(3,9)/.00007887526841582582D0/
1544  DATA BC(3,10)/.0001072081833420685D0/
1545  DATA BC(3,11)/.00003544595251288735D0/
1546  DATA BC(3,12)/-.00002201447920733824D0/
1547  DATA BC(3,13)/-.00002789336359620813D0/
1548  DATA BC(4,0)/.00058453330122076187255D0/
1549  DATA BC(4,1)/.0013854690422372401251D0/
1550  DATA BC(4,2)/.00086830374184946900245D0/
1551  DATA BC(4,3)/-.00093502904801345951693D0/
1552  DATA BC(4,4)/-.0019175486005525492059D0/
1553  DATA BC(4,5)/-.00090795047113308137941D0/
1554  DATA BC(4,6)/.00077050429806392235104D0/
1555  DATA BC(4,7)/.0012953100128255983488D0/
1556  DATA BC(4,8)/.00051933869471899550762D0/
1557  DATA BC(4,9)/-.00038482631948759834653D0/
1558  DATA BC(4,10)/-.00057335393099012476502D0/
1559  DATA BC(5,0)/-.0014301070053470410656D0/
1560  DATA BC(5,1)/-.0038637811942002539408D0/
1561  DATA BC(5,2)/-.0027322816261168328245D0/
1562  DATA BC(5,3)/.0033294980346743452748D0/
1563  DATA BC(5,4)/.0075972237660887795911D0/
1564  DATA BC(5,5)/.0039816655673062060620D0/
1565  DATA BC(5,6)/-.0037510180460986006595D0/
1566  DATA INDA/0,20,18,15,12,8/
1567  DATA INDB/20,20,15,13,10,6/
1568 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
1569  UNDER=D1MACH(1)*1.D+8
1570 CCCC CONSTANTS CCCCCCCCCCCCCC
1571  PI=ACOS(-1.D0)
1572  PIHALF=PI*0.5D0
1573  IERROK=0
1574 .EQ. IF (IFAC1) THEN
1575 .GE. IF (XA) THEN
1576  SINTH=A/X
1577  THET=ASIN(SINTH)
1578  COSTH=COS(THET)
1579  FDOMIN=X*(COSTH+THET*SINTH)
1580 .LE. IF (-FDOMINLOG(UNDER)) IERROK=1
1581  ELSE
1582 .LE. IF ((-PI*A*0.5D0)LOG(UNDER)) IERROK=1
1583  ENDIF
1584  ENDIF
1585 .EQ. IF (IERROK0) THEN
1586  DF=2.D0**(1.D0/3.D0)
1587 CCCC VARIABLES CCCCCCCCCCCCCCC
1588  Z=X/A
1589  A2=A*A
1590  A13=A**(1.D0/3.D0)
1591  A23=A13*A13
1592 .EQ. IF (IFAC1) THEN
1593  APIHAL=A*PIHALF
1594  EXPAPI=EXP(-APIHAL)
1595  FAC=PI*EXPAPI/A13
1596  FACD=2.D0*PI*EXPAPI/Z/A23
1597  ELSE
1598 .LT. IF (XA) THEN
1599  FAC=PI/A13
1600  FACD=2.D0*PI/Z/A23
1601  ELSE
1602  SINTH=A/X
1603  THET=ASIN(SINTH)
1604  COSTH=COS(THET)
1605  FDOMIN=X*(COSTH+THET*SINTH)
1606  EXPAM=EXP(-A*PIHALF+FDOMIN)
1607  FAC=PI*EXPAM/A13
1608  FACD=2.D0*PI*EXPAM/Z/A23
1609  ENDIF
1610  ENDIF
1611  F4=A13**4
1612 .LE. IF (Z0.9D0) THEN
1613  ZDS=SQRT((1.D0-Z)*(1.D0+Z))
1614  PSI=(1.5D0*(LOG((1.D0+ZDS)/Z)-ZDS))**(2.D0/3.D0)
1615 .GT. ELSEIF (Z1.1D0) THEN
1616  ZDS=SQRT((Z-1.D0)*(1.D0+Z))
1617  PSI=-(1.5D0*(ZDS-ACOS(1.D0/Z)))**(2.D0/3.D0)
1618  ELSE
1619  Y=Z-1.D0
1620  ZMASF=COZMAS(12)
1621  DO 8 K=0,11
1622  J=11-K
1623  ZMASF=COZMAS(J)+Y*ZMASF
1624  8 CONTINUE
1625  ZMASF=ABS(Y)**(1.5D0)*ZMASF
1626 .LT. IF (Z1.D0) THEN
1627  PSI=(1.5D0*ZMASF)**(2.D0/3.D0)
1628  ELSE
1629  PSI=-(1.5D0*ZMASF)**(2.D0/3.D0)
1630  ENDIF
1631  ENDIF
1632  ETA=PSI/DF
1633  ARG=-A23*PSI
1634  PSI2=PSI*PSI
1635  PSI3=PSI2*PSI
1636 .GT..AND..LT. IF ((Z0.8D0)(Z1.2D0)) THEN
1637  PHIS=0.D0
1638  CHI=0.D0
1639  SAS=0.D0
1640  SBS=0.D0
1641  SDS=0.D0
1642  SCS=1.D0
1643  BS=0.D0
1644  BSPO=0.D0
1645  DO 41 L=0,20
1646 .EQ. IF (L0) THEN
1647  ETAL=1.D0
1648  ELSE
1649  ETAL=ETAL*ETA
1650  ENDIF
1651  CHI=CHI+CHIN(L)*ETAL
1652  41 CONTINUE
1653  CHI=CHI/DF
1654  DO 10 K=0,20
1655  BSO=BS
1656  BSPO=BSP
1657  AS=0.D0
1658  BS=0.D0
1659  ASP=0.D0
1660  BSP=0.D0
1661 .EQ. IF (K0) THEN
1662  ETAK=1.D0
1663  A2K=1.D0
1664  SIG=1.D0
1665  ELSE
1666  ETAK=ETAK*ETA
1667  A2K=A2K*A2
1668  SIG=-1.D0*SIG
1669  ENDIF
1670  PHIS=PHIS+PHI(K)*ETAK
1671  F2=SIG/A2K
1672 .LE. IF (K5) THEN
1673  DO 20 J=0,20
1674 .EQ. IF (J0) THEN
1675  ETAJ=1.D0
1676  ELSE
1677  ETAJ=ETAJ*ETA
1678  ENDIF
1679 .LE. IF (JINDA(K)) THEN
1680  AS=AS+AC(K,J)*ETAJ
1681  ENDIF
1682 .LE. IF (JINDB(K)) THEN
1683  BS=BS+BC(K,J)*ETAJ
1684  ENDIF
1685 .LE. IF ((J+1)INDA(K)) THEN
1686  ASP=ASP+(J+1)*AC(K,J+1)*ETAJ
1687  ENDIF
1688 .LE. IF ((J+1)INDB(K)) THEN
1689  BSP=BSP+(J+1)*BC(K,J+1)*ETAJ
1690  ENDIF
1691  20 CONTINUE
1692  ASP=1.D0/DF*ASP
1693  BS=BS*DF
1694  SAS=SAS+AS*F2
1695  SBS=SBS+BS*F2/DF
1696  SDS=SDS-(CHI*AS+ASP+PSI*BS)*F2
1697  ENDIF
1698 .GT..AND..LE. IF ((K0)(K6)) THEN
1699  SCS=SCS+(AS+CHI*BSO+BSPO)*F2
1700  ENDIF
1701  10 CONTINUE
1702  PHIS=DF*PHIS
1703  SBS=DF*SBS
1704  ELSE
1705 CCCC EXACT VALUES CCCCCCCCCCCCCCCCCCCCC
1706  Z2=Z*Z
1707  Z21M=1.D0-Z2
1708  PHIEX=(4.D0*PSI/Z21M)**0.25D0
1709  PHIEX2=PHIEX*PHIEX
1710  A0EX=1.D0
1711  B0EX=-5.D0/48.D0/(PSI*PSI)+PHIEX2/48.D0/PSI*
1712  + (5.D0/Z21M-3.D0)
1713  CHIEX=0.25D0/PSI*(1.D0-Z2*PHIEX2**3*0.25D0)
1714  CHI=CHIEX
1715  D0EX=-(7.D0/48.D0/PSI+PHIEX2/48.D0*
1716  + (9.D0-7.D0/Z21M))
1717 .GT. IF (PSI0.D0) THEN
1718  PSI12=SQRT(PSI)
1719  DZZ=SQRT(Z21M)
1720  ELSE
1721  PSI12=SQRT(-PSI)
1722  DZZ=SQRT(ABS(Z21M))
1723  ENDIF
1724  B0PEX=5.D0/24.D0/PSI3+PHIEX2/48.D0*((2.D0*CHI*PSI-1.D0)/
1725  + PSI2*(5.D0/Z21M-3.D0)-10.D0*Z2*PSI12/
1726  + Z21M**2/DZZ/PSI)
1727  C0EX=1.D0
1728  SAS=A0EX
1729  SBS=B0EX
1730  SCS=C0EX
1731  SDS=D0EX
1732  BS=0.D0
1733  BSP=0.D0
1734  A2K=1.D0
1735  SIG=1.D0
1736  DO 30 K=1,6
1737  BSO=BS
1738  BSPO=BSP
1739  AS=0.D0
1740  BS=0.D0
1741  ASP=0.D0
1742  BSP=0.D0
1743  A2K=A2K*A2
1744  SIG=-1.D0*SIG
1745  F2=SIG/A2K
1746 .LE. IF (K5) THEN
1747  DO 40 J=0,20
1748 .EQ. IF (J0) THEN
1749  ETAJ=1.D0
1750  ELSE
1751  ETAJ=ETAJ*ETA
1752  ENDIF
1753 .LE. IF (JINDA(K)) THEN
1754  AS=AS+AC(K,J)*ETAJ
1755  ENDIF
1756 .LE. IF (JINDB(K)) THEN
1757  BS=BS+BC(K,J)*ETAJ
1758  ENDIF
1759 .LE. IF ((J+1)INDA(K)) THEN
1760  ASP=ASP+(J+1)*AC(K,J+1)*ETAJ
1761  ENDIF
1762 .LE. IF ((J+1)INDB(K)) THEN
1763  BSP=BSP+(J+1)*BC(K,J+1)*ETAJ
1764  ENDIF
1765  40 CONTINUE
1766  BS=DF*BS
1767  ASP=ASP/DF
1768  SAS=SAS+AS*F2
1769  SBS=SBS+BS*F2
1770  SDS=SDS-(CHI*AS+ASP+PSI*BS)*F2
1771  ENDIF
1772 .GE. IF (K2) THEN
1773  SCS=SCS+(AS+CHI*BSO+BSPO)*F2
1774  ELSE
1775  SCS=SCS+(AS+CHI*B0EX+B0PEX)*F2
1776  ENDIF
1777  30 CONTINUE
1778  PHIS=PHIEX
1779  ENDIF
1780 CCCCC CALL THE AIRY ROUTINE CCCCCCCCC
1781  IFUN=1
1782  IFACA=1
1783  CALL AIZ(IFUN,IFACA,ARG,0.D0,AIR,AII,IERRO)
1784  IFUN=2
1785  IFACA=1
1786  CALL AIZ(IFUN,IFACA,ARG,0.D0,DAIR,DAII,IERRO)
1787  DKAI=FAC*PHIS*(AIR*SAS+DAIR*SBS/F4)
1788  DKAID=FACD/PHIS*(DAIR*SCS+AIR*SDS/A23)
1789  ELSE
1790  DKAI=0.D0
1791  DKAID=0.D0
1792  ENDIF
1793  RETURN
1794  END
1795 
1796  SUBROUTINE DLIA(IFAC,X,A,DLI,DLID,IERRO)
1797 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1798 CCCC THE PURPOUSE OF THE ROUTINE IS THE CALCULATION OF THE
1799 CCCC MODIFIED BESSEL FUNCTIONS L_IA(X) AND L'_ia(x),
1800 cccc WHERE i is the imaginary unit and x is the argument of the
1801 cccc functions. we will refer to a as the order of the functions.
1802 cccc
1803 cccc the routine has the option of computing scaled functions.
1804 cccc this scaling can be used to enlarge the range of
1805 cccc computation.
1806 cccc the scaled functions are defined as follows:
1807 cccc(s stands for scaled and
1808 cccc l=sqrt{x**2-a**2} + a*arcsin(a/x)):
1809 cccc
1810 cccc exp(-l)*l_ia(x), IF x >=abs(a)
1811 cccc sl_ia(x) =
1812 cccc exp(-abs(a)*pi/2)*l_ia(x), IF x < abs(a)
1813 cccc
1814 cccc
1815 cccc exp(-l)*l'_IA(X), IF X >=ABS(A)
1816 CCCC SL'_ia(x) =
1817 cccc exp(-abs(a)*pi/2)*l'_IA(X), IF X <ABS(A)
1818 CCCC
1819 CCCC THE RANGE OF THE PARAMETERS (X,A) FOR THE COMPUTATION IS:
1820 CCCC 0 < X <= 1500, -1500 <= A <= 1500
1821 CCCC FOR FUNCTIONS WITHOUT SCALING, THE RANGE IS USUALLY LARGER
1822 CCCC THAN
1823 CCCC 0 < X <= 500, -400 <= A <= 400,
1824 CCCC DEPENDING ON THE MACHINE OVERFLOW/UNDERFLOW PARAMETERS, WHICH
1825 CCCC ARE SET UP BY THE ROUTINE D1MACH.
1826 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1827 CCCC METHODS OF COMPUTATION:
1828 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1829 CCCC 1) SERIES,
1830 CCCC IF ABS(A) > 0.044*ABS(X-3.1)**1.9+(X-3.1) OR
1831 CCCC ABS(A) <= 25 AND X <= 60
1832 CCCC 2) ASYMPTOTIC EXPANSION, IF ABS(A) < 0.7*X-10
1833 CCCC 3) AIRY-TYPE ASYMPTOTIC EXPANSION,
1834 CCCC IF ABS(A) < 3.7*X-10
1835 CCCC 4) QUADRATURES,
1836 CCCC IN THE REST OF THE PLANE (X,A)
1837 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1838 CCCC DESCRIPTION OF INPUT/OUTPUT VARIABLES:
1839 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1840 CCCC INPUTS:
1841 CCCC X : ARGUMENT OF THE FUNCTIONS
1842 CCC A : ORDER OF THE FUNCTIONS
1843 CCCC IFAC: INTEGER FLAG FOR THE SCALING
1844 CCCC * IFAC=1, THE CODE COMPUTES L_IA(X), L'_ia(x)
1845 cccc * otherwise, the code computes scaled l_ia(x), l'_IA(X)
1846 CCCC OUTPUTS:
1847 CCCC DLI : LIA(X) FUNCTION
1848 CCCC DLID : DERIVATIVE OF THE LIA(X) FUNCTION
1849 CCCC IERRO: ERROR FLAG
1850 CCCC * IF IERRO=0, COMPUTATION SUCCESSFUL.
1851 CCCC * IF IERRO=1, COMPUTATION OUT OF RANGE.
1852 CCCC * IF IERRO=2, ARGUMENT X AND/OR ORDER A OUT OF RANGE.
1853 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1854 CCCC ACCURACY:
1855 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1856 CCCC THE RELATIVE ACCURACY IS:
1857 CCCC BETTER THAN 10**(-13) FOR (X,A) IN (0,200]X[-200,200];
1858 CCCC BETTER THAN 5.10**(-13) FOR (X,A) IN (0,500]X[-500,500];
1859 CCCC CLOSE TO 10**(-12) FOR (X,A) IN (0,1500]X[-1500,1500].
1860 CCCC NEAR THE ZEROS OF THE FUNCTIONS (THERE ARE INFINITELY
1861 CCCC MANY OF THEM IN THE OSCILLATORY REGION) RELATIVE PRECISION
1862 CCCC LOOSES MEANING AND ONLY ABSOLUTE PRECISION MAKES SENSE.
1863 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1864 C AUTHORS:
1865 C AMPARO GIL (U. CANTABRIA, SANTANDER, SPAIN).
1866 C E-MAIL: AMPARO.GIL@UNICAN.ES
1867 C JAVIER SEGURA (U. CANTABRIA, SANTANDER, SPAIN).
1868 C E-MAIL: SEGURAJJ@UNICAN.ES
1869 C NICO M. TEMME (CWI, AMSTERDAM, THE NETHERLANDS).
1870 C E-MAIL: NICO.TEMME@CWI.NL
1871 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1872 C REFERENCES:
1873 C THIS IS THE COMPANION SOFTWARE OF THE ARTICLES
1874 C 1)'computing solutions of the modified bessel differential
1875 c equation for imaginary orders and positive arguments',
1876 C A. GIL, J. SEGURA, N.M. TEMME
1877 C ACM TRANS. MATH. SOFT. (2004)
1878 C 2)'modified bessel functions of imaginary order and
1879 c positive argument',
1880 C A. GIL, J. SEGURA, N.M. TEMME
1881 C ACM TRANS. MATH. SOFT. (2004)
1882 C ADDITIONAL REFERENCES:
1883 C - 'computation of the modified bessel FUNCTION of the
1884 c third kind for imaginary orders'
1885 C A. GIL, J. SEGURA, N.M. TEMME
1886 C J. COMPUT. PHYS. 175 (2002) 398-411
1887 C - 'computation of the modified bessel functions of the
1888 c third kind of imaginary orders:
1889 c uniform airy-TYPE asymptotic expansion'
1890 C A. GIL, J. SEGURA, N.M. TEMME, PROCEEDINGS OPSFA 2001
1891 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1892 CCC LOCAL VARIABLES:
1893 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1894 C D: X-3.1
1895 C DF1: 0.044*ABS(D)**1.9+(X-3.1)
1896 C DF2: 0.7*X-10
1897 C DF3: 3.7*X-10
1898 C PNU: ORDER OF THE FUNCTION
1899 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1900  DOUBLE PRECISION A,D,DF1,DF2,DF3,DLI,DLID,PNU,X
1901  INTEGER IERRO,IFAC
1902  IERRO=0
1903  PNU=A
1904 .GT..OR..LE. IF ((X1500.D0)(X0.D0)) THEN
1905  IERRO=2
1906  DLI=0.D0
1907  DLID=0.D0
1908  ENDIF
1909 .GT. IF (ABS(PNU)1500.D0) THEN
1910  IERRO=2
1911  DLI=0.D0
1912  DLID=0.D0
1913  ENDIF
1914 .EQ. IF (IERRO0) THEN
1915 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1916 CCC THESE FUNCTIONS ARE EVEN FUNCTIONS IN THE C
1917 CCC PARAMETER A (=PNU) C
1918 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1919 .LT. IF (PNU0.D0) PNU=-PNU
1920  D=X-3.1D0
1921  DF1=0.044D0*ABS(D)**1.9D0+D
1922  DF2=0.7D0*X-10.D0
1923  DF3=3.7D0*X-10.D0
1924 .GT..OR..LE..AND..LE. IF ((PNUDF1)(((PNU25.D0)(X60.D0))))
1925  + THEN
1926  CALL SERLIA(IFAC,X,PNU,DLI,DLID,IERRO)
1927 .LT. ELSEIF (PNUDF2) THEN
1928  CALL EXPLIA(IFAC,X,PNU,DLI,DLID,IERRO)
1929 .LT. ELSEIF (PNUDF3) THEN
1930  CALL AIEXLI(IFAC,X,PNU,DLI,DLID,IERRO)
1931  ELSE
1932  CALL DLINT(IFAC,X,PNU,DLI,DLID,IERRO)
1933  ENDIF
1934  ENDIF
1935  RETURN
1936  END
1937 
1938  SUBROUTINE DLINT(IFAC,XX,PNUA,DLINF,DLIND,IERRO)
1939 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1940 CCC CALCULATION OF L, L' using non-oscillating integral c
1941 ccc representations c
1942 ccc c
1943 ccc input: c
1944 ccc xx: argument of the functions c
1945 ccc pnua: order of the functions c
1946 ccc ifac: c
1947 ccc =1, non scaled functions c
1948 ccc otherwise, scaled functions c
1949 ccc output: c
1950 ccc dlinf: l FUNCTION c
1951 ccc dlind: derivative of the l FUNCTION c
1952 ccc ierro: error flag c
1953 ccc * IF ierro=0, computation successful. c
1954 ccc * IF ierro=1, computation out of range. c
1955 ccc * IF ierro=2, argument and/or order, out of range. c
1956 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1957 ccc method of computation:
1958 ccc * xx<pnua: the non-oscillating integral representations
1959 ccc for the oscillatory region are used
1960 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1961 ccc local variables:
1962 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1963 c chi : x*sinh(mu)-pnu*mu
1964 c contr1 : contribution of the semi-infinite integral
1965 c in the oscillatory CASE (including additional
1966 c factors: contr1=fos1*factors).
1967 c coschi : cos(chi)
1968 c cosh2m : cosh(2*mu)
1969 c coshm : cosh(mu)
1970 c df1 : factor(depending on ifac)
1971 c dff : (1-exp(-2*pi*pnu))*exp(-eta)
1972 c dmain : pi*pnu*0.5
1973 c dmu : solution mu of cosh(mu)=pnu/x
1974 c dmu2 : 2*mu
1975 c dmu3 : 3*mu
1976 c dmu5 : 5*mu
1977 c dmu7 : 7*mu
1978 c dmu9 : 9*mu
1979 c dmufac : mu*cosh(mu)-sinh(mu)
1980 c dmutan : mu-tanh(mu)
1981 c expo : exp(pi*pnu*0.5)
1982 c fos1 : contribution of the semi-infinite integral
1983 c in the oscillatory CASE (lia(x)).
1984 c fosd1 : contribution of the semi-infinite integral
1985 c in the oscillatory CASE (lia'(X)).
1986 C OVER : OVERFLOW NUMBER
1987 C PI : 3.1415...
1988 C PINU : PI*PNU
1989 C PNUA : SAME AS PNU (INPUT VARIABLE)
1990 C SINCHI : SIN(CHI)
1991 C SINH2M : SINH(2*MU)
1992 C SINHM : SINH(MU)
1993 C UNDER : UNDERFLOW NUMBER
1994 C X : ARGUMENT OF THE FUNCTION
1995 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
1996  DOUBLE PRECISION CHI,CONTR1,COSCHI,COSH2M,COSHM,
1997  + D1MACH,DF1,DFF,DLIND,DLINF,DMAIN,DMU,DMU2,DMU3,
1998  + DMU5,DMU7,DMU9,DMUFAC,DMUTAN,EXPO,FOS1,FOSD1,OVER,
1999  + PI,PINU,PNU,PNUA,SINCHI,SINH2M,SINHM,UNDER,X,XX
2000  INTEGER IERRO,IFAC
2001  COMMON/ARGU/X,PNU
2002  COMMON/PAROS1/DMU,COSHM,SINHM,DMUFAC,COSCHI,SINCHI
2003  COMMON/PAROS2/COSH2M,SINH2M,DMAIN
2004  COMMON/PAROS3/DMUTAN
2005 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
2006  OVER=D1MACH(2)*1.D-8
2007 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
2008  UNDER=D1MACH(1)*1.D+8
2009  X=XX
2010  PNU=PNUA
2011  PI=ACOS(-1.D0)
2012  IERRO=0
2013 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2014 CC THE USE OF INTEGRALS FOR THE LIA FUNCTION IS C
2015 CC RESTRICTED TO THE CASE X < A, WHICH IS THE C
2016 CC OSCILLATORY CASE C
2017 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2018  PI=ACOS(-1.D0)
2019 .EQ. IF (IFAC1) THEN
2020 .GE. IF ((PI*PNU*0.5D0)LOG(OVER)) IERRO=1
2021  ENDIF
2022 .EQ. IF (IERRO0) THEN
2023  COSHM=PNU/X
2024  DMU=LOG(COSHM+SQRT((COSHM-1.D0)*(COSHM+1.D0)))
2025  DMU2=2.D0*DMU
2026  COSH2M=COSH(DMU2)
2027  SINHM=SINH(DMU)
2028  SINH2M=SINH(DMU2)
2029 .GT. IF (DMU0.1D0) THEN
2030  CHI=X*SINHM-PNU*DMU
2031  DMUFAC=DMU*COSHM-SINHM
2032  ELSE
2033  DMU2=DMU*DMU
2034  DMU3=DMU2*DMU
2035  DMU5=DMU3*DMU2
2036  DMU7=DMU5*DMU2
2037  DMU9=DMU7*DMU2
2038  CHI=-2.D0*X*(1.D0/6.D0*DMU3+1.D0/60.D0*DMU5+
2039  + 3.D0/5040.D0*DMU7+4.D0/362880.D0*DMU9)
2040  DMUFAC=DMU3/3.D0+DMU5/30.D0+DMU7/840.D0+DMU9/45360.D0
2041  ENDIF
2042  DMUTAN=DMU-TANH(DMU)
2043  COSCHI=COS(CHI)
2044  SINCHI=SIN(CHI)
2045  PINU=PI*PNU
2046  DMAIN=PINU*0.5D0
2047 CCCCCCCCCCCCCCCCCC
2048 CCC INTEGRALS CC
2049 CCCCCCCCCCCCCCCCCC
2050  CALL TRAPRL(1,FOS1)
2051 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2052 CCC THEN, THE LIA FUNCTION IS GIVEN BY ... CC
2053 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2054 .EQ. IF (IFAC1) THEN
2055  EXPO=EXP(DMAIN)
2056  DFF=EXPO/PI
2057 .LE. IF ((-2.D0*PINU)LOG(UNDER)) THEN
2058  DF1=DFF*0.5D0
2059  ELSE
2060  DF1=DFF*(1.D0-EXP(-2.D0*PINU))*0.5D0
2061  ENDIF
2062  CONTR1=DF1*FOS1
2063  ELSE
2064  DFF=1.D0/PI
2065 .LE. IF ((-2.D0*PINU)LOG(UNDER)) THEN
2066  DF1=DFF*0.5D0
2067  ELSE
2068  DF1=DFF*(1.D0-EXP(-2.D0*PINU))*0.5D0
2069  ENDIF
2070  CONTR1=DF1*FOS1
2071  ENDIF
2072  DLINF=CONTR1
2073 CCCCCCCCCCCCCCCCCCCCCCCCCC
2074 CCC CALCULATION OF LIA' cc
2075 cccccccccccccccccccccccccc
2076 cccccccccccccccccc
2077 ccc integrals cc
2078 cccccccccccccccccc
2079  CALL traprl(10,fosd1)
2080 ccccccccccccccccccccccccccccccccccccccccccccc
2081 ccc THEN, the lia FUNCTION is given by ... cc
2082 ccccccccccccccccccccccccccccccccccccccccccccc
2083  IF (ifac.EQ.1) THEN
2084  expo=exp(dmain)
2085  dff=expo/pi
2086  IF ((-2.d0*pinu).LE.log(under)) THEN
2087  df1=dff*0.5d0
2088  ELSE
2089  df1=dff*(1.d0-exp(-2.d0*pinu))*0.5d0
2090  ENDIF
2091  contr1=df1*fosd1
2092  ELSE
2093  dff=1.d0/pi
2094  IF ((-2.d0*pinu).LE.log(under)) THEN
2095  df1=dff*0.5d0
2096  ELSE
2097  df1=dff*(1.d0-exp(-2.d0*pinu))*0.5d0
2098  ENDIF
2099  contr1=df1*fosd1
2100  ENDIF
2101  dlind=contr1
2102  ELSE
2103  dlinf=0.d0
2104  dlind=0.d0
2105  ENDIF
2106  RETURN
2107  END
2108 
2109  DOUBLE PRECISION FUNCTION fstal2(T)
2110 ccccccccccccccccccccccccccccccccccccccccccccccccc
2111 ccc computation of the integrand for the
2112 ccc integral for the l FUNCTION in eq.(37) OF
2113 ccc the reference:
2114 ccc 'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL
2115 CCC DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS
2116 CCC AND POSITIVE ARGUMENTS',
2117 ccc a. gil, j. segura, n.m. temme
2118 ccc acm trans. math. soft. (2004)
2119 cccccccccccccccccccccccccccccccccccccccccccccccccccc
2120 ccc input variable:
2121 cccccccccccccccccccccccccccccccccccccccccccccccccccc
2122 c t : argument of the FUNCTION
2123 c(variable of integration).
2124 ccccccccccccccccccccccccccccccccccccccccccccccccc
2125 ccc local variables:
2126 ccccccccccccccccccccccccccccccccccccccccccccccccc
2127 c argu : (cosh(mu)*u-dmufac)/sinh(u)
2128 c argu2 : argu**2
2129 c coschi : cosh(chi), chi=x*sinh(mu)-pnu*mu
2130 c cosh2m : cosh(2*mu)
2131 c coshm : cosh(mu)
2132 c d1 : cosh(mu)-dmufac
2133 c deri : (cosh(mu)/sinh(u)-d1*cosh(u)/sinh(u)**2)
2134 c /sqrt(1-argu2)
2135 c djaco : cosh(t)*exp(s)/sqrt(1+exp(s)**2)
2136 c dmain : pi*pnu*0.5
2137 c dmu : solution mu of cosh(mu)=pnu/x
2138 c dmu2 : 2*mu
2139 c dmufac : mu*cosh(mu)-sinh(mu)
2140 c expon : exp(-(phib(u)-pi*pnu*0.5))
2141 c f1 : sinh(u-mu)/(u-mu)
2142 c g1 : z/6+z3/120+z5/5040+z7/362880, z=2*y
2143 c phib(u) : x*cosh(u)*cos(sigma(u))+pnu*sigma(u),
2144 c WHERE x=argument of the function,
2145 c sigma(u)=arcsin((cosh(mu)*u-dmufac)/sinh(u)
2146 c resto : sin(chi)-cos(chi)*deri
2147 c s : sinh(t)
2148 c sinchi : sin(chi)
2149 c sinh2m : sinh(2*mu)
2150 c sinhm : sinh(mu)
2151 c sinhu : sinh(u)
2152 c sinhu2 : 2*sinh(u)
2153 c u : mu+log(x+sqrt(x**2+1))
2154 c under : underflow number
2155 c x : exp(s)
2156 c y : u-mu
2157 c z : 2*y
2158 c z2 : z**2
2159 c z3 : z**3
2160 c z5 : z**5
2161 c z7 : z**7
2162 cccccccccccccccccccccccccccccccccccccccccccccccccccccc
2163  DOUBLE PRECISION argu,argu2,coschi,cosh2m,coshm,
2164  + d1,d1mach,deri,djaco,dmain,dmu,dmu2,dmufac,dmutan,
2165  + expon,f1,g1,phib,resto,s,sinchi,sinh2m,sinhm,
2166  + sinhu,sinhu2,t,u,under,x,y,z,z2,z3,z5,z7
2167  common/paros1/dmu,coshm,sinhm,dmufac,coschi,sinchi
2168  common/paros2/cosh2m,sinh2m,dmain
2169  common/paros3/dmutan
2170 ccc machine dependent constant(underflow number)
2171  under=d1mach(1)*1.d+8
2172  s=sinh(t)
2173  x=exp(s)
2174  u=dmutan+log(x+sqrt(x*x+1.d0))
2175  y=u-dmu
2176  IF (abs(y).LE.0.1d0) THEN
2177  IF (abs(y).GT.under) THEN
2178  f1=sinh(y)/y
2179  ELSE
2180  f1=1.d0
2181  ENDIF
2182  z=y*2.d0
2183  z2=z*z
2184  z3=z2*z
2185  z5=z3*z2
2186  z7=z5*z2
2187  g1=z/6.d0+z3/120.d0+z5/5040.d0+z7/362880.d0
2188  dmu2=2.d0*dmu
2189  deri=sinh(u)/(f1-coshm*cosh(u))*
2190  + sqrt(cosh(dmu2)*f1*f1+2.d0*sinh(dmu2)*g1-
2191  + coshm*coshm)
2192  deri=1.d0/deri
2193  ELSE
2194  d1=coshm*u-dmufac
2195  sinhu=sinh(u)
2196  argu=d1/sinhu
2197  argu2=argu*argu
2198  sinhu2=sinhu*sinhu
2199  deri=1.d0/sqrt(1.d0-argu2)*
2200  + (coshm/sinhu-d1*cosh(u)/sinhu2)
2201  IF (u.LT.dmu) deri=-deri
2202  ENDIF
2203  djaco=cosh(t)*x/sqrt(1.d0+x*x)
2204  IF ((-(phib(u)-dmain)).LE.log(under)) THEN
2205  expon=0.d0
2206  fstal2=0.d0
2207  ELSE
2208  expon=exp(-(phib(u)-dmain))
2209  resto=(sinchi-coschi*deri)
2210  fstal2=expon*resto*djaco
2211  ENDIF
2212  RETURN
2213  END
2214 
2215  DOUBLE PRECISION FUNCTION fdtal2(T)
2216 cccccccccccccccccccccccccccccccccccccccccccccccccccc
2217 ccc computation of the integrand for the
2218 ccc integral for the derivative of the l function
2219 ccc in eq.(37) of the reference:
2220 ccc 'COMPUTING SOLUTIONS OF THE MODIFIED BESSEL
2221 CCC DIFFERENTIAL EQUATION FOR IMAGINARY ORDERS
2222 CCC AND POSITIVE ARGUMENTS',
2223 ccc a. gil, j. segura, n.m. temme
2224 ccc acm trans. math. soft. (2004)
2225 cccccccccccccccccccccccccccccccccccccccccccccccccccc
2226 ccc input variable:
2227 cccccccccccccccccccccccccccccccccccccccccccccccccccc
2228 c t : argument of the FUNCTION
2229 c(variable of integration).
2230 cccccccccccccccccccccccccccccccccccccccccccccccccccc
2231 ccc local variables:
2232 cccccccccccccccccccccccccccccccccccccccccccccccccccc
2233 c argu : (cosh(mu)*u-dmufac)/sinh(u)
2234 c argu2 : argu**2
2235 c coschi : cosh(chi), chi=x*sinh(mu)-pnu*mu
2236 c cosh2m : cosh(2*mu)
2237 c coshm : cosh(mu)
2238 c coshu : cosh(u)
2239 c coss : cos(sigma(u))
2240 c d1 : cosh(mu)-dmufac
2241 c delta : -sin(chi)*cos(sigma(u))*cosh(u)
2242 c +cos(chi)*sin(sigma(u))*sinh(u)
2243 c deri : (cosh(mu)/sinh(u)-d1*cosh(u)/sinh(u)**2)
2244 c /sqrt(1-argu2)
2245 c djaco : cosh(t)*exp(s)/sqrt(1+exp(s)**2)
2246 c dmain : pi*pnu*0.5
2247 c dmu : solution mu of cosh(mu)=pnu/x
2248 c dmu2 : 2*mu
2249 c dmufac : mu*cosh(mu)-sinh(mu)
2250 c expon : exp(-(phib(u)-pi*pnu*0.5))
2251 c f1 : sinh(u-mu)/(u-mu)
2252 c g1 : z/6+z3/120+z5/5040+z7/362880, z=2*y
2253 c gamma : sin(chi)*sin(sigma(u))*sinh(u)+
2254 c cos(chi)*cos(sigma(u))*cosh(u)
2255 c phib(u) : x*cosh(u)*cos(sigma(u))+pnu*sigma(u),
2256 c WHERE x=argument of the function,
2257 c sigma(u)=arcsin((cosh(mu)*u-dmufac)/sinh(u)
2258 c resto : -sin(chi)*cos(sigma(u))*cosh(u)
2259 c +cos(chi)*sin(sigma(u))*sinh(u)+
2260 c(sin(chi)*sin(sigma(u))*sinh(u)+
2261 c cos(chi)*cos(sigma(u))*cosh(u))*deri
2262 c s : sinh(t)
2263 c sigma(u): asin((cosh(mu)*u-dmufac)/sinh(u))
2264 c sigmau : sigma(u)
2265 c sinchi : sin(chi)
2266 c sinh2m : sinh(2*mu)
2267 c sinhm : sinh(mu)
2268 c sinhu : sinh(u)
2269 c sinhu2 : 2*sinh(u)
2270 c sins : sin(sigma(u))
2271 c u : mu+log(x+sqrt(x**2+1))
2272 c under : underflow number
2273 c x : exp(s)
2274 c y : u-mu
2275 c z : 2*y
2276 c z2 : z**2
2277 c z3 : z**3
2278 c z5 : z**5
2279 c z7 : z**7
2280 cccccccccccccccccccccccccccccccccccccccccccccccccccc
2281  DOUBLE PRECISION argu,argu2,coschi,cosh2m,coshm,
2282  + coshu,coss,d1,d1mach,delta,deri,djaco,dmain,dmu,
2283  + dmu2,dmufac,dmutan,expon,f1,g1,gamma,phib,resto,
2284  + s,sigma,sigmau,sinchi,sinh2m,sinhm,sinhu,sinhu2,
2285  + sins,t,u,under,x,y,z,z2,z3,z5,z7
2286  common/paros1/dmu,coshm,sinhm,dmufac,coschi,sinchi
2287  common/paros2/cosh2m,sinh2m,dmain
2288  common/paros3/dmutan
2289 ccc machine dependent constant(underflow number)
2290  under=d1mach(1)*1.d+8
2291  s=sinh(t)
2292  x=exp(s)
2293  u=dmutan+log(x+sqrt(x*x+1.d0))
2294  y=u-dmu
2295  coshu=cosh(u)
2296  sinhu=sinh(u)
2297  IF (abs(y).LE.0.1) THEN
2298  IF (abs(y).GT.under) THEN
2299  f1=sinh(y)/y
2300  ELSE
2301  f1=1.d0
2302  ENDIF
2303  z=y*2.d0
2304  z2=z*z
2305  z3=z2*z
2306  z5=z3*z2
2307  z7=z5*z2
2308  g1=z/6.d0+z3/120.d0+z5/5040.d0+z7/362880.d0
2309  dmu2=2.d0*dmu
2310  deri=sinhu/(f1-coshm*coshu)*sqrt(cosh(dmu2)*f1*f1+
2311  + 2.d0*sinh(dmu2)*g1-coshm*coshm)
2312  deri=1.d0/deri
2313  ELSE
2314  d1=coshm*u-dmufac
2315  argu=d1/sinhu
2316  argu2=argu*argu
2317  sinhu2=sinhu*sinhu
2318  deri=1.d0/sqrt(1.d0-argu2)*(coshm/sinhu-
2319  + d1*coshu/sinhu2)
2320  IF (u.LT.dmu) deri=-deri
2321  ENDIF
2322  djaco=cosh(t)*x/sqrt(1.d0+x*x)
2323  IF ((-(phib(u)-dmain)).LE.log(under)) THEN
2324  expon=0.d0
2325  fdtal2=0.d0
2326  ELSE
2327  expon=exp(-(phib(u)-dmain))
2328  sigmau=sigma(u)
2329  sins=sin(sigmau)
2330  coss=cos(sigmau)
2331  gamma=sinchi*sins*sinhu+coschi*coss*coshu
2332  delta=-sinchi*coss*coshu+coschi*sins*sinhu
2333  resto=delta+gamma*deri
2334  fdtal2=expon*resto*djaco
2335  ENDIF
2336  RETURN
2337  END
2338 
2339  SUBROUTINE traprl(IC,TI)
2340 cccccccccccccccccccccccccccccccccccccccccccccccccccc
2341 ccc implementation of an adaptive trapezoidal rule
2342 ccc for computing the integral representations of
2343 ccc the functions.
2344 ccc
2345 ccc input variable:
2346 ccc ic: depending on the values of ic,
2347 ccc different integrals are computed:
2348 ccc *ic=1, l function, oscillatory region
2349 ccc *ic=10, derivative of the l function,
2350 ccc oscillatory region.
2351 ccc output variable:
2352 ccc ti, result of the integral
2353 cccccccccccccccccccccccccccccccccccccccccccccccccccc
2354 ccc local variables
2355 cccccccccccccccccccccccccccccccccccccccccccccccccccc
2356 c a : lower integration limit
2357 c b : upper integration limit
2358 c delta : calculates the relative precision
2359 c eps : relative precision PARAMETER used in
2360 c the calculation
2361 c h : integration step
2362 c pnu : order of the FUNCTION
2363 c sum : accumulates the elementary
2364 c contributions
2365 c tin : evaluated integral
2366 c x : argument of the FUNCTION
2367 c xac : integration abcissa
2368 ccccccccccccccccccccccccccccccccccccccccccccccccc
2369  DOUBLE PRECISION a,b,d1mach,delta,eps,fintel,
2370  + h,pnu,sum,ti,tin,x,xac,z
2371  INTEGER i,ic,ifi,n
2372  common/argu/x,pnu
2373  eps=d1mach(3)
2374  IF (eps.LT.1.d-14) eps=1.d-14
2375  n=0
2376 ccccc integration limits: a,b
2377  z=x/pnu
2378  IF ((z.GT.0.999d0).AND.(z.LT.1.001d0)) THEN
2379  a=-4.5
2380  ELSE
2381  a=-2.5d0
2382  ENDIF
2383  b=-a
2384  h=b-a
2385  ti=0.5d0*h*(fintel(ic,a)+fintel(ic,b))
2386  delta=1.d0+eps
2387 11 n=n+1
2388  h=0.5d0*h
2389  IF (n.EQ.1) THEN
2390  ifi=1
2391  ELSE
2392  ifi=2*ifi
2393  ENDIF
2394  sum=0.d0
2395  DO 3 i=1,ifi
2396  xac=a+dble(2*i-1)*h
2397  sum=sum+fintel(ic,xac)
2398  3 CONTINUE
2399  tin=0.5d0*ti+h*sum
2400  IF ((tin.NE.0.d0).AND.(n.GT.4)) THEN
2401  delta=abs(1.d0-ti/tin)
2402  ENDIF
2403  ti=tin
2404  IF ((delta.GT.eps).AND.(n.LT.9)) goto 11
2405  RETURN
2406  END
2407 
2408  DOUBLE PRECISION FUNCTION fintel(IC,T)
2409 ccccccccccccccccccccccccccccccccccccccccccccccccccc
2410 ccc auxiliary FUNCTION for the subroutine traprl.
2411 ccc this FUNCTION calls the different functions
2412 ccc containing the integrands which are integrated
2413 ccc by traprl.
2414 ccccccccccccccccccccccccccccccccccccccccccccccccccc
2415 ccc input: c
2416 ccc ic: INTEGER parameter. its admissible c
2417 ccc values are the same as in the c
2418 ccc SUBROUTINE traprl. c
2419 ccc t: integration abscissa c
2420 ccccccccccccccccccccccccccccccccccccccccccccccccccc
2421 ccc local variables:
2422 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2423 c fdtal2 : oscillatory part: tau contribution routine
2424 c(lia'(X) FUNCTION). SEMI-INFINITE INTEGRAL.
2425 C USED FOR LARGE PNU
2426 C FSTAL2 : OSCILLATORY PART: TAU CONTRIBUTION ROUTINE
2427 C (LIA(X) FUNCTION). SEMI-INFINITE INTEGRAL.
2428 C USED FOR LARGE PNU
2429 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2430  DOUBLE PRECISION FDTAL2,FSTAL2,T
2431  INTEGER IC
2432 .EQ. IF (IC1) THEN
2433  FINTEL=FSTAL2(T)
2434  ENDIF
2435 .EQ. IF (IC10) THEN
2436  FINTEL=FDTAL2(T)
2437  ENDIF
2438  RETURN
2439  END
2440 
2441  SUBROUTINE SERLIA(IFAC,X,PNU,PSER,PSERD,IERRO)
2442 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2443 CCC CALCULATION OF POWER SERIES FOR L, L'
2444 ccccccccccccccccccccccccccccccccccccccccccccccccccccc
2445 ccc input variables: c
2446 ccc x: argument of the functions c
2447 ccc pnu: order of the functions c
2448 ccc ifac: c
2449 ccc =1, non scaled functions c
2450 ccc otherwise, scaled functions c
2451 ccc output variables: c
2452 ccc pser: l FUNCTION c
2453 ccc pserd: derivative of the l FUNCTION c
2454 ccc ierro: error flag c
2455 ccc * IF ierro=0, computation successful. c
2456 ccc * IF ierro=1, computation out of range. c
2457 ccc * IF ierro=2, argument and/or order, c
2458 ccc out of range. c
2459 ccccccccccccccccccccccccccccccccccccccccccccccccccccc
2460 ccc local variables:
2461 ccccccccccccccccccccccccccccccccccccccccccccccccccccc
2462 c a : order of the functions
2463 c a2 : a**2
2464 c a2h : a**2/2
2465 c a2n : a**(2*n)
2466 c accp : accumulates the p coefficients
2467 c accq : accumulates the q coefficients
2468 c argu : sig(0)-a*log(x/2)
2469 c c : x**(2*k)/k
2470 c coci : 1/(k**2+a**2)
2471 c costh : cos(thet)
2472 c deltak : accumulates the contribution for the
2473 c lia(x) function
2474 c deltkp : accumulates the contribution for the
2475 c lia'(X) FUNCTION
2476 C DF1 : FACTOR (DEPENDING ON IFAC)
2477 C ETA0 : PARAMETER FOR THE CALCULATION OF THE
2478 C COULOMB PHASE SHIFT
2479 C ETA02 : ETA0**2
2480 C F(K) : SIN(PHI(A,K)-A*LOG(X/2))
2481 C /(A**2*(1+A**2)...(K**2+A**2))**1/2,
2482 C WHERE PHI(A,K)=PHASE(GAMMA(1+K+IA))
2483 C FDOMIN : X*(COS(THET)+THET*SIN(THET))
2484 C K : CONTRIBUTION TO THE LIA(X) FUNCTION
2485 C KP : CONTRIBUTION TO THE LIA'(x) function
2486 c over : overflow number
2487 c p0 : parameters for the calculation of the coulomb
2488 c phase shift
2489 c p1 : parameters for the calculation of the coulomb
2490 c phase shift
2491 c p2 : parameters for the calculation of the coulomb
2492 c phase shift
2493 c pi : 3.1415...
2494 c pia : pi*a
2495 c pia2 : 2*pi*a
2496 c preci : relative precision used in the calculation
2497 c q0 : parameters for the calculation of the coulomb
2498 c phase shift
2499 c q1 : parameters for the calculation of the coulomb
2500 c phase shift
2501 c q2 : parameters for the calculation of the coulomb
2502 c phase shift
2503 c r(k) : f(k)*a/tan(phi(a,k)-a*log(x/2))
2504 c sig0 : coulomb phase shift
2505 c sinth : sin(thet)
2506 c thet : asin(a/x)
2507 c under : underflow number
2508 c x2 : x*x
2509 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2510  DOUBLE PRECISION a,a2,a2h,a2n,accp,accq,argu,c,coci,
2511  + costh,d1mach,dds,dee,deltak,deltkp,df1,dsmall,eta0,
2512  + eta02,euler,f(0:500),fdomin,k,kp,over,p0(0:9),p1(0:9),
2513  + p2(0:6),pi,pia,pia2,pnu,preci,pser,pserd,q0(0:9),
2514  + q1(0:9),q2(0:6),r(0:500),sig0,sinth,thet,under,x,x2
2515  INTEGER ierro,ifac,l,m,n
2516  SAVE p0,q0,p1,q1,p2,q2
2517  DATA p0/1.08871504904797411683d5,3.64707573081160914640d5,
2518  + 4.88801471582878013158d5,3.36275736298197324009d5,
2519  + 1.26899226277838479804d5,
2520  + 2.60795543527084582682d4,2.73352480554497990544d3,
2521  + 1.26447543569902963184d2,
2522  + 1.85446022125533909390d0,1.90716219990037648146d-3/
2523  DATA q0/6.14884786346071135090d5,2.29801588515708014282d6,
2524  + 3.50310844128424021934d6,
2525  + 2.81194990286041080264d6,1.28236441994358406742d6,
2526  + 3.35209348711803753154d5,
2527  + 4.84319580247948701171d4,3.54877039006873206531d3,
2528  + 1.11207201299804390166d2,1.d0/
2529  DATA p1/-1.044100987526487618670d10,-1.508574107180079913696d10,
2530  + -5.582652833355901160542d9,4.052529174369477275446d8,
2531  + 5.461712273118594275192d8,
2532  + 9.510404403068169395714d7,6.281126609997342119416d6,
2533  + 1.651178048950518520416d5,
2534  + 1.498824421329341285521d3,2.974686506595477984776d0/
2535  DATA q1/1.808868161493543887787d10,3.869142051704700267785d10,
2536  + 3.003264575147162634046d10,1.075554651494601843525d10,
2537  + 1.901298501823290694245d9,
2538  + 1.665999832151229472632d8,6.952188089169487375936d6,
2539  + 1.253235080625688652718d5,7.904420414560291396996d2,1.d0/
2540  DATA p2/7.08638611024520906826d-3,-6.54026368947801591128d-2,
2541  + 2.92684143106158043933d-1,4.66821392319665609167d0,
2542  + -3.43943790382690949054d0,
2543  + -7.72786486869252994370d0,-9.88841771200290647461d-01/
2544  DATA q2/-7.08638611024520908189d-3,6.59931690706339630254d-2,
2545  + -2.98754421632058618922d-1,-4.63752355513412248006d0,
2546  + 3.79700454098863541593d0,7.06184065426336718524d0,1.d0/
2547 ccc machine dependent constant(overflow number)
2548  over=d1mach(2)*1.d-8
2549 ccc machine dependent constant(underflow number)
2550  under=d1mach(1)*1.d+8
2551  pi=acos(-1.d0)
2552  ierro=0
2553  IF (ifac.EQ.1) THEN
2554  IF (x.GE.pnu) THEN
2555  sinth=pnu/x
2556  thet=asin(sinth)
2557  costh=cos(thet)
2558  fdomin=x*(costh+thet*sinth)
2559  IF (fdomin.GE.log(over)) ierro=1
2560  ELSE
2561  IF ((pi*pnu*0.5d0).GE.log(over)) ierro=1
2562  ENDIF
2563  ENDIF
2564  IF (ierro.EQ.0) THEN
2565  eta0=1.8055470716051069198764d0
2566  euler=0.577215664901532860606512d0
2567  preci=d1mach(3)*10
2568 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2569 ccc coefficients for the calculation of the coulomb phase shift
2570 ccc from cody & hillstrom, math. comput. 24(111) 1970
2571 cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2572  a=pnu
2573  eta02=eta0*eta0
2574  a2=a*a
2575  n=0
2576  accp=0.d0
2577  accq=0.d0
2578  IF (a.LE.2.d0) THEN
2579  33 a2n=a2**n
2580  accp=accp+p0(n)*a2n
2581  accq=accq+q0(n)*a2n
2582  n=n+1
2583  IF (n.LE.9) goto 33
2584  sig0=a*(a2-eta02)*accp/accq
2585  ELSE
2586  IF ((a.GT.2.d0).AND.(a.LE.4.d0)) THEN
2587  44 a2n=a2**n
2588  accp=accp+p1(n)*a2n
2589  accq=accq+q1(n)*a2n
2590  n=n+1
2591  IF (n.LE.9) goto 44
2592  sig0=a*accp/accq
2593  ELSE
2594  55 a2n=a2**n
2595  accp=accp+p2(n)/a2n
2596  accq=accq+q2(n)/a2n
2597  n=n+1
2598  IF (n.LE.6) goto 55
2599  sig0=atan(a)*0.5d0+a*(log(1.d0+a2)*0.5d0+accp/accq)
2600  ENDIF
2601  ENDIF
2602 ccccccccccccccccccccccccccccccccccccc
2603 ccc evaluation of f(0), r(0), r(1)
2604 ccccccccccccccccccccccccccccccccccccc
2605  pia=pi*a
2606  a2h=a2*0.5d0
2607  argu=sig0-a*log(x*0.5d0)
2608  r(0)=cos(argu)
2609  r(1)=1.d0/(1.d0+a2)*(r(0)-a*sin(argu))
2610  IF (a.LT.under) THEN
2611  f(0)=-(euler+log(x*0.5d0))
2612  ELSE
2613  f(0)=1.d0/a*sin(argu)
2614  ENDIF
2615  c=1.d0
2616  k=r(0)
2617  kp=a2h*f(0)
2618 cccccccccccccccccccccccccccccccccccc
2619 ccc recursion for f(k), r(k), c(k)
2620 cccccccccccccccccccccccccccccccccccc
2621  x2=0.25d0*x*x
2622  m=1
2623  coci=1.d0/(m*m+a2)
2624  deltak=k*10
2625  deltkp=kp*10
2626  66 f(m)=(m*f(m-1)+r(m-1))*coci
2627  c=x2*c/m
2628  deltak=r(m)*c
2629  k=k+deltak
2630  deltkp=(m*r(m)+a2h*f(m))*c
2631  kp=kp+deltkp
2632  m=m+1
2633  coci=1.d0/(m*m+a2)
2634  r(m)=((2.d0*m-1.d0)*r(m-1)-r(m-2))*coci
2635  IF (m.LT.500) THEN
2636  IF ((abs(deltak/k).GT.preci).OR.
2637  + (abs(deltkp/kp).GT.preci)) goto 66
2638  ENDIF
2639  pia2=2.d0*pia
2640  IF (ifac.EQ.1) THEN
2641  IF (-pia2.LE.log(under)) THEN
2642  dee=0.d0
2643  ELSE
2644  dee=exp(-pia2)
2645  ENDIF
2646  IF (a.LT.1.d-1) THEN
2647  IF (a.LT.under) THEN
2648  df1=1.d0
2649  ELSE
2650  l=0
2651  dds=1.d0
2652  dsmall=1.d0
2653  47 l=l+1
2654  dds=-dds*pia2/(l+1.d0)
2655  dsmall=dsmall+dds
2656  IF (abs(dds/dsmall).GT.preci) goto 47
2657  df1=exp(pia*0.5d0)*sqrt(dsmall)
2658  ENDIF
2659  ELSE
2660  df1=exp(pia*0.5d0)*sqrt((1.d0-dee)/pia2)
2661  ENDIF
2662  pser=k*df1
2663  pserd=df1*kp*2.d0/x
2664  ELSE
2665  IF (x.LT.a) THEN
2666  IF (-pia2.LE.log(under)) THEN
2667  dee=0.d0
2668  ELSE
2669  dee=exp(-pia2)
2670  ENDIF
2671  IF (a.LT.1.d-1) THEN
2672  IF (a.LT.under) THEN
2673  df1=1.d0
2674  ELSE
2675  l=0
2676  dds=1.d0
2677  dsmall=1.d0
2678  48 l=l+1
2679  dds=-dds*pia2/(l+1.d0)
2680  dsmall=dsmall+dds
2681  IF (abs(dds/dsmall).GT.preci) goto 48
2682  df1=sqrt(dsmall)
2683  ENDIF
2684  ELSE
2685  df1=sqrt((1.d0-dee)/pia2)
2686  ENDIF
2687  ELSE
2688  sinth=a/x
2689  thet=asin(sinth)
2690  costh=cos(thet)
2691  fdomin=x*(costh+thet*sinth)
2692  IF (-pia2.LE.log(under)) THEN
2693  dee=0.d0
2694  ELSE
2695  dee=exp(-pia2)
2696  ENDIF
2697  IF (a.LT.1.d-1) THEN
2698  IF (a.LT.under) THEN
2699  df1=1.d0
2700  ELSE
2701  l=0
2702  dds=1.d0
2703  dsmall=1.d0
2704  49 l=l+1
2705  dds=-dds*pia2/(l+1.d0)
2706  dsmall=dsmall+dds
2707  IF (abs(dds/dsmall).GT.preci) goto 49
2708  df1=exp(pia*0.5d0-fdomin)*sqrt(dsmall)
2709  ENDIF
2710  ELSE
2711  df1=exp(pia*0.5d0-fdomin)*sqrt((1.d0-dee)/pia2)
2712  ENDIF
2713  ENDIF
2714  pser=k*df1
2715  pserd=df1*kp*2.d0/x
2716  ENDIF
2717  ELSE
2718  pser=0.d0
2719  pserd=0.d0
2720  ENDIF
2721  RETURN
2722  END
2723 
2724  SUBROUTINE explia(IFAC,X,PNU,PEXP,PEXPP,IERRO)
2725 cccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2726 ccc calculation of the asymptotic expansions for c
2727 ccc the l, l' FUNCTIONS C
2728 CCC C
2729 CCC INPUT VARIABLES: C
2730 CCC X: ARGUMENT OF THE FUNCTIONS C
2731 CCC PNU: ORDER OF THE FUNCTIONS C
2732 CCC IFAC: C
2733 CCC =1, NON SCALED FUNCTIONS C
2734 CCC OTHERWISE, SCALED FUNCTIONS C
2735 CCC OUTPUT VARIABLES: C
2736 CCC PEXP: L FUNCTION C
2737 CCC PEXPP: DERIVATIVE OF THE L FUNCTION C
2738 CCC IERRO: ERROR FLAG C
2739 CCC * IF IERRO=0, COMPUTATION SUCCESSFUL. C
2740 CCC * IF IERRO=1, COMPUTATION OUT OF RANGE . C
2741 CCC * IF IERRO=2, ARGUMENT AND/OR ORDER, C
2742 CCC OUT OF RANGE. C
2743 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2744 CCC LOCAL VARIABLES:
2745 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2746 C A : ORDER OF THE FUNCTIONS
2747 C A2 : A**2
2748 C COSTH : COS(THET)
2749 C DELTAK : ACCUMULATES THE CONTRIBUTION FOR THE
2750 C LIA(X) FUNCTION
2751 C DELTKP : ACCUMULATES THE CONTRIBUTION FOR THE
2752 C LIA'(x) function
2753 c df : factor(depending on ifac)
2754 c fdomin : x*(cos(thet)+thet*sin(thet))
2755 c ind : (-1)**m
2756 c k : contribution to the lia(x) function
2757 c kp : contribution to the lia'(X) FUNCTION
2758 C OVER : OVERFLOW NUMBER
2759 C PI : 3.1415...
2760 C PIA : PI*A
2761 C PRECI : RELATIVE PRECISION USED IN THE CALCULATION
2762 C RAX(K+1) : -((K+1/2)**2+A**2)/(K+1)*RAX(K)
2763 C SINTH : SIN(THET)
2764 C THET : ASIN(PNU/X)
2765 C XI : 1/X
2766 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2767  DOUBLE PRECISION A,A2,COSTH,D1MACH,DELTAK,
2768  + DELTKP,DF,FDOMIN,K,KP,OVER,PEXP,PEXPP,PI,PIA,
2769  + PNU,PRECI,RAX,SINTH,THET,X,XI
2770  INTEGER IERRO,IFAC,IND,M
2771 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
2772  OVER=D1MACH(2)*1.D-8
2773  PI=ACOS(-1.D0)
2774  IERRO=0
2775 .EQ. IF (IFAC1) THEN
2776 .GE. IF (XPNU) THEN
2777  SINTH=PNU/X
2778  THET=ASIN(SINTH)
2779  COSTH=COS(THET)
2780  FDOMIN=X*(COSTH+THET*SINTH)
2781 .GE. IF (FDOMINLOG(OVER)) IERRO=1
2782  ELSE
2783 .GE. IF ((PI*PNU*0.5D0)LOG(OVER)) IERRO=1
2784  ENDIF
2785  ENDIF
2786 .EQ. IF (IERRO0) THEN
2787  PRECI=D1MACH(3)*10
2788  A=PNU
2789  XI=1.D0/X
2790  A2=A*A
2791  RAX=1.D0
2792  M=1
2793  DELTAK=RAX
2794  K=DELTAK
2795  DELTKP=RAX*(0.5D0*XI-1.D0)
2796  KP=DELTKP
2797  IND=-1
2798  44 RAX=-((M-0.5D0)*(1.D0-0.5D0/M)+A2/M)*0.5D0/X*RAX
2799  DELTAK=RAX*IND
2800  K=K+DELTAK
2801 .GE. IF (KOVER) M=2000
2802  M=M+1
2803  IND=-IND
2804 .GT..AND..LT. IF ((ABS(DELTAK/K)PRECI)(M1000)) GOTO 44
2805  RAX=1.D0
2806  M=1
2807  IND=-1
2808  45 RAX=-((M-0.5D0)*(1.D0-0.5D0/M)+A2/M)*0.5D0/X*RAX
2809  DELTAK=RAX*IND
2810  DELTKP=DELTAK*(XI*(0.5D0+M)-1.D0)
2811  KP=KP+DELTKP
2812 .GE. IF (KPOVER) M=2000
2813  M=M+1
2814  IND=-IND
2815 .GT..AND..LT. IF ((ABS(DELTKP/KP)PRECI)(M1000)) GOTO 45
2816 .EQ. IF (IFAC1) THEN
2817  DF=SQRT(0.5D0/(PI*X))*EXP(X)
2818  K=K*DF
2819  KP=-KP*DF
2820  ELSE
2821 .LT. IF (XA) THEN
2822  PIA=PI*A
2823  DF=SQRT(0.5D0/(PI*X))*EXP(X-0.5D0*PIA)
2824  K=K*DF
2825  KP=-KP*DF
2826  ELSE
2827  SINTH=A/X
2828  THET=ASIN(SINTH)
2829  COSTH=COS(THET)
2830  FDOMIN=X*(COSTH+THET*SINTH)
2831  DF=SQRT(0.5D0/(PI*X))*EXP(X-FDOMIN)
2832  K=K*DF
2833  KP=-KP*DF
2834  ENDIF
2835  ENDIF
2836  PEXP=K
2837  PEXPP=KP
2838  ELSE
2839  PEXP=0.D0
2840  PEXPP=0.D0
2841  ENDIF
2842  RETURN
2843  END
2844 
2845  SUBROUTINE AIEXLI(IFAC,X,A,DLAI,DLAID,IERROL)
2846 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2847 CCC AIRY-TYPE ASYMPTOTIC EXPANSION FOR THE K AND K' c
2848 ccc functions c
2849 ccc c
2850 ccc input variables: c
2851 ccc x: argument of the functions c
2852 ccc a: order of the functions c
2853 ccc ifac: c
2854 ccc =1, non scaled functions c
2855 ccc otherwise, scaled functions c
2856 ccc output variables: c
2857 ccc dlai: l FUNCTION c
2858 ccc dlaid: derivative of the l FUNCTION c
2859 ccc ierrol: error flag c
2860 ccc * IF ierrol=0, computation successful. c
2861 ccc * IF ierrol=1, computation out of range. c
2862 ccc * IF ierrol=2, argument and/or order, c
2863 ccc out of range. c
2864 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2865 ccc local variables:
2866 ccccccccccccccccccccccccccccccccccccccccccccccccccccccc
2867 c a0ex : exact value of the coefficient a_0
2868 c a13 : a**(1/3)
2869 c a2 : a**2
2870 c a23 : a**(2/3)
2871 c a2k : a**(2*k)
2872 c ac : coefficients a_s(psi) used in the expansions
2873 c apihal : a*pi/2
2874 c arg : -a**(2/3)*psi
2875 c as : accumulates the contribution of the ac coefficients
2876 c(for the lia(x) function)
2877 c asp : accumulates the contribution of the ac coefficients
2878 c(for the lia'(X) FUNCTION)
2879 C B0EX : EXACT VALUE OF THE COEFFICIENT B_0
2880 C B0PEX : EXACT VALUE OF THE COEFFICIENT B'_0
2881 c bc : coefficients b_s(psi) used in the expansions
2882 c bii : imaginary part of the airy FUNCTION bi(Z)
2883 c bir : REAL part of the airy function bi(z)
2884 c bs : accumulates the contribution of the bc coefficients
2885 c(for the lia(x) function)
2886 c bso : accumulates the old values of bs
2887 c bsp : accumulates the contribution of the bc coefficients
2888 c(for the lia'(X) FUNCTION)
2889 C BSPO : ACCUMULATES THE OLD VALUES OF BSP
2890 C C0EX : EXACT VALUE OF THE COEFFICIENT C_0
2891 C CHI : ACCUMULATES THE CONTRIBUTION OF THE CHIN COEFFICIENTS
2892 C CHIEX : EXACT VALUE OF THE FUNCTION CHI(PSI)
2893 C CHIN : COEFFICIENTS FOR THE EXPANSION OF THE FUNCTION
2894 C CHI(PSI)
2895 C COSTH : COS(THET)
2896 C D0EX : EXACT VALUE OF THE COEFFICIENT D_0
2897 C DBII : IMAGINARY PART OF THE AIRY FUNCTION BI'(z)
2898 c dbir : REAL part of the airy function bi'(Z)
2899 C DF : 2**(1/3)
2900 C DZZ : (ABS(1-Z**2))**(1/2)
2901 C ETA : PSI/2**(1/3)
2902 C ETAJ : ETA**J
2903 C ETAK : ETA**K
2904 C ETAL : ETA**L
2905 C EXPAM : EXP(A*PI/2-FDOMIN)
2906 C EXPAPI : EXP(A*PI/2)
2907 C F2 : (-)**K/A**(2*K)
2908 C F4 : (A**(1/3))**4
2909 C FAC : FACTOR FOR THE LIA(X) FUNCTION
2910 C FACD : FACTOR FOR THE LIA'(x) function
2911 c fdomin : x*(cos(thet)+thet*sin(thet))
2912 c ierro : error flag for the airy functions
2913 c ifaca :
2914 c *IF ifaca=1, computation of unscaled airy bi(z),bi'(Z)
2915 C *IF IFACA=2, COMPUTATION OF SCALED AIRY BI(Z),BI'(z)
2916 c ifun :
2917 c *IF ifun=1, computation of airy bi(z)
2918 c *IF ifun=2, computation of airy bi'(Z)
2919 C OVER : OVERFLOW NUMBER
2920 C PHI : COEFFICIENTS FOR THE EXPANSION OF THE FUNCTION
2921 C PHI(PSI)
2922 C PHIEX : EXACT VALUE OF THE FUNCTION PHI(PSI)
2923 C PHIEX2 : PHI(PSI)**2
2924 C PHIS : VALUE OF THE FUNCTION PHI(PSI)
2925 C PI : 3.1415...
2926 C PIHALF : PI/2
2927 C PSI :
2928 C *IF 0<=Z<=1, 2/3*PSI**3/2=LOG((1+SQRT(1-Z**2))/Z)-SQRT(1-Z**2)
2929 C *IF Z>=1, 2/3*(-PSI)**3/2=SQRT(Z**2-1)-ARCCOS(1/Z)
2930 C PSI12 : SQRT(PSI)
2931 C PSI2 : PSI*PSI
2932 C PSI3 : PSI**3
2933 C SAS : ACCUMULATES THE CONTRIBUTION OF A_S(PSI)
2934 C SBS : ACCUMULATES THE CONTRIBUTION OF B_S(PSI)
2935 C SCS : ACCUMULATES THE CONTRIBUTION OF C_S(PSI)
2936 C SDS : ACCUMULATES THE CONTRIBUTION OF D_S(PSI)
2937 C SIG : (-)**K
2938 C SINTH : SIN(THET)
2939 C THET : ASIN(A/X)
2940 C Y : Z-1
2941 C Z : X/A
2942 C Z2 : Z**2
2943 C Z21M : 1-Z**2
2944 C ZDS : SQRT(1-Z**2)
2945 C ZMASF : EXPANSION IN TERMS OF (Z-1) FOR THE
2946 C CALCULATION OF PSI(Z)
2947 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
2948  DOUBLE PRECISION A,A0EX,A13,A2,A23,A2K,AC(0:5,0:20),
2949  + APIHAL,ARG,AS,ASP,B0EX,B0PEX,BC(0:5,0:20),BII,BIR,
2950  + BS,BSO,BSP,BSPO,C0EX,CHI,CHIEX,CHIN(0:26),COSTH,
2951  + COZMAS(0:12),D0EX,D1MACH,DBII,DBIR,DF,DLAI,DLAID
2952  DOUBLE PRECISION DZZ,ETA,ETAJ,ETAK,ETAL,EXPAM,EXPAPI,
2953  + F2,F4,FAC,FACD,FDOMIN,OVER,PHI(0:35),PHIEX,PHIEX2,
2954  + PHIS,PI,PIHALF,PSI,PSI12,PSI2,PSI3,SAS,SBS,SCS,SDS,
2955  + SIG,SINTH,THET,X,Y,Z,Z2,Z21M,ZDS,ZMASF
2956  INTEGER IERRO,IERROL,IFAC,IFACA,IFUN,INDA(0:5),INDB(0:5),
2957  + J,K,L
2958  SAVE PHI,CHIN,AC,BC
2959 CCCC VALUES OF THE COEFFICIENTS
2960  DATA PHI/1.D0,0.2D0,.25714285714285714286D-1,
2961  + -.56507936507936507937D-2,-.39367965367965367965D-2,
2962  + -.5209362066504924D-3,.3708541264731741D-3,
2963  + .2123827840293627D-3,.2150629788753145D-4,
2964  + -.2636904062981799D-4,-.1405469826493129D-4,
2965  + -.1149328899029441D-5,.1972641193938624D-5,
2966  + .1014324305593329D-5,.7034331100192200D-7,
2967  + -.1525044777392676D-6,-.7677866256900572D-7,
2968  + -.4669842638693018D-8,.1206673645965329D-7,
2969  + .5990877668092344D-8,.3269102150077715D-9,
2970  + -.97138350244428335095D-9,-.47745934295232233834D-9,
2971  + -.23750318642839155779D-10,.79244598109106655567D-10,
2972  + .38653584230817865528D-10,.17734105846426807873D-11,
2973  + -.65332110030864751956D-11,-.31673512094772575686D-11,
2974  + -.13524195177121030660D-12,.54324103217903338951D-12,
2975  + .26204918647967626464D-12,.10488134973664584898D-13,
2976  + -.45490420121539001435D-13,-.21851238232690420878D-13,
2977  + -.82456080260379042800D-15/
2978  DATA CHIN/.2D0,.1142857142857142857142857D-1,
2979  + -.2438095238095238095238095D-1,-.1003471449185734900020614D-1,
2980  + .8811404468547325690182833D-3,.2318339655809043564145605D-2,
2981  + .7794494413564441575646057D-3,-.1373952504077228744949558D-3,
2982  + -.2162301322540308393022684D-3,-.6585004634375583559042795D-4,
2983  + .1502851367097217578058824D-4,.1991904617871647675455262D-4,
2984  + .5731972706719818525615427D-5,-.1496427612747891044606885D-5,
2985  + -.1821231830428939670133992D-5,-.5054254875930882534538217D-6,
2986  + .1433283859497625931203415D-6,.1657681319078678321113361D-6,
2987  + .4485592642702540575627044D-7,-.1346138242826094098161839D-7,
2988  + -.1504601862773535117708677D-7,-.3995660198654805921651406D-8,
2989  + .1250124931952495738882300D-8,.1363187174221864073749532D-8,
2990  + .3567608459777506132029204D-9,-.1152749290139859167119863D-9,
2991  + -.1233547289799408517691696D-9/
2992  DATA COZMAS/.9428090415820634D0,-.4242640687119285D0,
2993  + .2904188565587606D0,-.2234261009999160D0,
2994  + .1821754153944745D0,-.1539625154198624D0,
2995  + .1333737583085222D0,-.1176617834148007D0,
2996  + .1052687194772381D0,-.9524025714638822D-1,
2997  + .8695738197073783D-1,-.8000034897653656D-1,
2998  + .7407421797273086D-1/
2999  DATA AC(0,0)/1.D0/
3000  DATA AC(1,0)/-.44444444444444444445D-2/
3001  DATA AC(1,1)/-.18441558441558441558D-2/
3002  DATA AC(1,2)/.11213675213675213675D-2/
3003  DATA AC(1,3)/.13457752124418791086D-2/
3004  DATA AC(1,4)/.0003880626562979504D0/
3005  DATA AC(1,5)/-.0001830686723781799D0/
3006  DATA AC(1,6)/-.0001995460887806733D0/
3007  DATA AC(1,7)/-.00005256191234041587D0/
3008  DATA AC(1,8)/.00002460619652459158D0/
3009  DATA AC(1,9)/.00002519246780924541D0/
3010  DATA AC(1,10)/.6333157376533242D-5/
3011  DATA AC(1,11)/-.2957485733830202D-5/
3012  DATA AC(1,12)/-.2925255920564838D-5/
3013  DATA AC(1,13)/-.7159702610502009D-6/
3014  DATA AC(1,14)/.3331510720390949D-6/
3015  DATA AC(1,15)/.3227670475692310D-6/
3016  DATA AC(1,16)/.7767729381664199D-7/
3017  DATA AC(1,17)/-.3600954237921120D-7/
3018  DATA AC(1,18)/-.3441724449034226D-7/
3019  DATA AC(1,19)/-.8188194356398772D-8/
3020  DATA AC(1,20)/.3783148485152038D-8/
3021  DATA AC(2,0)/.69373554135458897365D-3/
3022  DATA AC(2,1)/.46448349036584330703D-3/
3023  DATA AC(2,2)/-.42838130171535112460D-3/
3024  DATA AC(2,3)/-.0007026702868771135D0/
3025  DATA AC(2,4)/-.0002632580046778811D0/
3026  DATA AC(2,5)/.0001663853666288703D0/
3027  DATA AC(2,6)/.0002212087687818584D0/
3028  DATA AC(2,7)/.00007020345615329662D0/
3029  DATA AC(2,8)/-.00004000421782540614D0/
3030  DATA AC(2,9)/-.00004786324966453962D0/
3031  DATA AC(2,10)/-.00001394600741473631D0/
3032  DATA AC(2,11)/.7536186591273727D-5/
3033  DATA AC(2,12)/.8478502161067667D-5/
3034  DATA AC(2,13)/.2345355228453912D-5/
3035  DATA AC(2,14)/-.1225943294710883D-5/
3036  DATA AC(2,15)/-.1325082343401027D-5/
3037  DATA AC(2,16)/-.3539954776569997D-6/
3038  DATA AC(2,17)/.1808291719376674D-6/
3039  DATA AC(2,18)/.1900383515233655D-6/
3040  DATA AC(3,0)/-.35421197145774384076D-3/
3041  DATA AC(3,1)/-.31232252789031883276D-3/
3042  DATA AC(3,2)/.3716442237502298D-3/
3043  DATA AC(3,3)/.0007539269155977733D0/
3044  DATA AC(3,4)/.0003408300059444739D0/
3045  DATA AC(3,5)/-.0002634968172069594D0/
3046  DATA AC(3,6)/-.0004089275726648432D0/
3047  DATA AC(3,7)/-.0001501108759563460D0/
3048  DATA AC(3,8)/.00009964015205538056D0/
3049  DATA AC(3,9)/.0001352492955751283D0/
3050  DATA AC(3,10)/.00004443117087272903D0/
3051  DATA AC(3,11)/-.00002713205071914117D0/
3052  DATA AC(3,12)/-.00003396796969771860D0/
3053  DATA AC(3,13)/-.00001040708865273043D0/
3054  DATA AC(3,14)/.6024639065414414D-5/
3055  DATA AC(3,15)/.7143919607846883D-5/
3056  DATA AC(4,0)/.378194199201773D-3/
3057  DATA AC(4,1)/.000404943905523634D0/
3058  DATA AC(4,2)/-.000579130526946453D0/
3059  DATA AC(4,3)/-.00138017901171011D0/
3060  DATA AC(4,4)/-.000722520056780091D0/
3061  DATA AC(4,5)/.000651265924036825D0/
3062  DATA AC(4,6)/.00114674563328389D0/
3063  DATA AC(4,7)/.000474423189340405D0/
3064  DATA AC(4,8)/-.000356495172735468D0/
3065  DATA AC(4,9)/-.000538157791035111D0/
3066  DATA AC(4,10)/-.000195687390661225D0/
3067  DATA AC(4,11)/.000132563525210293D0/
3068  DATA AC(4,12)/.000181949256267291D0/
3069  DATA AC(5,0)/-.69114139728829416760D-3/
3070  DATA AC(5,1)/-.00085995326611774383285D0/
3071  DATA AC(5,2)/.0014202335568143511489D0/
3072  DATA AC(5,3)/.0038535426995603052443D0/
3073  DATA AC(5,4)/.0022752811642901374595D0/
3074  DATA AC(5,5)/-.0023219572034556988366D0/
3075  DATA AC(5,6)/-.0045478643394434635622D0/
3076  DATA AC(5,7)/-.0020824431758272449919D0/
3077  DATA AC(5,8)/.0017370443573195808719D0/
3078  DATA BC(0,0)/.14285714285714285714D-1/
3079  DATA BC(0,1)/.88888888888888888889D-2/
3080  DATA BC(0,2)/.20482374768089053803D-2/
3081  DATA BC(0,3)/-.57826617826617826618D-3/
3082  DATA BC(0,4)/-.60412089799844901886D-3/
3083  DATA BC(0,5)/-.0001472685745626922D0/
3084  DATA BC(0,6)/.00005324102148009784D0/
3085  DATA BC(0,7)/.00005206561006583416D0/
3086  DATA BC(0,8)/.00001233115050894939D0/
3087  DATA BC(0,9)/-.4905932728531366D-5/
3088  DATA BC(0,10)/-.4632230987136350D-5/
3089  DATA BC(0,11)/-.1077174523455235D-5/
3090  DATA BC(0,12)/.4475963978932822D-6/
3091  DATA BC(0,13)/.4152586188464624D-6/
3092  DATA BC(0,14)/.9555819293589234D-7/
3093  DATA BC(0,15)/-.4060599208403059D-7/
3094  DATA BC(0,16)/-.3731367187988482D-7/
3095  DATA BC(0,17)/-.8532670645553778D-8/
3096  DATA BC(0,18)/.3673017245573624D-8/
3097  DATA BC(0,19)/.3355960460784536D-8/
3098  DATA BC(0,20)/.7643107095110475D-9/
3099  DATA BC(1,0)/-.11848595848595848596D-2/
3100  DATA BC(1,1)/-.13940630797773654917D-2/
3101  DATA BC(1,2)/-.48141005586383737645D-3/
3102  DATA BC(1,3)/.26841705366016142958D-3/
3103  DATA BC(1,4)/.0003419706982709903D0/
3104  DATA BC(1,5)/.0001034548234902078D0/
3105  DATA BC(1,6)/-.00005418191982095504D0/
3106  DATA BC(1,7)/-.00006202184830690167D0/
3107  DATA BC(1,8)/-.00001724885886056087D0/
3108  DATA BC(1,9)/.8744675992887053D-5/
3109  DATA BC(1,10)/.9420684216180929D-5/
3110  DATA BC(1,11)/.2494922112085850D-5/
3111  DATA BC(1,12)/-.1238458608836357D-5/
3112  DATA BC(1,13)/-.1285461713809769D-5/
3113  DATA BC(1,14)/-.3299710862537507D-6/
3114  DATA BC(1,15)/.1613441105788315D-6/
3115  DATA BC(1,16)/.1633623194402374D-6/
3116  DATA BC(1,17)/.4104252949605779D-7/
3117  DATA BC(1,18)/-.1984317042326989D-7/
3118  DATA BC(1,19)/-.1973948142769706D-7/
3119  DATA BC(1,20)/-.4882194808588752D-8/
3120  DATA BC(2,0)/.43829180944898810994D-3/
3121  DATA BC(2,1)/.71104865116708668943D-3/
3122  DATA BC(2,2)/.31858383945387580576D-3/
3123  DATA BC(2,3)/-.0002404809426804458D0/
3124  DATA BC(2,4)/-.0003722966038621536D0/
3125  DATA BC(2,5)/-.0001352752059595618D0/
3126  DATA BC(2,6)/.00008691694372704142D0/
3127  DATA BC(2,7)/.0001158750753591377D0/
3128  DATA BC(2,8)/.00003724965927846447D0/
3129  DATA BC(2,9)/-.00002198334949606935D0/
3130  DATA BC(2,10)/-.00002686449633870452D0/
3131  DATA BC(2,11)/-.8023061612032524D-5/
3132  DATA BC(2,12)/.4494756592180126D-5/
3133  DATA BC(2,13)/.5193504763856015D-5/
3134  DATA BC(2,14)/.1477156191529617D-5/
3135  DATA BC(2,15)/-.7988793826096784D-6/
3136  DATA BC(3,0)/-.37670439477105454219D-3/
3137  DATA BC(3,1)/-.75856271658798642365D-3/
3138  DATA BC(3,2)/-.0004103253968775048D0/
3139  DATA BC(3,3)/.0003791263310429010D0/
3140  DATA BC(3,4)/.0006850981673903450D0/
3141  DATA BC(3,5)/.0002878310571932216D0/
3142  DATA BC(3,6)/-.0002157010636115705D0/
3143  DATA BC(3,7)/-.0003260863991373500D0/
3144  DATA BC(3,8)/-.0001181317008748678D0/
3145  DATA BC(3,9)/.00007887526841582582D0/
3146  DATA BC(3,10)/.0001072081833420685D0/
3147  DATA BC(3,11)/.00003544595251288735D0/
3148  DATA BC(3,12)/-.00002201447920733824D0/
3149  DATA BC(3,13)/-.00002789336359620813D0/
3150  DATA BC(4,0)/.00058453330122076187255D0/
3151  DATA BC(4,1)/.0013854690422372401251D0/
3152  DATA BC(4,2)/.00086830374184946900245D0/
3153  DATA BC(4,3)/-.00093502904801345951693D0/
3154  DATA BC(4,4)/-.0019175486005525492059D0/
3155  DATA BC(4,5)/-.00090795047113308137941D0/
3156  DATA BC(4,6)/.00077050429806392235104D0/
3157  DATA BC(4,7)/.0012953100128255983488D0/
3158  DATA BC(4,8)/.00051933869471899550762D0/
3159  DATA BC(4,9)/-.00038482631948759834653D0/
3160  DATA BC(4,10)/-.00057335393099012476502D0/
3161  DATA BC(5,0)/-.0014301070053470410656D0/
3162  DATA BC(5,1)/-.0038637811942002539408D0/
3163  DATA BC(5,2)/-.0027322816261168328245D0/
3164  DATA BC(5,3)/.0033294980346743452748D0/
3165  DATA BC(5,4)/.0075972237660887795911D0/
3166  DATA BC(5,5)/.0039816655673062060620D0/
3167  DATA BC(5,6)/-.0037510180460986006595D0/
3168  DATA INDA/0,20,18,15,12,8/
3169  DATA INDB/20,20,15,13,10,6/
3170 CCCC CONSTANTS CCCCCCCCCCCCCC
3171  PI=ACOS(-1.D0)
3172 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
3173  OVER=D1MACH(2)*1.D-8
3174  IERROL=0
3175 .EQ. IF (IFAC1) THEN
3176 .GE. IF (XA) THEN
3177  SINTH=A/X
3178  THET=ASIN(SINTH)
3179  COSTH=COS(THET)
3180  FDOMIN=X*(COSTH+THET*SINTH)
3181 .GE. IF (FDOMINLOG(OVER)) IERROL=1
3182  ELSE
3183 .GE. IF ((PI*A*0.5D0)LOG(OVER)) IERROL=1
3184  ENDIF
3185  ENDIF
3186 .EQ. IF (IERROL0) THEN
3187  PIHALF=PI*0.5D0
3188  DF=2.D0**(1.D0/3.D0)
3189 CCCC VARIABLES CCCCCCCCCCCCCCC
3190  Z=X/A
3191  A2=A*A
3192  A13=A**(1.D0/3.D0)
3193  A23=A13*A13
3194 .EQ. IF (IFAC1) THEN
3195  APIHAL=A*PIHALF
3196  EXPAPI=EXP(APIHAL)
3197  FAC=EXPAPI*0.5D0/A13
3198  FACD=EXPAPI/Z/A23
3199  ELSE
3200 .LT. IF (XA) THEN
3201  FAC=0.5D0/A13
3202  FACD=1.D0/Z/A23
3203  ELSE
3204  SINTH=A/X
3205  THET=ASIN(SINTH)
3206  COSTH=COS(THET)
3207  FDOMIN=X*(COSTH+THET*SINTH)
3208  EXPAM=EXP(A*PIHALF-FDOMIN)
3209  FAC=EXPAM*0.5D0/A13
3210  FACD=EXPAM/Z/A23
3211  ENDIF
3212  ENDIF
3213  F4=A13**4
3214 .LE. IF (Z0.9D0) THEN
3215  ZDS=SQRT((1.D0-Z)*(1.D0+Z))
3216  PSI=(1.5D0*(LOG((1.D0+ZDS)/Z)-ZDS))**(2.D0/3.D0)
3217 .GT. ELSEIF (Z1.1D0) THEN
3218  ZDS=SQRT((Z-1.D0)*(1.D0+Z))
3219  PSI=-(1.5D0*(ZDS-ACOS(1.D0/Z)))**(2.D0/3.D0)
3220  ELSE
3221  Y=Z-1.D0
3222  ZMASF=COZMAS(12)
3223  DO 8 K=0,11
3224  J=11-K
3225  ZMASF=COZMAS(J)+Y*ZMASF
3226  8 CONTINUE
3227  ZMASF=ABS(Y)**(1.5D0)*ZMASF
3228 .LT. IF (Z1.D0) THEN
3229  PSI=(1.5D0*ZMASF)**(2.D0/3.D0)
3230  ELSE
3231  PSI=-(1.5D0*ZMASF)**(2.D0/3.D0)
3232  ENDIF
3233  ENDIF
3234  ETA=PSI/DF
3235  ARG=-A23*PSI
3236  PSI2=PSI*PSI
3237  PSI3=PSI2*PSI
3238 .GT..AND..LT. IF ((Z0.8D0)(Z1.2D0)) THEN
3239  PHIS=0.D0
3240  CHI=0.D0
3241  SAS=0.D0
3242  SBS=0.D0
3243  SDS=0.D0
3244  SCS=1.D0
3245  BS=0.D0
3246  BSPO=0.D0
3247  DO 41 L=0,20
3248 .EQ. IF (L0) THEN
3249  ETAL=1.D0
3250  ELSE
3251  ETAL=ETAL*ETA
3252  ENDIF
3253  CHI=CHI+CHIN(L)*ETAL
3254  41 CONTINUE
3255  CHI=CHI/DF
3256  DO 10 K=0,20
3257  BSO=BS
3258  BSPO=BSP
3259  AS=0.D0
3260  BS=0.D0
3261  ASP=0.D0
3262  BSP=0.D0
3263 .EQ. IF (K0) THEN
3264  ETAK=1.D0
3265  A2K=1.D0
3266  SIG=1.D0
3267  ELSE
3268  ETAK=ETAK*ETA
3269  A2K=A2K*A2
3270  SIG=-1.D0*SIG
3271  ENDIF
3272  PHIS=PHIS+PHI(K)*ETAK
3273  F2=SIG/A2K
3274 .LE. IF (K5) THEN
3275  DO 20 J=0,20
3276 .EQ. IF (J0) THEN
3277  ETAJ=1.D0
3278  ELSE
3279  ETAJ=ETAJ*ETA
3280  ENDIF
3281 .LE. IF (JINDA(K)) THEN
3282  AS=AS+AC(K,J)*ETAJ
3283  ENDIF
3284 .LE. IF (JINDB(K)) THEN
3285  BS=BS+BC(K,J)*ETAJ
3286  ENDIF
3287 .LE. IF ((J+1)INDA(K)) THEN
3288  ASP=ASP+(J+1)*AC(K,J+1)*ETAJ
3289  ENDIF
3290 .LE. IF ((J+1)INDB(K)) THEN
3291  BSP=BSP+(J+1)*BC(K,J+1)*ETAJ
3292  ENDIF
3293  20 CONTINUE
3294  ASP=1.D0/DF*ASP
3295  BS=BS*DF
3296  SAS=SAS+AS*F2
3297  SBS=SBS+BS*F2/DF
3298  SDS=SDS-(CHI*AS+ASP+PSI*BS)*F2
3299  ENDIF
3300 .GT..AND..LE. IF ((K0)(K6)) THEN
3301  SCS=SCS+(AS+CHI*BSO+BSPO)*F2
3302  ENDIF
3303  10 CONTINUE
3304  PHIS=DF*PHIS
3305  SBS=DF*SBS
3306  ELSE
3307 CCCC EXACT VALUES CCCCCCCCCCCCCCCCCCCCC
3308  Z2=Z*Z
3309  Z21M=1.D0-Z2
3310  PHIEX=(4.D0*PSI/Z21M)**0.25D0
3311  PHIEX2=PHIEX*PHIEX
3312  A0EX=1.D0
3313  B0EX=-5.D0/48.D0/(PSI*PSI)+PHIEX2/48.D0/PSI*
3314  + (5.D0/Z21M-3.D0)
3315  CHIEX=0.25D0/PSI*(1.D0-Z2*PHIEX2**3*0.25D0)
3316  CHI=CHIEX
3317  D0EX=-(7.D0/48.D0/PSI+PHIEX2/48.D0*
3318  + (9.D0-7.D0/Z21M))
3319 .GT. IF (PSI0.D0) THEN
3320  PSI12=SQRT(PSI)
3321  DZZ=SQRT(Z21M)
3322  ELSE
3323  PSI12=SQRT(-PSI)
3324  DZZ=SQRT(ABS(Z21M))
3325  ENDIF
3326  B0PEX=5.D0/24.D0/PSI3+PHIEX2/48.D0*((2.D0*CHI*PSI-1.D0)/
3327  + PSI2*(5.D0/Z21M-3.D0)-10.D0*Z2*PSI12/
3328  + Z21M**2/DZZ/PSI)
3329  C0EX=1.D0
3330  SAS=A0EX
3331  SBS=B0EX
3332  SCS=C0EX
3333  SDS=D0EX
3334  BS=0.D0
3335  BSP=0.D0
3336  A2K=1.D0
3337  SIG=1.D0
3338  DO 30 K=1,6
3339  BSO=BS
3340  BSPO=BSP
3341  AS=0.D0
3342  BS=0.D0
3343  ASP=0.D0
3344  BSP=0.D0
3345  A2K=A2K*A2
3346  SIG=-1.D0*SIG
3347  F2=SIG/A2K
3348 .LE. IF (K5) THEN
3349  DO 40 J=0,20
3350 .EQ. IF (J0) THEN
3351  ETAJ=1.D0
3352  ELSE
3353  ETAJ=ETAJ*ETA
3354  ENDIF
3355 .LE. IF (JINDA(K)) THEN
3356  AS=AS+AC(K,J)*ETAJ
3357  ENDIF
3358 .LE. IF (JINDB(K)) THEN
3359  BS=BS+BC(K,J)*ETAJ
3360  ENDIF
3361 .LE. IF ((J+1)INDA(K)) THEN
3362  ASP=ASP+(J+1)*AC(K,J+1)*ETAJ
3363  ENDIF
3364 .LE. IF ((J+1)INDB(K)) THEN
3365  BSP=BSP+(J+1)*BC(K,J+1)*ETAJ
3366  ENDIF
3367  40 CONTINUE
3368  BS=DF*BS
3369  ASP=ASP/DF
3370  SAS=SAS+AS*F2
3371  SBS=SBS+BS*F2
3372  SDS=SDS-(CHI*AS+ASP+PSI*BS)*F2
3373  ENDIF
3374 .GE. IF (K2) THEN
3375  SCS=SCS+(AS+CHI*BSO+BSPO)*F2
3376  ELSE
3377  SCS=SCS+(AS+CHI*B0EX+B0PEX)*F2
3378  ENDIF
3379  30 CONTINUE
3380  PHIS=PHIEX
3381  ENDIF
3382 CCCCC CALL THE AIRY ROUTINE CCCCCCCCC
3383  IFUN=1
3384  IFACA=1
3385  CALL BIZ(IFUN,IFACA,ARG,0.D0,BIR,BII,IERRO)
3386  IFUN=2
3387  IFACA=1
3388  CALL BIZ(IFUN,IFACA,ARG,0.D0,DBIR,DBII,IERRO)
3389  DLAI=FAC*PHIS*(BIR*SAS+DBIR*SBS/F4)
3390  DLAID=FACD/PHIS*(DBIR*SCS+BIR*SDS/A23)
3391  ELSE
3392  DLAI=0.D0
3393  DLAID=0.D0
3394  ENDIF
3395  RETURN
3396  END
3397 
3398 
3399  DOUBLE PRECISION FUNCTION V(U)
3400 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3401 CCC AUXILIARY FUNCTION:
3402 CCC THIS FUNCTION COMPUTES THE LAST EXPRESSION
3403 CCC IN EQ.(21) OF THE REFERENCE
3404 CCC 'computing solutions of the modified bessel
3405 ccc differential equation for imaginary orders
3406 ccc and positive arguments',
3407 CCC A. GIL, J. SEGURA, N.M. TEMME
3408 CCC ACM TRANS. MATH. SOFT. (2004)
3409 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3410 CCC INPUT VARIABLE:
3411 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3412 C U : ARGUMENT OF THE FUNCTION
3413 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3414 CCC LOCAL VARIABLES:
3415 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3416 C COSTH : COS(THET)
3417 C FDOMIN : X*(COS(THET)+THET*SIN(THET))
3418 C SINTH : SIN(THET)
3419 C THET : ASIN(PNU/X)
3420 C UNDER : UNDERFLOW NUMBER
3421 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3422  DOUBLE PRECISION COSTH,D1MACH,FDOMIN,SINTH,
3423  + THET,U,UNDER
3424  COMMON/PARMON/THET,SINTH,COSTH,FDOMIN
3425 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
3426  UNDER=D1MACH(1)*1.D+6
3427 .LT. IF (ABS(U)UNDER) THEN
3428  V=ASIN(SINTH)
3429  ELSE
3430  V=ASIN(U/SINH(U)*SINTH)
3431  ENDIF
3432  RETURN
3433  END
3434 
3435  DOUBLE PRECISION FUNCTION PHIR(U)
3436 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3437 CCC AUXILIARY FUNCTION:
3438 CCC THIS FUNCTION COMPUTES THE EXPONENT IN THE
3439 CCC INTEGRAND IN EQ.(20) OF THE REFERENCE
3440 CCC 'computing solutions of the modified bessel
3441 ccc differential equation for imaginary orders
3442 ccc and positive arguments',
3443 CCC A. GIL, J. SEGURA, N.M. TEMME
3444 CCC ACM TRANS. MATH. SOFT. (2004)
3445 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3446 CCC INPUT VARIABLE:
3447 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3448 C U : ARGUMENT OF THE FUNCTION
3449 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3450 CCC LOCAL VARIABLES:
3451 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3452 C COSTH : COS(THET)
3453 C COSVU : COS(V(U))
3454 C FDOMIN : X*(COS(THET)+THET*SIN(THET))
3455 C FU1 : 2*SINH(U/2)**2
3456 C FU2 : -2*SIN(FU3/2)*SIN(0.5*(V(U)+THET))
3457 C FU3 : ASIN(-SIN(THET)/(COS(THET)*U/SINH(U)
3458 C +COS(V(U))*(SINH(U)+U)*FUAC/(SINH(U)*SINH(U)))
3459 C FUAC : U**3/6+U**5/120+U**7/5040
3460 C OVER : OVERFLOW NUMBER
3461 C PNU : ORDER OF THE FUNCTION
3462 C SINHU : SINH(U)
3463 C SINHUH : SINH(U/2)
3464 C SINTH : SIN(THET)
3465 C THET : ASIN(PNU/X)
3466 C U2 : U**2
3467 C U3 : U**3
3468 C U5 : U**5
3469 C U7 : U**7
3470 C UH : U/2
3471 C UNDER : UNDERFLOW NUMBER
3472 C V(U) : ASIN(U/SINH(U)*SIN(THET))
3473 C VU : V(U)
3474 C X : ARGUMENT OF THE FUNCTIONS
3475 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3476  DOUBLE PRECISION COSTH,COSVU,D1MACH,FDOMIN,FU1,FU2,
3477  + FU3,FUAC,OVER,PNU,SINHU,SINHUH,SINTH,THET,U,U2,
3478  + U3,U5,U7,UH,UNDER,V,VU,X
3479  COMMON/ARGU/X,PNU
3480  COMMON/PARMON/THET,SINTH,COSTH,FDOMIN
3481 CCC MACHINE DEPENDENT CONSTANT (UNDERFLOW NUMBER)
3482  UNDER=D1MACH(1)*1.D+6
3483 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
3484  OVER=D1MACH(2)*1.D-6
3485 .LT. IF (U200.D0) THEN
3486 .LT. IF (ABS(U)0.1D0) THEN
3487 .LT. IF (ABS(U)UNDER) THEN
3488  PHIR=0.D0
3489  ELSE
3490  UH=U*0.5D0
3491  SINHUH=SINH(UH)
3492  FU1=2.D0*SINHUH*SINHUH
3493  U2=U*U
3494  U3=U2*U
3495  U5=U3*U2
3496  U7=U5*U2
3497  FUAC=U3/6.D0+U5/120.D0+U7/5040.D0
3498  SINHU=SINH(U)
3499  VU=V(U)
3500  COSVU=COS(VU)
3501  FU3=ASIN(-SINTH/(COSTH*U/SINHU+COSVU)*
3502  + (SINHU+U)*FUAC/(SINHU*SINHU))
3503  FU2=-2.D0*SIN(FU3*0.5D0)*SIN(0.5D0*(VU+THET))
3504  PHIR=X*(FU1*COSVU+FU2+SINTH*FU3)
3505  ENDIF
3506  ELSE
3507  VU=V(U)
3508  COSVU=COS(VU)
3509  PHIR=X*COSH(U)*COSVU+PNU*VU-FDOMIN
3510  ENDIF
3511  ELSE
3512  PHIR=OVER
3513  ENDIF
3514  RETURN
3515  END
3516 
3517  DOUBLE PRECISION FUNCTION SIGMA(U)
3518 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3519 CCC AUXILIARY FUNCTION:
3520 CCC THIS FUNCTION COMPUTES THE FUNCTION
3521 CCC SIGMA FROM EQ.(31) OF THE REFERENCE
3522 CCC 'computing solutions of the modified bessel
3523 ccc differential equation for imaginary orders
3524 ccc and positive arguments',
3525 CCC A. GIL, J. SEGURA, N.M. TEMME
3526 CCC ACM TRANS. MATH. SOFT. (2004)
3527 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3528 CCC INPUT VARIABLE:
3529 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3530 C U : ARGUMENT OF THE FUNCTION
3531 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3532 CCC LOCAL VARIABLES:
3533 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3534 C ARGU : (COSH(MU)*U-DMUFAC)/SINH(U)
3535 C COSCHI : COS(CHI)
3536 C COSHM : COSH(MU)
3537 C D1 : COSH(MU)*U-DMUFAC
3538 C DMU : SOLUTION MU OF COSH(MU)=PNU/X
3539 C DMUFAC : MU*COSH(MU)-SINH(MU)
3540 C SINCHI : SIN(CHI)
3541 C SINHM : SINH(MU)
3542 C SINHU : SINH(U)
3543 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3544  DOUBLE PRECISION ARGU,COSCHI,COSHM,D1,DMU,DMUFAC,
3545  + PI,SINCHI,SINHM,SINHU,U
3546  COMMON/PAROS1/DMU,COSHM,SINHM,DMUFAC,COSCHI,SINCHI
3547  PI=ACOS(-1.D0)
3548  D1=COSHM*U-DMUFAC
3549  SINHU=SINH(U)
3550  ARGU=D1/SINH(U)
3551 .GT. IF (ABS(ARGU)1.D0) THEN
3552  ARGU=1.D0
3553  ENDIF
3554  SIGMA=ASIN(ARGU)
3555 .LT. IF (UDMU) SIGMA=PI-SIGMA
3556  RETURN
3557  END
3558 
3559  DOUBLE PRECISION FUNCTION PHIB(U)
3560 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3561 CCC THIS FUNCTION COMPUTES THE FIRST EXPRESSION
3562 CCC OF EQ.(30) OF THE REFERENCE
3563 CCC 'computing solutions of the modified bessel
3564 ccc differential equation for imaginary orders
3565 ccc and positive arguments',
3566 CCC A. GIL, J. SEGURA, N.M. TEMME
3567 CCC ACM TRANS. MATH. SOFT. (2004)
3568 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3569 CCC INPUT VARIABLE:
3570 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3571 C U : ARGUMENT OF THE FUNCTION
3572 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3573 CCC LOCAL VARIABLES:
3574 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3575 C OVER : OVERFLOW NUMBER
3576 C PNU : ORDER OF THE FUNCTIONS
3577 C SIGMA(U) : ASIN((COSH(MU)*U-DMUFAC)/SINH(U))
3578 C SIGMAU : SIGMA(U)
3579 C X : ARGUMENT OF THE FUNCTIONS
3580 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
3581  DOUBLE PRECISION U,SIGMA,SIGMAU,OVER,D1MACH,
3582  + X,PNU
3583  COMMON/ARGU/X,PNU
3584 CCC MACHINE DEPENDENT CONSTANT (OVERFLOW NUMBER)
3585  OVER=D1MACH(2)*1.D-8
3586 .LT. IF (U200.D0) THEN
3587  SIGMAU=SIGMA(U)
3588  PHIB=X*COSH(U)*COS(SIGMAU)+PNU*SIGMAU
3589  ELSE
3590  PHIB=OVER
3591  ENDIF
3592  RETURN
3593  END
3594