TABLE OF CONTENTS
ABINIT/CALCK1 [ Functions ]
NAME
CALCK1
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 SUBROUTINE CALCK1(ARG,RESULT,JINT) 35 36 !-------------------------------------------------------------------- 37 ! 38 ! This packet computes modified Bessel functions of the second kind 39 ! and order one, K1(X) and EXP(X)*K1(X), for real arguments X. 40 ! It contains two function type subprograms, BESK1 and BESEK1, 41 ! and one subroutine type subprogram, CALCK1. The calling 42 ! statements for the primary entries are 43 ! 44 ! Y=BESK1(X) 45 ! and 46 ! Y=BESEK1(X) 47 ! 48 ! where the entry points correspond to the functions K1(X) and 49 ! EXP(X)*K1(X), respectively. The routine CALCK1 is intended 50 ! for internal packet use only, all computations within the 51 ! packet being concentrated in this routine. The function 52 ! subprograms invoke CALCK1 with the statement 53 ! CALL CALCK1(ARG,RESULT,JINT) 54 ! where the parameter usage is as follows 55 ! 56 ! Function Parameters for CALCK1 57 ! Call ARG RESULT JINT 58 ! 59 ! BESK1(ARG) XLEAST .LT. ARG .LT. XMAX K1(ARG) 1 60 ! BESEK1(ARG) XLEAST .LT. ARG EXP(ARG)*K1(ARG) 2 61 ! 62 ! The main computation evaluates slightly modified forms of near 63 ! minimax rational approximations generated by Russon and Blair, 64 ! Chalk River (Atomic Energy of Canada Limited) Report AECL-3461, 65 ! 1969. This transportable program is patterned after the 66 ! machine-dependent FUNPACK packet NATSK1, but cannot match that 67 ! version for efficiency or accuracy. This version uses rational 68 ! functions that theoretically approximate K-SUB-1(X) to at 69 ! least 18 significant decimal digits. The accuracy achieved 70 ! depends on the arithmetic system, the compiler, the intrinsic 71 ! functions, and proper selection of the machine-dependent 72 ! constants. 73 ! 74 !******************************************************************* 75 !******************************************************************* 76 ! 77 ! Explanation of machine-dependent constants 78 ! 79 ! beta = Radix for the floating-point system 80 ! minexp = Smallest representable power of beta 81 ! maxexp = Smallest power of beta that overflows 82 ! XLEAST = Smallest acceptable argument, i.e., smallest machine 83 ! number X such that 1/X is machine representable. 84 ! XSMALL = Argument below which BESK1(X) and BESEK1(X) may 85 ! each be represented by 1/X. A safe value is the 86 ! largest X such that 1.0 + X = 1.0 to machine 87 ! precision. 88 ! XINF = Largest positive machine number; approximately 89 ! beta**maxexp 90 ! XMAX = Largest argument acceptable to BESK1; Solution to 91 ! equation: 92 ! W(X) * (1+3/8X-15/128X**2) = beta**minexp 93 ! where W(X) = EXP(-X)*SQRT(PI/2X) 94 ! 95 ! 96 ! Approximate values for some important machines are: 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 ! XLEAST XSMALL XINF XMAX 113 ! 114 ! CRAY-1 1.84E-2466 3.55E-15 5.45E+2465 5674.858 115 ! Cyber 180/855 116 ! under NOS (S.P.) 3.14E-294 1.77E-15 1.26E+322 672.789 117 ! IEEE (IBM/XT, 118 ! SUN, etc.) (S.P.) 1.18E-38 5.95E-8 3.40E+38 85.343 119 ! IEEE (IBM/XT, 120 ! SUN, etc.) (D.P.) 2.23D-308 1.11D-16 1.79D+308 705.343 121 ! IBM 3033 (D.P.) 1.39D-76 1.11D-16 7.23D+75 177.855 122 ! VAX D-Format (D.P.) 5.88D-39 6.95D-18 1.70D+38 86.721 123 ! VAX G-Format (D.P.) 1.12D-308 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 ! BESK1 entry returns the value 0.0 for ARG .GT. XMAX. 132 ! 133 ! 134 ! Intrinsic functions required are: 135 ! 136 ! LOG, SQRT, EXP 137 ! 138 ! 139 ! Authors: W. J. Cody and Laura Stoltz 140 ! Mathematics and Computer Science Division 141 ! Argonne National Laboratory 142 ! Argonne, IL 60439 143 ! 144 ! Latest modification: January 28, 1988 145 ! Taken from http://www.netlib.org/specfun/k1 146 ! 147 !-------------------------------------------------------------------- 148 149 !This section has been created automatically by the script Abilint (TD). 150 !Do not modify the following lines by hand. 151 #undef ABI_FUNC 152 #define ABI_FUNC 'CALCK1' 153 !End of the abilint section 154 155 IMPLICIT NONE 156 INTEGER :: I,JINT 157 !CS REAL 158 DOUBLE PRECISION :: & 159 & ARG,F,G,ONE,P,PP,Q,QQ,RESULT,SUMF,SUMG,& 160 & SUMP,SUMQ,X,XINF,XMAX,XLEAST,XSMALL,XX,ZERO 161 DIMENSION P(5),Q(3),PP(11),QQ(9),F(5),G(3) 162 !-------------------------------------------------------------------- 163 ! Mathematical constants 164 !-------------------------------------------------------------------- 165 !CS DATA ONE/1.0E0/,ZERO/0.0E0/ 166 DATA ONE/1.0D0/,ZERO/0.0D0/ 167 !-------------------------------------------------------------------- 168 ! Machine-dependent constants 169 !-------------------------------------------------------------------- 170 !CS DATA XLEAST/1.18E-38/,XSMALL/5.95E-8/,XINF/3.40E+38/, 171 !CS 1 XMAX/85.343E+0/ 172 DATA XLEAST/2.23D-308/,XSMALL/1.11D-16/,XINF/1.79D+308/,& 173 & XMAX/705.343D+0/ 174 !-------------------------------------------------------------------- 175 ! Coefficients for XLEAST .LE. ARG .LE. 1.0 176 !-------------------------------------------------------------------- 177 !CS DATA P/ 4.8127070456878442310E-1, 9.9991373567429309922E+1, 178 !CS 1 7.1885382604084798576E+3, 1.7733324035147015630E+5, 179 !CS 2 7.1938920065420586101E+5/ 180 !CS DATA Q/-2.8143915754538725829E+2, 3.7264298672067697862E+4, 181 !CS 1 -2.2149374878243304548E+6/ 182 !CS DATA F/-2.2795590826955002390E-1,-5.3103913335180275253E+1, 183 !CS 1 -4.5051623763436087023E+3,-1.4758069205414222471E+5, 184 !CS 2 -1.3531161492785421328E+6/ 185 !CS DATA G/-3.0507151578787595807E+2, 4.3117653211351080007E+4, 186 !CS 2 -2.7062322985570842656E+6/ 187 DATA P/ 4.8127070456878442310D-1, 9.9991373567429309922D+1,& 188 & 7.1885382604084798576D+3, 1.7733324035147015630D+5,& 189 & 7.1938920065420586101D+5/ 190 DATA Q/-2.8143915754538725829D+2, 3.7264298672067697862D+4,& 191 & -2.2149374878243304548D+6/ 192 DATA F/-2.2795590826955002390D-1,-5.3103913335180275253D+1,& 193 & -4.5051623763436087023D+3,-1.4758069205414222471D+5,& 194 & -1.3531161492785421328D+6/ 195 DATA G/-3.0507151578787595807D+2, 4.3117653211351080007D+4,& 196 & -2.7062322985570842656D+6/ 197 !-------------------------------------------------------------------- 198 ! Coefficients for 1.0 .LT. ARG 199 !-------------------------------------------------------------------- 200 !CS DATA PP/ 6.4257745859173138767E-2, 7.5584584631176030810E+0, 201 !CS 1 1.3182609918569941308E+2, 8.1094256146537402173E+2, 202 !CS 2 2.3123742209168871550E+3, 3.4540675585544584407E+3, 203 !CS 3 2.8590657697910288226E+3, 1.3319486433183221990E+3, 204 !CS 4 3.4122953486801312910E+2, 4.4137176114230414036E+1, 205 !CS 5 2.2196792496874548962E+0/ 206 !CS DATA QQ/ 3.6001069306861518855E+1, 3.3031020088765390854E+2, 207 !CS 1 1.2082692316002348638E+3, 2.1181000487171943810E+3, 208 !CS 2 1.9448440788918006154E+3, 9.6929165726802648634E+2, 209 !CS 3 2.5951223655579051357E+2, 3.4552228452758912848E+1, 210 !CS 4 1.7710478032601086579E+0/ 211 DATA PP/ 6.4257745859173138767D-2, 7.5584584631176030810D+0,& 212 & 1.3182609918569941308D+2, 8.1094256146537402173D+2,& 213 & 2.3123742209168871550D+3, 3.4540675585544584407D+3,& 214 & 2.8590657697910288226D+3, 1.3319486433183221990D+3,& 215 & 3.4122953486801312910D+2, 4.4137176114230414036D+1,& 216 & 2.2196792496874548962D+0/ 217 DATA QQ/ 3.6001069306861518855D+1, 3.3031020088765390854D+2,& 218 & 1.2082692316002348638D+3, 2.1181000487171943810D+3,& 219 & 1.9448440788918006154D+3, 9.6929165726802648634D+2,& 220 & 2.5951223655579051357D+2, 3.4552228452758912848D+1,& 221 & 1.7710478032601086579D+0/ 222 !-------------------------------------------------------------------- 223 X = ARG 224 IF (X .LT. XLEAST) THEN 225 !-------------------------------------------------------------------- 226 ! Error return for ARG .LT. XLEAST 227 !-------------------------------------------------------------------- 228 RESULT = XINF 229 ELSE IF (X .LE. ONE) THEN 230 !-------------------------------------------------------------------- 231 ! XLEAST .LE. ARG .LE. 1.0 232 !-------------------------------------------------------------------- 233 IF (X .LT. XSMALL) THEN 234 !-------------------------------------------------------------------- 235 ! Return for small ARG 236 !-------------------------------------------------------------------- 237 RESULT = ONE / X 238 ELSE 239 XX = X * X 240 SUMP = ((((P(1)*XX + P(2))*XX + P(3))*XX + P(4))*XX & 241 & + P(5))*XX + Q(3) 242 SUMQ = ((XX + Q(1))*XX + Q(2))*XX + Q(3) 243 SUMF = (((F(1)*XX + F(2))*XX + F(3))*XX + F(4))*XX & 244 & + F(5) 245 SUMG = ((XX + G(1))*XX + G(2))*XX + G(3) 246 RESULT = (XX * LOG(X) * SUMF/SUMG + SUMP/SUMQ) / X 247 IF (JINT .EQ. 2) RESULT = RESULT * EXP(X) 248 END IF 249 ELSE IF ((JINT .EQ. 1) .AND. (X .GT. XMAX)) THEN 250 !-------------------------------------------------------------------- 251 ! Error return for ARG .GT. XMAX 252 !-------------------------------------------------------------------- 253 RESULT = ZERO 254 ELSE 255 !-------------------------------------------------------------------- 256 ! 1.0 .LT. ARG 257 !-------------------------------------------------------------------- 258 XX = ONE / X 259 SUMP = PP(1) 260 DO 120 I = 2, 11 261 SUMP = SUMP * XX + PP(I) 262 120 CONTINUE 263 SUMQ = XX 264 DO 140 I = 1, 8 265 SUMQ = (SUMQ + QQ(I)) * XX 266 140 CONTINUE 267 SUMQ = SUMQ + QQ(9) 268 RESULT = SUMP / SUMQ / SQRT(X) 269 IF (JINT .EQ. 1) RESULT = RESULT * EXP(-X) 270 END IF 271 RETURN 272 !---------- Last line of CALCK1 ---------- 273 END 274 !CS REAL 275 DOUBLE PRECISION FUNCTION BESK1(X) 276 !-------------------------------------------------------------------- 277 ! 278 ! This function program computes approximate values for the 279 ! modified Bessel function of the second kind of order one 280 ! for arguments XLEAST .LE. ARG .LE. XMAX. 281 ! 282 !-------------------------------------------------------------------- 283 284 !This section has been created automatically by the script Abilint (TD). 285 !Do not modify the following lines by hand. 286 #undef ABI_FUNC 287 #define ABI_FUNC 'BESK1' 288 use interfaces_28_numeric_noabirule, except_this_one => BESK1 289 !End of the abilint section 290 291 IMPLICIT NONE 292 INTEGER :: JINT 293 !CS REAL 294 DOUBLE PRECISION :: & 295 & X, RESULT 296 !-------------------------------------------------------------------- 297 JINT = 1 298 CALL CALCK1(X,RESULT,JINT) 299 BESK1 = RESULT 300 RETURN 301 !---------- Last line of BESK1 ---------- 302 END 303 !CS REAL 304 DOUBLE PRECISION FUNCTION BESEK1(X) 305 !-------------------------------------------------------------------- 306 ! 307 ! This function program computes approximate values for the 308 ! modified Bessel function of the second kind of order one 309 ! multiplied by the exponential function, for arguments 310 ! XLEAST .LE. ARG .LE. XMAX. 311 ! 312 !-------------------------------------------------------------------- 313 314 !This section has been created automatically by the script Abilint (TD). 315 !Do not modify the following lines by hand. 316 #undef ABI_FUNC 317 #define ABI_FUNC 'BESEK1' 318 use interfaces_28_numeric_noabirule, except_this_one => BESEK1 319 !End of the abilint section 320 321 IMPLICIT NONE 322 INTEGER JINT 323 !CS REAL 324 DOUBLE PRECISION :: & 325 & X, RESULT 326 !-------------------------------------------------------------------- 327 JINT = 2 328 CALL CALCK1(X,RESULT,JINT) 329 BESEK1 = RESULT 330 RETURN 331 !---------- Last line of BESEK1 ---------- 332 END