Branch data 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 : : !!***
|
