Line data Source code
1 : !!****m* ABINIT/m_bessel
2 : !! NAME
3 : !! m_bessel
4 : !!
5 : !! FUNCTION
6 : !! Bessel functions
7 : !!
8 : !! COPYRIGHT
9 : !! Copyright (C) 2008-2026 ABINIT group (MG)
10 : !! This file is distributed under the terms of the
11 : !! GNU General Public License, see ~abinit/COPYING
12 : !! or http://www.gnu.org/copyleft/gpl.txt .
13 : !!
14 : !! SOURCE
15 :
16 : #if defined HAVE_CONFIG_H
17 : #include "config.h"
18 : #endif
19 :
20 : #include "abi_common.h"
21 :
22 : module m_bessel
23 :
24 : implicit none
25 :
26 : private
27 :
28 : public :: CALJY0
29 : public :: CALJY1
30 : public :: CALCK0
31 : public :: CALCK1
32 :
33 : contains
34 : !!***
35 :
36 0 : SUBROUTINE CALJY0(ARG,RESULT,JINT)
37 :
38 : !---------------------------------------------------------------------
39 : !
40 : ! This packet computes zero-order Bessel functions of the first and
41 : ! second kind (J0 and Y0), for real arguments X, where 0 < X <= XMAX
42 : ! for Y0, and |X| <= XMAX for J0. It contains two function-type
43 : ! subprograms, BESJ0 and BESY0, and one subroutine-type
44 : ! subprogram, CALJY0. The calling statements for the primary
45 : ! entries are:
46 : !
47 : ! Y = BESJ0(X)
48 : ! and
49 : ! Y = BESY0(X),
50 : !
51 : ! where the entry points correspond to the functions J0(X) and Y0(X),
52 : ! respectively. The routine CALJY0 is intended for internal packet
53 : ! use only, all computations within the packet being concentrated in
54 : ! this one routine. The function subprograms invoke CALJY0 with
55 : ! the statement
56 : ! CALL CALJY0(ARG,RESULT,JINT),
57 : ! where the parameter usage is as follows:
58 : !
59 : ! Function Parameters for CALJY0
60 : ! call ARG RESULT JINT
61 : !
62 : ! BESJ0(ARG) |ARG| .LE. XMAX J0(ARG) 0
63 : ! BESY0(ARG) 0 .LT. ARG .LE. XMAX Y0(ARG) 1
64 : !
65 : ! The main computation uses unpublished minimax rational
66 : ! approximations for X .LE. 8.0, and an approximation from the
67 : ! book Computer Approximations by Hart, et. al., Wiley and Sons,
68 : ! New York, 1968, for arguments larger than 8.0 Part of this
69 : ! transportable packet is patterned after the machine-dependent
70 : ! FUNPACK program BESJ0(X), but cannot match that version for
71 : ! efficiency or accuracy. This version uses rational functions
72 : ! that are theoretically accurate to at least 18 significant decimal
73 : ! digits for X <= 8, and at least 18 decimal places for X > 8. The
74 : ! accuracy achieved depends on the arithmetic system, the compiler,
75 : ! the intrinsic functions, and proper selection of the machine-
76 : ! dependent constants.
77 : !
78 : !*******************************************************************
79 : !
80 : ! Explanation of machine-dependent constants
81 : !
82 : ! XINF = largest positive machine number
83 : ! XMAX = largest acceptable argument. The functions AINT, SIN
84 : ! and COS must perform properly for ABS(X) .LE. XMAX.
85 : ! We recommend that XMAX be a small integer multiple of
86 : ! sqrt(1/eps), where eps is the smallest positive number
87 : ! such that 1+eps > 1.
88 : ! XSMALL = positive argument such that 1.0-(X/2)**2 = 1.0
89 : ! to machine precision for all ABS(X) .LE. XSMALL.
90 : ! We recommend that XSMALL < sqrt(eps)/beta, where beta
91 : ! is the floating-point radix (usually 2 or 16).
92 : !
93 : ! Approximate values for some important machines are
94 : !
95 : ! eps XMAX XSMALL XINF
96 : !
97 : ! CDC 7600 (S.P.) 7.11E-15 1.34E+08 2.98E-08 1.26E+322
98 : ! CRAY-1 (S.P.) 7.11E-15 1.34E+08 2.98E-08 5.45E+2465
99 : ! IBM PC (8087) (S.P.) 5.96E-08 8.19E+03 1.22E-04 3.40E+38
100 : ! IBM PC (8087) (D.P.) 1.11D-16 2.68D+08 3.72D-09 1.79D+308
101 : ! IBM 195 (D.P.) 2.22D-16 6.87D+09 9.09D-13 7.23D+75
102 : ! UNIVAC 1108 (D.P.) 1.73D-18 4.30D+09 2.33D-10 8.98D+307
103 : ! VAX 11/780 (D.P.) 1.39D-17 1.07D+09 9.31D-10 1.70D+38
104 : !
105 : !*******************************************************************
106 : !*******************************************************************
107 : !
108 : ! Error Returns
109 : !
110 : ! The program returns the value zero for X .GT. XMAX, and returns
111 : ! -XINF when BESLY0 is called with a negative or zero argument.
112 : !
113 : !
114 : ! Intrinsic functions required are:
115 : !
116 : ! ABS, AINT, COS, LOG, SIN, SQRT
117 : !
118 : !
119 : ! Latest modification: June 2, 1989
120 : !
121 : ! Author: W. J. Cody
122 : ! Mathematics and Computer Science Division
123 : ! Argonne National Laboratory
124 : ! Argonne, IL 60439
125 : !
126 : ! Taken from http://www.netlib.org/specfun/j0y0
127 : !
128 : !--------------------------------------------------------------------
129 : !IMPLICIT NONE
130 : INTEGER :: I,JINT
131 : !CS REAL
132 : DOUBLE PRECISION :: &
133 : & ARG,AX,CONS,DOWN,EIGHT,FIVE5,FOUR,ONE,ONEOV8,PI2,PJ0,&
134 : & PJ1,PLG,PROD,PY0,PY1,PY2,P0,P1,P17,QJ0,QJ1,QLG,QY0,QY1,&
135 : & QY2,Q0,Q1,RESJ,RESULT,R0,R1,SIXTY4,THREE,TWOPI,TWOPI1,&
136 : & TWOPI2,TWO56,UP,W,WSQ,XDEN,XINF,XMAX,XNUM,XSMALL,XJ0,&
137 : & XJ1,XJ01,XJ02,XJ11,XJ12,XY,XY0,XY01,XY02,XY1,XY11,XY12,&
138 : & XY2,XY21,XY22,Z,ZERO,ZSQ
139 : DIMENSION :: PJ0(7),PJ1(8),PLG(4),PY0(6),PY1(7),PY2(8),P0(6),P1(6),&
140 : & QJ0(5),QJ1(7),QLG(4),QY0(5),QY1(6),QY2(7),Q0(5),Q1(5)
141 : !-------------------------------------------------------------------
142 : ! Mathematical constants
143 : ! CONS = ln(.5) + Euler's gamma
144 : !-------------------------------------------------------------------
145 : !CS DATA ZERO,ONE,THREE,FOUR,EIGHT/0.0E0,1.0E0,3.0E0,4.0E0,8.0E0/,
146 : !CS 1 FIVE5,SIXTY4,ONEOV8,P17/5.5E0,64.0E0,0.125E0,1.716E-1/,
147 : !CS 2 TWO56,CONS/256.0E0,-1.1593151565841244881E-1/,
148 : !CS 3 PI2,TWOPI/6.3661977236758134308E-1,6.2831853071795864769E0/,
149 : !CS 4 TWOPI1,TWOPI2/6.28125E0,1.9353071795864769253E-3/
150 : DATA ZERO,ONE,THREE,FOUR,EIGHT/0.0D0,1.0D0,3.0D0,4.0D0,8.0D0/,&
151 : & FIVE5,SIXTY4,ONEOV8,P17/5.5D0,64.0D0,0.125D0,1.716D-1/,&
152 : & TWO56,CONS/256.0D0,-1.1593151565841244881D-1/,&
153 : & PI2,TWOPI/6.3661977236758134308D-1,6.2831853071795864769D0/,&
154 : & TWOPI1,TWOPI2/6.28125D0,1.9353071795864769253D-3/
155 : !-------------------------------------------------------------------
156 : ! Machine-dependent constants
157 : !-------------------------------------------------------------------
158 : !CS DATA XMAX/8.19E+03/,XSMALL/1.22E-09/,XINF/1.7E+38/
159 : DATA XMAX/1.07D+09/,XSMALL/9.31D-10/,XINF/1.7D+38/
160 : !-------------------------------------------------------------------
161 : ! Zeroes of Bessel functions
162 : !-------------------------------------------------------------------
163 : !CS DATA XJ0/2.4048255576957727686E+0/,XJ1/5.5200781102863106496E+0/,
164 : !CS 1 XY0/8.9357696627916752158E-1/,XY1/3.9576784193148578684E+0/,
165 : !CS 2 XY2/7.0860510603017726976E+0/,
166 : !CS 3 XJ01/ 616.0E+0/, XJ02/-1.4244423042272313784E-03/,
167 : !CS 4 XJ11/1413.0E+0/, XJ12/ 5.4686028631064959660E-04/,
168 : !CS 5 XY01/ 228.0E+0/, XY02/ 2.9519662791675215849E-03/,
169 : !CS 6 XY11/1013.0E+0/, XY12/ 6.4716931485786837568E-04/,
170 : !CS 7 XY21/1814.0E+0/, XY22/ 1.1356030177269762362E-04/
171 : DATA XJ0/2.4048255576957727686D+0/,XJ1/5.5200781102863106496D+0/,&
172 : & XY0/8.9357696627916752158D-1/,XY1/3.9576784193148578684D+0/,&
173 : & XY2/7.0860510603017726976D+0/,&
174 : & XJ01/ 616.0D+0/, XJ02/-1.4244423042272313784D-03/,&
175 : & XJ11/1413.0D+0/, XJ12/ 5.4686028631064959660D-04/,&
176 : & XY01/ 228.0D+0/, XY02/ 2.9519662791675215849D-03/,&
177 : & XY11/1013.0D+0/, XY12/ 6.4716931485786837568D-04/,&
178 : & XY21/1814.0D+0/, XY22/ 1.1356030177269762362D-04/
179 : !C-------------------------------------------------------------------
180 : !C Coefficients for rational approximation to ln(x/a)
181 : !C--------------------------------------------------------------------
182 : !CS DATA PLG/-2.4562334077563243311E+01,2.3642701335621505212E+02,
183 : !CS 1 -5.4989956895857911039E+02,3.5687548468071500413E+02/
184 : !CS DATA QLG/-3.5553900764052419184E+01,1.9400230218539473193E+02,
185 : !CS 1 -3.3442903192607538956E+02,1.7843774234035750207E+02/
186 : DATA PLG/-2.4562334077563243311D+01,2.3642701335621505212D+02,&
187 : & -5.4989956895857911039D+02,3.5687548468071500413D+02/
188 : DATA QLG/-3.5553900764052419184D+01,1.9400230218539473193D+02,&
189 : & -3.3442903192607538956D+02,1.7843774234035750207D+02/
190 : !C-------------------------------------------------------------------
191 : !C Coefficients for rational approximation of
192 : !C J0(X) / (X**2 - XJ0**2), XSMALL < |X| <= 4.0
193 : !C--------------------------------------------------------------------
194 : !CS DATA PJ0/6.6302997904833794242E+06,-6.2140700423540120665E+08,
195 : !CS 1 2.7282507878605942706E+10,-4.1298668500990866786E+11,
196 : !CS 2 -1.2117036164593528341E-01, 1.0344222815443188943E+02,
197 : !CS 3 -3.6629814655107086448E+04/
198 : !CS DATA QJ0/4.5612696224219938200E+05, 1.3985097372263433271E+08,
199 : !CS 1 2.6328198300859648632E+10, 2.3883787996332290397E+12,
200 : !CS 2 9.3614022392337710626E+02/
201 : DATA PJ0/6.6302997904833794242D+06,-6.2140700423540120665D+08,&
202 : & 2.7282507878605942706D+10,-4.1298668500990866786D+11,&
203 : & -1.2117036164593528341D-01, 1.0344222815443188943D+02,&
204 : & -3.6629814655107086448D+04/
205 : DATA QJ0/4.5612696224219938200D+05, 1.3985097372263433271D+08,&
206 : & 2.6328198300859648632D+10, 2.3883787996332290397D+12,&
207 : & 9.3614022392337710626D+02/
208 : !C-------------------------------------------------------------------
209 : !C Coefficients for rational approximation of
210 : !C J0(X) / (X**2 - XJ1**2), 4.0 < |X| <= 8.0
211 : !C-------------------------------------------------------------------
212 : !CS DATA PJ1/4.4176707025325087628E+03, 1.1725046279757103576E+04,
213 : !CS 1 1.0341910641583726701E+04,-7.2879702464464618998E+03,
214 : !CS 2 -1.2254078161378989535E+04,-1.8319397969392084011E+03,
215 : !CS 3 4.8591703355916499363E+01, 7.4321196680624245801E+02/
216 : !CS DATA QJ1/3.3307310774649071172E+02,-2.9458766545509337327E+03,
217 : !CS 1 1.8680990008359188352E+04,-8.4055062591169562211E+04,
218 : !CS 2 2.4599102262586308984E+05,-3.5783478026152301072E+05,
219 : !CS 3 -2.5258076240801555057E+01/
220 : DATA PJ1/4.4176707025325087628D+03, 1.1725046279757103576D+04,&
221 : & 1.0341910641583726701D+04,-7.2879702464464618998D+03,&
222 : & -1.2254078161378989535D+04,-1.8319397969392084011D+03,&
223 : & 4.8591703355916499363D+01, 7.4321196680624245801D+02/
224 : DATA QJ1/3.3307310774649071172D+02,-2.9458766545509337327D+03,&
225 : & 1.8680990008359188352D+04,-8.4055062591169562211D+04,&
226 : & 2.4599102262586308984D+05,-3.5783478026152301072D+05,&
227 : & -2.5258076240801555057D+01/
228 : !C-------------------------------------------------------------------
229 : !C Coefficients for rational approximation of
230 : !C (Y0(X) - 2 LN(X/XY0) J0(X)) / (X**2 - XY0**2),
231 : !C XSMALL < |X| <= 3.0
232 : !C--------------------------------------------------------------------
233 : !CS DATA PY0/1.0102532948020907590E+04,-2.1287548474401797963E+06,
234 : !CS 1 2.0422274357376619816E+08,-8.3716255451260504098E+09,
235 : !CS 2 1.0723538782003176831E+11,-1.8402381979244993524E+01/
236 : !CS DATA QY0/6.6475986689240190091E+02, 2.3889393209447253406E+05,
237 : !CS 1 5.5662956624278251596E+07, 8.1617187777290363573E+09,
238 : !CS 2 5.8873865738997033405E+11/
239 : DATA PY0/1.0102532948020907590D+04,-2.1287548474401797963D+06,&
240 : & 2.0422274357376619816D+08,-8.3716255451260504098D+09,&
241 : & 1.0723538782003176831D+11,-1.8402381979244993524D+01/
242 : DATA QY0/6.6475986689240190091D+02, 2.3889393209447253406D+05,&
243 : & 5.5662956624278251596D+07, 8.1617187777290363573D+09,&
244 : & 5.8873865738997033405D+11/
245 : !C-------------------------------------------------------------------
246 : !C Coefficients for rational approximation of
247 : !C (Y0(X) - 2 LN(X/XY1) J0(X)) / (X**2 - XY1**2),
248 : !C 3.0 < |X| <= 5.5
249 : !C--------------------------------------------------------------------
250 : !CS DATA PY1/-1.4566865832663635920E+04, 4.6905288611678631510E+06,
251 : !CS 1 -6.9590439394619619534E+08, 4.3600098638603061642E+10,
252 : !CS 2 -5.5107435206722644429E+11,-2.2213976967566192242E+13,
253 : !CS 3 1.7427031242901594547E+01/
254 : !CS DATA QY1/ 8.3030857612070288823E+02, 4.0669982352539552018E+05,
255 : !CS 1 1.3960202770986831075E+08, 3.4015103849971240096E+10,
256 : !CS 2 5.4266824419412347550E+12, 4.3386146580707264428E+14/
257 : DATA PY1/-1.4566865832663635920D+04, 4.6905288611678631510D+06,&
258 : & -6.9590439394619619534D+08, 4.3600098638603061642D+10,&
259 : & -5.5107435206722644429D+11,-2.2213976967566192242D+13,&
260 : & 1.7427031242901594547D+01/
261 : DATA QY1/ 8.3030857612070288823D+02, 4.0669982352539552018D+05,&
262 : & 1.3960202770986831075D+08, 3.4015103849971240096D+10,&
263 : & 5.4266824419412347550D+12, 4.3386146580707264428D+14/
264 : !C-------------------------------------------------------------------
265 : !C Coefficients for rational approximation of
266 : !C (Y0(X) - 2 LN(X/XY2) J0(X)) / (X**2 - XY2**2),
267 : !C 5.5 < |X| <= 8.0
268 : !C--------------------------------------------------------------------
269 : !CS DATA PY2/ 2.1363534169313901632E+04,-1.0085539923498211426E+07,
270 : !CS 1 2.1958827170518100757E+09,-1.9363051266772083678E+11,
271 : !CS 2 -1.2829912364088687306E+11, 6.7016641869173237784E+14,
272 : !CS 3 -8.0728726905150210443E+15,-1.7439661319197499338E+01/
273 : !CS DATA QY2/ 8.7903362168128450017E+02, 5.3924739209768057030E+05,
274 : !CS 1 2.4727219475672302327E+08, 8.6926121104209825246E+10,
275 : !CS 2 2.2598377924042897629E+13, 3.9272425569640309819E+15,
276 : !CS 3 3.4563724628846457519E+17/
277 : DATA PY2/ 2.1363534169313901632D+04,-1.0085539923498211426D+07,&
278 : & 2.1958827170518100757D+09,-1.9363051266772083678D+11,&
279 : & -1.2829912364088687306D+11, 6.7016641869173237784D+14,&
280 : & -8.0728726905150210443D+15,-1.7439661319197499338D+01/
281 : DATA QY2/ 8.7903362168128450017D+02, 5.3924739209768057030D+05,&
282 : & 2.4727219475672302327D+08, 8.6926121104209825246D+10,&
283 : & 2.2598377924042897629D+13, 3.9272425569640309819D+15,&
284 : & 3.4563724628846457519D+17/
285 : !C-------------------------------------------------------------------
286 : !C Coefficients for Hart,s approximation, |X| > 8.0
287 : !C-------------------------------------------------------------------
288 : !CS DATA P0/3.4806486443249270347E+03, 2.1170523380864944322E+04,
289 : !CS 1 4.1345386639580765797E+04, 2.2779090197304684302E+04,
290 : !CS 2 8.8961548424210455236E-01, 1.5376201909008354296E+02/
291 : !CS DATA Q0/3.5028735138235608207E+03, 2.1215350561880115730E+04,
292 : !CS 1 4.1370412495510416640E+04, 2.2779090197304684318E+04,
293 : !CS 2 1.5711159858080893649E+02/
294 : !CS DATA P1/-2.2300261666214198472E+01,-1.1183429920482737611E+02,
295 : !CS 1 -1.8591953644342993800E+02,-8.9226600200800094098E+01,
296 : !CS 2 -8.8033303048680751817E-03,-1.2441026745835638459E+00/
297 : !CS DATA Q1/1.4887231232283756582E+03, 7.2642780169211018836E+03,
298 : !CS 1 1.1951131543434613647E+04, 5.7105024128512061905E+03,
299 : !CS 2 9.0593769594993125859E+01/
300 : DATA P0/3.4806486443249270347D+03, 2.1170523380864944322D+04,&
301 : & 4.1345386639580765797D+04, 2.2779090197304684302D+04,&
302 : & 8.8961548424210455236D-01, 1.5376201909008354296D+02/
303 : DATA Q0/3.5028735138235608207D+03, 2.1215350561880115730D+04,&
304 : & 4.1370412495510416640D+04, 2.2779090197304684318D+04,&
305 : & 1.5711159858080893649D+02/
306 : DATA P1/-2.2300261666214198472D+01,-1.1183429920482737611D+02,&
307 : & -1.8591953644342993800D+02,-8.9226600200800094098D+01,&
308 : & -8.8033303048680751817D-03,-1.2441026745835638459D+00/
309 : DATA Q1/1.4887231232283756582D+03, 7.2642780169211018836D+03,&
310 : & 1.1951131543434613647D+04, 5.7105024128512061905D+03,&
311 : & 9.0593769594993125859D+01/
312 : !C-------------------------------------------------------------------
313 : !C Check for error conditions
314 : !C-------------------------------------------------------------------
315 0 : AX = ABS(ARG)
316 0 : IF ((JINT .EQ. 1) .AND. (ARG .LE. ZERO)) THEN
317 0 : RESULT = -XINF
318 0 : GO TO 2000
319 0 : ELSE IF (AX .GT. XMAX) THEN
320 0 : RESULT = ZERO
321 0 : GO TO 2000
322 : END IF
323 0 : IF (AX .GT. EIGHT) GO TO 800
324 0 : IF (AX .LE. XSMALL) THEN
325 0 : IF (JINT .EQ. 0) THEN
326 0 : RESULT = ONE
327 : ELSE
328 0 : RESULT = PI2 * (LOG(AX) + CONS)
329 : END IF
330 : GO TO 2000
331 : END IF
332 : !C-------------------------------------------------------------------
333 : !C Calculate J0 for appropriate interval, preserving
334 : !C accuracy near the zero of J0
335 : !C-------------------------------------------------------------------
336 0 : ZSQ = AX * AX
337 0 : IF (AX .LE. FOUR) THEN
338 0 : XNUM = (PJ0(5) * ZSQ + PJ0(6)) * ZSQ + PJ0(7)
339 0 : XDEN = ZSQ + QJ0(5)
340 0 : DO 50 I = 1, 4
341 0 : XNUM = XNUM * ZSQ + PJ0(I)
342 0 : XDEN = XDEN * ZSQ + QJ0(I)
343 0 : 50 CONTINUE
344 0 : PROD = ((AX - XJ01/TWO56) - XJ02) * (AX + XJ0)
345 : ELSE
346 0 : WSQ = ONE - ZSQ / SIXTY4
347 0 : XNUM = PJ1(7) * WSQ + PJ1(8)
348 0 : XDEN = WSQ + QJ1(7)
349 0 : DO 220 I = 1, 6
350 0 : XNUM = XNUM * WSQ + PJ1(I)
351 0 : XDEN = XDEN * WSQ + QJ1(I)
352 0 : 220 CONTINUE
353 0 : PROD = (AX + XJ1) * ((AX - XJ11/TWO56) - XJ12)
354 : END IF
355 0 : RESULT = PROD * XNUM / XDEN
356 0 : IF (JINT .EQ. 0) GO TO 2000
357 : !C-------------------------------------------------------------------
358 : !C Calculate Y0. First find RESJ = pi/2 ln(x/xn) J0(x),
359 : !C where xn is a zero of Y0
360 : !C-------------------------------------------------------------------
361 0 : IF (AX .LE. THREE) THEN
362 0 : UP = (AX-XY01/TWO56)-XY02
363 0 : XY = XY0
364 0 : ELSE IF (AX .LE. FIVE5) THEN
365 0 : UP = (AX-XY11/TWO56)-XY12
366 0 : XY = XY1
367 : ELSE
368 0 : UP = (AX-XY21/TWO56)-XY22
369 0 : XY = XY2
370 : END IF
371 0 : DOWN = AX + XY
372 0 : IF (ABS(UP) .LT. P17*DOWN) THEN
373 0 : W = UP/DOWN
374 0 : WSQ = W*W
375 0 : XNUM = PLG(1)
376 0 : XDEN = WSQ + QLG(1)
377 0 : DO 320 I = 2, 4
378 0 : XNUM = XNUM*WSQ + PLG(I)
379 0 : XDEN = XDEN*WSQ + QLG(I)
380 0 : 320 CONTINUE
381 0 : RESJ = PI2 * RESULT * W * XNUM/XDEN
382 : ELSE
383 0 : RESJ = PI2 * RESULT * LOG(AX/XY)
384 : END IF
385 : !C-------------------------------------------------------------------
386 : !C Now calculate Y0 for appropriate interval, preserving
387 : !C accuracy near the zero of Y0
388 : !C-------------------------------------------------------------------
389 0 : IF (AX .LE. THREE) THEN
390 0 : XNUM = PY0(6) * ZSQ + PY0(1)
391 0 : XDEN = ZSQ + QY0(1)
392 0 : DO 340 I = 2, 5
393 0 : XNUM = XNUM * ZSQ + PY0(I)
394 0 : XDEN = XDEN * ZSQ + QY0(I)
395 0 : 340 CONTINUE
396 0 : ELSE IF (AX .LE. FIVE5) THEN
397 0 : XNUM = PY1(7) * ZSQ + PY1(1)
398 0 : XDEN = ZSQ + QY1(1)
399 0 : DO 360 I = 2, 6
400 0 : XNUM = XNUM * ZSQ + PY1(I)
401 0 : XDEN = XDEN * ZSQ + QY1(I)
402 0 : 360 CONTINUE
403 : ELSE
404 0 : XNUM = PY2(8) * ZSQ + PY2(1)
405 0 : XDEN = ZSQ + QY2(1)
406 0 : DO 380 I = 2, 7
407 0 : XNUM = XNUM * ZSQ + PY2(I)
408 0 : XDEN = XDEN * ZSQ + QY2(I)
409 0 : 380 CONTINUE
410 : END IF
411 0 : RESULT = RESJ + UP * DOWN * XNUM / XDEN
412 0 : GO TO 2000
413 : !C-------------------------------------------------------------------
414 : !C Calculate J0 or Y0 for |ARG| > 8.0
415 : !C-------------------------------------------------------------------
416 0 : 800 Z = EIGHT / AX
417 0 : W = AX / TWOPI
418 0 : W = AINT(W) + ONEOV8
419 0 : W = (AX - W * TWOPI1) - W * TWOPI2
420 0 : ZSQ = Z * Z
421 0 : XNUM = P0(5) * ZSQ + P0(6)
422 0 : XDEN = ZSQ + Q0(5)
423 0 : UP = P1(5) * ZSQ + P1(6)
424 0 : DOWN = ZSQ + Q1(5)
425 0 : DO 850 I = 1, 4
426 0 : XNUM = XNUM * ZSQ + P0(I)
427 0 : XDEN = XDEN * ZSQ + Q0(I)
428 0 : UP = UP * ZSQ + P1(I)
429 0 : DOWN = DOWN * ZSQ + Q1(I)
430 0 : 850 CONTINUE
431 0 : R0 = XNUM / XDEN
432 0 : R1 = UP / DOWN
433 0 : IF (JINT .EQ. 0) THEN
434 0 : RESULT = SQRT(PI2/AX) * (R0*COS(W) - Z*R1*SIN(W))
435 : ELSE
436 0 : RESULT = SQRT(PI2/AX) * (R0*SIN(W) + Z*R1*COS(W))
437 : END IF
438 0 : 2000 RETURN
439 : !C---------- Last line of CALJY0 ----------
440 : END subroutine caljy0
441 :
442 : DOUBLE PRECISION FUNCTION BESJ0(X)
443 : !CS REAL FUNCTION BESJ0(X)
444 : !--------------------------------------------------------------------
445 : !
446 : ! This subprogram computes approximate values for Bessel functions
447 : ! of the first kind of order zero for arguments |X| <= XMAX
448 : ! (see comments heading CALJY0).
449 : !
450 : !--------------------------------------------------------------------
451 : !IMPLICIT NONE
452 : INTEGER :: JINT
453 : !S REAL X, RESULT
454 : DOUBLE PRECISION :: X, RESULT
455 : !--------------------------------------------------------------------
456 : JINT=0
457 : CALL CALJY0(X,RESULT,JINT)
458 : BESJ0 = RESULT
459 : RETURN
460 : !---------- Last line of BESJ0 ----------
461 : END function besj0
462 : DOUBLE PRECISION FUNCTION BESY0(X)
463 : !CS REAL FUNCTION BESY0(X)
464 : !C--------------------------------------------------------------------
465 : !C
466 : !C This subprogram computes approximate values for Bessel functions
467 : !C of the second kind of order zero for arguments 0 < X <= XMAX
468 : !C (see comments heading CALJY0).
469 : !C
470 : !C--------------------------------------------------------------------
471 : !IMPLICIT NONE
472 : INTEGER JINT
473 : !CS REAL X, RESULT
474 : DOUBLE PRECISION :: X, RESULT
475 : !--------------------------------------------------------------------
476 : JINT=1
477 : CALL CALJY0(X,RESULT,JINT)
478 : BESY0 = RESULT
479 : RETURN
480 : !---------- Last line of BESY0 ----------
481 : END function besy0
482 :
483 0 : SUBROUTINE CALJY1(ARG,RESULT,JINT)
484 :
485 : !---------------------------------------------------------------------
486 : !
487 : ! This packet computes first-order Bessel functions of the first and
488 : ! second kind (J1 and Y1), for real arguments X, where 0 < X <= XMAX
489 : ! for Y1, and |X| <= XMAX for J1. It contains two function-type
490 : ! subprograms, BESJ1 and BESY1, and one subroutine-type
491 : ! subprogram, CALJY1. The calling statements for the primary
492 : ! entries are:
493 : !
494 : ! Y = BESJ1(X)
495 : ! and
496 : ! Y = BESY1(X),
497 : !
498 : ! where the entry points correspond to the functions J1(X) and Y1(X),
499 : ! respectively. The routine CALJY1 is intended for internal packet
500 : ! use only, all computations within the packet being concentrated in
501 : ! this one routine. The function subprograms invoke CALJY1 with
502 : ! the statement
503 : ! CALL CALJY1(ARG,RESULT,JINT),
504 : ! where the parameter usage is as follows:
505 : !
506 : ! Function Parameters for CALJY1
507 : ! call ARG RESULT JINT
508 : !
509 : ! BESJ1(ARG) |ARG| .LE. XMAX J1(ARG) 0
510 : ! BESY1(ARG) 0 .LT. ARG .LE. XMAX Y1(ARG) 1
511 : !
512 : ! The main computation uses unpublished minimax rational
513 : ! approximations for X .LE. 8.0, and an approximation from the
514 : ! book Computer Approximations by Hart, et. al., Wiley and Sons,
515 : ! New York, 1968, for arguments larger than 8.0 Part of this
516 : ! transportable packet is patterned after the machine-dependent
517 : ! FUNPACK program BESJ1(X), but cannot match that version for
518 : ! efficiency or accuracy. This version uses rational functions
519 : ! that are theoretically accurate to at least 18 significant decimal
520 : ! digits for X <= 8, and at least 18 decimal places for X > 8. The
521 : ! accuracy achieved depends on the arithmetic system, the compiler,
522 : ! the intrinsic functions, and proper selection of the machine-
523 : ! dependent constants.
524 : !
525 : !*******************************************************************
526 : !
527 : ! Explanation of machine-dependent constants
528 : !
529 : ! XINF = largest positive machine number
530 : ! XMAX = largest acceptable argument. The functions AINT, SIN
531 : ! and COS must perform properly for ABS(X) .LE. XMAX.
532 : ! We recommend that XMAX be a small integer multiple of
533 : ! sqrt(1/eps), where eps is the smallest positive number
534 : ! such that 1+eps > 1.
535 : ! XSMALL = positive argument such that 1.0-(1/2)(X/2)**2 = 1.0
536 : ! to machine precision for all ABS(X) .LE. XSMALL.
537 : ! We recommend that XSMALL < sqrt(eps)/beta, where beta
538 : ! is the floating-point radix (usually 2 or 16).
539 : !
540 : ! Approximate values for some important machines are
541 : !
542 : ! eps XMAX XSMALL XINF
543 : !
544 : ! CDC 7600 (S.P.) 7.11E-15 1.34E+08 2.98E-08 1.26E+322
545 : ! CRAY-1 (S.P.) 7.11E-15 1.34E+08 2.98E-08 5.45E+2465
546 : ! IBM PC (8087) (S.P.) 5.96E-08 8.19E+03 1.22E-04 3.40E+38
547 : ! IBM PC (8087) (D.P.) 1.11D-16 2.68D+08 3.72D-09 1.79D+308
548 : ! IBM 195 (D.P.) 2.22D-16 6.87D+09 9.09D-13 7.23D+75
549 : ! UNIVAC 1108 (D.P.) 1.73D-18 4.30D+09 2.33D-10 8.98D+307
550 : ! VAX 11/780 (D.P.) 1.39D-17 1.07D+09 9.31D-10 1.70D+38
551 : !
552 : !*******************************************************************
553 : !*******************************************************************
554 : !
555 : ! Error Returns
556 : !
557 : ! The program returns the value zero for X .GT. XMAX, and returns
558 : ! -XINF when BESLY1 is called with a negative or zero argument.
559 : !
560 : !
561 : ! Intrinsic functions required are:
562 : !
563 : ! ABS, AINT, COS, LOG, SIN, SQRT
564 : !
565 : !
566 : ! Author: W. J. Cody
567 : ! Mathematics and Computer Science Division
568 : ! Argonne National Laboratory
569 : ! Argonne, IL 60439
570 : !
571 : ! Latest modification: November 10, 1987
572 : !
573 : ! Taken from http://www.netlib.org/specfun/j1y1
574 : !
575 : !--------------------------------------------------------------------
576 : !IMPLICIT NONE
577 : INTEGER :: I,JINT
578 : DIMENSION :: PJ0(7),PJ1(8),PLG(4),PY0(7),PY1(9),P0(6),P1(6),&
579 : & QJ0(5),QJ1(7),QLG(4),QY0(6),QY1(8),Q0(6),Q1(6)
580 : !CS REAL
581 : DOUBLE PRECISION :: &
582 : & ARG,AX,DOWN,EIGHT,FOUR,HALF,PI2,PJ0,PJ1,PLG,PROD,PY0,&
583 : & PY1,P0,P1,P17,QJ0,QJ1,QLG,QY0,QY1,Q0,Q1,RESJ,RESULT,&
584 : & RTPI2,R0,R1,THROV8,TWOPI,TWOPI1,TWOPI2,TWO56,UP,W,WSQ,&
585 : & XDEN,XINF,XMAX,XNUM,XSMALL,XJ0,XJ1,XJ01,XJ02,XJ11,XJ12,&
586 : & XY,XY0,XY01,XY02,XY1,XY11,XY12,Z,ZERO,ZSQ
587 : !-------------------------------------------------------------------
588 : ! Mathematical constants
589 : !-------------------------------------------------------------------
590 : !CS DATA EIGHT/8.0E0/,
591 : !CS 1 FOUR/4.0E0/,HALF/0.5E0/,THROV8/0.375E0/,
592 : !CS 2 PI2/6.3661977236758134308E-1/,P17/1.716E-1/
593 : !CS 3 TWOPI/6.2831853071795864769E+0/,ZERO/0.0E0/,
594 : !CS 4 TWOPI1/6.28125E0/,TWOPI2/1.9353071795864769253E-03/
595 : !CS 5 TWO56/256.0E+0/,RTPI2/7.9788456080286535588E-1/
596 : DATA EIGHT/8.0D0/,&
597 : & FOUR/4.0D0/,HALF/0.5D0/,THROV8/0.375D0/,&
598 : & PI2/6.3661977236758134308D-1/,P17/1.716D-1/&
599 : & TWOPI/6.2831853071795864769D+0/,ZERO/0.0D0/,&
600 : & TWOPI1/6.28125D0/,TWOPI2/1.9353071795864769253D-03/&
601 : & TWO56/256.0D+0/,RTPI2/7.9788456080286535588D-1/
602 : !-------------------------------------------------------------------
603 : ! Machine-dependent constants
604 : !-------------------------------------------------------------------
605 : !CS DATA XMAX/8.19E+03/,XSMALL/1.22E-09/,XINF/1.7E+38/
606 : DATA XMAX/1.07D+09/,XSMALL/9.31D-10/,XINF/1.7D+38/
607 : !-------------------------------------------------------------------
608 : ! Zeroes of Bessel functions
609 : !-------------------------------------------------------------------
610 : !CS DATA XJ0/3.8317059702075123156E+0/,XJ1/7.0155866698156187535E+0/,
611 : !CS 1 XY0/2.1971413260310170351E+0/,XY1/5.4296810407941351328E+0/,
612 : !CS 2 XJ01/ 981.0E+0/, XJ02/-3.2527979248768438556E-04/,
613 : !CS 3 XJ11/1796.0E+0/, XJ12/-3.8330184381246462950E-05/,
614 : !CS 4 XY01/ 562.0E+0/, XY02/ 1.8288260310170351490E-03/,
615 : !CS 5 XY11/1390.0E+0/, XY12/-6.4592058648672279948E-06/
616 : DATA XJ0/3.8317059702075123156D+0/,XJ1/7.0155866698156187535D+0/,&
617 : & XY0/2.1971413260310170351D+0/,XY1/5.4296810407941351328D+0/,&
618 : & XJ01/ 981.0D+0/, XJ02/-3.2527979248768438556D-04/,&
619 : & XJ11/1796.0D+0/, XJ12/-3.8330184381246462950D-05/,&
620 : & XY01/ 562.0D+0/, XY02/ 1.8288260310170351490D-03/,&
621 : & XY11/1390.0D+0/, XY12/-6.4592058648672279948D-06/
622 : !-------------------------------------------------------------------
623 : ! Coefficients for rational approximation to ln(x/a)
624 : !--------------------------------------------------------------------
625 : !CS DATA PLG/-2.4562334077563243311E+01,2.3642701335621505212E+02,
626 : !CS 1 -5.4989956895857911039E+02,3.5687548468071500413E+02/
627 : !CS DATA QLG/-3.5553900764052419184E+01,1.9400230218539473193E+02,
628 : !CS 1 -3.3442903192607538956E+02,1.7843774234035750207E+02/
629 : DATA PLG/-2.4562334077563243311D+01,2.3642701335621505212D+02,&
630 : & -5.4989956895857911039D+02,3.5687548468071500413D+02/
631 : DATA QLG/-3.5553900764052419184D+01,1.9400230218539473193D+02,&
632 : & -3.3442903192607538956D+02,1.7843774234035750207D+02/
633 : !-------------------------------------------------------------------
634 : ! Coefficients for rational approximation of
635 : ! J1(X) / (X * (X**2 - XJ0**2)), XSMALL < |X| <= 4.0
636 : !--------------------------------------------------------------------
637 : !CS DATA PJ0/9.8062904098958257677E+05,-1.1548696764841276794E+08,
638 : !CS 1 6.6781041261492395835E+09,-1.4258509801366645672E+11,
639 : !CS 2 -4.4615792982775076130E+03, 1.0650724020080236441E+01,
640 : !CS 3 -1.0767857011487300348E-02/
641 : !CS DATA QJ0/5.9117614494174794095E+05, 2.0228375140097033958E+08,
642 : !CS 1 4.2091902282580133541E+10, 4.1868604460820175290E+12,
643 : !CS 2 1.0742272239517380498E+03/
644 : DATA PJ0/9.8062904098958257677D+05,-1.1548696764841276794D+08,&
645 : & 6.6781041261492395835D+09,-1.4258509801366645672D+11,&
646 : & -4.4615792982775076130D+03, 1.0650724020080236441D+01,&
647 : & -1.0767857011487300348D-02/
648 : DATA QJ0/5.9117614494174794095D+05, 2.0228375140097033958D+08,&
649 : & 4.2091902282580133541D+10, 4.1868604460820175290D+12,&
650 : & 1.0742272239517380498D+03/
651 : !-------------------------------------------------------------------
652 : ! Coefficients for rational approximation of
653 : ! J1(X) / (X * (X**2 - XJ1**2)), 4.0 < |X| <= 8.0
654 : !-------------------------------------------------------------------
655 : !CS DATA PJ1/4.6179191852758252280E+00,-7.1329006872560947377E+03,
656 : !CS 1 4.5039658105749078904E+06,-1.4437717718363239107E+09,
657 : !CS 2 2.3569285397217157313E+11,-1.6324168293282543629E+13,
658 : !CS 3 1.1357022719979468624E+14, 1.0051899717115285432E+15/
659 : !CS DATA QJ1/1.1267125065029138050E+06, 6.4872502899596389593E+08,
660 : !CS 1 2.7622777286244082666E+11, 8.4899346165481429307E+13,
661 : !CS 2 1.7128800897135812012E+16, 1.7253905888447681194E+18,
662 : !CS 3 1.3886978985861357615E+03/
663 : DATA PJ1/4.6179191852758252280D+00,-7.1329006872560947377D+03,&
664 : & 4.5039658105749078904D+06,-1.4437717718363239107D+09,&
665 : & 2.3569285397217157313D+11,-1.6324168293282543629D+13,&
666 : & 1.1357022719979468624D+14, 1.0051899717115285432D+15/
667 : DATA QJ1/1.1267125065029138050D+06, 6.4872502899596389593D+08,&
668 : & 2.7622777286244082666D+11, 8.4899346165481429307D+13,&
669 : & 1.7128800897135812012D+16, 1.7253905888447681194D+18,&
670 : & 1.3886978985861357615D+03/
671 : !-------------------------------------------------------------------
672 : ! Coefficients for rational approximation of
673 : ! (Y1(X) - 2 LN(X/XY0) J1(X)) / (X**2 - XY0**2),
674 : ! XSMALL < |X| <= 4.0
675 : !--------------------------------------------------------------------
676 : !CS DATA PY0/2.2157953222280260820E+05,-5.9157479997408395984E+07,
677 : !CS 1 7.2144548214502560419E+09,-3.7595974497819597599E+11,
678 : !CS 2 5.4708611716525426053E+12, 4.0535726612579544093E+13,
679 : !CS 3 -3.1714424660046133456E+02/
680 : !CS DATA QY0/8.2079908168393867438E+02, 3.8136470753052572164E+05,
681 : !CS 1 1.2250435122182963220E+08, 2.7800352738690585613E+10,
682 : !CS 2 4.1272286200406461981E+12, 3.0737873921079286084E+14/
683 : DATA PY0/2.2157953222280260820D+05,-5.9157479997408395984D+07,&
684 : & 7.2144548214502560419D+09,-3.7595974497819597599D+11,&
685 : & 5.4708611716525426053D+12, 4.0535726612579544093D+13,&
686 : & -3.1714424660046133456D+02/
687 : DATA QY0/8.2079908168393867438D+02, 3.8136470753052572164D+05,&
688 : & 1.2250435122182963220D+08, 2.7800352738690585613D+10,&
689 : & 4.1272286200406461981D+12, 3.0737873921079286084D+14/
690 : !--------------------------------------------------------------------
691 : ! Coefficients for rational approximation of
692 : ! (Y1(X) - 2 LN(X/XY1) J1(X)) / (X**2 - XY1**2),
693 : ! 4.0 < |X| <= 8.0
694 : !--------------------------------------------------------------------
695 : !CS DATA PY1/ 1.9153806858264202986E+06,-1.1957961912070617006E+09,
696 : !CS 1 3.7453673962438488783E+11,-5.9530713129741981618E+13,
697 : !CS 2 4.0686275289804744814E+15,-2.3638408497043134724E+16,
698 : !CS 3 -5.6808094574724204577E+18, 1.1514276357909013326E+19,
699 : !CS 4 -1.2337180442012953128E+03/
700 : !CS DATA QY1/ 1.2855164849321609336E+03, 1.0453748201934079734E+06,
701 : !CS 1 6.3550318087088919566E+08, 3.0221766852960403645E+11,
702 : !CS 2 1.1187010065856971027E+14, 3.0837179548112881950E+16,
703 : !CS 3 5.6968198822857178911E+18, 5.3321844313316185697E+20/
704 : DATA PY1/ 1.9153806858264202986D+06,-1.1957961912070617006D+09,&
705 : & 3.7453673962438488783D+11,-5.9530713129741981618D+13,&
706 : & 4.0686275289804744814D+15,-2.3638408497043134724D+16,&
707 : & -5.6808094574724204577D+18, 1.1514276357909013326D+19,&
708 : & -1.2337180442012953128D+03/
709 : DATA QY1/ 1.2855164849321609336D+03, 1.0453748201934079734D+06,&
710 : & 6.3550318087088919566D+08, 3.0221766852960403645D+11,&
711 : & 1.1187010065856971027D+14, 3.0837179548112881950D+16,&
712 : & 5.6968198822857178911D+18, 5.3321844313316185697D+20/
713 : !-------------------------------------------------------------------
714 : ! Coefficients for Hart,s approximation, |X| > 8.0
715 : !-------------------------------------------------------------------
716 : !CS DATA P0/-1.0982405543459346727E+05,-1.5235293511811373833E+06,
717 : !CS 1 -6.6033732483649391093E+06,-9.9422465050776411957E+06,
718 : !CS 2 -4.4357578167941278571E+06,-1.6116166443246101165E+03/
719 : !CS DATA Q0/-1.0726385991103820119E+05,-1.5118095066341608816E+06,
720 : !CS 1 -6.5853394797230870728E+06,-9.9341243899345856590E+06,
721 : !CS 2 -4.4357578167941278568E+06,-1.4550094401904961825E+03/
722 : !CS DATA P1/ 1.7063754290207680021E+03, 1.8494262873223866797E+04,
723 : !CS 1 6.6178836581270835179E+04, 8.5145160675335701966E+04,
724 : !CS 2 3.3220913409857223519E+04, 3.5265133846636032186E+01/
725 : !CS DATA Q1/ 3.7890229745772202641E+04, 4.0029443582266975117E+05,
726 : !CS 1 1.4194606696037208929E+06, 1.8194580422439972989E+06,
727 : !CS 2 7.0871281941028743574E+05, 8.6383677696049909675E+02/
728 : DATA P0/-1.0982405543459346727D+05,-1.5235293511811373833D+06,&
729 : & -6.6033732483649391093D+06,-9.9422465050776411957D+06,&
730 : & -4.4357578167941278571D+06,-1.6116166443246101165D+03/
731 : DATA Q0/-1.0726385991103820119D+05,-1.5118095066341608816D+06,&
732 : & -6.5853394797230870728D+06,-9.9341243899345856590D+06,&
733 : & -4.4357578167941278568D+06,-1.4550094401904961825D+03/
734 : DATA P1/ 1.7063754290207680021D+03, 1.8494262873223866797D+04,&
735 : & 6.6178836581270835179D+04, 8.5145160675335701966D+04,&
736 : & 3.3220913409857223519D+04, 3.5265133846636032186D+01/
737 : DATA Q1/ 3.7890229745772202641D+04, 4.0029443582266975117D+05,&
738 : & 1.4194606696037208929D+06, 1.8194580422439972989D+06,&
739 : & 7.0871281941028743574D+05, 8.6383677696049909675D+02/
740 : !-------------------------------------------------------------------
741 : ! Check for error conditions
742 : !-------------------------------------------------------------------
743 0 : AX = ABS(ARG)
744 0 : IF ((JINT .EQ. 1) .AND. ((ARG .LE. ZERO) .OR.&
745 : & ((ARG .LT. HALF) .AND. (AX*XINF .LT. PI2)))) THEN
746 0 : RESULT = -XINF
747 0 : GO TO 2000
748 0 : ELSE IF (AX .GT. XMAX) THEN
749 0 : RESULT = ZERO
750 0 : GO TO 2000
751 : END IF
752 0 : IF (AX .GT. EIGHT) THEN
753 : GO TO 800
754 0 : ELSE IF (AX .LE. XSMALL) THEN
755 0 : IF (JINT .EQ. 0) THEN
756 0 : RESULT = ARG * HALF
757 : ELSE
758 0 : RESULT = -PI2 / AX
759 : END IF
760 : GO TO 2000
761 : END IF
762 : !-------------------------------------------------------------------
763 : ! Calculate J1 for appropriate interval, preserving
764 : ! accuracy near the zero of J1
765 : !-------------------------------------------------------------------
766 0 : ZSQ = AX * AX
767 0 : IF (AX .LE. FOUR) THEN
768 0 : XNUM = (PJ0(7) * ZSQ + PJ0(6)) * ZSQ + PJ0(5)
769 0 : XDEN = ZSQ + QJ0(5)
770 0 : DO 50 I = 1, 4
771 0 : XNUM = XNUM * ZSQ + PJ0(I)
772 0 : XDEN = XDEN * ZSQ + QJ0(I)
773 0 : 50 CONTINUE
774 0 : PROD = ARG * ((AX - XJ01/TWO56) - XJ02) * (AX + XJ0)
775 : ELSE
776 0 : XNUM = PJ1(1)
777 0 : XDEN = (ZSQ + QJ1(7)) * ZSQ + QJ1(1)
778 0 : DO 220 I = 2, 6
779 0 : XNUM = XNUM * ZSQ + PJ1(I)
780 0 : XDEN = XDEN * ZSQ + QJ1(I)
781 0 : 220 CONTINUE
782 0 : XNUM = XNUM * (AX - EIGHT) * (AX + EIGHT) + PJ1(7)
783 0 : XNUM = XNUM * (AX - FOUR) * (AX + FOUR) + PJ1(8)
784 0 : PROD = ARG * ((AX - XJ11/TWO56) - XJ12) * (AX + XJ1)
785 : END IF
786 0 : RESULT = PROD * (XNUM / XDEN)
787 0 : IF (JINT .EQ. 0) GO TO 2000
788 : !-------------------------------------------------------------------
789 : ! Calculate Y1. First find RESJ = pi/2 ln(x/xn) J1(x),
790 : ! where xn is a zero of Y1
791 : !-------------------------------------------------------------------
792 0 : IF (AX .LE. FOUR) THEN
793 0 : UP = (AX-XY01/TWO56)-XY02
794 0 : XY = XY0
795 : ELSE
796 0 : UP = (AX-XY11/TWO56)-XY12
797 0 : XY = XY1
798 : END IF
799 0 : DOWN = AX + XY
800 0 : IF (ABS(UP) .LT. P17*DOWN) THEN
801 0 : W = UP/DOWN
802 0 : WSQ = W*W
803 0 : XNUM = PLG(1)
804 0 : XDEN = WSQ + QLG(1)
805 0 : DO 320 I = 2, 4
806 0 : XNUM = XNUM*WSQ + PLG(I)
807 0 : XDEN = XDEN*WSQ + QLG(I)
808 0 : 320 CONTINUE
809 0 : RESJ = PI2 * RESULT * W * XNUM/XDEN
810 : ELSE
811 0 : RESJ = PI2 * RESULT * LOG(AX/XY)
812 : END IF
813 : !-------------------------------------------------------------------
814 : ! Now calculate Y1 for appropriate interval, preserving
815 : ! accuracy near the zero of Y1
816 : !-------------------------------------------------------------------
817 0 : IF (AX .LE. FOUR) THEN
818 0 : XNUM = PY0(7) * ZSQ + PY0(1)
819 0 : XDEN = ZSQ + QY0(1)
820 0 : DO 340 I = 2, 6
821 0 : XNUM = XNUM * ZSQ + PY0(I)
822 0 : XDEN = XDEN * ZSQ + QY0(I)
823 0 : 340 CONTINUE
824 : ELSE
825 0 : XNUM = PY1(9) * ZSQ + PY1(1)
826 0 : XDEN = ZSQ + QY1(1)
827 0 : DO 360 I = 2, 8
828 0 : XNUM = XNUM * ZSQ + PY1(I)
829 0 : XDEN = XDEN * ZSQ + QY1(I)
830 0 : 360 CONTINUE
831 : END IF
832 0 : RESULT = RESJ + (UP*DOWN/AX) * XNUM / XDEN
833 0 : GO TO 2000
834 : !-------------------------------------------------------------------
835 : ! Calculate J1 or Y1 for |ARG| > 8.0
836 : !-------------------------------------------------------------------
837 0 : 800 Z = EIGHT / AX
838 0 : W = AINT(AX/TWOPI) + THROV8
839 0 : W = (AX - W * TWOPI1) - W * TWOPI2
840 0 : ZSQ = Z * Z
841 0 : XNUM = P0(6)
842 0 : XDEN = ZSQ + Q0(6)
843 0 : UP = P1(6)
844 0 : DOWN = ZSQ + Q1(6)
845 0 : DO 850 I = 1, 5
846 0 : XNUM = XNUM * ZSQ + P0(I)
847 0 : XDEN = XDEN * ZSQ + Q0(I)
848 0 : UP = UP * ZSQ + P1(I)
849 0 : DOWN = DOWN * ZSQ + Q1(I)
850 0 : 850 CONTINUE
851 0 : R0 = XNUM / XDEN
852 0 : R1 = UP / DOWN
853 0 : IF (JINT .EQ. 0) THEN
854 0 : RESULT = (RTPI2/SQRT(AX)) * (R0*COS(W) - Z*R1*SIN(W))
855 : ELSE
856 0 : RESULT = (RTPI2/SQRT(AX)) * (R0*SIN(W) + Z*R1*COS(W))
857 : END IF
858 0 : IF ((JINT .EQ. 0) .AND. (ARG .LT. ZERO)) RESULT = -RESULT
859 0 : 2000 RETURN
860 : !---------- Last card of CALJY1 ----------
861 : END subroutine caljy1
862 :
863 : DOUBLE PRECISION FUNCTION BESJ1(X)
864 : !--------------------------------------------------------------------
865 : !
866 : ! This subprogram computes approximate values for Bessel functions
867 : ! of the first kind of order zero for arguments |X| <= XMAX
868 : ! (see comments heading CALJY1).
869 : !
870 : !--------------------------------------------------------------------
871 : !IMPLICIT NONE
872 : INTEGER :: JINT
873 : !CS REAL
874 : DOUBLE PRECISION :: &
875 : & RESULT,X
876 : !--------------------------------------------------------------------
877 : JINT=0
878 : CALL CALJY1(X,RESULT,JINT)
879 : BESJ1 = RESULT
880 : RETURN
881 : !---------- Last card of BESJ1 ----------
882 : END function besj1
883 :
884 :
885 : DOUBLE PRECISION FUNCTION BESY1(X)
886 : !--------------------------------------------------------------------
887 : !
888 : ! This subprogram computes approximate values for Bessel functions
889 : ! of the second kind of order zero for arguments 0 < X <= XMAX
890 : ! (see comments heading CALJY1).
891 : !
892 : !--------------------------------------------------------------------
893 : !IMPLICIT NONE
894 : INTEGER :: JINT
895 : !CS REAL
896 : DOUBLE PRECISION :: &
897 : & RESULT,X
898 : !--------------------------------------------------------------------
899 : JINT=1
900 : CALL CALJY1(X,RESULT,JINT)
901 : BESY1 = RESULT
902 : RETURN
903 : !---------- Last card of BESY1 ----------
904 : END function besy1
905 : !!***
906 :
907 : !!****f* ABINIT/CALCK0
908 : !! NAME
909 : !!
910 : !! FUNCTION
911 : !!
912 : !! INPUTS
913 : !!
914 : !! OUTPUT
915 : !!
916 : !! SIDE EFFECTS
917 : !!
918 : !! NOTES
919 : !!
920 : !! SOURCE
921 :
922 5345190 : SUBROUTINE CALCK0(ARG,RESULT,JINT)
923 :
924 : !--------------------------------------------------------------------
925 : !
926 : ! This packet computes modified Bessel functions of the second kind
927 : ! and order zero, K0(X) and EXP(X)*K0(X), for real
928 : ! arguments X. It contains two function type subprograms, BESK0
929 : ! and BESEK0, and one subroutine type subprogram, CALCK0.
930 : ! the calling statements for the primary entries are
931 : !
932 : ! Y=BESK0(X)
933 : ! and
934 : ! Y=BESEK0(X)
935 : !
936 : ! where the entry points correspond to the functions K0(X) and
937 : ! EXP(X)*K0(X), respectively. The routine CALCK0 is
938 : ! intended for internal packet use only, all computations within
939 : ! the packet being concentrated in this routine. The function
940 : ! subprograms invoke CALCK0 with the statement
941 : ! CALL CALCK0(ARG,RESULT,JINT)
942 : ! where the parameter usage is as follows
943 : !
944 : ! Function Parameters for CALCK0
945 : ! Call ARG RESULT JINT
946 : !
947 : ! BESK0(ARG) 0 .LT. ARG .LE. XMAX K0(ARG) 1
948 : ! BESEK0(ARG) 0 .LT. ARG EXP(ARG)*K0(ARG) 2
949 : !
950 : ! The main computation evaluates slightly modified forms of near
951 : ! minimax rational approximations generated by Russon and Blair,
952 : ! Chalk River (Atomic Energy of Canada Limited) Report AECL-3461,
953 : ! 1969. This transportable program is patterned after the
954 : ! machine-dependent FUNPACK packet NATSK0, but cannot match that
955 : ! version for efficiency or accuracy. This version uses rational
956 : ! functions that theoretically approximate K-SUB-0(X) to at
957 : ! least 18 significant decimal digits. The accuracy achieved
958 : ! depends on the arithmetic system, the compiler, the intrinsic
959 : ! functions, and proper selection of the machine-dependent
960 : ! constants.
961 : !
962 : !*******************************************************************
963 : !*******************************************************************
964 : !
965 : ! Explanation of machine-dependent constants
966 : !
967 : ! beta = Radix for the floating-point system
968 : ! minexp = Smallest representable power of beta
969 : ! maxexp = Smallest power of beta that overflows
970 : ! XSMALL = Argument below which BESK0 and BESEK0 may
971 : ! each be represented by a constant and a log.
972 : ! largest X such that 1.0 + X = 1.0 to machine
973 : ! precision.
974 : ! XINF = Largest positive machine number; approximately
975 : ! beta**maxexp
976 : ! XMAX = Largest argument acceptable to BESK0; Solution to
977 : ! equation:
978 : ! W(X) * (1-1/8X+9/128X**2) = beta**minexp
979 : ! where W(X) = EXP(-X)*SQRT(PI/2X)
980 : !
981 : !
982 : ! Approximate values for some important machines are:
983 : !
984 : !
985 : ! beta minexp maxexp
986 : !
987 : ! CRAY-1 (S.P.) 2 -8193 8191
988 : ! Cyber 180/185
989 : ! under NOS (S.P.) 2 -975 1070
990 : ! IEEE (IBM/XT,
991 : ! SUN, etc.) (S.P.) 2 -126 128
992 : ! IEEE (IBM/XT,
993 : ! SUN, etc.) (D.P.) 2 -1022 1024
994 : ! IBM 3033 (D.P.) 16 -65 63
995 : ! VAX D-Format (D.P.) 2 -128 127
996 : ! VAX G-Format (D.P.) 2 -1024 1023
997 : !
998 : !
999 : ! XSMALL XINF XMAX
1000 : !
1001 : ! CRAY-1 (S.P.) 3.55E-15 5.45E+2465 5674.858
1002 : ! Cyber 180/855
1003 : ! under NOS (S.P.) 1.77E-15 1.26E+322 672.788
1004 : ! IEEE (IBM/XT,
1005 : ! SUN, etc.) (S.P.) 5.95E-8 3.40E+38 85.337
1006 : ! IEEE (IBM/XT,
1007 : ! SUN, etc.) (D.P.) 1.11D-16 1.79D+308 705.342
1008 : ! IBM 3033 (D.P.) 1.11D-16 7.23D+75 177.852
1009 : ! VAX D-Format (D.P.) 6.95D-18 1.70D+38 86.715
1010 : ! VAX G-Format (D.P.) 5.55D-17 8.98D+307 706.728
1011 : !
1012 : !*******************************************************************
1013 : !*******************************************************************
1014 : !
1015 : ! Error returns
1016 : !
1017 : ! The program returns the value XINF for ARG .LE. 0.0, and the
1018 : ! BESK0 entry returns the value 0.0 for ARG .GT. XMAX.
1019 : !
1020 : !
1021 : ! Intrinsic functions required are:
1022 : !
1023 : ! EXP, LOG, SQRT
1024 : !
1025 : ! Latest modification: March 19, 1990
1026 : !
1027 : ! Authors: W. J. Cody and Laura Stoltz
1028 : ! Mathematics and Computer Science Division
1029 : ! Argonne National Laboratory
1030 : ! Argonne, IL 60439
1031 : !
1032 : ! Original subroutine from netlib http://www.netlib.org/specfun/k0
1033 : ! Slightly modified by MG to follow f90 rules and double precision arithmetic
1034 : !
1035 : !--------------------------------------------------------------------
1036 : !IMPLICIT NONE
1037 : INTEGER :: I,JINT
1038 : !CS REAL
1039 : DOUBLE PRECISION :: ARG,RESULT,SUMF,SUMG,SUMP,SUMQ,TEMP
1040 : !CS REAL
1041 : DOUBLE PRECISION :: X,XX
1042 : !CS REAL
1043 : DOUBLE PRECISION :: P(6),Q(2),PP(10),QQ(10),F(4),G(3)
1044 : !C--------------------------------------------------------------------
1045 : !C Mathematical constants
1046 : !C--------------------------------------------------------------------
1047 : !CS REAL, PARAMETER :: ONE=1.0E0,ZERO=0.0E0
1048 : DOUBLE PRECISION,PARAMETER :: ONE=1.0D0,ZERO=0.0D0
1049 : !C--------------------------------------------------------------------
1050 : !C Machine-dependent constants
1051 : !C--------------------------------------------------------------------
1052 : !CS REAL.PARAMETER :: XSMALL=5.95E-8, XINF=3.40E+38 ,XMAX=85.337E0
1053 : DOUBLE PRECISION,PARAMETER :: XSMALL=1.11D-16,XINF=1.79D+308,XMAX=705.342D0
1054 : !--------------------------------------------------------------------
1055 : !
1056 : ! Coefficients for XSMALL .LE. ARG .LE. 1.0
1057 : !
1058 : !--------------------------------------------------------------------
1059 : !S DATA P/ 5.8599221412826100000E-04, 1.3166052564989571850E-01,
1060 : !S 1 1.1999463724910714109E+01, 4.6850901201934832188E+02,
1061 : !S 2 5.9169059852270512312E+03, 2.4708152720399552679E+03/
1062 : !S DATA Q/-2.4994418972832303646E+02, 2.1312714303849120380E+04/
1063 : !S DATA F/-1.6414452837299064100E+00,-2.9601657892958843866E+02,
1064 : !S 1 -1.7733784684952985886E+04,-4.0320340761145482298E+05/
1065 : !S DATA G/-2.5064972445877992730E+02, 2.9865713163054025489E+04,
1066 : !S 1 -1.6128136304458193998E+06/
1067 : DATA P/5.8599221412826100000D-04,1.3166052564989571850D-01,&
1068 : & 1.1999463724910714109D+01,4.6850901201934832188D+02,&
1069 : & 5.9169059852270512312D+03,2.4708152720399552679D+03/
1070 : DATA Q/-2.4994418972832303646D+02, 2.1312714303849120380D+04/
1071 : DATA F/-1.6414452837299064100D+00,-2.9601657892958843866D+02,&
1072 : & -1.7733784684952985886D+04,-4.0320340761145482298D+05/
1073 : DATA G/-2.5064972445877992730D+02, 2.9865713163054025489D+04,&
1074 : & -1.6128136304458193998D+06/
1075 : !--------------------------------------------------------------------
1076 : !
1077 : ! Coefficients for 1.0 .LT. ARG
1078 : !
1079 : !--------------------------------------------------------------------
1080 : !S DATA PP/ 1.1394980557384778174E+02, 3.6832589957340267940E+03,
1081 : !S 1 3.1075408980684392399E+04, 1.0577068948034021957E+05,
1082 : !S 2 1.7398867902565686251E+05, 1.5097646353289914539E+05,
1083 : !S 3 7.1557062783764037541E+04, 1.8321525870183537725E+04,
1084 : !S 4 2.3444738764199315021E+03, 1.1600249425076035558E+02/
1085 : !S DATA QQ/ 2.0013443064949242491E+02, 4.4329628889746408858E+03,
1086 : !S 1 3.1474655750295278825E+04, 9.7418829762268075784E+04,
1087 : !S 2 1.5144644673520157801E+05, 1.2689839587977598727E+05,
1088 : !S 3 5.8824616785857027752E+04, 1.4847228371802360957E+04,
1089 : !S 4 1.8821890840982713696E+03, 9.2556599177304839811E+01/
1090 : DATA PP/ 1.1394980557384778174D+02, 3.6832589957340267940D+03,&
1091 : & 3.1075408980684392399D+04, 1.0577068948034021957D+05,&
1092 : & 1.7398867902565686251D+05, 1.5097646353289914539D+05,&
1093 : & 7.1557062783764037541D+04, 1.8321525870183537725D+04,&
1094 : & 2.3444738764199315021D+03, 1.1600249425076035558D+02/
1095 : DATA QQ/ 2.0013443064949242491D+02, 4.4329628889746408858D+03, &
1096 : & 3.1474655750295278825D+04, 9.7418829762268075784D+04,&
1097 : & 1.5144644673520157801D+05, 1.2689839587977598727D+05,&
1098 : & 5.8824616785857027752D+04, 1.4847228371802360957D+04,&
1099 : & 1.8821890840982713696D+03, 9.2556599177304839811D+01/
1100 : !--------------------------------------------------------------------
1101 5345190 : X = ARG
1102 5345190 : IF (X .GT. ZERO) THEN
1103 5345190 : IF (X .LE. ONE) THEN
1104 : !--------------------------------------------------------------------
1105 : ! 0.0 .LT. ARG .LE. 1.0
1106 : !--------------------------------------------------------------------
1107 1708190 : TEMP = LOG(X)
1108 1708190 : IF (X .LT. XSMALL) THEN
1109 : !--------------------------------------------------------------------
1110 : ! Return for small ARG
1111 : !--------------------------------------------------------------------
1112 0 : RESULT = P(6)/Q(2) - TEMP
1113 : ELSE
1114 1708190 : XX = X * X
1115 : SUMP = ((((P(1)*XX + P(2))*XX + P(3))*XX +&
1116 1708190 : P(4))*XX + P(5))*XX + P(6)
1117 1708190 : SUMQ = (XX + Q(1))*XX + Q(2)
1118 1708190 : SUMF = ((F(1)*XX + F(2))*XX + F(3))*XX + F(4)
1119 1708190 : SUMG = ((XX + G(1))*XX + G(2))*XX + G(3)
1120 1708190 : RESULT = SUMP/SUMQ - XX*SUMF*TEMP/SUMG - TEMP
1121 1708190 : IF (JINT .EQ. 2) RESULT = RESULT * EXP(X)
1122 : END IF
1123 3637000 : ELSE IF ((JINT .EQ. 1) .AND. (X .GT. XMAX)) THEN
1124 : !--------------------------------------------------------------------
1125 : ! Error return for ARG .GT. XMAX
1126 : !--------------------------------------------------------------------
1127 0 : RESULT = ZERO
1128 : ELSE
1129 : !--------------------------------------------------------------------
1130 : ! 1.0 .LT. ARG
1131 : !--------------------------------------------------------------------
1132 3637000 : XX = ONE / X
1133 3637000 : SUMP = PP(1)
1134 36370000 : DO 120 I = 2, 10
1135 32733000 : SUMP = SUMP*XX + PP(I)
1136 3637000 : 120 CONTINUE
1137 : SUMQ = XX
1138 36370000 : DO 140 I = 1, 9
1139 32733000 : SUMQ = (SUMQ + QQ(I))*XX
1140 3637000 : 140 CONTINUE
1141 3637000 : SUMQ = SUMQ + QQ(10)
1142 3637000 : RESULT = SUMP / SUMQ / SQRT(X)
1143 3637000 : IF (JINT .EQ. 1) RESULT = RESULT * EXP(-X)
1144 : END IF
1145 : ELSE
1146 : !--------------------------------------------------------------------
1147 : ! Error return for ARG .LE. 0.0
1148 : !--------------------------------------------------------------------
1149 0 : RESULT = XINF
1150 : END IF
1151 : !--------------------------------------------------------------------
1152 : ! Update error counts, etc.
1153 : !--------------------------------------------------------------------
1154 5345190 : RETURN
1155 : !---------- Last line of CALCK0 ----------
1156 : END subroutine calck0
1157 : !!***
1158 :
1159 : !S REAL
1160 : DOUBLE PRECISION FUNCTION BESK0(X)
1161 : !--------------------------------------------------------------------
1162 : !
1163 : ! This function program computes approximate values for the
1164 : ! modified Bessel function of the second kind of order zero
1165 : ! for arguments 0.0 .LT. ARG .LE. XMAX (see comments heading
1166 : ! CALCK0).
1167 : !
1168 : ! Authors: W. J. Cody and Laura Stoltz
1169 : !
1170 : ! Latest Modification: January 19, 1988
1171 : !
1172 : !--------------------------------------------------------------------
1173 : !IMPLICIT NONE
1174 : INTEGER :: JINT
1175 : !S REAL
1176 : DOUBLE PRECISION :: X, RESULT
1177 : !--------------------------------------------------------------------
1178 : JINT = 1
1179 : CALL CALCK0(X,RESULT,JINT)
1180 : BESK0 = RESULT
1181 : RETURN
1182 : !---------- Last line of BESK0 ----------
1183 : END function besk0
1184 : !!***
1185 :
1186 : !S REAL
1187 : DOUBLE PRECISION FUNCTION BESEK0(X)
1188 : !--------------------------------------------------------------------
1189 : !
1190 : ! This function program computes approximate values for the
1191 : ! modified Bessel function of the second kind of order zero
1192 : ! multiplied by the Exponential function, for arguments
1193 : ! 0.0 .LT. ARG.
1194 : !
1195 : ! Authors: W. J. Cody and Laura Stoltz
1196 : !
1197 : ! Latest Modification: January 19, 1988
1198 : !
1199 : !--------------------------------------------------------------------
1200 : !IMPLICIT NONE
1201 : INTEGER JINT
1202 : !S REAL
1203 : DOUBLE PRECISION :: X,RESULT
1204 : !--------------------------------------------------------------------
1205 : JINT = 2
1206 : CALL CALCK0(X,RESULT,JINT)
1207 : BESEK0 = RESULT
1208 : RETURN
1209 : !---------- Last line of BESEK0 ----------
1210 : END function BESEK0
1211 : !!***
1212 :
1213 : !!****f* ABINIT/CALCK1
1214 : !! NAME
1215 : !! CALCK1
1216 : !!
1217 : !! FUNCTION
1218 : !!
1219 : !! INPUTS
1220 : !!
1221 : !! OUTPUT
1222 : !!
1223 : !! SIDE EFFECTS
1224 : !!
1225 : !! NOTES
1226 : !!
1227 : !! SOURCE
1228 :
1229 0 : SUBROUTINE CALCK1(ARG,RESULT,JINT)
1230 :
1231 : !--------------------------------------------------------------------
1232 : !
1233 : ! This packet computes modified Bessel functions of the second kind
1234 : ! and order one, K1(X) and EXP(X)*K1(X), for real arguments X.
1235 : ! It contains two function type subprograms, BESK1 and BESEK1,
1236 : ! and one subroutine type subprogram, CALCK1. The calling
1237 : ! statements for the primary entries are
1238 : !
1239 : ! Y=BESK1(X)
1240 : ! and
1241 : ! Y=BESEK1(X)
1242 : !
1243 : ! where the entry points correspond to the functions K1(X) and
1244 : ! EXP(X)*K1(X), respectively. The routine CALCK1 is intended
1245 : ! for internal packet use only, all computations within the
1246 : ! packet being concentrated in this routine. The function
1247 : ! subprograms invoke CALCK1 with the statement
1248 : ! CALL CALCK1(ARG,RESULT,JINT)
1249 : ! where the parameter usage is as follows
1250 : !
1251 : ! Function Parameters for CALCK1
1252 : ! Call ARG RESULT JINT
1253 : !
1254 : ! BESK1(ARG) XLEAST .LT. ARG .LT. XMAX K1(ARG) 1
1255 : ! BESEK1(ARG) XLEAST .LT. ARG EXP(ARG)*K1(ARG) 2
1256 : !
1257 : ! The main computation evaluates slightly modified forms of near
1258 : ! minimax rational approximations generated by Russon and Blair,
1259 : ! Chalk River (Atomic Energy of Canada Limited) Report AECL-3461,
1260 : ! 1969. This transportable program is patterned after the
1261 : ! machine-dependent FUNPACK packet NATSK1, but cannot match that
1262 : ! version for efficiency or accuracy. This version uses rational
1263 : ! functions that theoretically approximate K-SUB-1(X) to at
1264 : ! least 18 significant decimal digits. The accuracy achieved
1265 : ! depends on the arithmetic system, the compiler, the intrinsic
1266 : ! functions, and proper selection of the machine-dependent
1267 : ! constants.
1268 : !
1269 : !*******************************************************************
1270 : !*******************************************************************
1271 : !
1272 : ! Explanation of machine-dependent constants
1273 : !
1274 : ! beta = Radix for the floating-point system
1275 : ! minexp = Smallest representable power of beta
1276 : ! maxexp = Smallest power of beta that overflows
1277 : ! XLEAST = Smallest acceptable argument, i.e., smallest machine
1278 : ! number X such that 1/X is machine representable.
1279 : ! XSMALL = Argument below which BESK1(X) and BESEK1(X) may
1280 : ! each be represented by 1/X. A safe value is the
1281 : ! largest X such that 1.0 + X = 1.0 to machine
1282 : ! precision.
1283 : ! XINF = Largest positive machine number; approximately
1284 : ! beta**maxexp
1285 : ! XMAX = Largest argument acceptable to BESK1; Solution to
1286 : ! equation:
1287 : ! W(X) * (1+3/8X-15/128X**2) = beta**minexp
1288 : ! where W(X) = EXP(-X)*SQRT(PI/2X)
1289 : !
1290 : !
1291 : ! Approximate values for some important machines are:
1292 : !
1293 : ! beta minexp maxexp
1294 : !
1295 : ! CRAY-1 (S.P.) 2 -8193 8191
1296 : ! Cyber 180/185
1297 : ! under NOS (S.P.) 2 -975 1070
1298 : ! IEEE (IBM/XT,
1299 : ! SUN, etc.) (S.P.) 2 -126 128
1300 : ! IEEE (IBM/XT,
1301 : ! SUN, etc.) (D.P.) 2 -1022 1024
1302 : ! IBM 3033 (D.P.) 16 -65 63
1303 : ! VAX D-Format (D.P.) 2 -128 127
1304 : ! VAX G-Format (D.P.) 2 -1024 1023
1305 : !
1306 : !
1307 : ! XLEAST XSMALL XINF XMAX
1308 : !
1309 : ! CRAY-1 1.84E-2466 3.55E-15 5.45E+2465 5674.858
1310 : ! Cyber 180/855
1311 : ! under NOS (S.P.) 3.14E-294 1.77E-15 1.26E+322 672.789
1312 : ! IEEE (IBM/XT,
1313 : ! SUN, etc.) (S.P.) 1.18E-38 5.95E-8 3.40E+38 85.343
1314 : ! IEEE (IBM/XT,
1315 : ! SUN, etc.) (D.P.) 2.23D-308 1.11D-16 1.79D+308 705.343
1316 : ! IBM 3033 (D.P.) 1.39D-76 1.11D-16 7.23D+75 177.855
1317 : ! VAX D-Format (D.P.) 5.88D-39 6.95D-18 1.70D+38 86.721
1318 : ! VAX G-Format (D.P.) 1.12D-308 5.55D-17 8.98D+307 706.728
1319 : !
1320 : !*******************************************************************
1321 : !*******************************************************************
1322 : !
1323 : ! Error returns
1324 : !
1325 : ! The program returns the value XINF for ARG .LE. 0.0 and the
1326 : ! BESK1 entry returns the value 0.0 for ARG .GT. XMAX.
1327 : !
1328 : !
1329 : ! Intrinsic functions required are:
1330 : !
1331 : ! LOG, SQRT, EXP
1332 : !
1333 : !
1334 : ! Authors: W. J. Cody and Laura Stoltz
1335 : ! Mathematics and Computer Science Division
1336 : ! Argonne National Laboratory
1337 : ! Argonne, IL 60439
1338 : !
1339 : ! Latest modification: January 28, 1988
1340 : ! Taken from http://www.netlib.org/specfun/k1
1341 : !
1342 : !--------------------------------------------------------------------
1343 : !IMPLICIT NONE
1344 : INTEGER :: I,JINT
1345 : !CS REAL
1346 : DOUBLE PRECISION :: &
1347 : & ARG,F,G,ONE,P,PP,Q,QQ,RESULT,SUMF,SUMG,&
1348 : & SUMP,SUMQ,X,XINF,XMAX,XLEAST,XSMALL,XX,ZERO
1349 : DIMENSION P(5),Q(3),PP(11),QQ(9),F(5),G(3)
1350 : !--------------------------------------------------------------------
1351 : ! Mathematical constants
1352 : !--------------------------------------------------------------------
1353 : !CS DATA ONE/1.0E0/,ZERO/0.0E0/
1354 : DATA ONE/1.0D0/,ZERO/0.0D0/
1355 : !--------------------------------------------------------------------
1356 : ! Machine-dependent constants
1357 : !--------------------------------------------------------------------
1358 : !CS DATA XLEAST/1.18E-38/,XSMALL/5.95E-8/,XINF/3.40E+38/,
1359 : !CS 1 XMAX/85.343E+0/
1360 : DATA XLEAST/2.23D-308/,XSMALL/1.11D-16/,XINF/1.79D+308/,&
1361 : & XMAX/705.343D+0/
1362 : !--------------------------------------------------------------------
1363 : ! Coefficients for XLEAST .LE. ARG .LE. 1.0
1364 : !--------------------------------------------------------------------
1365 : !CS DATA P/ 4.8127070456878442310E-1, 9.9991373567429309922E+1,
1366 : !CS 1 7.1885382604084798576E+3, 1.7733324035147015630E+5,
1367 : !CS 2 7.1938920065420586101E+5/
1368 : !CS DATA Q/-2.8143915754538725829E+2, 3.7264298672067697862E+4,
1369 : !CS 1 -2.2149374878243304548E+6/
1370 : !CS DATA F/-2.2795590826955002390E-1,-5.3103913335180275253E+1,
1371 : !CS 1 -4.5051623763436087023E+3,-1.4758069205414222471E+5,
1372 : !CS 2 -1.3531161492785421328E+6/
1373 : !CS DATA G/-3.0507151578787595807E+2, 4.3117653211351080007E+4,
1374 : !CS 2 -2.7062322985570842656E+6/
1375 : DATA P/ 4.8127070456878442310D-1, 9.9991373567429309922D+1,&
1376 : & 7.1885382604084798576D+3, 1.7733324035147015630D+5,&
1377 : & 7.1938920065420586101D+5/
1378 : DATA Q/-2.8143915754538725829D+2, 3.7264298672067697862D+4,&
1379 : & -2.2149374878243304548D+6/
1380 : DATA F/-2.2795590826955002390D-1,-5.3103913335180275253D+1,&
1381 : & -4.5051623763436087023D+3,-1.4758069205414222471D+5,&
1382 : & -1.3531161492785421328D+6/
1383 : DATA G/-3.0507151578787595807D+2, 4.3117653211351080007D+4,&
1384 : & -2.7062322985570842656D+6/
1385 : !--------------------------------------------------------------------
1386 : ! Coefficients for 1.0 .LT. ARG
1387 : !--------------------------------------------------------------------
1388 : !CS DATA PP/ 6.4257745859173138767E-2, 7.5584584631176030810E+0,
1389 : !CS 1 1.3182609918569941308E+2, 8.1094256146537402173E+2,
1390 : !CS 2 2.3123742209168871550E+3, 3.4540675585544584407E+3,
1391 : !CS 3 2.8590657697910288226E+3, 1.3319486433183221990E+3,
1392 : !CS 4 3.4122953486801312910E+2, 4.4137176114230414036E+1,
1393 : !CS 5 2.2196792496874548962E+0/
1394 : !CS DATA QQ/ 3.6001069306861518855E+1, 3.3031020088765390854E+2,
1395 : !CS 1 1.2082692316002348638E+3, 2.1181000487171943810E+3,
1396 : !CS 2 1.9448440788918006154E+3, 9.6929165726802648634E+2,
1397 : !CS 3 2.5951223655579051357E+2, 3.4552228452758912848E+1,
1398 : !CS 4 1.7710478032601086579E+0/
1399 : DATA PP/ 6.4257745859173138767D-2, 7.5584584631176030810D+0,&
1400 : & 1.3182609918569941308D+2, 8.1094256146537402173D+2,&
1401 : & 2.3123742209168871550D+3, 3.4540675585544584407D+3,&
1402 : & 2.8590657697910288226D+3, 1.3319486433183221990D+3,&
1403 : & 3.4122953486801312910D+2, 4.4137176114230414036D+1,&
1404 : & 2.2196792496874548962D+0/
1405 : DATA QQ/ 3.6001069306861518855D+1, 3.3031020088765390854D+2,&
1406 : & 1.2082692316002348638D+3, 2.1181000487171943810D+3,&
1407 : & 1.9448440788918006154D+3, 9.6929165726802648634D+2,&
1408 : & 2.5951223655579051357D+2, 3.4552228452758912848D+1,&
1409 : & 1.7710478032601086579D+0/
1410 : !--------------------------------------------------------------------
1411 0 : X = ARG
1412 0 : IF (X .LT. XLEAST) THEN
1413 : !--------------------------------------------------------------------
1414 : ! Error return for ARG .LT. XLEAST
1415 : !--------------------------------------------------------------------
1416 0 : RESULT = XINF
1417 0 : ELSE IF (X .LE. ONE) THEN
1418 : !--------------------------------------------------------------------
1419 : ! XLEAST .LE. ARG .LE. 1.0
1420 : !--------------------------------------------------------------------
1421 0 : IF (X .LT. XSMALL) THEN
1422 : !--------------------------------------------------------------------
1423 : ! Return for small ARG
1424 : !--------------------------------------------------------------------
1425 0 : RESULT = ONE / X
1426 : ELSE
1427 0 : XX = X * X
1428 : SUMP = ((((P(1)*XX + P(2))*XX + P(3))*XX + P(4))*XX &
1429 0 : & + P(5))*XX + Q(3)
1430 0 : SUMQ = ((XX + Q(1))*XX + Q(2))*XX + Q(3)
1431 : SUMF = (((F(1)*XX + F(2))*XX + F(3))*XX + F(4))*XX &
1432 0 : & + F(5)
1433 0 : SUMG = ((XX + G(1))*XX + G(2))*XX + G(3)
1434 0 : RESULT = (XX * LOG(X) * SUMF/SUMG + SUMP/SUMQ) / X
1435 0 : IF (JINT .EQ. 2) RESULT = RESULT * EXP(X)
1436 : END IF
1437 0 : ELSE IF ((JINT .EQ. 1) .AND. (X .GT. XMAX)) THEN
1438 : !--------------------------------------------------------------------
1439 : ! Error return for ARG .GT. XMAX
1440 : !--------------------------------------------------------------------
1441 0 : RESULT = ZERO
1442 : ELSE
1443 : !--------------------------------------------------------------------
1444 : ! 1.0 .LT. ARG
1445 : !--------------------------------------------------------------------
1446 0 : XX = ONE / X
1447 0 : SUMP = PP(1)
1448 0 : DO 120 I = 2, 11
1449 0 : SUMP = SUMP * XX + PP(I)
1450 0 : 120 CONTINUE
1451 : SUMQ = XX
1452 0 : DO 140 I = 1, 8
1453 0 : SUMQ = (SUMQ + QQ(I)) * XX
1454 0 : 140 CONTINUE
1455 0 : SUMQ = SUMQ + QQ(9)
1456 0 : RESULT = SUMP / SUMQ / SQRT(X)
1457 0 : IF (JINT .EQ. 1) RESULT = RESULT * EXP(-X)
1458 : END IF
1459 0 : RETURN
1460 : !---------- Last line of CALCK1 ----------
1461 : END subroutine calck1
1462 : !!***
1463 :
1464 : !CS REAL
1465 : DOUBLE PRECISION FUNCTION BESK1(X)
1466 : !--------------------------------------------------------------------
1467 : !
1468 : ! This function program computes approximate values for the
1469 : ! modified Bessel function of the second kind of order one
1470 : ! for arguments XLEAST .LE. ARG .LE. XMAX.
1471 : !
1472 : !--------------------------------------------------------------------
1473 : !IMPLICIT NONE
1474 : INTEGER :: JINT
1475 : !CS REAL
1476 : DOUBLE PRECISION :: &
1477 : & X, RESULT
1478 : !--------------------------------------------------------------------
1479 : JINT = 1
1480 : CALL CALCK1(X,RESULT,JINT)
1481 : BESK1 = RESULT
1482 : RETURN
1483 : !---------- Last line of BESK1 ----------
1484 : END function besk1
1485 : !!***
1486 :
1487 : !CS REAL
1488 : DOUBLE PRECISION FUNCTION BESEK1(X)
1489 : !--------------------------------------------------------------------
1490 : !
1491 : ! This function program computes approximate values for the
1492 : ! modified Bessel function of the second kind of order one
1493 : ! multiplied by the exponential function, for arguments
1494 : ! XLEAST .LE. ARG .LE. XMAX.
1495 : !
1496 : !--------------------------------------------------------------------
1497 : !IMPLICIT NONE
1498 : INTEGER JINT
1499 : !CS REAL
1500 : DOUBLE PRECISION :: &
1501 : & X, RESULT
1502 : !--------------------------------------------------------------------
1503 : JINT = 2
1504 : CALL CALCK1(X,RESULT,JINT)
1505 : BESEK1 = RESULT
1506 : RETURN
1507 : !---------- Last line of BESEK1 ----------
1508 : END function besek1
1509 : !!***
1510 :
1511 : end module m_bessel
1512 : !!***
|
