TABLE OF CONTENTS


ABINIT/CALCK0 [ Functions ]

[ Top ] [ Functions ]

NAME

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

CHILDREN

SOURCE

 969 SUBROUTINE CALCK0(ARG,RESULT,JINT)
 970 
 971 !--------------------------------------------------------------------
 972 !
 973 ! This packet computes modified Bessel functions of the second kind
 974 !   and order zero, K0(X) and EXP(X)*K0(X), for real
 975 !   arguments X.  It contains two function type subprograms, BESK0
 976 !   and BESEK0, and one subroutine type subprogram, CALCK0.
 977 !   the calling statements for the primary entries are
 978 !
 979 !                   Y=BESK0(X)
 980 !   and
 981 !                   Y=BESEK0(X)
 982 !
 983 !   where the entry points correspond to the functions K0(X) and
 984 !   EXP(X)*K0(X), respectively.  The routine CALCK0 is
 985 !   intended for internal packet use only, all computations within
 986 !   the packet being concentrated in this routine.  The function
 987 !   subprograms invoke CALCK0 with the statement
 988 !          CALL CALCK0(ARG,RESULT,JINT)
 989 !   where the parameter usage is as follows
 990 !
 991 !      Function                     Parameters for CALCK0
 992 !       Call              ARG                  RESULT          JINT
 993 !
 994 !     BESK0(ARG)   0 .LT. ARG .LE. XMAX       K0(ARG)           1
 995 !     BESEK0(ARG)     0 .LT. ARG           EXP(ARG)*K0(ARG)     2
 996 !
 997 !   The main computation evaluates slightly modified forms of near
 998 !   minimax rational approximations generated by Russon and Blair,
 999 !   Chalk River (Atomic Energy of Canada Limited) Report AECL-3461,
1000 !   1969.  This transportable program is patterned after the
1001 !   machine-dependent FUNPACK packet NATSK0, but cannot match that
1002 !   version for efficiency or accuracy.  This version uses rational
1003 !   functions that theoretically approximate K-SUB-0(X) to at
1004 !   least 18 significant decimal digits.  The accuracy achieved
1005 !   depends on the arithmetic system, the compiler, the intrinsic
1006 !   functions, and proper selection of the machine-dependent
1007 !   constants.
1008 !
1009 !*******************************************************************
1010 !*******************************************************************
1011 !
1012 ! Explanation of machine-dependent constants
1013 !
1014 !   beta   = Radix for the floating-point system
1015 !   minexp = Smallest representable power of beta
1016 !   maxexp = Smallest power of beta that overflows
1017 !   XSMALL = Argument below which BESK0 and BESEK0 may
1018 !            each be represented by a constant and a log.
1019 !            largest X such that  1.0 + X = 1.0  to machine
1020 !            precision.
1021 !   XINF   = Largest positive machine number; approximately
1022 !            beta**maxexp
1023 !   XMAX   = Largest argument acceptable to BESK0;  Solution to
1024 !            equation:
1025 !               W(X) * (1-1/8X+9/128X**2) = beta**minexp
1026 !            where  W(X) = EXP(-X)*SQRT(PI/2X)
1027 !
1028 !
1029 !     Approximate values for some important machines are:
1030 !
1031 !
1032 !                           beta       minexp       maxexp
1033 !
1034 !  CRAY-1        (S.P.)       2        -8193         8191
1035 !  Cyber 180/185
1036 !    under NOS   (S.P.)       2         -975         1070
1037 !  IEEE (IBM/XT,
1038 !    SUN, etc.)  (S.P.)       2         -126          128
1039 !  IEEE (IBM/XT,
1040 !    SUN, etc.)  (D.P.)       2        -1022         1024
1041 !  IBM 3033      (D.P.)      16          -65           63
1042 !  VAX D-Format  (D.P.)       2         -128          127
1043 !  VAX G-Format  (D.P.)       2        -1024         1023
1044 !
1045 !
1046 !                          XSMALL       XINF         XMAX
1047 !
1048 ! CRAY-1        (S.P.)    3.55E-15   5.45E+2465    5674.858
1049 ! Cyber 180/855
1050 !   under NOS   (S.P.)    1.77E-15   1.26E+322      672.788
1051 ! IEEE (IBM/XT,
1052 !   SUN, etc.)  (S.P.)    5.95E-8    3.40E+38        85.337
1053 ! IEEE (IBM/XT,
1054 !   SUN, etc.)  (D.P.)    1.11D-16   1.79D+308      705.342
1055 ! IBM 3033      (D.P.)    1.11D-16   7.23D+75       177.852
1056 ! VAX D-Format  (D.P.)    6.95D-18   1.70D+38        86.715
1057 ! VAX G-Format  (D.P.)    5.55D-17   8.98D+307      706.728
1058 !
1059 !*******************************************************************
1060 !*******************************************************************
1061 !
1062 ! Error returns
1063 !
1064 !  The program returns the value XINF for ARG .LE. 0.0, and the
1065 !  BESK0 entry returns the value 0.0 for ARG .GT. XMAX.
1066 !
1067 !
1068 !  Intrinsic functions required are:
1069 !
1070 !     EXP, LOG, SQRT
1071 !
1072 !  Latest modification: March 19, 1990
1073 !
1074 !  Authors: W. J. Cody and Laura Stoltz
1075 !           Mathematics and Computer Science Division
1076 !           Argonne National Laboratory
1077 !           Argonne, IL 60439
1078 !
1079 !  Original subroutine from netlib http://www.netlib.org/specfun/k0
1080 !  Slightly modified by MG to follow f90 rules and double precision arithmetic
1081 !
1082 !--------------------------------------------------------------------
1083 
1084 !This section has been created automatically by the script Abilint (TD).
1085 !Do not modify the following lines by hand.
1086 #undef ABI_FUNC
1087 #define ABI_FUNC 'CALCK0'
1088 !End of the abilint section
1089 
1090       IMPLICIT NONE
1091       INTEGER :: I,JINT
1092 !CS    REAL
1093       DOUBLE PRECISION :: ARG,RESULT,SUMF,SUMG,SUMP,SUMQ,TEMP
1094 !CS    REAL
1095       DOUBLE PRECISION :: X,XX
1096 !CS    REAL
1097       DOUBLE PRECISION :: P(6),Q(2),PP(10),QQ(10),F(4),G(3)
1098 !C--------------------------------------------------------------------
1099 !C  Mathematical constants
1100 !C--------------------------------------------------------------------
1101 !CS    REAL, PARAMETER  ::             ONE=1.0E0,ZERO=0.0E0
1102       DOUBLE PRECISION,PARAMETER ::  ONE=1.0D0,ZERO=0.0D0
1103 !C--------------------------------------------------------------------
1104 !C  Machine-dependent constants
1105 !C--------------------------------------------------------------------
1106 !CS    REAL.PARAMETER ::             XSMALL=5.95E-8, XINF=3.40E+38 ,XMAX=85.337E0
1107       DOUBLE PRECISION,PARAMETER :: XSMALL=1.11D-16,XINF=1.79D+308,XMAX=705.342D0
1108 !--------------------------------------------------------------------
1109 !
1110 !     Coefficients for XSMALL .LE.  ARG  .LE. 1.0
1111 !
1112 !--------------------------------------------------------------------
1113 !S    DATA   P/ 5.8599221412826100000E-04, 1.3166052564989571850E-01,
1114 !S   1          1.1999463724910714109E+01, 4.6850901201934832188E+02,
1115 !S   2          5.9169059852270512312E+03, 2.4708152720399552679E+03/
1116 !S    DATA   Q/-2.4994418972832303646E+02, 2.1312714303849120380E+04/
1117 !S    DATA   F/-1.6414452837299064100E+00,-2.9601657892958843866E+02,
1118 !S   1         -1.7733784684952985886E+04,-4.0320340761145482298E+05/
1119 !S    DATA   G/-2.5064972445877992730E+02, 2.9865713163054025489E+04,
1120 !S   1         -1.6128136304458193998E+06/
1121      DATA    P/5.8599221412826100000D-04,1.3166052564989571850D-01,&
1122 &              1.1999463724910714109D+01,4.6850901201934832188D+02,&
1123 &              5.9169059852270512312D+03,2.4708152720399552679D+03/
1124      DATA    Q/-2.4994418972832303646D+02, 2.1312714303849120380D+04/
1125      DATA    F/-1.6414452837299064100D+00,-2.9601657892958843866D+02,&
1126 &              -1.7733784684952985886D+04,-4.0320340761145482298D+05/
1127      DATA    G/-2.5064972445877992730D+02, 2.9865713163054025489D+04,&
1128 &              -1.6128136304458193998D+06/
1129 !--------------------------------------------------------------------
1130 !
1131 !     Coefficients for  1.0 .LT. ARG
1132 !
1133 !--------------------------------------------------------------------
1134 !S    DATA  PP/ 1.1394980557384778174E+02, 3.6832589957340267940E+03,
1135 !S   1          3.1075408980684392399E+04, 1.0577068948034021957E+05,
1136 !S   2          1.7398867902565686251E+05, 1.5097646353289914539E+05,
1137 !S   3          7.1557062783764037541E+04, 1.8321525870183537725E+04,
1138 !S   4          2.3444738764199315021E+03, 1.1600249425076035558E+02/
1139 !S    DATA  QQ/ 2.0013443064949242491E+02, 4.4329628889746408858E+03,
1140 !S   1          3.1474655750295278825E+04, 9.7418829762268075784E+04,
1141 !S   2          1.5144644673520157801E+05, 1.2689839587977598727E+05,
1142 !S   3          5.8824616785857027752E+04, 1.4847228371802360957E+04,
1143 !S   4          1.8821890840982713696E+03, 9.2556599177304839811E+01/
1144      DATA  PP/  1.1394980557384778174D+02, 3.6832589957340267940D+03,&
1145 &               3.1075408980684392399D+04, 1.0577068948034021957D+05,&
1146 &               1.7398867902565686251D+05, 1.5097646353289914539D+05,&
1147 &               7.1557062783764037541D+04, 1.8321525870183537725D+04,&
1148 &               2.3444738764199315021D+03, 1.1600249425076035558D+02/
1149      DATA  QQ/  2.0013443064949242491D+02, 4.4329628889746408858D+03, &
1150 &                3.1474655750295278825D+04, 9.7418829762268075784D+04,&
1151 &                1.5144644673520157801D+05, 1.2689839587977598727D+05,&
1152 &                5.8824616785857027752D+04, 1.4847228371802360957D+04,&
1153 &                1.8821890840982713696D+03, 9.2556599177304839811D+01/
1154 !--------------------------------------------------------------------
1155       X = ARG
1156       IF (X .GT. ZERO) THEN
1157             IF (X .LE. ONE) THEN
1158 !--------------------------------------------------------------------
1159 !     0.0 .LT.  ARG  .LE. 1.0
1160 !--------------------------------------------------------------------
1161                   TEMP = LOG(X)
1162                   IF (X .LT. XSMALL) THEN
1163 !--------------------------------------------------------------------
1164 !     Return for small ARG
1165 !--------------------------------------------------------------------
1166                         RESULT = P(6)/Q(2) - TEMP
1167                      ELSE
1168                         XX = X * X
1169                         SUMP = ((((P(1)*XX + P(2))*XX + P(3))*XX +&
1170                                P(4))*XX + P(5))*XX + P(6)
1171                         SUMQ = (XX + Q(1))*XX + Q(2)
1172                         SUMF = ((F(1)*XX + F(2))*XX + F(3))*XX + F(4)
1173                         SUMG = ((XX + G(1))*XX + G(2))*XX + G(3)
1174                         RESULT = SUMP/SUMQ - XX*SUMF*TEMP/SUMG - TEMP
1175                         IF (JINT .EQ. 2) RESULT = RESULT * EXP(X)
1176                   END IF
1177                ELSE IF ((JINT .EQ. 1) .AND. (X .GT. XMAX)) THEN
1178 !--------------------------------------------------------------------
1179 !     Error return for ARG .GT. XMAX
1180 !--------------------------------------------------------------------
1181                   RESULT = ZERO
1182                ELSE
1183 !--------------------------------------------------------------------
1184 !     1.0 .LT. ARG
1185 !--------------------------------------------------------------------
1186                   XX = ONE / X
1187                   SUMP = PP(1)
1188                   DO 120 I = 2, 10
1189                      SUMP = SUMP*XX + PP(I)
1190   120             CONTINUE
1191                   SUMQ = XX
1192                   DO 140 I = 1, 9
1193                      SUMQ = (SUMQ + QQ(I))*XX
1194   140             CONTINUE
1195                   SUMQ = SUMQ + QQ(10)
1196                   RESULT = SUMP / SUMQ / SQRT(X)
1197                   IF (JINT .EQ. 1) RESULT = RESULT * EXP(-X)
1198             END IF
1199          ELSE
1200 !--------------------------------------------------------------------
1201 !     Error return for ARG .LE. 0.0
1202 !--------------------------------------------------------------------
1203             RESULT = XINF
1204       END IF
1205 !--------------------------------------------------------------------
1206 !     Update error counts, etc.
1207 !--------------------------------------------------------------------
1208       RETURN
1209 !---------- Last line of CALCK0 ----------
1210       END subroutine calck0

ABINIT/CALCK1 [ Functions ]

[ Top ] [ Functions ]

NAME

  CALCK1

FUNCTION

INPUTS

OUTPUT

SIDE EFFECTS

NOTES

PARENTS

CHILDREN

SOURCE

1301 SUBROUTINE CALCK1(ARG,RESULT,JINT)
1302 
1303 !--------------------------------------------------------------------
1304 !
1305 ! This packet computes modified Bessel functions of the second kind
1306 !   and order one,  K1(X)  and  EXP(X)*K1(X), for real arguments X.
1307 !   It contains two function type subprograms, BESK1  and  BESEK1,
1308 !   and one subroutine type subprogram, CALCK1.  The calling
1309 !   statements for the primary entries are
1310 !
1311 !                   Y=BESK1(X)
1312 !   and
1313 !                   Y=BESEK1(X)
1314 !
1315 !   where the entry points correspond to the functions K1(X) and
1316 !   EXP(X)*K1(X), respectively.  The routine CALCK1 is intended
1317 !   for internal packet use only, all computations within the
1318 !   packet being concentrated in this routine.  The function
1319 !   subprograms invoke CALCK1 with the statement
1320 !          CALL CALCK1(ARG,RESULT,JINT)
1321 !   where the parameter usage is as follows
1322 !
1323 !      Function                      Parameters for CALCK1
1324 !        Call             ARG                  RESULT          JINT
1325 !
1326 !     BESK1(ARG)  XLEAST .LT. ARG .LT. XMAX    K1(ARG)          1
1327 !     BESEK1(ARG)     XLEAST .LT. ARG       EXP(ARG)*K1(ARG)    2
1328 !
1329 !   The main computation evaluates slightly modified forms of near
1330 !   minimax rational approximations generated by Russon and Blair,
1331 !   Chalk River (Atomic Energy of Canada Limited) Report AECL-3461,
1332 !   1969.  This transportable program is patterned after the
1333 !   machine-dependent FUNPACK packet NATSK1, but cannot match that
1334 !   version for efficiency or accuracy.  This version uses rational
1335 !   functions that theoretically approximate K-SUB-1(X) to at
1336 !   least 18 significant decimal digits.  The accuracy achieved
1337 !   depends on the arithmetic system, the compiler, the intrinsic
1338 !   functions, and proper selection of the machine-dependent
1339 !   constants.
1340 !
1341 !*******************************************************************
1342 !*******************************************************************
1343 !
1344 ! Explanation of machine-dependent constants
1345 !
1346 !   beta   = Radix for the floating-point system
1347 !   minexp = Smallest representable power of beta
1348 !   maxexp = Smallest power of beta that overflows
1349 !   XLEAST = Smallest acceptable argument, i.e., smallest machine
1350 !            number X such that 1/X is machine representable.
1351 !   XSMALL = Argument below which BESK1(X) and BESEK1(X) may
1352 !            each be represented by 1/X.  A safe value is the
1353 !            largest X such that  1.0 + X = 1.0  to machine
1354 !            precision.
1355 !   XINF   = Largest positive machine number; approximately
1356 !            beta**maxexp
1357 !   XMAX   = Largest argument acceptable to BESK1;  Solution to
1358 !            equation:
1359 !               W(X) * (1+3/8X-15/128X**2) = beta**minexp
1360 !            where  W(X) = EXP(-X)*SQRT(PI/2X)
1361 !
1362 !
1363 !     Approximate values for some important machines are:
1364 !
1365 !                           beta       minexp       maxexp
1366 !
1367 !  CRAY-1        (S.P.)       2        -8193         8191
1368 !  Cyber 180/185
1369 !    under NOS   (S.P.)       2         -975         1070
1370 !  IEEE (IBM/XT,
1371 !    SUN, etc.)  (S.P.)       2         -126          128
1372 !  IEEE (IBM/XT,
1373 !    SUN, etc.)  (D.P.)       2        -1022         1024
1374 !  IBM 3033      (D.P.)      16          -65           63
1375 !  VAX D-Format  (D.P.)       2         -128          127
1376 !  VAX G-Format  (D.P.)       2        -1024         1023
1377 !
1378 !
1379 !                         XLEAST     XSMALL      XINF       XMAX
1380 !
1381 ! CRAY-1                1.84E-2466  3.55E-15  5.45E+2465  5674.858
1382 ! Cyber 180/855
1383 !   under NOS   (S.P.)  3.14E-294   1.77E-15  1.26E+322    672.789
1384 ! IEEE (IBM/XT,
1385 !   SUN, etc.)  (S.P.)  1.18E-38    5.95E-8   3.40E+38      85.343
1386 ! IEEE (IBM/XT,
1387 !   SUN, etc.)  (D.P.)  2.23D-308   1.11D-16  1.79D+308    705.343
1388 ! IBM 3033      (D.P.)  1.39D-76    1.11D-16  7.23D+75     177.855
1389 ! VAX D-Format  (D.P.)  5.88D-39    6.95D-18  1.70D+38      86.721
1390 ! VAX G-Format  (D.P.)  1.12D-308   5.55D-17  8.98D+307    706.728
1391 !
1392 !*******************************************************************
1393 !*******************************************************************
1394 !
1395 ! Error returns
1396 !
1397 !  The program returns the value XINF for ARG .LE. 0.0 and the
1398 !   BESK1 entry returns the value 0.0 for ARG .GT. XMAX.
1399 !
1400 !
1401 !  Intrinsic functions required are:
1402 !
1403 !     LOG, SQRT, EXP
1404 !
1405 !
1406 !  Authors: W. J. Cody and Laura Stoltz
1407 !           Mathematics and Computer Science Division
1408 !           Argonne National Laboratory
1409 !           Argonne, IL 60439
1410 !
1411 !  Latest modification: January 28, 1988
1412 !  Taken from http://www.netlib.org/specfun/k1
1413 !
1414 !--------------------------------------------------------------------
1415 
1416 !This section has been created automatically by the script Abilint (TD).
1417 !Do not modify the following lines by hand.
1418 #undef ABI_FUNC
1419 #define ABI_FUNC 'CALCK1'
1420 !End of the abilint section
1421 
1422       IMPLICIT NONE
1423       INTEGER :: I,JINT
1424 !CS    REAL
1425       DOUBLE PRECISION :: &
1426 &         ARG,F,G,ONE,P,PP,Q,QQ,RESULT,SUMF,SUMG,&
1427 &         SUMP,SUMQ,X,XINF,XMAX,XLEAST,XSMALL,XX,ZERO
1428       DIMENSION P(5),Q(3),PP(11),QQ(9),F(5),G(3)
1429 !--------------------------------------------------------------------
1430 !  Mathematical constants
1431 !--------------------------------------------------------------------
1432 !CS    DATA ONE/1.0E0/,ZERO/0.0E0/
1433        DATA ONE/1.0D0/,ZERO/0.0D0/
1434 !--------------------------------------------------------------------
1435 !  Machine-dependent constants
1436 !--------------------------------------------------------------------
1437 !CS    DATA XLEAST/1.18E-38/,XSMALL/5.95E-8/,XINF/3.40E+38/,
1438 !CS   1     XMAX/85.343E+0/
1439      DATA XLEAST/2.23D-308/,XSMALL/1.11D-16/,XINF/1.79D+308/,&
1440 &         XMAX/705.343D+0/
1441 !--------------------------------------------------------------------
1442 !  Coefficients for  XLEAST .LE.  ARG  .LE. 1.0
1443 !--------------------------------------------------------------------
1444 !CS    DATA   P/ 4.8127070456878442310E-1, 9.9991373567429309922E+1,
1445 !CS   1          7.1885382604084798576E+3, 1.7733324035147015630E+5,
1446 !CS   2          7.1938920065420586101E+5/
1447 !CS    DATA   Q/-2.8143915754538725829E+2, 3.7264298672067697862E+4,
1448 !CS   1         -2.2149374878243304548E+6/
1449 !CS    DATA   F/-2.2795590826955002390E-1,-5.3103913335180275253E+1,
1450 !CS   1         -4.5051623763436087023E+3,-1.4758069205414222471E+5,
1451 !CS   2         -1.3531161492785421328E+6/
1452 !CS    DATA   G/-3.0507151578787595807E+2, 4.3117653211351080007E+4,
1453 !CS   2         -2.7062322985570842656E+6/
1454       DATA   P/ 4.8127070456878442310D-1, 9.9991373567429309922D+1,&
1455 &             7.1885382604084798576D+3, 1.7733324035147015630D+5,&
1456 &             7.1938920065420586101D+5/
1457     DATA   Q/-2.8143915754538725829D+2, 3.7264298672067697862D+4,&
1458 &            -2.2149374878243304548D+6/
1459     DATA   F/-2.2795590826955002390D-1,-5.3103913335180275253D+1,&
1460 &            -4.5051623763436087023D+3,-1.4758069205414222471D+5,&
1461 &            -1.3531161492785421328D+6/
1462     DATA   G/-3.0507151578787595807D+2, 4.3117653211351080007D+4,&
1463 &            -2.7062322985570842656D+6/
1464 !--------------------------------------------------------------------
1465 !  Coefficients for  1.0 .LT.  ARG
1466 !--------------------------------------------------------------------
1467 !CS    DATA  PP/ 6.4257745859173138767E-2, 7.5584584631176030810E+0,
1468 !CS   1          1.3182609918569941308E+2, 8.1094256146537402173E+2,
1469 !CS   2          2.3123742209168871550E+3, 3.4540675585544584407E+3,
1470 !CS   3          2.8590657697910288226E+3, 1.3319486433183221990E+3,
1471 !CS   4          3.4122953486801312910E+2, 4.4137176114230414036E+1,
1472 !CS   5          2.2196792496874548962E+0/
1473 !CS    DATA  QQ/ 3.6001069306861518855E+1, 3.3031020088765390854E+2,
1474 !CS   1          1.2082692316002348638E+3, 2.1181000487171943810E+3,
1475 !CS   2          1.9448440788918006154E+3, 9.6929165726802648634E+2,
1476 !CS   3          2.5951223655579051357E+2, 3.4552228452758912848E+1,
1477 !CS   4          1.7710478032601086579E+0/
1478     DATA  PP/ 6.4257745859173138767D-2, 7.5584584631176030810D+0,&
1479 &             1.3182609918569941308D+2, 8.1094256146537402173D+2,&
1480 &             2.3123742209168871550D+3, 3.4540675585544584407D+3,&
1481 &             2.8590657697910288226D+3, 1.3319486433183221990D+3,&
1482 &             3.4122953486801312910D+2, 4.4137176114230414036D+1,&
1483 &             2.2196792496874548962D+0/
1484     DATA  QQ/ 3.6001069306861518855D+1, 3.3031020088765390854D+2,&
1485 &             1.2082692316002348638D+3, 2.1181000487171943810D+3,&
1486 &             1.9448440788918006154D+3, 9.6929165726802648634D+2,&
1487 &             2.5951223655579051357D+2, 3.4552228452758912848D+1,&
1488 &             1.7710478032601086579D+0/
1489 !--------------------------------------------------------------------
1490       X = ARG
1491       IF (X .LT. XLEAST) THEN
1492 !--------------------------------------------------------------------
1493 !  Error return for  ARG  .LT. XLEAST
1494 !--------------------------------------------------------------------
1495             RESULT = XINF
1496          ELSE IF (X .LE. ONE) THEN
1497 !--------------------------------------------------------------------
1498 !  XLEAST .LE.  ARG  .LE. 1.0
1499 !--------------------------------------------------------------------
1500             IF (X .LT. XSMALL) THEN
1501 !--------------------------------------------------------------------
1502 !  Return for small ARG
1503 !--------------------------------------------------------------------
1504                   RESULT = ONE / X
1505                ELSE
1506                   XX = X * X
1507                   SUMP = ((((P(1)*XX + P(2))*XX + P(3))*XX + P(4))*XX &
1508 &                        + P(5))*XX + Q(3)
1509                   SUMQ = ((XX + Q(1))*XX + Q(2))*XX + Q(3)
1510                   SUMF = (((F(1)*XX + F(2))*XX + F(3))*XX + F(4))*XX &
1511 &                        + F(5)
1512                   SUMG = ((XX + G(1))*XX + G(2))*XX + G(3)
1513                   RESULT = (XX * LOG(X) * SUMF/SUMG + SUMP/SUMQ) / X
1514                   IF (JINT .EQ. 2) RESULT = RESULT * EXP(X)
1515             END IF
1516          ELSE IF ((JINT .EQ. 1) .AND. (X .GT. XMAX)) THEN
1517 !--------------------------------------------------------------------
1518 !  Error return for  ARG  .GT. XMAX
1519 !--------------------------------------------------------------------
1520             RESULT = ZERO
1521          ELSE
1522 !--------------------------------------------------------------------
1523 !  1.0 .LT.  ARG
1524 !--------------------------------------------------------------------
1525             XX = ONE / X
1526             SUMP = PP(1)
1527             DO 120 I = 2, 11
1528                SUMP = SUMP * XX + PP(I)
1529   120       CONTINUE
1530             SUMQ = XX
1531             DO 140 I = 1, 8
1532                SUMQ = (SUMQ + QQ(I)) * XX
1533   140       CONTINUE
1534             SUMQ = SUMQ + QQ(9)
1535             RESULT = SUMP / SUMQ / SQRT(X)
1536             IF (JINT .EQ. 1) RESULT = RESULT * EXP(-X)
1537       END IF
1538       RETURN
1539 !---------- Last line of CALCK1 ----------
1540       END subroutine calck1

ABINIT/m_bessel [ Modules ]

[ Top ] [ Modules ]

NAME

 m_bessel

FUNCTION

 Bessel functions

COPYRIGHT

  Copyright (C) 2008-2018 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

17 #if defined HAVE_CONFIG_H
18 #include "config.h"
19 #endif
20 
21 #include "abi_common.h"
22 
23 module m_bessel
24 
25  implicit none
26 
27  private
28 
29  public :: CALJY0
30  public :: CALJY1
31  public :: CALCK0
32  public :: CALCK1
33 
34 contains