TABLE OF CONTENTS
ABINIT/CALCK0 [ Functions ]
NAME
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
922 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 X = ARG 1102 IF (X .GT. ZERO) THEN 1103 IF (X .LE. ONE) THEN 1104 !-------------------------------------------------------------------- 1105 ! 0.0 .LT. ARG .LE. 1.0 1106 !-------------------------------------------------------------------- 1107 TEMP = LOG(X) 1108 IF (X .LT. XSMALL) THEN 1109 !-------------------------------------------------------------------- 1110 ! Return for small ARG 1111 !-------------------------------------------------------------------- 1112 RESULT = P(6)/Q(2) - TEMP 1113 ELSE 1114 XX = X * X 1115 SUMP = ((((P(1)*XX + P(2))*XX + P(3))*XX +& 1116 P(4))*XX + P(5))*XX + P(6) 1117 SUMQ = (XX + Q(1))*XX + Q(2) 1118 SUMF = ((F(1)*XX + F(2))*XX + F(3))*XX + F(4) 1119 SUMG = ((XX + G(1))*XX + G(2))*XX + G(3) 1120 RESULT = SUMP/SUMQ - XX*SUMF*TEMP/SUMG - TEMP 1121 IF (JINT .EQ. 2) RESULT = RESULT * EXP(X) 1122 END IF 1123 ELSE IF ((JINT .EQ. 1) .AND. (X .GT. XMAX)) THEN 1124 !-------------------------------------------------------------------- 1125 ! Error return for ARG .GT. XMAX 1126 !-------------------------------------------------------------------- 1127 RESULT = ZERO 1128 ELSE 1129 !-------------------------------------------------------------------- 1130 ! 1.0 .LT. ARG 1131 !-------------------------------------------------------------------- 1132 XX = ONE / X 1133 SUMP = PP(1) 1134 DO 120 I = 2, 10 1135 SUMP = SUMP*XX + PP(I) 1136 120 CONTINUE 1137 SUMQ = XX 1138 DO 140 I = 1, 9 1139 SUMQ = (SUMQ + QQ(I))*XX 1140 140 CONTINUE 1141 SUMQ = SUMQ + QQ(10) 1142 RESULT = SUMP / SUMQ / SQRT(X) 1143 IF (JINT .EQ. 1) RESULT = RESULT * EXP(-X) 1144 END IF 1145 ELSE 1146 !-------------------------------------------------------------------- 1147 ! Error return for ARG .LE. 0.0 1148 !-------------------------------------------------------------------- 1149 RESULT = XINF 1150 END IF 1151 !-------------------------------------------------------------------- 1152 ! Update error counts, etc. 1153 !-------------------------------------------------------------------- 1154 RETURN 1155 !---------- Last line of CALCK0 ---------- 1156 END subroutine calck0
ABINIT/CALCK1 [ Functions ]
NAME
CALCK1
FUNCTION
INPUTS
OUTPUT
SIDE EFFECTS
NOTES
SOURCE
1229 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 X = ARG 1412 IF (X .LT. XLEAST) THEN 1413 !-------------------------------------------------------------------- 1414 ! Error return for ARG .LT. XLEAST 1415 !-------------------------------------------------------------------- 1416 RESULT = XINF 1417 ELSE IF (X .LE. ONE) THEN 1418 !-------------------------------------------------------------------- 1419 ! XLEAST .LE. ARG .LE. 1.0 1420 !-------------------------------------------------------------------- 1421 IF (X .LT. XSMALL) THEN 1422 !-------------------------------------------------------------------- 1423 ! Return for small ARG 1424 !-------------------------------------------------------------------- 1425 RESULT = ONE / X 1426 ELSE 1427 XX = X * X 1428 SUMP = ((((P(1)*XX + P(2))*XX + P(3))*XX + P(4))*XX & 1429 & + P(5))*XX + Q(3) 1430 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 & + F(5) 1433 SUMG = ((XX + G(1))*XX + G(2))*XX + G(3) 1434 RESULT = (XX * LOG(X) * SUMF/SUMG + SUMP/SUMQ) / X 1435 IF (JINT .EQ. 2) RESULT = RESULT * EXP(X) 1436 END IF 1437 ELSE IF ((JINT .EQ. 1) .AND. (X .GT. XMAX)) THEN 1438 !-------------------------------------------------------------------- 1439 ! Error return for ARG .GT. XMAX 1440 !-------------------------------------------------------------------- 1441 RESULT = ZERO 1442 ELSE 1443 !-------------------------------------------------------------------- 1444 ! 1.0 .LT. ARG 1445 !-------------------------------------------------------------------- 1446 XX = ONE / X 1447 SUMP = PP(1) 1448 DO 120 I = 2, 11 1449 SUMP = SUMP * XX + PP(I) 1450 120 CONTINUE 1451 SUMQ = XX 1452 DO 140 I = 1, 8 1453 SUMQ = (SUMQ + QQ(I)) * XX 1454 140 CONTINUE 1455 SUMQ = SUMQ + QQ(9) 1456 RESULT = SUMP / SUMQ / SQRT(X) 1457 IF (JINT .EQ. 1) RESULT = RESULT * EXP(-X) 1458 END IF 1459 RETURN 1460 !---------- Last line of CALCK1 ---------- 1461 END subroutine calck1
ABINIT/m_bessel [ Modules ]
NAME
m_bessel
FUNCTION
Bessel functions
COPYRIGHT
Copyright (C) 2008-2024 ABINIT group (MG) This file is distributed under the terms of the GNU General Public License, see ~abinit/COPYING or http://www.gnu.org/copyleft/gpl.txt .
SOURCE
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