TABLE OF CONTENTS
ABINIT/CALCK0 [ Functions ]
NAME
FUNCTION
COPYRIGHT
Copyright (C) 2014-2017 ABINIT group (XG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
PARENTS
CHILDREN
SOURCE
28 #if defined HAVE_CONFIG_H 29 #include "config.h" 30 #endif 31 32 #include "abi_common.h" 33 34 35 SUBROUTINE CALCK0(ARG,RESULT,JINT) 36 37 !-------------------------------------------------------------------- 38 ! 39 ! This packet computes modified Bessel functions of the second kind 40 ! and order zero, K0(X) and EXP(X)*K0(X), for real 41 ! arguments X. It contains two function type subprograms, BESK0 42 ! and BESEK0, and one subroutine type subprogram, CALCK0. 43 ! the calling statements for the primary entries are 44 ! 45 ! Y=BESK0(X) 46 ! and 47 ! Y=BESEK0(X) 48 ! 49 ! where the entry points correspond to the functions K0(X) and 50 ! EXP(X)*K0(X), respectively. The routine CALCK0 is 51 ! intended for internal packet use only, all computations within 52 ! the packet being concentrated in this routine. The function 53 ! subprograms invoke CALCK0 with the statement 54 ! CALL CALCK0(ARG,RESULT,JINT) 55 ! where the parameter usage is as follows 56 ! 57 ! Function Parameters for CALCK0 58 ! Call ARG RESULT JINT 59 ! 60 ! BESK0(ARG) 0 .LT. ARG .LE. XMAX K0(ARG) 1 61 ! BESEK0(ARG) 0 .LT. ARG EXP(ARG)*K0(ARG) 2 62 ! 63 ! The main computation evaluates slightly modified forms of near 64 ! minimax rational approximations generated by Russon and Blair, 65 ! Chalk River (Atomic Energy of Canada Limited) Report AECL-3461, 66 ! 1969. This transportable program is patterned after the 67 ! machine-dependent FUNPACK packet NATSK0, but cannot match that 68 ! version for efficiency or accuracy. This version uses rational 69 ! functions that theoretically approximate K-SUB-0(X) to at 70 ! least 18 significant decimal digits. The accuracy achieved 71 ! depends on the arithmetic system, the compiler, the intrinsic 72 ! functions, and proper selection of the machine-dependent 73 ! constants. 74 ! 75 !******************************************************************* 76 !******************************************************************* 77 ! 78 ! Explanation of machine-dependent constants 79 ! 80 ! beta = Radix for the floating-point system 81 ! minexp = Smallest representable power of beta 82 ! maxexp = Smallest power of beta that overflows 83 ! XSMALL = Argument below which BESK0 and BESEK0 may 84 ! each be represented by a constant and a log. 85 ! largest X such that 1.0 + X = 1.0 to machine 86 ! precision. 87 ! XINF = Largest positive machine number; approximately 88 ! beta**maxexp 89 ! XMAX = Largest argument acceptable to BESK0; Solution to 90 ! equation: 91 ! W(X) * (1-1/8X+9/128X**2) = beta**minexp 92 ! where W(X) = EXP(-X)*SQRT(PI/2X) 93 ! 94 ! 95 ! Approximate values for some important machines are: 96 ! 97 ! 98 ! beta minexp maxexp 99 ! 100 ! CRAY-1 (S.P.) 2 -8193 8191 101 ! Cyber 180/185 102 ! under NOS (S.P.) 2 -975 1070 103 ! IEEE (IBM/XT, 104 ! SUN, etc.) (S.P.) 2 -126 128 105 ! IEEE (IBM/XT, 106 ! SUN, etc.) (D.P.) 2 -1022 1024 107 ! IBM 3033 (D.P.) 16 -65 63 108 ! VAX D-Format (D.P.) 2 -128 127 109 ! VAX G-Format (D.P.) 2 -1024 1023 110 ! 111 ! 112 ! XSMALL XINF XMAX 113 ! 114 ! CRAY-1 (S.P.) 3.55E-15 5.45E+2465 5674.858 115 ! Cyber 180/855 116 ! under NOS (S.P.) 1.77E-15 1.26E+322 672.788 117 ! IEEE (IBM/XT, 118 ! SUN, etc.) (S.P.) 5.95E-8 3.40E+38 85.337 119 ! IEEE (IBM/XT, 120 ! SUN, etc.) (D.P.) 1.11D-16 1.79D+308 705.342 121 ! IBM 3033 (D.P.) 1.11D-16 7.23D+75 177.852 122 ! VAX D-Format (D.P.) 6.95D-18 1.70D+38 86.715 123 ! VAX G-Format (D.P.) 5.55D-17 8.98D+307 706.728 124 ! 125 !******************************************************************* 126 !******************************************************************* 127 ! 128 ! Error returns 129 ! 130 ! The program returns the value XINF for ARG .LE. 0.0, and the 131 ! BESK0 entry returns the value 0.0 for ARG .GT. XMAX. 132 ! 133 ! 134 ! Intrinsic functions required are: 135 ! 136 ! EXP, LOG, SQRT 137 ! 138 ! Latest modification: March 19, 1990 139 ! 140 ! Authors: W. J. Cody and Laura Stoltz 141 ! Mathematics and Computer Science Division 142 ! Argonne National Laboratory 143 ! Argonne, IL 60439 144 ! 145 ! Original subroutine from netlib http://www.netlib.org/specfun/k0 146 ! Slightly modified by MG to follow f90 rules and double precision arithmetic 147 ! 148 !-------------------------------------------------------------------- 149 150 !This section has been created automatically by the script Abilint (TD). 151 !Do not modify the following lines by hand. 152 #undef ABI_FUNC 153 #define ABI_FUNC 'CALCK0' 154 !End of the abilint section 155 156 IMPLICIT NONE 157 INTEGER :: I,JINT 158 !CS REAL 159 DOUBLE PRECISION :: ARG,RESULT,SUMF,SUMG,SUMP,SUMQ,TEMP 160 !CS REAL 161 DOUBLE PRECISION :: X,XX 162 !CS REAL 163 DOUBLE PRECISION :: P(6),Q(2),PP(10),QQ(10),F(4),G(3) 164 !C-------------------------------------------------------------------- 165 !C Mathematical constants 166 !C-------------------------------------------------------------------- 167 !CS REAL, PARAMETER :: ONE=1.0E0,ZERO=0.0E0 168 DOUBLE PRECISION,PARAMETER :: ONE=1.0D0,ZERO=0.0D0 169 !C-------------------------------------------------------------------- 170 !C Machine-dependent constants 171 !C-------------------------------------------------------------------- 172 !CS REAL.PARAMETER :: XSMALL=5.95E-8, XINF=3.40E+38 ,XMAX=85.337E0 173 DOUBLE PRECISION,PARAMETER :: XSMALL=1.11D-16,XINF=1.79D+308,XMAX=705.342D0 174 !-------------------------------------------------------------------- 175 ! 176 ! Coefficients for XSMALL .LE. ARG .LE. 1.0 177 ! 178 !-------------------------------------------------------------------- 179 !S DATA P/ 5.8599221412826100000E-04, 1.3166052564989571850E-01, 180 !S 1 1.1999463724910714109E+01, 4.6850901201934832188E+02, 181 !S 2 5.9169059852270512312E+03, 2.4708152720399552679E+03/ 182 !S DATA Q/-2.4994418972832303646E+02, 2.1312714303849120380E+04/ 183 !S DATA F/-1.6414452837299064100E+00,-2.9601657892958843866E+02, 184 !S 1 -1.7733784684952985886E+04,-4.0320340761145482298E+05/ 185 !S DATA G/-2.5064972445877992730E+02, 2.9865713163054025489E+04, 186 !S 1 -1.6128136304458193998E+06/ 187 DATA P/5.8599221412826100000D-04,1.3166052564989571850D-01,& 188 & 1.1999463724910714109D+01,4.6850901201934832188D+02,& 189 & 5.9169059852270512312D+03,2.4708152720399552679D+03/ 190 DATA Q/-2.4994418972832303646D+02, 2.1312714303849120380D+04/ 191 DATA F/-1.6414452837299064100D+00,-2.9601657892958843866D+02,& 192 & -1.7733784684952985886D+04,-4.0320340761145482298D+05/ 193 DATA G/-2.5064972445877992730D+02, 2.9865713163054025489D+04,& 194 & -1.6128136304458193998D+06/ 195 !-------------------------------------------------------------------- 196 ! 197 ! Coefficients for 1.0 .LT. ARG 198 ! 199 !-------------------------------------------------------------------- 200 !S DATA PP/ 1.1394980557384778174E+02, 3.6832589957340267940E+03, 201 !S 1 3.1075408980684392399E+04, 1.0577068948034021957E+05, 202 !S 2 1.7398867902565686251E+05, 1.5097646353289914539E+05, 203 !S 3 7.1557062783764037541E+04, 1.8321525870183537725E+04, 204 !S 4 2.3444738764199315021E+03, 1.1600249425076035558E+02/ 205 !S DATA QQ/ 2.0013443064949242491E+02, 4.4329628889746408858E+03, 206 !S 1 3.1474655750295278825E+04, 9.7418829762268075784E+04, 207 !S 2 1.5144644673520157801E+05, 1.2689839587977598727E+05, 208 !S 3 5.8824616785857027752E+04, 1.4847228371802360957E+04, 209 !S 4 1.8821890840982713696E+03, 9.2556599177304839811E+01/ 210 DATA PP/ 1.1394980557384778174D+02, 3.6832589957340267940D+03,& 211 & 3.1075408980684392399D+04, 1.0577068948034021957D+05,& 212 & 1.7398867902565686251D+05, 1.5097646353289914539D+05,& 213 & 7.1557062783764037541D+04, 1.8321525870183537725D+04,& 214 & 2.3444738764199315021D+03, 1.1600249425076035558D+02/ 215 DATA QQ/ 2.0013443064949242491D+02, 4.4329628889746408858D+03, & 216 & 3.1474655750295278825D+04, 9.7418829762268075784D+04,& 217 & 1.5144644673520157801D+05, 1.2689839587977598727D+05,& 218 & 5.8824616785857027752D+04, 1.4847228371802360957D+04,& 219 & 1.8821890840982713696D+03, 9.2556599177304839811D+01/ 220 !-------------------------------------------------------------------- 221 X = ARG 222 IF (X .GT. ZERO) THEN 223 IF (X .LE. ONE) THEN 224 !-------------------------------------------------------------------- 225 ! 0.0 .LT. ARG .LE. 1.0 226 !-------------------------------------------------------------------- 227 TEMP = LOG(X) 228 IF (X .LT. XSMALL) THEN 229 !-------------------------------------------------------------------- 230 ! Return for small ARG 231 !-------------------------------------------------------------------- 232 RESULT = P(6)/Q(2) - TEMP 233 ELSE 234 XX = X * X 235 SUMP = ((((P(1)*XX + P(2))*XX + P(3))*XX +& 236 P(4))*XX + P(5))*XX + P(6) 237 SUMQ = (XX + Q(1))*XX + Q(2) 238 SUMF = ((F(1)*XX + F(2))*XX + F(3))*XX + F(4) 239 SUMG = ((XX + G(1))*XX + G(2))*XX + G(3) 240 RESULT = SUMP/SUMQ - XX*SUMF*TEMP/SUMG - TEMP 241 IF (JINT .EQ. 2) RESULT = RESULT * EXP(X) 242 END IF 243 ELSE IF ((JINT .EQ. 1) .AND. (X .GT. XMAX)) THEN 244 !-------------------------------------------------------------------- 245 ! Error return for ARG .GT. XMAX 246 !-------------------------------------------------------------------- 247 RESULT = ZERO 248 ELSE 249 !-------------------------------------------------------------------- 250 ! 1.0 .LT. ARG 251 !-------------------------------------------------------------------- 252 XX = ONE / X 253 SUMP = PP(1) 254 DO 120 I = 2, 10 255 SUMP = SUMP*XX + PP(I) 256 120 CONTINUE 257 SUMQ = XX 258 DO 140 I = 1, 9 259 SUMQ = (SUMQ + QQ(I))*XX 260 140 CONTINUE 261 SUMQ = SUMQ + QQ(10) 262 RESULT = SUMP / SUMQ / SQRT(X) 263 IF (JINT .EQ. 1) RESULT = RESULT * EXP(-X) 264 END IF 265 ELSE 266 !-------------------------------------------------------------------- 267 ! Error return for ARG .LE. 0.0 268 !-------------------------------------------------------------------- 269 RESULT = XINF 270 END IF 271 !-------------------------------------------------------------------- 272 ! Update error counts, etc. 273 !-------------------------------------------------------------------- 274 RETURN 275 !---------- Last line of CALCK0 ---------- 276 END 277 278 !S REAL 279 DOUBLE PRECISION FUNCTION BESK0(X) 280 !-------------------------------------------------------------------- 281 ! 282 ! This function program computes approximate values for the 283 ! modified Bessel function of the second kind of order zero 284 ! for arguments 0.0 .LT. ARG .LE. XMAX (see comments heading 285 ! CALCK0). 286 ! 287 ! Authors: W. J. Cody and Laura Stoltz 288 ! 289 ! Latest Modification: January 19, 1988 290 ! 291 !-------------------------------------------------------------------- 292 293 !This section has been created automatically by the script Abilint (TD). 294 !Do not modify the following lines by hand. 295 #undef ABI_FUNC 296 #define ABI_FUNC 'BESK0' 297 use interfaces_28_numeric_noabirule, except_this_one => BESK0 298 !End of the abilint section 299 300 IMPLICIT NONE 301 INTEGER :: JINT 302 !S REAL 303 DOUBLE PRECISION :: X, RESULT 304 !-------------------------------------------------------------------- 305 JINT = 1 306 CALL CALCK0(X,RESULT,JINT) 307 BESK0 = RESULT 308 RETURN 309 !---------- Last line of BESK0 ---------- 310 END 311 !S REAL 312 DOUBLE PRECISION FUNCTION BESEK0(X) 313 !-------------------------------------------------------------------- 314 ! 315 ! This function program computes approximate values for the 316 ! modified Bessel function of the second kind of order zero 317 ! multiplied by the Exponential function, for arguments 318 ! 0.0 .LT. ARG. 319 ! 320 ! Authors: W. J. Cody and Laura Stoltz 321 ! 322 ! Latest Modification: January 19, 1988 323 ! 324 !-------------------------------------------------------------------- 325 326 !This section has been created automatically by the script Abilint (TD). 327 !Do not modify the following lines by hand. 328 #undef ABI_FUNC 329 #define ABI_FUNC 'BESEK0' 330 use interfaces_28_numeric_noabirule, except_this_one => BESEK0 331 !End of the abilint section 332 333 IMPLICIT NONE 334 INTEGER JINT 335 !S REAL 336 DOUBLE PRECISION :: X,RESULT 337 !-------------------------------------------------------------------- 338 JINT = 2 339 CALL CALCK0(X,RESULT,JINT) 340 BESEK0 = RESULT 341 RETURN 342 !---------- Last line of BESEK0 ---------- 343 END