C ALGORITHM 779, COLLECTED ALGORITHMS FROM ACM. C THIS WORK PUBLISHED IN TRANSACTIONS ON MATHEMATICAL SOFTWARE, C VOL. 24, NO. 1, March, 1998, P. 1--12. #! /bin/sh # This is a shell archive, meaning: # 1. Remove everything above the #! /bin/sh line. # 2. Save the resulting text in a file. # 3. Execute the file with /bin/sh (not csh) to create the files: # Drivers # Src # This archive created: Tue Jun 10 09:22:46 1997 export PATH; PATH=/bin:$PATH if test ! -d 'Drivers' then mkdir 'Drivers' fi cd 'Drivers' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'RES' then echo shar: will not over-write existing file "'RES'" else cat << \SHAR_EOF > 'RES' TESTING FDM0P5 X-VALUE EXACT VALUE COMP. VALUE REL. ERROR -70.0000 0.3975450D-30 0.3975450D-30 0.0000000D+00 -50.0000 0.1928750D-21 0.1928750D-21 0.1218918D-15 -25.0000 0.1388794D-10 0.1388794D-10 0.0000000D+00 -16.0000 0.1125352D-06 0.1125352D-06 0.2352134D-15 -10.0000 0.4539847D-04 0.4539847D-04 0.0000000D+00 -7.5000 0.5528682D-03 0.5528682D-03 0.0000000D+00 -5.0000 0.6706020D-02 0.6706020D-02 0.2586815D-15 -4.0000 0.1808192D-01 0.1808192D-01 0.3837476D-15 -3.0000 0.4810264D-01 0.4810264D-01 0.2885037D-15 -2.5000 0.7761872D-01 0.7761872D-01 0.1787943D-15 -2.0000 0.1236656D+00 0.1236656D+00 0.5611013D-15 -1.5000 0.1933054D+00 0.1933054D+00 0.1435841D-15 -1.0000 0.2940276D+00 0.2940276D+00 0.0000000D+00 -0.5000 0.4312314D+00 0.4312314D+00 0.1287270D-15 -0.2500 0.5137933D+00 0.5137933D+00 0.0000000D+00 0.0000 0.6048986D+00 0.6048986D+00 0.1835387D-15 0.2500 0.7033904D+00 0.7033904D+00 0.4735164D-15 0.7500 0.9162281D+00 0.9162281D+00 0.1211732D-15 1.2500 0.1138568D+01 0.1138568D+01 0.1950210D-15 1.5000 0.1249323D+01 0.1249323D+01 0.3554638D-15 2.0000 0.1464295D+01 0.1464295D+01 0.0000000D+00 2.5000 0.1666313D+01 0.1666313D+01 0.2665101D-15 3.0000 0.1853485D+01 0.1853485D+01 0.2395969D-15 4.0000 0.2185870D+01 0.2185870D+01 0.2031636D-15 5.0000 0.2472988D+01 0.2472988D+01 0.1795760D-15 7.5000 0.3065171D+01 0.3065171D+01 0.0000000D+00 10.0000 0.3552779D+01 0.3552779D+01 0.1249977D-15 20.0000 0.5041019D+01 0.5041019D+01 0.0000000D+00 50.0000 0.7977531D+01 0.7977531D+01 0.0000000D+00 100.0000 0.1128333D+02 0.1128333D+02 0.0000000D+00 TESTING FDP0P5 X-VALUE EXACT VALUE COMP. VALUE REL. ERROR -70.0000 0.3975450D-30 0.3975450D-30 0.0000000D+00 -50.0000 0.1928750D-21 0.1928750D-21 0.1218918D-15 -25.0000 0.1388794D-10 0.1388794D-10 0.0000000D+00 -16.0000 0.1125352D-06 0.1125352D-06 0.2352134D-15 -10.0000 0.4539920D-04 0.4539920D-04 0.1492595D-15 -7.5000 0.5529762D-03 0.5529762D-03 0.3921334D-15 -5.0000 0.6721954D-02 0.6721954D-02 0.1290342D-15 -4.0000 0.1819820D-01 0.1819820D-01 0.1906478D-15 -3.0000 0.4893371D-01 0.4893371D-01 0.1418019D-15 -2.5000 0.7980385D-01 0.7980385D-01 0.0000000D+00 -2.0000 0.1292985D+00 0.1292985D+00 0.2146628D-15 -1.5000 0.2073982D+00 0.2073982D+00 0.0000000D+00 -1.0000 0.3277952D+00 0.3277952D+00 0.1693471D-15 -0.5000 0.5075371D+00 0.5075371D+00 0.2187472D-15 -0.2500 0.6254781D+00 0.6254781D+00 0.0000000D+00 0.0000 0.7651470D+00 0.7651470D+00 0.1450993D-15 0.2500 0.9285440D+00 0.9285440D+00 0.2391320D-15 0.7500 0.1332761D+01 0.1332761D+01 0.1666050D-15 1.2500 0.1846346D+01 0.1846346D+01 0.1202617D-15 1.5000 0.2144861D+01 0.2144861D+01 0.2070480D-15 2.0000 0.2823721D+01 0.2823721D+01 0.0000000D+00 2.5000 0.3606975D+01 0.3606975D+01 0.2462391D-15 3.0000 0.4487547D+01 0.4487547D+01 0.0000000D+00 4.0000 0.6511568D+01 0.6511568D+01 0.0000000D+00 5.0000 0.8844209D+01 0.8844209D+01 0.2008497D-15 7.5000 0.1579629D+02 0.1579629D+02 0.2249081D-15 10.0000 0.2408466D+02 0.2408466D+02 0.2950188D-15 20.0000 0.6749151D+02 0.6749151D+02 0.0000000D+00 50.0000 0.2660928D+03 0.2660928D+03 0.0000000D+00 100.0000 0.7523456D+03 0.7523456D+03 0.1511099D-15 TESTING FDP1P5 X-VALUE EXACT VALUE COMP. VALUE REL. ERROR -70.0000 0.3975450D-30 0.3975450D-30 0.0000000D+00 -50.0000 0.1928750D-21 0.1928750D-21 0.1218918D-15 -25.0000 0.1388794D-10 0.1388794D-10 0.4653208D-15 -16.0000 0.1125352D-06 0.1125352D-06 0.3528201D-15 -10.0000 0.4539957D-04 0.4539957D-04 0.1492583D-15 -7.5000 0.5530303D-03 0.5530303D-03 0.0000000D+00 -5.0000 0.6729941D-02 0.6729941D-02 0.2577621D-15 -4.0000 0.1825673D-01 0.1825673D-01 0.1900366D-15 -3.0000 0.4935661D-01 0.4935661D-01 0.1405869D-15 -2.5000 0.8092801D-01 0.8092801D-01 0.1714831D-15 -2.0000 0.1322468D+00 0.1322468D+00 0.0000000D+00 -1.5000 0.2149728D+00 0.2149728D+00 0.0000000D+00 -1.0000 0.3466748D+00 0.3466748D+00 0.3202491D-15 -0.5000 0.5526495D+00 0.5526495D+00 0.4017820D-15 -0.2500 0.6938464D+00 0.6938464D+00 0.1600099D-15 0.0000 0.8671999D+00 0.8671999D+00 0.1280239D-15 0.2500 0.1078398D+01 0.1078398D+01 0.2059022D-15 0.7500 0.1639285D+01 0.1639285D+01 0.1354521D-15 1.2500 0.2429426D+01 0.2429426D+01 0.0000000D+00 1.5000 0.2927749D+01 0.2927749D+01 0.1516828D-15 2.0000 0.4165414D+01 0.4165414D+01 0.2132269D-15 2.5000 0.5768879D+01 0.5768879D+01 0.1539603D-15 3.0000 0.7788611D+01 0.7788611D+01 0.4561422D-15 4.0000 0.1326049D+02 0.1326049D+02 0.1339586D-15 5.0000 0.2091447D+02 0.2091447D+02 0.0000000D+00 7.5000 0.5140755D+02 0.5140755D+02 0.1382176D-15 10.0000 0.1010051D+03 0.1010051D+03 0.1406944D-15 20.0000 0.5465630D+03 0.5465630D+03 0.0000000D+00 50.0000 0.5332354D+04 0.5332354D+04 0.1705616D-15 100.0000 0.3010867D+05 0.3010867D+05 0.1208283D-15 TESTING FDP2P5 X-VALUE EXACT VALUE COMP. VALUE REL. ERROR -70.0000 0.3975450D-30 0.3975450D-30 0.0000000D+00 -50.0000 0.1928750D-21 0.1928750D-21 0.1218918D-15 -25.0000 0.1388794D-10 0.1388794D-10 0.4653208D-15 -16.0000 0.1125352D-06 0.1125352D-06 0.0000000D+00 -10.0000 0.4539975D-04 0.4539975D-04 0.2985155D-15 -7.5000 0.5530573D-03 0.5530573D-03 0.1960379D-15 -5.0000 0.6733941D-02 0.6733941D-02 0.0000000D+00 -4.0000 0.1828612D-01 0.1828612D-01 0.0000000D+00 -3.0000 0.4957057D-01 0.4957057D-01 0.2799602D-15 -2.5000 0.8150093D-01 0.8150093D-01 0.1702777D-15 -2.0000 0.1337669D+00 0.1337669D+00 0.4149841D-15 -1.5000 0.2189495D+00 0.2189495D+00 0.1267670D-15 -1.0000 0.3568591D+00 0.3568591D+00 0.1244438D-14 -0.5000 0.5779522D+00 0.5779522D+00 0.0000000D+00 -0.2500 0.7331501D+00 0.7331501D+00 0.6057276D-15 0.0000 0.9275536D+00 0.9275536D+00 0.1196937D-15 0.2500 0.1169902D+01 0.1169902D+01 0.1897975D-15 0.7500 0.1840905D+01 0.1840905D+01 0.2412342D-15 1.2500 0.2847384D+01 0.2847384D+01 0.1559640D-15 1.5000 0.3515476D+01 0.3515476D+01 0.1263241D-15 2.0000 0.5274622D+01 0.5274622D+01 0.0000000D+00 2.5000 0.7741875D+01 0.7741875D+01 0.0000000D+00 3.0000 0.1111290D+02 0.1111290D+02 0.1598464D-15 4.0000 0.2146871D+02 0.2146871D+02 0.3309667D-15 5.0000 0.3836175D+02 0.3836175D+02 0.1852217D-15 7.5000 0.1251403D+03 0.1251403D+03 0.1135593D-15 10.0000 0.3113376D+03 0.3113376D+03 0.0000000D+00 20.0000 0.3186735D+04 0.3186735D+04 0.0000000D+00 50.0000 0.7642665D+05 0.7642665D+05 0.0000000D+00 100.0000 0.8609550D+06 0.8609550D+06 0.1352165D-15 SHAR_EOF fi # end of overwriting check if test -f 'fdtestdp.f' then echo shar: will not over-write existing file "'fdtestdp.f'" else cat << \SHAR_EOF > 'fdtestdp.f' PROGRAM TESTFD C C This program runs a very simple test procedure to assess C the accuracy of the 4 Fermi-Dirac functions, in double C precision form. C C To run the test program load the 4 Fermi-Dirac functions C FDM0P5, FDP0P5, FDP1P5, FDP2P5, and the Chebyshev series C evaluation function CHEVAL. These should have been put C into double-precision form according to the instructions C in the codes. C C The INTEGER variable IOUT defines the output stream for the C results. It is currently set to 6, but can be easily changed. C C The complete program can then be compiled and run. C INTEGER IOUT,NTEST,NUMFUN, N1, N2, N3 DOUBLE PRECISION ARMP5M(45),ARMP5P(45),AR0P5M(45),AR0P5P(45), 1 AR1P5M(45),AR1P5P(45),AR2P5M(45),AR2P5P(45), 2 EXVAL,FDM0P5,FDP0P5,FDP1P5,FDP2P5,FDVAL,RELERR,XVAL C C Define the arguments and exact function values. C DATA ARMP5M/-70.0 D 0, 1.0 D 0, 0.3975 4497 3590 8647 D -30, & -50.0 D 0, 1.0 D 0, 0.1928 7498 4796 3918 D -21, & -25.0 D 0, 1.0 D 0, 0.1388 7943 8648 2764 D -10, & -16.0 D 0, 1.0 D 0, 0.1125 3516 5764 3426 D -6, & -10.0 D 0, 1.0 D 0, 0.4539 8472 3608 0550 D -4, & -15.0 D 0, 2.0 D 0, 0.5528 6816 2177 6332 D -3, & -5.0 D 0, 1.0 D 0, 0.6706 0199 8926 8209 D -2, & -4.0 D 0, 1.0 D 0, 0.1808 1922 9913 9939 D -1, & -3.0 D 0, 1.0 D 0, 0.4810 2635 3322 0408 D -1, & -5.0 D 0, 2.0 D 0, 0.7761 8724 5948 2337 D -1, & -2.0 D 0, 1.0 D 0, 0.1236 6562 1801 2099 D 0, & -3.0 D 0, 2.0 D 0, 0.1933 0537 0255 0534 D 0, & -1.0 D 0, 1.0 D 0, 0.2940 2761 7611 4512 D 0, & -1.0 D 0, 2.0 D 0, 0.4312 3144 1926 3971 D 0, & -1.0 D 0, 4.0 D 0, 0.5137 9331 0787 9051 D 0/ DATA ARMP5P/0.0 D 0, 1.0 D 0, 0.6048 9864 3421 6304 D 0, & 1.0 D 0, 4.0 D 0, 0.7033 9042 6255 4945 D 0, & 3.0 D 0, 4.0 D 0, 0.9162 2808 8973 4550 D 0, & 5.0 D 0, 4.0 D 0, 0.1138 5676 6277 3306 D 1, & 3.0 D 0, 2.0 D 0, 0.1249 3233 4785 2712 D 1, & 2.0 D 0, 1.0 D 0, 0.1464 2945 8908 7629 D 1, & 5.0 D 0, 2.0 D 0, 0.1666 3128 8343 5831 D 1, & 3.0 D 0, 1.0 D 0, 0.1853 4850 8860 1518 D 1, & 4.0 D 0, 1.0 D 0, 0.2185 8696 9062 3812 D 1, & 5.0 D 0, 1.0 D 0, 0.2472 9876 2248 2944 D 1, & 15.0 D 0, 2.0 D 0, 0.3065 1707 1013 6853 D 1, & 10.0 D 0, 1.0 D 0, 0.3552 7792 3953 6617 D 1, & 20.0 D 0, 1.0 D 0, 0.5041 0185 0753 5329 D 1, & 50.0 D 0, 1.0 D 0, 0.7977 5308 5858 1870 D 1, & 100.0 D 0, 1.0 D 0, 0.1128 3327 4429 2768 D 2/ DATA AR0P5M/-70.0 D 0, 1.0 D 0, 0.3975 4497 3590 8647 D -30, & -50.0 D 0, 1.0 D 0, 0.1928 7498 4796 3918 D -21, & -25.0 D 0, 1.0 D 0, 0.1388 7943 8648 9583 D -10, & -16.0 D 0, 1.0 D 0, 0.1125 3517 0241 8007 D -6, & -10.0 D 0, 1.0 D 0, 0.4539 9201 0526 4133 D -4, & -15.0 D 0, 2.0 D 0, 0.5529 7624 9894 1281 D -3, & -5.0 D 0, 1.0 D 0, 0.6721 9543 1450 5913 D -2, & -4.0 D 0, 1.0 D 0, 0.1819 8203 5083 6711 D -1, & -3.0 D 0, 1.0 D 0, 0.4893 3705 6964 9578 D -1, & -5.0 D 0, 2.0 D 0, 0.7980 3854 5407 3086 D -1, & -2.0 D 0, 1.0 D 0, 0.1292 9851 3320 0756 D 0, & -3.0 D 0, 2.0 D 0, 0.2073 9818 7032 0298 D 0, & -1.0 D 0, 1.0 D 0, 0.3277 9515 9260 7115 D 0, & -1.0 D 0, 2.0 D 0, 0.5075 3710 3554 6378 D 0, & -1.0 D 0, 4.0 D 0, 0.6254 7810 2300 7525 D 0/ DATA AR0P5P/0.0 D 0, 1.0 D 0, 0.7651 4702 4625 4079 D 0, & 1.0 D 0, 4.0 D 0, 0.9285 4399 6038 6521 D 0, & 3.0 D 0, 4.0 D 0, 0.1332 7609 7454 5192 D 1, & 5.0 D 0, 4.0 D 0, 0.1846 3456 7453 9936 D 1, & 3.0 D 0, 2.0 D 0, 0.2144 8608 7758 3114 D 1, & 2.0 D 0, 1.0 D 0, 0.2823 7212 7740 1584 D 1, & 5.0 D 0, 2.0 D 0, 0.3606 9753 7795 0373 D 1, & 3.0 D 0, 1.0 D 0, 0.4487 5474 2135 1709 D 1, & 4.0 D 0, 1.0 D 0, 0.6511 5675 9275 4791 D 1, & 5.0 D 0, 1.0 D 0, 0.8844 2088 9524 2954 D 1, & 15.0 D 0, 2.0 D 0, 0.1579 6289 9021 8740 D 2, & 10.0 D 0, 1.0 D 0, 0.2408 4656 9646 3765 D 2, & 20.0 D 0, 1.0 D 0, 0.6749 1512 2216 5892 D 2, & 50.0 D 0, 1.0 D 0, 0.2660 9281 2521 3626 D 3, & 100.0 D 0, 1.0 D 0, 0.7523 4559 1552 1961 D 3/ DATA AR1P5M/-70.0 D 0, 1.0 D 0, 0.3975 4497 3590 8647 D -30, & -50.0 D 0, 1.0 D 0, 0.1928 7498 4796 3918 D -21, & -25.0 D 0, 1.0 D 0, 0.1388 7943 8649 2992 D -10, & -16.0 D 0, 1.0 D 0, 0.1125 3517 2480 5299 D -6, & -10.0 D 0, 1.0 D 0, 0.4539 9565 4045 6176 D -4, & -15.0 D 0, 2.0 D 0, 0.5530 3030 4597 1386 D -3, & -5.0 D 0, 1.0 D 0, 0.6729 9409 0901 4693 D -2, & -4.0 D 0, 1.0 D 0, 0.1825 6727 5851 0827 D -1, & -3.0 D 0, 1.0 D 0, 0.4935 6612 7906 8416 D -1, & -5.0 D 0, 2.0 D 0, 0.8092 8011 6297 5109 D -1, & -2.0 D 0, 1.0 D 0, 0.1322 4678 2251 7724 D 0, & -3.0 D 0, 2.0 D 0, 0.2149 7282 5847 4116 D 0, & -1.0 D 0, 1.0 D 0, 0.3466 7479 4799 0574 D 0, & -1.0 D 0, 2.0 D 0, 0.5526 4952 5947 3541 D 0, & -1.0 D 0, 4.0 D 0, 0.6938 4635 6450 9715 D 0/ DATA AR1P5P/0.0 D 0, 1.0 D 0, 0.8671 9988 9012 1841 D 0, & 1.0 D 0, 4.0 D 0, 0.1078 3981 5492 7105 D 1, & 3.0 D 0, 4.0 D 0, 0.1639 2853 1015 5971 D 1, & 5.0 D 0, 4.0 D 0, 0.2429 4255 4581 1454 D 1, & 3.0 D 0, 2.0 D 0, 0.2927 7494 1279 3217 D 1, & 2.0 D 0, 1.0 D 0, 0.4165 4144 5986 8322 D 1, & 5.0 D 0, 2.0 D 0, 0.5768 8793 1221 0517 D 1, & 3.0 D 0, 1.0 D 0, 0.7788 6107 7029 5970 D 1, & 4.0 D 0, 1.0 D 0, 0.1326 0488 1772 9107 D 2, & 5.0 D 0, 1.0 D 0, 0.2091 4467 4027 6263 D 2, & 15.0 D 0, 2.0 D 0, 0.5140 7547 0579 3215 D 2, & 10.0 D 0, 1.0 D 0, 0.1010 0510 0843 3260 D 3, & 20.0 D 0, 1.0 D 0, 0.5465 6301 0065 7602 D 3, & 50.0 D 0, 1.0 D 0, 0.5332 3535 6668 7146 D 4, & 100.0 D 0, 1.0 D 0, 0.3010 8671 6813 5487 D 5/ DATA AR2P5M/-70.0 D 0, 1.0 D 0, 0.3975 4497 3590 8647 D -30, & -50.0 D 0, 1.0 D 0, 0.1928 7498 4796 3918 D -21, & -25.0 D 0, 1.0 D 0, 0.1388 7943 8649 4697 D -10, & -16.0 D 0, 1.0 D 0, 0.1125 3517 3599 8945 D -6, & -10.0 D 0, 1.0 D 0, 0.4539 9747 5825 2285 D -4, & -15.0 D 0, 2.0 D 0, 0.5530 5733 5564 2952 D -3, & -5.0 D 0, 1.0 D 0, 0.6733 9406 9947 1496 D -2, & -4.0 D 0, 1.0 D 0, 0.1828 6118 4132 8075 D -1, & -3.0 D 0, 1.0 D 0, 0.4957 0567 5377 4198 D -1, & -5.0 D 0, 2.0 D 0, 0.8150 0927 5091 5113 D -1, & -2.0 D 0, 1.0 D 0, 0.1337 6692 9045 9733 D 0, & -3.0 D 0, 2.0 D 0, 0.2189 4951 7846 0930 D 0, & -1.0 D 0, 1.0 D 0, 0.3568 5914 1868 9825 D 0, & -1.0 D 0, 2.0 D 0, 0.5779 5216 0541 0087 D 0, & -1.0 D 0, 4.0 D 0, 0.7331 5006 4773 5718 D 0/ DATA AR2P5P/0.0 D 0, 1.0 D 0, 0.9275 5357 7773 9480 D 0, & 1.0 D 0, 4.0 D 0, 0.1169 9024 5238 0993 D 1, & 3.0 D 0, 4.0 D 0, 0.1840 9049 9761 4609 D 1, & 5.0 D 0, 4.0 D 0, 0.2847 3835 0168 1993 D 1, & 3.0 D 0, 2.0 D 0, 0.3515 4755 7474 9159 D 1, & 2.0 D 0, 1.0 D 0, 0.5274 6217 1262 1000 D 1, & 5.0 D 0, 2.0 D 0, 0.7741 8748 4883 8761 D 1, & 3.0 D 0, 1.0 D 0, 0.1111 2899 5227 1169 D 2, & 4.0 D 0, 1.0 D 0, 0.2146 8708 2280 1143 D 2, & 5.0 D 0, 1.0 D 0, 0.3836 1745 7340 1834 D 2, & 15.0 D 0, 2.0 D 0, 0.1251 4034 4040 0588 D 3, & 10.0 D 0, 1.0 D 0, 0.3113 3764 1850 7202 D 3, & 20.0 D 0, 1.0 D 0, 0.3186 7350 9683 1266 D 4, & 50.0 D 0, 1.0 D 0, 0.7642 6646 1278 3266 D 5, & 100.0 D 0, 1.0 D 0, 0.8609 5497 3735 2777 D 6/ C C Define the DATA output stream. C DATA IOUT/6/ C C Run the tests C DO 200 NUMFUN = 1 , 8 IF ( NUMFUN .EQ. 1 ) THEN WRITE(IOUT,1000) WRITE(IOUT,1040) ENDIF IF ( NUMFUN .EQ. 3 ) THEN WRITE(IOUT,1010) WRITE(IOUT,1040) ENDIF IF ( NUMFUN .EQ. 5 ) THEN WRITE(IOUT,1020) WRITE(IOUT,1040) ENDIF IF ( NUMFUN .EQ. 7 ) THEN WRITE(IOUT,1030) WRITE(IOUT,1040) ENDIF DO 100 NTEST = 1 , 15 N1 = 3 * NTEST - 2 N2 = N1 + 1 N3 = N2 + 1 IF ( NUMFUN .EQ. 1 ) THEN XVAL = ARMP5M(N1) / ARMP5M(N2) EXVAL = ARMP5M(N3) FDVAL = FDM0P5 ( XVAL ) ENDIF IF ( NUMFUN .EQ. 2 ) THEN XVAL = ARMP5P(N1) / ARMP5P(N2) EXVAL = ARMP5P(N3) FDVAL = FDM0P5 ( XVAL ) ENDIF IF ( NUMFUN .EQ. 3 ) THEN XVAL = AR0P5M(N1) / AR0P5M(N2) EXVAL = AR0P5M(N3) FDVAL = FDP0P5 ( XVAL ) ENDIF IF ( NUMFUN .EQ. 4 ) THEN XVAL = AR0P5P(N1) / AR0P5P(N2) EXVAL = AR0P5P(N3) FDVAL = FDP0P5 ( XVAL ) ENDIF IF ( NUMFUN .EQ. 5 ) THEN XVAL = AR1P5M(N1) / AR1P5M(N2) EXVAL = AR1P5M(N3) FDVAL = FDP1P5 ( XVAL ) ENDIF IF ( NUMFUN .EQ. 6 ) THEN XVAL = AR1P5P(N1) / AR1P5P(N2) EXVAL = AR1P5P(N3) FDVAL = FDP1P5 ( XVAL ) ENDIF IF ( NUMFUN .EQ. 7 ) THEN XVAL = AR2P5M(N1) / AR2P5M(N2) EXVAL = AR2P5M(N3) FDVAL = FDP2P5 ( XVAL ) ENDIF IF ( NUMFUN .EQ. 8 ) THEN XVAL = AR2P5P(N1) / AR2P5P(N2) EXVAL = AR2P5P(N3) FDVAL = FDP2P5 ( XVAL ) ENDIF RELERR = ABS ( ( FDVAL - EXVAL ) /EXVAL ) WRITE(IOUT,1100)XVAL,EXVAL,FDVAL,RELERR 100 CONTINUE 200 CONTINUE 1000 FORMAT(///25X,'TESTING FDM0P5') 1010 FORMAT(///25X,'TESTING FDP0P5') 1020 FORMAT(///25X,'TESTING FDP1P5') 1030 FORMAT(///25X,'TESTING FDP2P5') 1040 FORMAT(/5X,'X-VALUE',7X,'EXACT VALUE',8X,'COMP. VALUE', 1 9X,'REL. ERROR') 1100 FORMAT(2X,F10.4,3X,D16.7,3X,D16.7,3X,D16.7) END SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'RES' then echo shar: will not over-write existing file "'RES'" else cat << \SHAR_EOF > 'RES' TESTING FDM0P5 X-VALUE EXACT VALUE COMP. VALUE REL. ERROR -70.0000 0.3975450E-30 0.3975449E-30 0.1182754E-06 -50.0000 0.1928750E-21 0.1928750E-21 0.0000000E+00 -25.0000 0.1388794E-10 0.1388794E-10 0.0000000E+00 -16.0000 0.1125352E-06 0.1125352E-06 0.6313962E-07 -10.0000 0.4539847E-04 0.4539847E-04 0.0000000E+00 -7.5000 0.5528682E-03 0.5528682E-03 0.1052831E-06 -5.0000 0.6706020E-02 0.6706020E-02 0.0000000E+00 -4.0000 0.1808192E-01 0.1808192E-01 0.0000000E+00 -3.0000 0.4810264E-01 0.4810263E-01 0.7744462E-07 -2.5000 0.7761873E-01 0.7761873E-01 0.9598947E-07 -2.0000 0.1236656E+00 0.1236656E+00 0.1204956E-06 -1.5000 0.1933054E+00 0.1933054E+00 0.7708612E-07 -1.0000 0.2940276E+00 0.2940276E+00 0.0000000E+00 -0.5000 0.4312314E+00 0.4312315E+00 0.6910981E-07 -0.2500 0.5137933E+00 0.5137933E+00 0.0000000E+00 0.0000 0.6048986E+00 0.6048986E+00 0.0000000E+00 0.2500 0.7033904E+00 0.7033904E+00 0.0000000E+00 0.7500 0.9162281E+00 0.9162281E+00 0.6505437E-07 1.2500 0.1138568E+01 0.1138568E+01 0.1047011E-06 1.5000 0.1249323E+01 0.1249323E+01 0.0000000E+00 2.0000 0.1464295E+01 0.1464295E+01 0.0000000E+00 2.5000 0.1666313E+01 0.1666313E+01 0.0000000E+00 3.0000 0.1853485E+01 0.1853485E+01 0.6431629E-07 4.0000 0.2185870E+01 0.2185870E+01 0.0000000E+00 5.0000 0.2472988E+01 0.2472988E+01 0.0000000E+00 7.5000 0.3065171E+01 0.3065171E+01 0.0000000E+00 10.0000 0.3552779E+01 0.3552779E+01 0.6710763E-07 20.0000 0.5041018E+01 0.5041018E+01 0.0000000E+00 50.0000 0.7977531E+01 0.7977530E+01 0.5977252E-07 100.0000 0.1128333E+02 0.1128333E+02 0.0000000E+00 TESTING FDP0P5 X-VALUE EXACT VALUE COMP. VALUE REL. ERROR -70.0000 0.3975450E-30 0.3975449E-30 0.1182754E-06 -50.0000 0.1928750E-21 0.1928750E-21 0.0000000E+00 -25.0000 0.1388794E-10 0.1388794E-10 0.0000000E+00 -16.0000 0.1125352E-06 0.1125352E-06 0.0000000E+00 -10.0000 0.4539920E-04 0.4539920E-04 0.0000000E+00 -7.5000 0.5529763E-03 0.5529763E-03 0.1052625E-06 -5.0000 0.6721954E-02 0.6721953E-02 0.1385494E-06 -4.0000 0.1819820E-01 0.1819820E-01 0.0000000E+00 -3.0000 0.4893371E-01 0.4893370E-01 0.7612933E-07 -2.5000 0.7980385E-01 0.7980385E-01 0.0000000E+00 -2.0000 0.1292985E+00 0.1292985E+00 0.1152462E-06 -1.5000 0.2073982E+00 0.2073982E+00 0.7184808E-07 -1.0000 0.3277951E+00 0.3277953E+00 0.3636701E-06 -0.5000 0.5075371E+00 0.5075372E+00 0.1174390E-06 -0.2500 0.6254781E+00 0.6254781E+00 0.9529454E-07 0.0000 0.7651470E+00 0.7651471E+00 0.7789960E-07 0.2500 0.9285440E+00 0.9285440E+00 0.0000000E+00 0.7500 0.1332761E+01 0.1332761E+01 0.0000000E+00 1.2500 0.1846346E+01 0.1846346E+01 0.1291300E-06 1.5000 0.2144861E+01 0.2144861E+01 0.0000000E+00 2.0000 0.2823721E+01 0.2823721E+01 0.0000000E+00 2.5000 0.3606975E+01 0.3606976E+01 0.6609931E-07 3.0000 0.4487547E+01 0.4487547E+01 0.0000000E+00 4.0000 0.6511568E+01 0.6511568E+01 0.7322924E-07 5.0000 0.8844209E+01 0.8844210E+01 0.1078304E-06 7.5000 0.1579629E+02 0.1579629E+02 0.0000000E+00 10.0000 0.2408466E+02 0.2408466E+02 0.0000000E+00 20.0000 0.6749151E+02 0.6749152E+02 0.1130423E-06 50.0000 0.2660928E+03 0.2660928E+03 0.1146877E-06 100.0000 0.7523456E+03 0.7523456E+03 0.8112649E-07 TESTING FDP1P5 X-VALUE EXACT VALUE COMP. VALUE REL. ERROR -70.0000 0.3975450E-30 0.3975449E-30 0.1182754E-06 -50.0000 0.1928750E-21 0.1928750E-21 0.0000000E+00 -25.0000 0.1388794E-10 0.1388794E-10 0.0000000E+00 -16.0000 0.1125352E-06 0.1125352E-06 0.0000000E+00 -10.0000 0.4539957E-04 0.4539956E-04 0.1602649E-06 -7.5000 0.5530303E-03 0.5530303E-03 0.1052522E-06 -5.0000 0.6729941E-02 0.6729940E-02 0.6919248E-07 -4.0000 0.1825673E-01 0.1825673E-01 0.1020251E-06 -3.0000 0.4935661E-01 0.4935661E-01 0.0000000E+00 -2.5000 0.8092801E-01 0.8092801E-01 0.0000000E+00 -2.0000 0.1322468E+00 0.1322468E+00 0.0000000E+00 -1.5000 0.2149728E+00 0.2149728E+00 0.6931648E-07 -1.0000 0.3466748E+00 0.3466748E+00 0.0000000E+00 -0.5000 0.5526495E+00 0.5526495E+00 0.0000000E+00 -0.2500 0.6938463E+00 0.6938465E+00 0.1718093E-06 0.0000 0.8671999E+00 0.8672000E+00 0.6873230E-07 0.2500 0.1078398E+01 0.1078398E+01 0.1105429E-06 0.7500 0.1639285E+01 0.1639286E+01 0.1454406E-06 1.2500 0.2429425E+01 0.2429426E+01 0.9813785E-07 1.5000 0.2927749E+01 0.2927750E+01 0.8143408E-07 2.0000 0.4165414E+01 0.4165415E+01 0.1144753E-06 2.5000 0.5768879E+01 0.5768879E+01 0.0000000E+00 3.0000 0.7788611E+01 0.7788611E+01 0.6122236E-07 4.0000 0.1326049E+02 0.1326049E+02 0.0000000E+00 5.0000 0.2091447E+02 0.2091446E+02 0.9119757E-07 7.5000 0.5140755E+02 0.5140755E+02 0.7420500E-07 10.0000 0.1010051E+03 0.1010051E+03 0.0000000E+00 20.0000 0.5465630E+03 0.5465630E+03 0.0000000E+00 50.0000 0.5332354E+04 0.5332354E+04 0.9156955E-07 100.0000 0.3010867E+05 0.3010867E+05 0.6486918E-07 TESTING FDP2P5 X-VALUE EXACT VALUE COMP. VALUE REL. ERROR -70.0000 0.3975450E-30 0.3975449E-30 0.1182754E-06 -50.0000 0.1928750E-21 0.1928750E-21 0.0000000E+00 -25.0000 0.1388794E-10 0.1388794E-10 0.0000000E+00 -16.0000 0.1125352E-06 0.1125352E-06 0.6313961E-07 -10.0000 0.4539975E-04 0.4539975E-04 0.8013213E-07 -7.5000 0.5530574E-03 0.5530574E-03 0.0000000E+00 -5.0000 0.6733941E-02 0.6733940E-02 0.6915138E-07 -4.0000 0.1828612E-01 0.1828612E-01 0.1018612E-06 -3.0000 0.4957057E-01 0.4957057E-01 0.0000000E+00 -2.5000 0.8150093E-01 0.8150093E-01 0.9141713E-07 -2.0000 0.1337669E+00 0.1337669E+00 0.0000000E+00 -1.5000 0.2189495E+00 0.2189495E+00 0.6805752E-07 -1.0000 0.3568591E+00 0.3568590E+00 0.3340514E-06 -0.5000 0.5779521E+00 0.5779521E+00 0.0000000E+00 -0.2500 0.7331501E+00 0.7331500E+00 0.8129938E-07 0.0000 0.9275536E+00 0.9275536E+00 0.0000000E+00 0.2500 0.1169902E+01 0.1169902E+01 0.0000000E+00 0.7500 0.1840905E+01 0.1840905E+01 0.0000000E+00 1.2500 0.2847383E+01 0.2847383E+01 0.0000000E+00 1.5000 0.3515476E+01 0.3515476E+01 0.0000000E+00 2.0000 0.5274621E+01 0.5274621E+01 0.0000000E+00 2.5000 0.7741875E+01 0.7741874E+01 0.6159195E-07 3.0000 0.1111290E+02 0.1111290E+02 0.0000000E+00 4.0000 0.2146871E+02 0.2146871E+02 0.8884319E-07 5.0000 0.3836174E+02 0.3836175E+02 0.9944014E-07 7.5000 0.1251403E+03 0.1251404E+03 0.1219334E-06 10.0000 0.3113376E+03 0.3113376E+03 0.9802084E-07 20.0000 0.3186735E+04 0.3186735E+04 0.7661152E-07 50.0000 0.7642665E+05 0.7642665E+05 0.0000000E+00 100.0000 0.8609550E+06 0.8609550E+06 0.0000000E+00 SHAR_EOF fi # end of overwriting check if test -f 'fdtestsp.f' then echo shar: will not over-write existing file "'fdtestsp.f'" else cat << \SHAR_EOF > 'fdtestsp.f' PROGRAM TESTFD C C This program runs a very simple test procedure to assess C the accuracy of the 4 Fermi-Dirac functions, in single C precision form. C C To run the test program load the 4 Fermi-Dirac functions C FDM0P5, FDP0P5, FDP1P5, FDP2P5, and the Chebyshev series C evaluation function CHEVAL. These should have been put C into single-precision form according to the instructions C in the codes. C C The INTEGER variable IOUT defines the output stream for the C results. It is currently set to 6, but can be easily changed. C C The complete program can then be compiled and run. C INTEGER IOUT,NTEST,NUMFUN, N1, N2, N3 REAL ARMP5M(45),ARMP5P(45),AR0P5M(45),AR0P5P(45), 1 AR1P5M(45),AR1P5P(45),AR2P5M(45),AR2P5P(45), 2 EXVAL,FDM0P5,FDP0P5,FDP1P5,FDP2P5,FDVAL,RELERR,XVAL C C Define the arguments and exact function values. C DATA ARMP5M/-70.0 E 0, 1.0 E 0, 0.3975 4497 3590 8647 E -30, & -50.0 E 0, 1.0 E 0, 0.1928 7498 4796 3918 E -21, & -25.0 E 0, 1.0 E 0, 0.1388 7943 8648 2764 E -10, & -16.0 E 0, 1.0 E 0, 0.1125 3516 5764 3426 E -6, & -10.0 E 0, 1.0 E 0, 0.4539 8472 3608 0550 E -4, & -15.0 E 0, 2.0 E 0, 0.5528 6816 2177 6332 E -3, & -5.0 E 0, 1.0 E 0, 0.6706 0199 8926 8209 E -2, & -4.0 E 0, 1.0 E 0, 0.1808 1922 9913 9939 E -1, & -3.0 E 0, 1.0 E 0, 0.4810 2635 3322 0408 E -1, & -5.0 E 0, 2.0 E 0, 0.7761 8724 5948 2337 E -1, & -2.0 E 0, 1.0 E 0, 0.1236 6562 1801 2099 E 0, & -3.0 E 0, 2.0 E 0, 0.1933 0537 0255 0534 E 0, & -1.0 E 0, 1.0 E 0, 0.2940 2761 7611 4512 E 0, & -1.0 E 0, 2.0 E 0, 0.4312 3144 1926 3971 E 0, & -1.0 E 0, 4.0 E 0, 0.5137 9331 0787 9051 E 0/ DATA ARMP5P/0.0 E 0, 1.0 E 0, 0.6048 9864 3421 6304 E 0, & 1.0 E 0, 4.0 E 0, 0.7033 9042 6255 4945 E 0, & 3.0 E 0, 4.0 E 0, 0.9162 2808 8973 4550 E 0, & 5.0 E 0, 4.0 E 0, 0.1138 5676 6277 3306 E 1, & 3.0 E 0, 2.0 E 0, 0.1249 3233 4785 2712 E 1, & 2.0 E 0, 1.0 E 0, 0.1464 2945 8908 7629 E 1, & 5.0 E 0, 2.0 E 0, 0.1666 3128 8343 5831 E 1, & 3.0 E 0, 1.0 E 0, 0.1853 4850 8860 1518 E 1, & 4.0 E 0, 1.0 E 0, 0.2185 8696 9062 3812 E 1, & 5.0 E 0, 1.0 E 0, 0.2472 9876 2248 2944 E 1, & 15.0 E 0, 2.0 E 0, 0.3065 1707 1013 6853 E 1, & 10.0 E 0, 1.0 E 0, 0.3552 7792 3953 6617 E 1, & 20.0 E 0, 1.0 E 0, 0.5041 0185 0753 5329 E 1, & 50.0 E 0, 1.0 E 0, 0.7977 5308 5858 1870 E 1, & 100.0 E 0, 1.0 E 0, 0.1128 3327 4429 2768 E 2/ DATA AR0P5M/-70.0 E 0, 1.0 E 0, 0.3975 4497 3590 8647 E -30, & -50.0 E 0, 1.0 E 0, 0.1928 7498 4796 3918 E -21, & -25.0 E 0, 1.0 E 0, 0.1388 7943 8648 9583 E -10, & -16.0 E 0, 1.0 E 0, 0.1125 3517 0241 8007 E -6, & -10.0 E 0, 1.0 E 0, 0.4539 9201 0526 4133 E -4, & -15.0 E 0, 2.0 E 0, 0.5529 7624 9894 1281 E -3, & -5.0 E 0, 1.0 E 0, 0.6721 9543 1450 5913 E -2, & -4.0 E 0, 1.0 E 0, 0.1819 8203 5083 6711 E -1, & -3.0 E 0, 1.0 E 0, 0.4893 3705 6964 9578 E -1, & -5.0 E 0, 2.0 E 0, 0.7980 3854 5407 3086 E -1, & -2.0 E 0, 1.0 E 0, 0.1292 9851 3320 0756 E 0, & -3.0 E 0, 2.0 E 0, 0.2073 9818 7032 0298 E 0, & -1.0 E 0, 1.0 E 0, 0.3277 9515 9260 7115 E 0, & -1.0 E 0, 2.0 E 0, 0.5075 3710 3554 6378 E 0, & -1.0 E 0, 4.0 E 0, 0.6254 7810 2300 7525 E 0/ DATA AR0P5P/0.0 E 0, 1.0 E 0, 0.7651 4702 4625 4079 E 0, & 1.0 E 0, 4.0 E 0, 0.9285 4399 6038 6521 E 0, & 3.0 E 0, 4.0 E 0, 0.1332 7609 7454 5192 E 1, & 5.0 E 0, 4.0 E 0, 0.1846 3456 7453 9936 E 1, & 3.0 E 0, 2.0 E 0, 0.2144 8608 7758 3114 E 1, & 2.0 E 0, 1.0 E 0, 0.2823 7212 7740 1584 E 1, & 5.0 E 0, 2.0 E 0, 0.3606 9753 7795 0373 E 1, & 3.0 E 0, 1.0 E 0, 0.4487 5474 2135 1709 E 1, & 4.0 E 0, 1.0 E 0, 0.6511 5675 9275 4791 E 1, & 5.0 E 0, 1.0 E 0, 0.8844 2088 9524 2954 E 1, & 15.0 E 0, 2.0 E 0, 0.1579 6289 9021 8740 E 2, & 10.0 E 0, 1.0 E 0, 0.2408 4656 9646 3765 E 2, & 20.0 E 0, 1.0 E 0, 0.6749 1512 2216 5892 E 2, & 50.0 E 0, 1.0 E 0, 0.2660 9281 2521 3626 E 3, & 100.0 E 0, 1.0 E 0, 0.7523 4559 1552 1961 E 3/ DATA AR1P5M/-70.0 E 0, 1.0 E 0, 0.3975 4497 3590 8647 E -30, & -50.0 E 0, 1.0 E 0, 0.1928 7498 4796 3918 E -21, & -25.0 E 0, 1.0 E 0, 0.1388 7943 8649 2992 E -10, & -16.0 E 0, 1.0 E 0, 0.1125 3517 2480 5299 E -6, & -10.0 E 0, 1.0 E 0, 0.4539 9565 4045 6176 E -4, & -15.0 E 0, 2.0 E 0, 0.5530 3030 4597 1386 E -3, & -5.0 E 0, 1.0 E 0, 0.6729 9409 0901 4693 E -2, & -4.0 E 0, 1.0 E 0, 0.1825 6727 5851 0827 E -1, & -3.0 E 0, 1.0 E 0, 0.4935 6612 7906 8416 E -1, & -5.0 E 0, 2.0 E 0, 0.8092 8011 6297 5109 E -1, & -2.0 E 0, 1.0 E 0, 0.1322 4678 2251 7724 E 0, & -3.0 E 0, 2.0 E 0, 0.2149 7282 5847 4116 E 0, & -1.0 E 0, 1.0 E 0, 0.3466 7479 4799 0574 E 0, & -1.0 E 0, 2.0 E 0, 0.5526 4952 5947 3541 E 0, & -1.0 E 0, 4.0 E 0, 0.6938 4635 6450 9715 E 0/ DATA AR1P5P/0.0 E 0, 1.0 E 0, 0.8671 9988 9012 1841 E 0, & 1.0 E 0, 4.0 E 0, 0.1078 3981 5492 7105 E 1, & 3.0 E 0, 4.0 E 0, 0.1639 2853 1015 5971 E 1, & 5.0 E 0, 4.0 E 0, 0.2429 4255 4581 1454 E 1, & 3.0 E 0, 2.0 E 0, 0.2927 7494 1279 3217 E 1, & 2.0 E 0, 1.0 E 0, 0.4165 4144 5986 8322 E 1, & 5.0 E 0, 2.0 E 0, 0.5768 8793 1221 0517 E 1, & 3.0 E 0, 1.0 E 0, 0.7788 6107 7029 5970 E 1, & 4.0 E 0, 1.0 E 0, 0.1326 0488 1772 9107 E 2, & 5.0 E 0, 1.0 E 0, 0.2091 4467 4027 6263 E 2, & 15.0 E 0, 2.0 E 0, 0.5140 7547 0579 3215 E 2, & 10.0 E 0, 1.0 E 0, 0.1010 0510 0843 3260 E 3, & 20.0 E 0, 1.0 E 0, 0.5465 6301 0065 7602 E 3, & 50.0 E 0, 1.0 E 0, 0.5332 3535 6668 7146 E 4, & 100.0 E 0, 1.0 E 0, 0.3010 8671 6813 5487 E 5/ DATA AR2P5M/-70.0 E 0, 1.0 E 0, 0.3975 4497 3590 8647 E -30, & -50.0 E 0, 1.0 E 0, 0.1928 7498 4796 3918 E -21, & -25.0 E 0, 1.0 E 0, 0.1388 7943 8649 4697 E -10, & -16.0 E 0, 1.0 E 0, 0.1125 3517 3599 8945 E -6, & -10.0 E 0, 1.0 E 0, 0.4539 9747 5825 2285 E -4, & -15.0 E 0, 2.0 E 0, 0.5530 5733 5564 2952 E -3, & -5.0 E 0, 1.0 E 0, 0.6733 9406 9947 1496 E -2, & -4.0 E 0, 1.0 E 0, 0.1828 6118 4132 8075 E -1, & -3.0 E 0, 1.0 E 0, 0.4957 0567 5377 4198 E -1, & -5.0 E 0, 2.0 E 0, 0.8150 0927 5091 5113 E -1, & -2.0 E 0, 1.0 E 0, 0.1337 6692 9045 9733 E 0, & -3.0 E 0, 2.0 E 0, 0.2189 4951 7846 0930 E 0, & -1.0 E 0, 1.0 E 0, 0.3568 5914 1868 9825 E 0, & -1.0 E 0, 2.0 E 0, 0.5779 5216 0541 0087 E 0, & -1.0 E 0, 4.0 E 0, 0.7331 5006 4773 5718 E 0/ DATA AR2P5P/0.0 E 0, 1.0 E 0, 0.9275 5357 7773 9480 E 0, & 1.0 E 0, 4.0 E 0, 0.1169 9024 5238 0993 E 1, & 3.0 E 0, 4.0 E 0, 0.1840 9049 9761 4609 E 1, & 5.0 E 0, 4.0 E 0, 0.2847 3835 0168 1993 E 1, & 3.0 E 0, 2.0 E 0, 0.3515 4755 7474 9159 E 1, & 2.0 E 0, 1.0 E 0, 0.5274 6217 1262 1000 E 1, & 5.0 E 0, 2.0 E 0, 0.7741 8748 4883 8761 E 1, & 3.0 E 0, 1.0 E 0, 0.1111 2899 5227 1169 E 2, & 4.0 E 0, 1.0 E 0, 0.2146 8708 2280 1143 E 2, & 5.0 E 0, 1.0 E 0, 0.3836 1745 7340 1834 E 2, & 15.0 E 0, 2.0 E 0, 0.1251 4034 4040 0588 E 3, & 10.0 E 0, 1.0 E 0, 0.3113 3764 1850 7202 E 3, & 20.0 E 0, 1.0 E 0, 0.3186 7350 9683 1266 E 4, & 50.0 E 0, 1.0 E 0, 0.7642 6646 1278 3266 E 5, & 100.0 E 0, 1.0 E 0, 0.8609 5497 3735 2777 E 6/ C C Define the DATA output stream. C DATA IOUT/6/ C C Run the tests C DO 200 NUMFUN = 1 , 8 IF ( NUMFUN .EQ. 1 ) THEN WRITE(IOUT,1000) WRITE(IOUT,1040) ENDIF IF ( NUMFUN .EQ. 3 ) THEN WRITE(IOUT,1010) WRITE(IOUT,1040) ENDIF IF ( NUMFUN .EQ. 5 ) THEN WRITE(IOUT,1020) WRITE(IOUT,1040) ENDIF IF ( NUMFUN .EQ. 7 ) THEN WRITE(IOUT,1030) WRITE(IOUT,1040) ENDIF DO 100 NTEST = 1 , 15 N1 = 3 * NTEST - 2 N2 = N1 + 1 N3 = N2 + 1 IF ( NUMFUN .EQ. 1 ) THEN XVAL = ARMP5M(N1) / ARMP5M(N2) EXVAL = ARMP5M(N3) FDVAL = FDM0P5 ( XVAL ) ENDIF IF ( NUMFUN .EQ. 2 ) THEN XVAL = ARMP5P(N1) / ARMP5P(N2) EXVAL = ARMP5P(N3) FDVAL = FDM0P5 ( XVAL ) ENDIF IF ( NUMFUN .EQ. 3 ) THEN XVAL = AR0P5M(N1) / AR0P5M(N2) EXVAL = AR0P5M(N3) FDVAL = FDP0P5 ( XVAL ) ENDIF IF ( NUMFUN .EQ. 4 ) THEN XVAL = AR0P5P(N1) / AR0P5P(N2) EXVAL = AR0P5P(N3) FDVAL = FDP0P5 ( XVAL ) ENDIF IF ( NUMFUN .EQ. 5 ) THEN XVAL = AR1P5M(N1) / AR1P5M(N2) EXVAL = AR1P5M(N3) FDVAL = FDP1P5 ( XVAL ) ENDIF IF ( NUMFUN .EQ. 6 ) THEN XVAL = AR1P5P(N1) / AR1P5P(N2) EXVAL = AR1P5P(N3) FDVAL = FDP1P5 ( XVAL ) ENDIF IF ( NUMFUN .EQ. 7 ) THEN XVAL = AR2P5M(N1) / AR2P5M(N2) EXVAL = AR2P5M(N3) FDVAL = FDP2P5 ( XVAL ) ENDIF IF ( NUMFUN .EQ. 8 ) THEN XVAL = AR2P5P(N1) / AR2P5P(N2) EXVAL = AR2P5P(N3) FDVAL = FDP2P5 ( XVAL ) ENDIF RELERR = ABS ( ( FDVAL - EXVAL ) /EXVAL ) WRITE(IOUT,1100)XVAL,EXVAL,FDVAL,RELERR 100 CONTINUE 200 CONTINUE 1000 FORMAT(///25X,'TESTING FDM0P5') 1010 FORMAT(///25X,'TESTING FDP0P5') 1020 FORMAT(///25X,'TESTING FDP1P5') 1030 FORMAT(///25X,'TESTING FDP2P5') 1040 FORMAT(/5X,'X-VALUE',7X,'EXACT VALUE',8X,'COMP. VALUE', 1 9X,'REL. ERROR') 1100 FORMAT(2X,F10.4,3X,E16.7,3X,E16.7,3X,E16.7) END SHAR_EOF fi # end of overwriting check cd .. cd .. if test ! -d 'Src' then mkdir 'Src' fi cd 'Src' if test ! -d 'Dp' then mkdir 'Dp' fi cd 'Dp' if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << \SHAR_EOF > 'src.f' DOUBLE PRECISION FUNCTION FDM0P5(XVALUE) C C DESCRIPTION: C C This function computes the Fermi-Dirac function of C order -1/2, defined as C C Int{0 to inf} t**(-1/2) / (1+exp(t-x)) dt C FDM0P5(x) = ----------------------------------------- C Gamma(1/2) C C The function uses Chebyshev expansions which are given to C 16 decimal places for x <= 2, but only 10 decimal places C for x > 2. C C C ERROR RETURNS: C C None. C C C MACHINE-DEPENDENT CONSTANTS: C C NTERMS1 - INTEGER - The number of terms used from the array C ARRFD1. The recommended value is such that C ABS(ARRFD1(NTERMS1)) < EPS/10 C subject to 1 <= NTERMS1 <= 14. C C NTERMS2 - INTEGER - The number of terms used from the array C ARRFD2. The recommended value is such that C ABS(ARRFD2(NTERMS2)) < EPS/10 C subject to 1 <= NTERMS1 <= 23. C C NTERMS3 - INTEGER - The number of terms used from the array C ARRFD3. The recommended value is such that C ABS(ARRFD3(NTERMS3)) < EPS/10 C subject to 1 <= NTERMS3 <= 28. C C XMIN1 - REAL - The value of x below which C FDM0P5(x) = exp(x) C to machine precision. The recommended value C is LN ( SQRT(2) * EPSNEG ) C C XMIN2 - REAL - The value of x below which C FDM0P5(x) = 0.0 C to machine precision. The recommended value C is LN ( XMIN ) C C XHIGH - REAL - The value of x above which C FDM0P5(x) = 2 sqrt (x/pi) C to machine precision. The recommended value C is 1 / sqrt( 2 * EPSNEG ) C C For values of EPS, EPSNEG, and XMIN the user should refer to the C paper by Cody in ACM. Trans. Math. Soft. Vol. 14 (1988) p303-311. C C This code is provided with single and double precision values C of the machine-dependent parameters, suitable for machines C which satisfy the IEEE floating-point standard. C C C AUTHOR: C DR. ALLAN MACLEOD, C DEPT. OF MATHEMATICS AND STATISTICS, C UNIVERSITY OF PAISLEY, C HIGH ST., C PAISLEY, C SCOTLAND C PA1 2BE C C (e-mail: macl-ms0@paisley.ac.uk ) C C C LATEST UPDATE: C 20 NOVEMBER, 1996 C INTEGER NTERM1,NTERM2,NTERM3 DOUBLE PRECISION 1 ARRFD1(0:14),ARRFD2(0:23),ARRFD3(0:58), 2 CHEVAL,CHV,EXPX,FIFTY,FORTY2, 3 GAM1P5,ONE,T,THREE,TWO,TWOE, 4 X,XHIGH,XMIN1,XMIN2,XSQ,XVALUE,ZERO DATA ARRFD1/1.7863 5963 8510 2264 D 0, 1 -0.9993 7200 7632 333 D -1, 2 0.6414 4652 2160 54 D -2, 3 -0.4356 4153 7134 5 D -3, 4 0.3052 1670 0310 D -4, 5 -0.2181 0648 110 D -5, 6 0.1580 0507 81 D -6, 7 -0.1156 2057 0 D -7, 8 0.8525 860 D -9, 9 -0.6325 29 D -10, X 0.4715 9 D -11, 1 -0.3530 D -12, 2 0.265 D -13, 3 -0.20 D -14, 4 0.2 D -15/ DATA ARRFD2( 0)/ 1.6877 1115 2605 2352 D 0/ DATA ARRFD2( 1)/ 0.5978 3602 2633 6983 D 0/ DATA ARRFD2( 2)/ 0.3572 2600 4541 669 D -1/ DATA ARRFD2( 3)/-0.1321 4478 6506 426 D -1/ DATA ARRFD2( 4)/-0.4040 1342 0744 7 D -3/ DATA ARRFD2( 5)/ 0.5330 0118 4688 7 D -3/ DATA ARRFD2( 6)/-0.1489 2350 4863 D -4/ DATA ARRFD2( 7)/-0.2188 6382 2916 D -4/ DATA ARRFD2( 8)/ 0.1965 2084 277 D -5/ DATA ARRFD2( 9)/ 0.8565 8304 66 D -6/ DATA ARRFD2(10)/-0.1407 7231 33 D -6/ DATA ARRFD2(11)/-0.3051 7580 3 D -7/ DATA ARRFD2(12)/ 0.8352 4532 D -8/ DATA ARRFD2(13)/ 0.9025 750 D -9/ DATA ARRFD2(14)/-0.4455 471 D -9/ DATA ARRFD2(15)/-0.1483 42 D -10/ DATA ARRFD2(16)/ 0.2192 66 D -10/ DATA ARRFD2(17)/-0.6579 D -12/ DATA ARRFD2(18)/-0.1000 9 D -11/ DATA ARRFD2(19)/ 0.936 D -13/ DATA ARRFD2(20)/ 0.420 D -13/ DATA ARRFD2(21)/-0.71 D -14/ DATA ARRFD2(22)/-0.16 D -14/ DATA ARRFD2(23)/ 0.4 D -15/ DATA ARRFD3(0)/ 0.8707 1950 2959 0563 D 0/ DATA ARRFD3(1)/ 0.5983 3110 2317 33 D -2/ DATA ARRFD3(2)/ -0.4326 7047 0895 746 D -1/ DATA ARRFD3(3)/ -0.3930 8368 1608 590 D -1/ DATA ARRFD3(4)/ -0.1914 8268 8045 932 D -1/ DATA ARRFD3(5)/ -0.6558 2880 9801 58 D -2/ DATA ARRFD3(6)/ -0.2227 6691 5163 12 D -2/ DATA ARRFD3(7)/ -0.8466 7869 3617 8 D -3/ DATA ARRFD3(8)/ -0.2807 4594 8921 9 D -3/ DATA ARRFD3(9)/ -0.9555 7502 4348 D -4/ DATA ARRFD3(10)/-0.3623 6766 2803 D -4/ DATA ARRFD3(11)/-0.1091 5846 8869 D -4/ DATA ARRFD3(12)/-0.3935 6701 000 D -5/ DATA ARRFD3(13)/-0.1310 8192 725 D -5/ DATA ARRFD3(14)/-0.2468 8163 88 D -6/ DATA ARRFD3(15)/-0.1048 3803 11 D -6/ DATA ARRFD3(16)/ 0.2361 8148 7 D -7/ DATA ARRFD3(17)/ 0.2271 4535 9 D -7/ DATA ARRFD3(18)/ 0.1457 7517 4 D -7/ DATA ARRFD3(19)/ 0.1539 2676 7 D -7/ DATA ARRFD3(20)/ 0.5692 4772 D -8/ DATA ARRFD3(21)/ 0.5062 3068 D -8/ DATA ARRFD3(22)/ 0.2342 6075 D -8/ DATA ARRFD3(23)/ 0.1265 2275 D -8/ DATA ARRFD3(24)/ 0.8927 773 D -9/ DATA ARRFD3(25)/ 0.2994 501 D -9/ DATA ARRFD3(26)/ 0.2822 785 D -9/ DATA ARRFD3(27)/ 0.9106 85 D -10/ DATA ARRFD3(28)/ 0.6962 85 D -10/ DATA ARRFD3(29)/ 0.3662 25 D -10/ DATA ARRFD3(30)/ 0.1243 51 D -10/ DATA ARRFD3(31)/ 0.1450 19 D -10/ DATA ARRFD3(32)/ 0.1664 5 D -11/ DATA ARRFD3(33)/ 0.4585 6 D -11/ DATA ARRFD3(34)/ 0.6092 D -12/ DATA ARRFD3(35)/ 0.9331 D -12/ DATA ARRFD3(36)/ 0.5238 D -12/ DATA ARRFD3(37)/-0.56 D -14/ DATA ARRFD3(38)/ 0.3170 D -12/ DATA ARRFD3(39)/-0.926 D -13/ DATA ARRFD3(40)/ 0.1265 D -12/ DATA ARRFD3(41)/-0.327 D -13/ DATA ARRFD3(42)/ 0.276 D -13/ DATA ARRFD3(43)/ 0.33 D -14/ DATA ARRFD3(44)/-0.42 D -14/ DATA ARRFD3(45)/ 0.101 D -13/ DATA ARRFD3(46)/-0.73 D -14/ DATA ARRFD3(47)/ 0.64 D -14/ DATA ARRFD3(48)/-0.37 D -14/ DATA ARRFD3(49)/ 0.23 D -14/ DATA ARRFD3(50)/-0.9 D -15/ DATA ARRFD3(51)/ 0.2 D -15/ DATA ARRFD3(52)/ 0.2 D -15/ DATA ARRFD3(53)/-0.3 D -15/ DATA ARRFD3(54)/ 0.4 D -15/ DATA ARRFD3(55)/-0.3 D -15/ DATA ARRFD3(56)/ 0.2 D -15/ DATA ARRFD3(57)/-0.1 D -15/ DATA ARRFD3(58)/ 0.1 D -15/ DATA ZERO,ONE,TWO/ 0.0 D 0 , 1.0 D 0 , 2.0 D 0/ DATA THREE,FORTY2,FIFTY/ 3.0 D 0 , 42.0 D 0 , 50.0 D 0/ DATA GAM1P5/0.8862 2692 5452 7580 D 0/ DATA TWOE/5.4365 6365 6918 0905 D 0/ C C Machine-dependent constants C DATA NTERM1,NTERM2,NTERM3/14,23,58/ DATA XMIN1,XMIN2,XHIGH/-36.39023D0,-708.39641D0,67108864.0D0/ C C Start calculation C X=XVALUE C C Code for x < -1 C IF ( X .LT. -ONE ) THEN IF ( X .GT. XMIN1 ) THEN EXPX = EXP(X) T = TWOE * EXPX - ONE FDM0P5 = EXPX * CHEVAL ( NTERM1 , ARRFD1 , T ) ELSE IF ( X .LT. XMIN2 ) THEN FDM0P5 = ZERO ELSE FDM0P5 = EXP(X) ENDIF ENDIF ELSE C C Code for -1 <= x <= 2 C IF ( X .LE. TWO ) THEN T = ( TWO * X - ONE ) / THREE FDM0P5 = CHEVAL ( NTERM2 , ARRFD2 , T ) ELSE C C Code for x > 2 C FDM0P5 = SQRT(X) / GAM1P5 IF ( X .LE. XHIGH ) THEN XSQ = X * X T = ( FIFTY - XSQ ) / ( FORTY2 + XSQ ) CHV = CHEVAL ( NTERM3 , ARRFD3 , T ) FDM0P5 = FDM0P5 * ( ONE - CHV / XSQ ) ENDIF ENDIF ENDIF RETURN END DOUBLE PRECISION FUNCTION FDP0P5(XVALUE) C C DESCRIPTION: C C This function computes the Fermi-Dirac function of C order 1/2, defined as C C Int{0 to inf} t**(1/2) / (1+exp(t-x)) dt C FDP0P5(x) = ----------------------------------------- C Gamma(3/2) C C The function uses Chebyshev expansions which are given to C 16 decimal places for x <= 2, but only 10 decimal places C for x > 2. C C C ERROR RETURNS: C C If XVALUE too large and positive, the function value C will overflow. An error message is printed and the function C returns the value 0.0. C C C MACHINE-DEPENDENT CONSTANTS: C C NTERMS1 - INTEGER - The number of terms used from the array C ARRFD1. The recommended value is such that C ABS(ARRFD1(NTERMS1)) < EPS/10 C subject to 1 <= NTERMS1 <= 13. C C NTERMS2 - INTEGER - The number of terms used from the array C ARRFD2. The recommended value is such that C ABS(ARRFD2(NTERMS2)) < EPS/10 C subject to 1 <= NTERMS1 <= 23. C C NTERMS3 - INTEGER - The number of terms used from the array C ARRFD3. The recommended value is such that C ABS(ARRFD3(NTERMS3)) < EPS/10 C subject to 1 <= NTERMS3 <= 32. C C XMIN1 - REAL - The value of x below which C FDP0P5(x) = exp(x) C to machine precision. The recommended value C is 1.5*LN(2) + LN(EPSNEG) C C XMIN2 - REAL - The value of x below which C FDP0P5(x) = 0.0 C to machine precision. The recommended value C is LN ( XMIN ) C C XHIGH1 - REAL - The value of x above which C FDP0P5(x) = x**(3/2)/GAMMA(5/2) C to machine precision. The recommended value C is pi / SQRT(8*EPS) C C XHIGH2 - REAL - The value of x above which FDP0P5 would C overflow. The reommended value is C (1.329*XMAX)**(2/3) C C For values of EPS, EPSNEG, and XMIN the user should refer to the C paper by Cody in ACM. Trans. Math. Soft. Vol. 14 (1988) p303-311. C C This code is provided with single and double precision values C of the machine-dependent parameters, suitable for machines C which satisfy the IEEE floating-point standard. C C C AUTHOR: C DR. ALLAN MACLEOD, C DEPT. OF MATHEMATICS AND STATISTICS, C UNIVERSITY OF PAISLEY, C HIGH ST., C PAISLEY, C SCOTLAND C PA1 2BE C C (e-mail: macl-ms0@paisley.ac.uk ) C C C LATEST UPDATE: C 20 NOVEMBER, 1996 C INTEGER NTERM1,NTERM2,NTERM3 DOUBLE PRECISION 1 ARRFD1(0:13),ARRFD2(0:23),ARRFD3(0:53), 2 CHEVAL,CHV,EXPX,FIFTY,FORTY2, 3 GAM2P5,ONE,T,THREE,TWO,TWOE,X,XHIGH1, 4 XHIGH2,XMIN1,XMIN2,XSQ,XVALUE,ZERO DATA ARRFD1/1.8862 9683 9273 4597 D 0, 1 -0.5435 8081 7644 053 D -1, 2 0.2364 4975 4397 20 D -2, 3 -0.1216 9293 6588 0 D -3, 4 0.6869 5130 622 D -5, 5 -0.4112 0761 72 D -6, 6 0.2563 5162 8 D -7, 7 -0.1646 5008 D -8, 8 0.1081 948 D -9, 9 -0.7239 2 D -11, X 0.4915 D -12, 1 -0.338 D -13, 2 0.23 D -14, 3 -0.2 D -15/ DATA ARRFD2( 0)/ 2.6982 4927 8817 0612 D 0/ DATA ARRFD2( 1)/ 1.2389 9141 4113 3012 D 0/ DATA ARRFD2( 2)/ 0.2291 4393 7981 6278 D 0/ DATA ARRFD2( 3)/ 0.9031 6534 6872 79 D -2/ DATA ARRFD2( 4)/-0.2577 6524 6912 46 D -2/ DATA ARRFD2( 5)/-0.5836 8160 5388 D -4/ DATA ARRFD2( 6)/ 0.6936 0945 8725 D -4/ DATA ARRFD2( 7)/-0.1806 1670 265 D -5/ DATA ARRFD2( 8)/-0.2132 1530 005 D -5/ DATA ARRFD2( 9)/ 0.1754 9839 51 D -6/ DATA ARRFD2(10)/ 0.6653 2547 0 D -7/ DATA ARRFD2(11)/-0.1016 7597 7 D -7/ DATA ARRFD2(12)/-0.1963 7597 D -8/ DATA ARRFD2(13)/ 0.5075 769 D -9/ DATA ARRFD2(14)/ 0.4914 69 D -10/ DATA ARRFD2(15)/-0.2337 37 D -10/ DATA ARRFD2(16)/-0.6645 D -12/ DATA ARRFD2(17)/ 0.1011 5 D -11/ DATA ARRFD2(18)/-0.313 D -13/ DATA ARRFD2(19)/-0.412 D -13/ DATA ARRFD2(20)/ 0.38 D -14/ DATA ARRFD2(21)/ 0.16 D -14/ DATA ARRFD2(22)/-0.3 D -15/ DATA ARRFD2(23)/-0.1 D -15/ DATA ARRFD3(0)/ 2.5484 3841 9800 9122 D 0/ DATA ARRFD3(1)/ 0.5104 3940 8960 652 D -1/ DATA ARRFD3(2)/ 0.7749 3527 6282 94 D -2/ DATA ARRFD3(3)/ -0.7504 1656 5849 53 D -2/ DATA ARRFD3(4)/ -0.7754 0826 3202 96 D -2/ DATA ARRFD3(5)/ -0.4581 0844 5399 77 D -2/ DATA ARRFD3(6)/ -0.2343 1641 5873 63 D -2/ DATA ARRFD3(7)/ -0.1178 8049 5135 91 D -2/ DATA ARRFD3(8)/ -0.5802 7393 5970 2 D -3/ DATA ARRFD3(9)/ -0.2825 3507 0053 7 D -3/ DATA ARRFD3(10)/-0.1388 1366 5179 9 D -3/ DATA ARRFD3(11)/-0.6806 9508 4875 D -4/ DATA ARRFD3(12)/-0.3353 5635 0608 D -4/ DATA ARRFD3(13)/-0.1665 3301 8734 D -4/ DATA ARRFD3(14)/-0.8271 4908 266 D -5/ DATA ARRFD3(15)/-0.4142 5714 409 D -5/ DATA ARRFD3(16)/-0.2080 5255 294 D -5/ DATA ARRFD3(17)/-0.1047 9767 478 D -5/ DATA ARRFD3(18)/-0.5315 2738 02 D -6/ DATA ARRFD3(19)/-0.2694 0611 78 D -6/ DATA ARRFD3(20)/-0.1374 8787 49 D -6/ DATA ARRFD3(21)/-0.7023 0888 7 D -7/ DATA ARRFD3(22)/-0.3595 4394 2 D -7/ DATA ARRFD3(23)/-0.1851 0612 6 D -7/ DATA ARRFD3(24)/-0.9502 3937 D -8/ DATA ARRFD3(25)/-0.4918 4811 D -8/ DATA ARRFD3(26)/-0.2537 1950 D -8/ DATA ARRFD3(27)/-0.1315 1532 D -8/ DATA ARRFD3(28)/-0.6835 168 D -9/ DATA ARRFD3(29)/-0.3538 244 D -9/ DATA ARRFD3(30)/-0.1853 182 D -9/ DATA ARRFD3(31)/-0.9589 83 D -10/ DATA ARRFD3(32)/-0.5040 83 D -10/ DATA ARRFD3(33)/-0.2622 38 D -10/ DATA ARRFD3(34)/-0.1372 55 D -10/ DATA ARRFD3(35)/-0.7234 0 D -11/ DATA ARRFD3(36)/-0.3742 9 D -11/ DATA ARRFD3(37)/-0.2005 9 D -11/ DATA ARRFD3(38)/-0.1026 9 D -11/ DATA ARRFD3(39)/-0.5551 D -12/ DATA ARRFD3(40)/-0.2857 D -12/ DATA ARRFD3(41)/-0.1520 D -12/ DATA ARRFD3(42)/-0.811 D -13/ DATA ARRFD3(43)/-0.410 D -13/ DATA ARRFD3(44)/-0.234 D -13/ DATA ARRFD3(45)/-0.110 D -13/ DATA ARRFD3(46)/-0.67 D -14/ DATA ARRFD3(47)/-0.30 D -14/ DATA ARRFD3(48)/-0.19 D -14/ DATA ARRFD3(49)/-0.9 D -15/ DATA ARRFD3(50)/-0.5 D -15/ DATA ARRFD3(51)/-0.3 D -15/ DATA ARRFD3(52)/-0.1 D -15/ DATA ARRFD3(53)/-0.1 D -15/ DATA ZERO,ONE,TWO/ 0.0 D 0 , 1.0 D 0 , 2.0 D 0/ DATA THREE,FORTY2,FIFTY/ 3.0 D 0 , 42.0 D 0 , 50.0 D 0/ DATA GAM2P5/0.1329 3403 8817 9137 D 1/ DATA TWOE/5.4365 6365 6918 0905 D 0/ C C Machine-dependent constants (suitable for IEEE machines) C DATA NTERM1,NTERM2,NTERM3/13,23,53/ DATA XMIN1,XMIN2/-35.7D00,-708.394D00/ DATA XHIGH1,XHIGH2/7.45467D7,3.8392996D205/ C C Start calculation C X=XVALUE C C Test for error condition C IF ( X .GT. XHIGH2 ) THEN PRINT*,'** ERROR ** - X TOO LARGE FOR FDP0P5' FDP0P5 = ZERO RETURN ENDIF C C Code for x < -1 C IF ( X .LT. -ONE ) THEN IF ( X .GT. XMIN1 ) THEN EXPX = EXP(X) T = TWOE * EXPX - ONE FDP0P5 = EXPX * CHEVAL ( NTERM1 , ARRFD1 , T ) ELSE IF ( X .LT. XMIN2 ) THEN FDP0P5 = ZERO ELSE FDP0P5 = EXP(X) ENDIF ENDIF ELSE C C Code for -1 <= x <= 2 C IF ( X .LE. TWO ) THEN T = ( TWO * X - ONE ) / THREE FDP0P5 = CHEVAL ( NTERM2 , ARRFD2 , T ) ELSE C C Code for x > 2 C FDP0P5 = X * SQRT(X) / GAM2P5 IF ( X .LE. XHIGH1 ) THEN XSQ = X * X T = ( FIFTY - XSQ ) / ( FORTY2 + XSQ ) CHV = CHEVAL ( NTERM3 , ARRFD3 , T ) FDP0P5 = FDP0P5 * ( ONE + CHV / XSQ ) ENDIF ENDIF ENDIF RETURN END DOUBLE PRECISION FUNCTION FDP1P5(XVALUE) C C DESCRIPTION: C C This function computes the Fermi-Dirac function of C order 3/2, defined as C C Int{0 to inf} t**(3/2) / (1+exp(t-x)) dt C FDP1P5(x) = ----------------------------------------- C Gamma(5/2) C C The function uses Chebyshev expansions which are given to C 16 decimal places for x <= 2, but only 10 decimal places C for x > 2. C C C ERROR RETURNS: C C If XVALUE too large and positive, the function value C will overflow. An error message is printed and the function C returns the value 0.0. C C C MACHINE-DEPENDENT CONSTANTS: C C NTERMS1 - INTEGER - The number of terms used from the array C ARRFD1. The recommended value is such that C ABS(ARRFD1(NTERMS1)) < EPS/10 C subject to 1 <= NTERMS1 <= 12. C C NTERMS2 - INTEGER - The number of terms used from the array C ARRFD2. The recommended value is such that C ABS(ARRFD2(NTERMS2)) < EPS/10 C subject to 1 <= NTERMS1 <= 22. C C NTERMS3 - INTEGER - The number of terms used from the array C ARRFD3. The recommended value is such that C ABS(ARRFD3(NTERMS3)) < EPS/10 C subject to 1 <= NTERMS3 <= 33. C C XMIN1 - REAL - The value of x below which C FDP1P5(x) = exp(x) C to machine precision. The recommended value C is 2.5*LN(2) + LN(EPSNEG) C C XMIN2 - REAL - The value of x below which C FDP1P5(x) = 0.0 C to machine precision. The recommended value C is LN ( XMIN ) C C XHIGH1 - REAL - The value of x above which C FDP1P5(x) = x**(5/2)/GAMMA(7/2) C to machine precision. The recommended value C is pi * SQRT(1.6/EPS) C C XHIGH2 - REAL - The value of x above which FDP1P5 would C overflow. The reommended value is C (3.233509*XMAX)**(2/5) C C For values of EPS, EPSNEG, and XMIN the user should refer to the C paper by Cody in ACM. Trans. Math. Soft. Vol. 14 (1988) p303-311. C C This code is provided with single and double precision values C of the machine-dependent parameters, suitable for machines C which satisfy the IEEE floating-point standard. C C C AUTHOR: C DR. ALLAN MACLEOD, C DEPT. OF MATHEMATICS AND STATISTICS, C UNIVERSITY OF PAISLEY, C HIGH ST., C PAISLEY, C SCOTLAND C PA1 2BE C C (e-mail: macl_ms0@paisley.ac.uk ) C C C LATEST UPDATE: C 21 NOVEMBER, 1996 C INTEGER NTERM1,NTERM2,NTERM3 DOUBLE PRECISION 1 ARRFD1(0:12),ARRFD2(0:22),ARRFD3(0:55), 2 CHEVAL,CHV,EXPX,FIFTY,FORTY2, 3 GAM3P5,ONE,T,THREE,TWO,TWOE,X,XHIGH1, 4 XHIGH2,XMIN1,XMIN2,XSQ,XVALUE,ZERO DATA ARRFD1/1.9406 5492 1037 8650 D 0, 1 -0.2878 6747 5518 043 D -1, 2 0.8509 1579 5231 3 D -3, 3 -0.3327 8452 5669 D -4, 4 0.1517 1202 058 D -5, 5 -0.7622 0087 4 D -7, 6 0.4095 5489 D -8, 7 -0.2311 964 D -9, 8 0.1355 37 D -10, 9 -0.8187 D -12, X 0.507 D -13, 1 -0.32 D -14, 2 0.2 D -15/ DATA ARRFD2( 0)/ 3.5862 2516 1563 4306 D 0/ DATA ARRFD2( 1)/ 1.8518 2900 5626 5751 D 0/ DATA ARRFD2( 2)/ 0.4612 3491 0241 7150 D 0/ DATA ARRFD2( 3)/ 0.5793 0397 6126 881 D -1/ DATA ARRFD2( 4)/ 0.1704 3790 5548 75 D -2/ DATA ARRFD2( 5)/-0.3970 5201 2249 6 D -3/ DATA ARRFD2( 6)/-0.7070 2491 890 D -5/ DATA ARRFD2( 7)/ 0.7659 9748 792 D -5/ DATA ARRFD2( 8)/-0.1857 8113 33 D -6/ DATA ARRFD2( 9)/-0.1832 2379 56 D -6/ DATA ARRFD2(10)/ 0.1392 4949 5 D -7/ DATA ARRFD2(11)/ 0.4670 2027 D -8/ DATA ARRFD2(12)/-0.6671 984 D -9/ DATA ARRFD2(13)/-0.1161 292 D -9/ DATA ARRFD2(14)/ 0.2844 38 D -10/ DATA ARRFD2(15)/ 0.2490 6 D -11/ DATA ARRFD2(16)/-0.1143 1 D -11/ DATA ARRFD2(17)/-0.279 D -13/ DATA ARRFD2(18)/ 0.439 D -13/ DATA ARRFD2(19)/-0.14 D -14/ DATA ARRFD2(20)/-0.16 D -14/ DATA ARRFD2(21)/ 0.1 D -15/ DATA ARRFD2(22)/ 0.1 D -15/ DATA ARRFD3( 0)/12.1307 5817 3688 4627 D 0/ DATA ARRFD3( 1)/-0.1547 5011 1128 7255 D 0/ DATA ARRFD3( 2)/-0.7390 0738 8850 999 D -1/ DATA ARRFD3( 3)/-0.3072 3537 7959 258 D -1/ DATA ARRFD3( 4)/-0.1145 4857 9330 328 D -1/ DATA ARRFD3( 5)/-0.4056 7636 8095 39 D -2/ DATA ARRFD3( 6)/-0.1398 0158 3732 27 D -2/ DATA ARRFD3( 7)/-0.4454 9018 1015 3 D -3/ DATA ARRFD3( 8)/-0.1173 9461 1270 4 D -3/ DATA ARRFD3( 9)/-0.1484 0898 0093 D -4/ DATA ARRFD3(10)/ 0.1188 9515 4223 D -4/ DATA ARRFD3(11)/ 0.1464 7695 8178 D -4/ DATA ARRFD3(12)/ 0.1132 2874 1730 D -4/ DATA ARRFD3(13)/ 0.7576 2292 948 D -5/ DATA ARRFD3(14)/ 0.4712 0400 466 D -5/ DATA ARRFD3(15)/ 0.2813 2628 202 D -5/ DATA ARRFD3(16)/ 0.1637 0517 341 D -5/ DATA ARRFD3(17)/ 0.9351 0762 72 D -6/ DATA ARRFD3(18)/ 0.5278 6892 10 D -6/ DATA ARRFD3(19)/ 0.2951 0798 70 D -6/ DATA ARRFD3(20)/ 0.1638 6001 90 D -6/ DATA ARRFD3(21)/ 0.9052 0540 9 D -7/ DATA ARRFD3(22)/ 0.4977 5697 5 D -7/ DATA ARRFD3(23)/ 0.2729 5586 3 D -7/ DATA ARRFD3(24)/ 0.1492 1458 5 D -7/ DATA ARRFD3(25)/ 0.8142 0359 D -8/ DATA ARRFD3(26)/ 0.4434 9200 D -8/ DATA ARRFD3(27)/ 0.2411 6032 D -8/ DATA ARRFD3(28)/ 0.1310 5018 D -8/ DATA ARRFD3(29)/ 0.7109 736 D -9/ DATA ARRFD3(30)/ 0.3856 721 D -9/ DATA ARRFD3(31)/ 0.2089 529 D -9/ DATA ARRFD3(32)/ 0.1131 735 D -9/ DATA ARRFD3(33)/ 0.6127 85 D -10/ DATA ARRFD3(34)/ 0.3314 48 D -10/ DATA ARRFD3(35)/ 0.1794 19 D -10/ DATA ARRFD3(36)/ 0.9695 3 D -11/ DATA ARRFD3(37)/ 0.5246 3 D -11/ DATA ARRFD3(38)/ 0.2834 3 D -11/ DATA ARRFD3(39)/ 0.1532 3 D -11/ DATA ARRFD3(40)/ 0.8284 D -12/ DATA ARRFD3(41)/ 0.4472 D -12/ DATA ARRFD3(42)/ 0.2421 D -12/ DATA ARRFD3(43)/ 0.1304 D -12/ DATA ARRFD3(44)/ 0.707 D -13/ DATA ARRFD3(45)/ 0.381 D -13/ DATA ARRFD3(46)/ 0.206 D -13/ DATA ARRFD3(47)/ 0.111 D -13/ DATA ARRFD3(48)/ 0.60 D -14/ DATA ARRFD3(49)/ 0.33 D -14/ DATA ARRFD3(50)/ 0.17 D -14/ DATA ARRFD3(51)/ 0.11 D -14/ DATA ARRFD3(52)/ 0.5 D -15/ DATA ARRFD3(53)/ 0.3 D -15/ DATA ARRFD3(54)/ 0.1 D -15/ DATA ARRFD3(55)/ 0.1 D -15/ DATA ZERO,ONE,TWO/ 0.0 D 0 , 1.0 D 0 , 2.0 D 0/ DATA THREE,FORTY2,FIFTY/ 3.0 D 0 , 42.0 D 0 , 50.0 D 0/ DATA GAM3P5/0.3323 3509 7044 7843 D 1/ DATA TWOE/5.4365 6365 6918 0905 D 0/ C C Machine-dependent constants (suitable for IEEE machines) C DATA NTERM1,NTERM2,NTERM3/12,22,55/ DATA XMIN1,XMIN2/-35.004D0,-708.396418D0/ DATA XHIGH1,XHIGH2/166674733.2D0,3.204467D123/ C C Start calculation C X=XVALUE C C Test for error condition C IF ( X .GT. XHIGH2 ) THEN PRINT*,'** ERROR ** - X TOO LARGE FOR FDP1P5' FDP1P5 = ZERO RETURN ENDIF C C Code for x < -1 C IF ( X .LT. -ONE ) THEN IF ( X .GT. XMIN1 ) THEN EXPX = EXP(X) T = TWOE * EXPX - ONE FDP1P5 = EXPX * CHEVAL ( NTERM1 , ARRFD1 , T ) ELSE IF ( X .LT. XMIN2 ) THEN FDP1P5 = ZERO ELSE FDP1P5 = EXP(X) ENDIF ENDIF ELSE C C Code for -1 <= x <= 2 C IF ( X .LE. TWO ) THEN T = ( TWO * X - ONE ) / THREE FDP1P5 = CHEVAL ( NTERM2 , ARRFD2 , T ) ELSE C C Code for x > 2 C FDP1P5 = X * X * SQRT(X) / GAM3P5 IF ( X .LE. XHIGH1 ) THEN XSQ = X * X T = ( FIFTY - XSQ ) / ( FORTY2 + XSQ ) CHV = CHEVAL ( NTERM3 , ARRFD3 , T ) FDP1P5 = FDP1P5 * ( ONE + CHV / XSQ ) ENDIF ENDIF ENDIF RETURN END DOUBLE PRECISION FUNCTION FDP2P5(XVALUE) C C DESCRIPTION: C C This function computes the Fermi-Dirac function of C order 5/2, defined as C C Int{0 to inf} t**(5/2) / (1+exp(t-x)) dt C FDP2P5(x) = ----------------------------------------- C Gamma(7/2) C C The function uses Chebyshev expansions which are given to C 16 decimal places for x <= 2, but only 10 decimal places C for x > 2. C C C ERROR RETURNS: C C If XVALUE too large and positive, the function value C will overflow. An error message is printed and the function C returns the value 0.0. C C C MACHINE-DEPENDENT CONSTANTS: C C NTERMS1 - INTEGER - The number of terms used from the array C ARRFD1. The recommended value is such that C ABS(ARRFD1(NTERMS1)) < EPS/10 C subject to 1 <= NTERMS1 <= 11. C C NTERMS2 - INTEGER - The number of terms used from the array C ARRFD2. The recommended value is such that C ABS(ARRFD2(NTERMS2)) < EPS/10 C subject to 1 <= NTERMS1 <= 21. C C NTERMS3 - INTEGER - The number of terms used from the array C ARRFD3. The recommended value is such that C ABS(ARRFD3(NTERMS3)) < EPS/10 C subject to 1 <= NTERMS3 <= 39. C C XMIN1 - REAL - The value of x below which C FDP2P5(x) = exp(x) C to machine precision. The recommended value C is 3.5*LN(2) + LN(EPSNEG) C C XMIN2 - REAL - The value of x below which C FDP2P5(x) = 0.0 C to machine precision. The recommended value C is LN ( XMIN ) C C XHIGH1 - REAL - The value of x above which C FDP2P5(x) = x**(7/2)/GAMMA(9/2) C to machine precision. The recommended value C is pi * SQRT(35/(12*EPS)) C C XHIGH2 - REAL - The value of x above which FDP2P5 would C overflow. The reommended value is C (11.6317*XMAX)**(2/7) C C For values of EPS, EPSNEG, and XMIN the user should refer to the C paper by Cody in ACM. Trans. Math. Soft. Vol. 14 (1988) p303-311. C C This code is provided with single and double precision values C of the machine-dependent parameters, suitable for machines C which satisfy the IEEE floating-point standard. C C C AUTHOR: C DR. ALLAN MACLEOD, C DEPT. OF MATHEMATICS AND STATISTICS, C UNIVERSITY OF PAISLEY, C HIGH ST., C PAISLEY, C SCOTLAND C PA1 2BE C C (e-mail: macl-ms0@paisley.ac.uk ) C C C LATEST UPDATE: C 21 NOVEMBER, 1996 C INTEGER NTERM1,NTERM2,NTERM3 DOUBLE PRECISION 1 ARRFD1(0:11),ARRFD2(0:21),ARRFD3(0:61), 2 CHEVAL,CHV,EXPX,FIFTY,FORTY2, 3 GAM4P5,ONE,T,THREE,TWO,TWOE,X,XHIGH1, 4 XHIGH2,XMIN1,XMIN2,XSQ,XVALUE,ZERO DATA ARRFD1/1.9694 4166 8589 6693 D 0, 1 -0.1496 9179 4643 492 D -1, 2 0.3006 9558 1662 7 D -3, 3 -0.8946 2485 950 D -5, 4 0.3298 0720 25 D -6, 5 -0.1392 3929 8 D -7, 6 0.6455 885 D -9, 7 -0.3206 23 D -10, 8 0.1678 3 D -11, 9 -0.916 D -13, X 0.52 D -14, 1 -0.3 D -15/ DATA ARRFD2( 0)/ 4.2642 8383 9865 5301 D 0/ DATA ARRFD2( 1)/ 2.3437 4268 8491 2867 D 0/ DATA ARRFD2( 2)/ 0.6727 1197 8005 2076 D 0/ DATA ARRFD2( 3)/ 0.1148 8263 2796 5569 D 0/ DATA ARRFD2( 4)/ 0.1093 6396 8046 758 D -1/ DATA ARRFD2( 5)/ 0.2567 1739 5701 5 D -3/ DATA ARRFD2( 6)/-0.5058 8998 3911 D -4/ DATA ARRFD2( 7)/-0.7376 2157 74 D -6/ DATA ARRFD2( 8)/ 0.7352 9987 58 D -6/ DATA ARRFD2( 9)/-0.1664 2173 6 D -7/ DATA ARRFD2(10)/-0.1409 2049 9 D -7/ DATA ARRFD2(11)/ 0.9949 192 D -9/ DATA ARRFD2(12)/ 0.2991 457 D -9/ DATA ARRFD2(13)/-0.4013 32 D -10/ DATA ARRFD2(14)/-0.6354 6 D -11/ DATA ARRFD2(15)/ 0.1479 3 D -11/ DATA ARRFD2(16)/ 0.1181 D -12/ DATA ARRFD2(17)/-0.524 D -13/ DATA ARRFD2(18)/-0.11 D -14/ DATA ARRFD2(19)/ 0.18 D -14/ DATA ARRFD2(20)/-0.1 D -15/ DATA ARRFD2(21)/-0.1 D -15/ DATA ARRFD3( 0)/30.2895 6768 5980 2579 D 0/ DATA ARRFD3( 1)/ 1.1678 9766 4206 0562 D 0/ DATA ARRFD3( 2)/ 0.6420 5918 0082 1472 D 0/ DATA ARRFD3( 3)/ 0.3461 7238 6840 7417 D 0/ DATA ARRFD3( 4)/ 0.1840 8167 9078 1889 D 0/ DATA ARRFD3( 5)/ 0.9730 9243 5354 509 D -1/ DATA ARRFD3( 6)/ 0.5139 7329 2675 393 D -1/ DATA ARRFD3( 7)/ 0.2717 0980 1041 757 D -1/ DATA ARRFD3( 8)/ 0.1438 3327 1401 165 D -1/ DATA ARRFD3( 9)/ 0.7626 4863 9521 55 D -2/ DATA ARRFD3(10)/ 0.4050 3695 7672 02 D -2/ DATA ARRFD3(11)/ 0.2154 3961 4641 49 D -2/ DATA ARRFD3(12)/ 0.1147 5689 9017 77 D -2/ DATA ARRFD3(13)/ 0.6120 6223 6928 2 D -3/ DATA ARRFD3(14)/ 0.3268 3403 3785 9 D -3/ DATA ARRFD3(15)/ 0.1747 1455 2274 2 D -3/ DATA ARRFD3(16)/ 0.9348 7845 7860 D -4/ DATA ARRFD3(17)/ 0.5006 9221 2553 D -4/ DATA ARRFD3(18)/ 0.2683 7382 1846 D -4/ DATA ARRFD3(19)/ 0.1439 5719 1251 D -4/ DATA ARRFD3(20)/ 0.7727 2440 700 D -5/ DATA ARRFD3(21)/ 0.4150 3820 336 D -5/ DATA ARRFD3(22)/ 0.2230 5118 261 D -5/ DATA ARRFD3(23)/ 0.1199 3697 093 D -5/ DATA ARRFD3(24)/ 0.6452 3443 69 D -6/ DATA ARRFD3(25)/ 0.3472 8228 81 D -6/ DATA ARRFD3(26)/ 0.1869 9642 15 D -6/ DATA ARRFD3(27)/ 0.1007 3002 72 D -6/ DATA ARRFD3(28)/ 0.5428 0756 1 D -7/ DATA ARRFD3(29)/ 0.2926 0782 9 D -7/ DATA ARRFD3(30)/ 0.1577 8591 8 D -7/ DATA ARRFD3(31)/ 0.8511 0768 D -8/ DATA ARRFD3(32)/ 0.4592 2760 D -8/ DATA ARRFD3(33)/ 0.2478 5001 D -8/ DATA ARRFD3(34)/ 0.1338 0255 D -8/ DATA ARRFD3(35)/ 0.7225 103 D -9/ DATA ARRFD3(36)/ 0.3902 350 D -9/ DATA ARRFD3(37)/ 0.2108 157 D -9/ DATA ARRFD3(38)/ 0.1139 122 D -9/ DATA ARRFD3(39)/ 0.6156 38 D -10/ DATA ARRFD3(40)/ 0.3327 81 D -10/ DATA ARRFD3(41)/ 0.1799 19 D -10/ DATA ARRFD3(42)/ 0.9728 8 D -11/ DATA ARRFD3(43)/ 0.5261 7 D -11/ DATA ARRFD3(44)/ 0.2846 1 D -11/ DATA ARRFD3(45)/ 0.1539 7 D -11/ DATA ARRFD3(46)/ 0.8331 D -12/ DATA ARRFD3(47)/ 0.4508 D -12/ DATA ARRFD3(48)/ 0.2440 D -12/ DATA ARRFD3(49)/ 0.1321 D -12/ DATA ARRFD3(50)/ 0.715 D -13/ DATA ARRFD3(51)/ 0.387 D -13/ DATA ARRFD3(52)/ 0.210 D -13/ DATA ARRFD3(53)/ 0.114 D -13/ DATA ARRFD3(54)/ 0.61 D -14/ DATA ARRFD3(55)/ 0.33 D -14/ DATA ARRFD3(56)/ 0.18 D -14/ DATA ARRFD3(57)/ 0.11 D -14/ DATA ARRFD3(58)/ 0.5 D -15/ DATA ARRFD3(59)/ 0.3 D -15/ DATA ARRFD3(60)/ 0.2 D -15/ DATA ARRFD3(61)/ 0.1 D -15/ DATA ZERO,ONE,TWO/ 0.0 D 0 , 1.0 D 0 , 2.0 D 0/ DATA THREE,FORTY2,FIFTY/ 3.0 D 0 , 42.0 D 0 , 50.0 D 0/ DATA GAM4P5/0.1163 1728 3965 6745 D 2/ DATA TWOE/5.4365 6365 6918 0905 D 0/ C C Machine-dependent constants (suitable for IEEE machines) C DATA NTERM1,NTERM2,NTERM3/11,21,61/ DATA XMIN1,XMIN2/-34.3107854D0,-708.396418D0/ DATA XHIGH1,XHIGH2/254599860.5D0,2.383665D88/ C C Start calculation C X=XVALUE C C Test for error condition C IF ( X .GT. XHIGH2 ) THEN PRINT*,'** ERROR ** - X TOO LARGE FOR FDP2P5' FDP2P5 = ZERO RETURN ENDIF C C Code for x < -1 C IF ( X .LT. -ONE ) THEN IF ( X .GT. XMIN1 ) THEN EXPX = EXP(X) T = TWOE * EXPX - ONE FDP2P5 = EXPX * CHEVAL ( NTERM1 , ARRFD1 , T ) ELSE IF ( X .LT. XMIN2 ) THEN FDP2P5 = ZERO ELSE FDP2P5 = EXP(X) ENDIF ENDIF ELSE C C Code for -1 <= x <= 2 C IF ( X .LE. TWO ) THEN T = ( TWO * X - ONE ) / THREE FDP2P5 = CHEVAL ( NTERM2 , ARRFD2 , T ) ELSE C C Code for x > 2 C FDP2P5 = X * X * X * SQRT(X) / GAM4P5 IF ( X .LE. XHIGH1 ) THEN XSQ = X * X T = ( FIFTY - XSQ ) / ( FORTY2 + XSQ ) CHV = CHEVAL ( NTERM3 , ARRFD3 , T ) FDP2P5 = FDP2P5 * ( ONE + CHV / XSQ ) ENDIF ENDIF ENDIF RETURN END DOUBLE PRECISION FUNCTION CHEVAL(N,A,T) C C DESCRIPTION: C C This function evaluates a Chebyshev series, using the C Clenshaw method with Reinsch modification, as analysed C in the paper by Oliver. C C C INPUT PARAMETERS C C N - INTEGER - The no. of terms in the sequence C C A - REAL ARRAY, dimension 0 to N - The coefficients of C the Chebyshev series C C T - REAL - The value at which the series is to be C evaluated C C C REFERENCES C C "An error analysis of the modified Clenshaw method for C evaluating Chebyshev and Fourier series" J. Oliver, C J.I.M.A., vol. 20, 1977, pp379-391 C C C MACHINE-DEPENDENT CONSTANTS: NONE C C C INTRINSIC FUNCTIONS USED; C C ABS C C C AUTHOR: Dr. Allan J. MacLeod, C Dept. of Mathematics and Statistics, C University of Paisley , C High St., C PAISLEY, C SCOTLAND C ( e-mail: macl-ms0@paisley.ac.uk ) C C C LATEST MODIFICATION: C 21 September , 1995 C C INTEGER I,N DOUBLE PRECISION 1 A(0:N),D1,D2,HALF,T,TEST,TT,TWO,U0,U1,U2,ZERO DATA ZERO,HALF/ 0.0 D 0 , 0.5 D 0 / DATA TEST,TWO/ 0.6 D 0 , 2.0 D 0 / U1 = ZERO C C If ABS ( T ) < 0.6 use the standard Clenshaw method C IF ( ABS( T ) .LT. TEST ) THEN U0 = ZERO TT = T + T DO 100 I = N , 0 , -1 U2 = U1 U1 = U0 U0 = TT * U1 + A( I ) - U2 100 CONTINUE CHEVAL = ( U0 - U2 ) / TWO ELSE C C If ABS ( T ) > = 0.6 use the Reinsch modification C D1 = ZERO C C T > = 0.6 code C IF ( T .GT. ZERO ) THEN TT = ( T - HALF ) - HALF TT = TT + TT DO 200 I = N , 0 , -1 D2 = D1 U2 = U1 D1 = TT * U2 + A( I ) + D2 U1 = D1 + U2 200 CONTINUE CHEVAL = ( D1 + D2 ) / TWO ELSE C C T < = -0.6 code C TT = ( T + HALF ) + HALF TT = TT + TT DO 300 I = N , 0 , -1 D2 = D1 U2 = U1 D1 = TT * U2 + A( I ) - D2 U1 = D1 - U2 300 CONTINUE CHEVAL = ( D1 - D2 ) / TWO ENDIF ENDIF RETURN END SHAR_EOF fi # end of overwriting check cd .. if test ! -d 'Sp' then mkdir 'Sp' fi cd 'Sp' if test -f 'src.f' then echo shar: will not over-write existing file "'src.f'" else cat << \SHAR_EOF > 'src.f' REAL FUNCTION FDM0P5(XVALUE) C C DESCRIPTION: C C This function computes the Fermi-Dirac function of C order -1/2, defined as C C Int{0 to inf} t**(-1/2) / (1+exp(t-x)) dt C FDM0P5(x) = ----------------------------------------- C Gamma(1/2) C C The function uses Chebyshev expansions which are given to C 16 decimal places for x <= 2, but only 10 decimal places C for x > 2. C C C ERROR RETURNS: C C None. C C C MACHINE-DEPENDENT CONSTANTS: C C NTERMS1 - INTEGER - The number of terms used from the array C ARRFD1. The recommended value is such that C ABS(ARRFD1(NTERMS1)) < EPS/10 C subject to 1 <= NTERMS1 <= 14. C C NTERMS2 - INTEGER - The number of terms used from the array C ARRFD2. The recommended value is such that C ABS(ARRFD2(NTERMS2)) < EPS/10 C subject to 1 <= NTERMS1 <= 23. C C NTERMS3 - INTEGER - The number of terms used from the array C ARRFD3. The recommended value is such that C ABS(ARRFD3(NTERMS3)) < EPS/10 C subject to 1 <= NTERMS3 <= 28. C C XMIN1 - REAL - The value of x below which C FDM0P5(x) = exp(x) C to machine precision. The recommended value C is LN ( SQRT(2) * EPSNEG ) C C XMIN2 - REAL - The value of x below which C FDM0P5(x) = 0.0 C to machine precision. The recommended value C is LN ( XMIN ) C C XHIGH - REAL - The value of x above which C FDM0P5(x) = 2 sqrt (x/pi) C to machine precision. The recommended value C is 1 / sqrt( 2 * EPSNEG ) C C For values of EPS, EPSNEG, and XMIN the user should refer to the C paper by Cody in ACM. Trans. Math. Soft. Vol. 14 (1988) p303-311. C C This code is provided with single and double precision values C of the machine-dependent parameters, suitable for machines C which satisfy the IEEE floating-point standard. C C C AUTHOR: C DR. ALLAN MACLEOD, C DEPT. OF MATHEMATICS AND STATISTICS, C UNIVERSITY OF PAISLEY, C HIGH ST., C PAISLEY, C SCOTLAND C PA1 2BE C C (e-mail: macl-ms0@paisley.ac.uk ) C C C LATEST UPDATE: C 20 NOVEMBER, 1996 C INTEGER NTERM1,NTERM2,NTERM3 REAL 1 ARRFD1(0:14),ARRFD2(0:23),ARRFD3(0:58), 2 CHEVAL,CHV,EXPX,FIFTY,FORTY2, 3 GAM1P5,ONE,T,THREE,TWO,TWOE, 4 X,XHIGH,XMIN1,XMIN2,XSQ,XVALUE,ZERO DATA ARRFD1/1.7863 5963 8510 2264 E 0, 1 -0.9993 7200 7632 333 E -1, 2 0.6414 4652 2160 54 E -2, 3 -0.4356 4153 7134 5 E -3, 4 0.3052 1670 0310 E -4, 5 -0.2181 0648 110 E -5, 6 0.1580 0507 81 E -6, 7 -0.1156 2057 0 E -7, 8 0.8525 860 E -9, 9 -0.6325 29 E -10, X 0.4715 9 E -11, 1 -0.3530 E -12, 2 0.265 E -13, 3 -0.20 E -14, 4 0.2 E -15/ DATA ARRFD2( 0)/ 1.6877 1115 2605 2352 E 0/ DATA ARRFD2( 1)/ 0.5978 3602 2633 6983 E 0/ DATA ARRFD2( 2)/ 0.3572 2600 4541 669 E -1/ DATA ARRFD2( 3)/-0.1321 4478 6506 426 E -1/ DATA ARRFD2( 4)/-0.4040 1342 0744 7 E -3/ DATA ARRFD2( 5)/ 0.5330 0118 4688 7 E -3/ DATA ARRFD2( 6)/-0.1489 2350 4863 E -4/ DATA ARRFD2( 7)/-0.2188 6382 2916 E -4/ DATA ARRFD2( 8)/ 0.1965 2084 277 E -5/ DATA ARRFD2( 9)/ 0.8565 8304 66 E -6/ DATA ARRFD2(10)/-0.1407 7231 33 E -6/ DATA ARRFD2(11)/-0.3051 7580 3 E -7/ DATA ARRFD2(12)/ 0.8352 4532 E -8/ DATA ARRFD2(13)/ 0.9025 750 E -9/ DATA ARRFD2(14)/-0.4455 471 E -9/ DATA ARRFD2(15)/-0.1483 42 E -10/ DATA ARRFD2(16)/ 0.2192 66 E -10/ DATA ARRFD2(17)/-0.6579 E -12/ DATA ARRFD2(18)/-0.1000 9 E -11/ DATA ARRFD2(19)/ 0.936 E -13/ DATA ARRFD2(20)/ 0.420 E -13/ DATA ARRFD2(21)/-0.71 E -14/ DATA ARRFD2(22)/-0.16 E -14/ DATA ARRFD2(23)/ 0.4 E -15/ DATA ARRFD3(0)/ 0.8707 1950 2959 0563 E 0/ DATA ARRFD3(1)/ 0.5983 3110 2317 33 E -2/ DATA ARRFD3(2)/ -0.4326 7047 0895 746 E -1/ DATA ARRFD3(3)/ -0.3930 8368 1608 590 E -1/ DATA ARRFD3(4)/ -0.1914 8268 8045 932 E -1/ DATA ARRFD3(5)/ -0.6558 2880 9801 58 E -2/ DATA ARRFD3(6)/ -0.2227 6691 5163 12 E -2/ DATA ARRFD3(7)/ -0.8466 7869 3617 8 E -3/ DATA ARRFD3(8)/ -0.2807 4594 8921 9 E -3/ DATA ARRFD3(9)/ -0.9555 7502 4348 E -4/ DATA ARRFD3(10)/-0.3623 6766 2803 E -4/ DATA ARRFD3(11)/-0.1091 5846 8869 E -4/ DATA ARRFD3(12)/-0.3935 6701 000 E -5/ DATA ARRFD3(13)/-0.1310 8192 725 E -5/ DATA ARRFD3(14)/-0.2468 8163 88 E -6/ DATA ARRFD3(15)/-0.1048 3803 11 E -6/ DATA ARRFD3(16)/ 0.2361 8148 7 E -7/ DATA ARRFD3(17)/ 0.2271 4535 9 E -7/ DATA ARRFD3(18)/ 0.1457 7517 4 E -7/ DATA ARRFD3(19)/ 0.1539 2676 7 E -7/ DATA ARRFD3(20)/ 0.5692 4772 E -8/ DATA ARRFD3(21)/ 0.5062 3068 E -8/ DATA ARRFD3(22)/ 0.2342 6075 E -8/ DATA ARRFD3(23)/ 0.1265 2275 E -8/ DATA ARRFD3(24)/ 0.8927 773 E -9/ DATA ARRFD3(25)/ 0.2994 501 E -9/ DATA ARRFD3(26)/ 0.2822 785 E -9/ DATA ARRFD3(27)/ 0.9106 85 E -10/ DATA ARRFD3(28)/ 0.6962 85 E -10/ DATA ARRFD3(29)/ 0.3662 25 E -10/ DATA ARRFD3(30)/ 0.1243 51 E -10/ DATA ARRFD3(31)/ 0.1450 19 E -10/ DATA ARRFD3(32)/ 0.1664 5 E -11/ DATA ARRFD3(33)/ 0.4585 6 E -11/ DATA ARRFD3(34)/ 0.6092 E -12/ DATA ARRFD3(35)/ 0.9331 E -12/ DATA ARRFD3(36)/ 0.5238 E -12/ DATA ARRFD3(37)/-0.56 E -14/ DATA ARRFD3(38)/ 0.3170 E -12/ DATA ARRFD3(39)/-0.926 E -13/ DATA ARRFD3(40)/ 0.1265 E -12/ DATA ARRFD3(41)/-0.327 E -13/ DATA ARRFD3(42)/ 0.276 E -13/ DATA ARRFD3(43)/ 0.33 E -14/ DATA ARRFD3(44)/-0.42 E -14/ DATA ARRFD3(45)/ 0.101 E -13/ DATA ARRFD3(46)/-0.73 E -14/ DATA ARRFD3(47)/ 0.64 E -14/ DATA ARRFD3(48)/-0.37 E -14/ DATA ARRFD3(49)/ 0.23 E -14/ DATA ARRFD3(50)/-0.9 E -15/ DATA ARRFD3(51)/ 0.2 E -15/ DATA ARRFD3(52)/ 0.2 E -15/ DATA ARRFD3(53)/-0.3 E -15/ DATA ARRFD3(54)/ 0.4 E -15/ DATA ARRFD3(55)/-0.3 E -15/ DATA ARRFD3(56)/ 0.2 E -15/ DATA ARRFD3(57)/-0.1 E -15/ DATA ARRFD3(58)/ 0.1 E -15/ DATA ZERO,ONE,TWO/ 0.0 E 0 , 1.0 E 0 , 2.0 E 0/ DATA THREE,FORTY2,FIFTY/ 3.0 E 0 , 42.0 E 0 , 50.0 E 0/ DATA GAM1P5/0.8862 2692 5452 7580 E 0/ DATA TWOE/5.4365 6365 6918 0905 E 0/ C C Machine-dependent constants C DATA NTERM1,NTERM2,NTERM3/8,13,23/ DATA XMIN1,XMIN2,XHIGH/-16.2890E0,-87.3365E0,2896.309E0/ C C Start calculation C X=XVALUE C C Code for x < -1 C IF ( X .LT. -ONE ) THEN IF ( X .GT. XMIN1 ) THEN EXPX = EXP(X) T = TWOE * EXPX - ONE FDM0P5 = EXPX * CHEVAL ( NTERM1 , ARRFD1 , T ) ELSE IF ( X .LT. XMIN2 ) THEN FDM0P5 = ZERO ELSE FDM0P5 = EXP(X) ENDIF ENDIF ELSE C C Code for -1 <= x <= 2 C IF ( X .LE. TWO ) THEN T = ( TWO * X - ONE ) / THREE FDM0P5 = CHEVAL ( NTERM2 , ARRFD2 , T ) ELSE C C Code for x > 2 C FDM0P5 = SQRT(X) / GAM1P5 IF ( X .LE. XHIGH ) THEN XSQ = X * X T = ( FIFTY - XSQ ) / ( FORTY2 + XSQ ) CHV = CHEVAL ( NTERM3 , ARRFD3 , T ) FDM0P5 = FDM0P5 * ( ONE - CHV / XSQ ) ENDIF ENDIF ENDIF RETURN END REAL FUNCTION FDP0P5(XVALUE) C C DESCRIPTION: C C This function computes the Fermi-Dirac function of C order 1/2, defined as C C Int{0 to inf} t**(1/2) / (1+exp(t-x)) dt C FDP0P5(x) = ----------------------------------------- C Gamma(3/2) C C The function uses Chebyshev expansions which are given to C 16 decimal places for x <= 2, but only 10 decimal places C for x > 2. C C C ERROR RETURNS: C C If XVALUE too large and positive, the function value C will overflow. An error message is printed and the function C returns the value 0.0. C C C MACHINE-DEPENDENT CONSTANTS: C C NTERMS1 - INTEGER - The number of terms used from the array C ARRFD1. The recommended value is such that C ABS(ARRFD1(NTERMS1)) < EPS/10 C subject to 1 <= NTERMS1 <= 13. C C NTERMS2 - INTEGER - The number of terms used from the array C ARRFD2. The recommended value is such that C ABS(ARRFD2(NTERMS2)) < EPS/10 C subject to 1 <= NTERMS1 <= 23. C C NTERMS3 - INTEGER - The number of terms used from the array C ARRFD3. The recommended value is such that C ABS(ARRFD3(NTERMS3)) < EPS/10 C subject to 1 <= NTERMS3 <= 32. C C XMIN1 - REAL - The value of x below which C FDP0P5(x) = exp(x) C to machine precision. The recommended value C is 1.5*LN(2) + LN(EPSNEG) C C XMIN2 - REAL - The value of x below which C FDP0P5(x) = 0.0 C to machine precision. The recommended value C is LN ( XMIN ) C C XHIGH1 - REAL - The value of x above which C FDP0P5(x) = x**(3/2)/GAMMA(5/2) C to machine precision. The recommended value C is pi / SQRT(8*EPS) C C XHIGH2 - REAL - The value of x above which FDP0P5 would C overflow. The reommended value is C (1.329*XMAX)**(2/3) C C For values of EPS, EPSNEG, and XMIN the user should refer to the C paper by Cody in ACM. Trans. Math. Soft. Vol. 14 (1988) p303-311. C C This code is provided with single and double precision values C of the machine-dependent parameters, suitable for machines C which satisfy the IEEE floating-point standard. C C C C AUTHOR: C DR. ALLAN MACLEOD, C DEPT. OF MATHEMATICS AND STATISTICS, C UNIVERSITY OF PAISLEY, C HIGH ST., C PAISLEY, C SCOTLAND C PA1 2BE C C (e-mail: macl-ms0@paisley.ac.uk ) C C C LATEST UPDATE: C 20 NOVEMBER, 1996 C INTEGER NTERM1,NTERM2,NTERM3 REAL 1 ARRFD1(0:13),ARRFD2(0:23),ARRFD3(0:53), 2 CHEVAL,CHV,EXPX,FIFTY,FORTY2, 3 GAM2P5,ONE,T,THREE,TWO,TWOE,X,XHIGH1, 4 XHIGH2,XMIN1,XMIN2,XSQ,XVALUE,ZERO DATA ARRFD1/1.8862 9683 9273 4597 E 0, 1 -0.5435 8081 7644 053 E -1, 2 0.2364 4975 4397 20 E -2, 3 -0.1216 9293 6588 0 E -3, 4 0.6869 5130 622 E -5, 5 -0.4112 0761 72 E -6, 6 0.2563 5162 8 E -7, 7 -0.1646 5008 E -8, 8 0.1081 948 E -9, 9 -0.7239 2 E -11, X 0.4915 E -12, 1 -0.338 E -13, 2 0.23 E -14, 3 -0.2 E -15/ DATA ARRFD2( 0)/ 2.6982 4927 8817 0612 E 0/ DATA ARRFD2( 1)/ 1.2389 9141 4113 3012 E 0/ DATA ARRFD2( 2)/ 0.2291 4393 7981 6278 E 0/ DATA ARRFD2( 3)/ 0.9031 6534 6872 79 E -2/ DATA ARRFD2( 4)/-0.2577 6524 6912 46 E -2/ DATA ARRFD2( 5)/-0.5836 8160 5388 E -4/ DATA ARRFD2( 6)/ 0.6936 0945 8725 E -4/ DATA ARRFD2( 7)/-0.1806 1670 265 E -5/ DATA ARRFD2( 8)/-0.2132 1530 005 E -5/ DATA ARRFD2( 9)/ 0.1754 9839 51 E -6/ DATA ARRFD2(10)/ 0.6653 2547 0 E -7/ DATA ARRFD2(11)/-0.1016 7597 7 E -7/ DATA ARRFD2(12)/-0.1963 7597 E -8/ DATA ARRFD2(13)/ 0.5075 769 E -9/ DATA ARRFD2(14)/ 0.4914 69 E -10/ DATA ARRFD2(15)/-0.2337 37 E -10/ DATA ARRFD2(16)/-0.6645 E -12/ DATA ARRFD2(17)/ 0.1011 5 E -11/ DATA ARRFD2(18)/-0.313 E -13/ DATA ARRFD2(19)/-0.412 E -13/ DATA ARRFD2(20)/ 0.38 E -14/ DATA ARRFD2(21)/ 0.16 E -14/ DATA ARRFD2(22)/-0.3 E -15/ DATA ARRFD2(23)/-0.1 E -15/ DATA ARRFD3(0)/ 2.5484 3841 9800 9122 E 0/ DATA ARRFD3(1)/ 0.5104 3940 8960 652 E -1/ DATA ARRFD3(2)/ 0.7749 3527 6282 94 E -2/ DATA ARRFD3(3)/ -0.7504 1656 5849 53 E -2/ DATA ARRFD3(4)/ -0.7754 0826 3202 96 E -2/ DATA ARRFD3(5)/ -0.4581 0844 5399 77 E -2/ DATA ARRFD3(6)/ -0.2343 1641 5873 63 E -2/ DATA ARRFD3(7)/ -0.1178 8049 5135 91 E -2/ DATA ARRFD3(8)/ -0.5802 7393 5970 2 E -3/ DATA ARRFD3(9)/ -0.2825 3507 0053 7 E -3/ DATA ARRFD3(10)/-0.1388 1366 5179 9 E -3/ DATA ARRFD3(11)/-0.6806 9508 4875 E -4/ DATA ARRFD3(12)/-0.3353 5635 0608 E -4/ DATA ARRFD3(13)/-0.1665 3301 8734 E -4/ DATA ARRFD3(14)/-0.8271 4908 266 E -5/ DATA ARRFD3(15)/-0.4142 5714 409 E -5/ DATA ARRFD3(16)/-0.2080 5255 294 E -5/ DATA ARRFD3(17)/-0.1047 9767 478 E -5/ DATA ARRFD3(18)/-0.5315 2738 02 E -6/ DATA ARRFD3(19)/-0.2694 0611 78 E -6/ DATA ARRFD3(20)/-0.1374 8787 49 E -6/ DATA ARRFD3(21)/-0.7023 0888 7 E -7/ DATA ARRFD3(22)/-0.3595 4394 2 E -7/ DATA ARRFD3(23)/-0.1851 0612 6 E -7/ DATA ARRFD3(24)/-0.9502 3937 E -8/ DATA ARRFD3(25)/-0.4918 4811 E -8/ DATA ARRFD3(26)/-0.2537 1950 E -8/ DATA ARRFD3(27)/-0.1315 1532 E -8/ DATA ARRFD3(28)/-0.6835 168 E -9/ DATA ARRFD3(29)/-0.3538 244 E -9/ DATA ARRFD3(30)/-0.1853 182 E -9/ DATA ARRFD3(31)/-0.9589 83 E -10/ DATA ARRFD3(32)/-0.5040 83 E -10/ DATA ARRFD3(33)/-0.2622 38 E -10/ DATA ARRFD3(34)/-0.1372 55 E -10/ DATA ARRFD3(35)/-0.7234 0 E -11/ DATA ARRFD3(36)/-0.3742 9 E -11/ DATA ARRFD3(37)/-0.2005 9 E -11/ DATA ARRFD3(38)/-0.1026 9 E -11/ DATA ARRFD3(39)/-0.5551 E -12/ DATA ARRFD3(40)/-0.2857 E -12/ DATA ARRFD3(41)/-0.1520 E -12/ DATA ARRFD3(42)/-0.811 E -13/ DATA ARRFD3(43)/-0.410 E -13/ DATA ARRFD3(44)/-0.234 E -13/ DATA ARRFD3(45)/-0.110 E -13/ DATA ARRFD3(46)/-0.67 E -14/ DATA ARRFD3(47)/-0.30 E -14/ DATA ARRFD3(48)/-0.19 E -14/ DATA ARRFD3(49)/-0.9 E -15/ DATA ARRFD3(50)/-0.5 E -15/ DATA ARRFD3(51)/-0.3 E -15/ DATA ARRFD3(52)/-0.1 E -15/ DATA ARRFD3(53)/-0.1 E -15/ DATA ZERO,ONE,TWO/ 0.0 E 0 , 1.0 E 0 , 2.0 E 0/ DATA THREE,FORTY2,FIFTY/ 3.0 E 0 , 42.0 E 0 , 50.0 E 0/ DATA GAM2P5/0.1329 3403 8817 9137 E 1/ DATA TWOE/5.4365 6365 6918 0905 E 0/ C C Machine-dependent constants (suitable for IEEE machines) C DATA NTERM1,NTERM2,NTERM3/8,13,27/ DATA XMIN1,XMIN2/-15.6E0,-87.3365E0/ DATA XHIGH1,XHIGH2/3219.82E0,5.88847E25/ C C Start calculation C X=XVALUE C C Test for error condition C IF ( X .GT. XHIGH2 ) THEN PRINT*,'** ERROR ** - X TOO LARGE FOR FDP0P5' FDP0P5 = ZERO RETURN ENDIF C C Code for x < -1 C IF ( X .LT. -ONE ) THEN IF ( X .GT. XMIN1 ) THEN EXPX = EXP(X) T = TWOE * EXPX - ONE FDP0P5 = EXPX * CHEVAL ( NTERM1 , ARRFD1 , T ) ELSE IF ( X .LT. XMIN2 ) THEN FDP0P5 = ZERO ELSE FDP0P5 = EXP(X) ENDIF ENDIF ELSE C C Code for -1 <= x <= 2 C IF ( X .LE. TWO ) THEN T = ( TWO * X - ONE ) / THREE FDP0P5 = CHEVAL ( NTERM2 , ARRFD2 , T ) ELSE C C Code for x > 2 C FDP0P5 = X * SQRT(X) / GAM2P5 IF ( X .LE. XHIGH1 ) THEN XSQ = X * X T = ( FIFTY - XSQ ) / ( FORTY2 + XSQ ) CHV = CHEVAL ( NTERM3 , ARRFD3 , T ) FDP0P5 = FDP0P5 * ( ONE + CHV / XSQ ) ENDIF ENDIF ENDIF RETURN END REAL FUNCTION FDP1P5(XVALUE) C C DESCRIPTION: C C This function computes the Fermi-Dirac function of C order 3/2, defined as C C Int{0 to inf} t**(3/2) / (1+exp(t-x)) dt C FDP1P5(x) = ----------------------------------------- C Gamma(5/2) C C The function uses Chebyshev expansions which are given to C 16 decimal places for x <= 2, but only 10 decimal places C for x > 2. C C C ERROR RETURNS: C C If XVALUE too large and positive, the function value C will overflow. An error message is printed and the function C returns the value 0.0. C C C MACHINE-DEPENDENT CONSTANTS: C C NTERMS1 - INTEGER - The number of terms used from the array C ARRFD1. The recommended value is such that C ABS(ARRFD1(NTERMS1)) < EPS/10 C subject to 1 <= NTERMS1 <= 12. C C NTERMS2 - INTEGER - The number of terms used from the array C ARRFD2. The recommended value is such that C ABS(ARRFD2(NTERMS2)) < EPS/10 C subject to 1 <= NTERMS1 <= 22. C C NTERMS3 - INTEGER - The number of terms used from the array C ARRFD3. The recommended value is such that C ABS(ARRFD3(NTERMS3)) < EPS/10 C subject to 1 <= NTERMS3 <= 33. C C XMIN1 - REAL - The value of x below which C FDP1P5(x) = exp(x) C to machine precision. The recommended value C is 2.5*LN(2) + LN(EPSNEG) C C XMIN2 - REAL - The value of x below which C FDP1P5(x) = 0.0 C to machine precision. The recommended value C is LN ( XMIN ) C C XHIGH1 - REAL - The value of x above which C FDP1P5(x) = x**(5/2)/GAMMA(7/2) C to machine precision. The recommended value C is pi * SQRT(1.6/EPS) C C XHIGH2 - REAL - The value of x above which FDP1P5 would C overflow. The reommended value is C (3.233509*XMAX)**(2/5) C C For values of EPS, EPSNEG, and XMIN the user should refer to the C paper by Cody in ACM. Trans. Math. Soft. Vol. 14 (1988) p303-311. C C This code is provided with single and double precision values C of the machine-dependent parameters, suitable for machines C which satisfy the IEEE floating-point standard. C C C C AUTHOR: C DR. ALLAN MACLEOD, C DEPT. OF MATHEMATICS AND STATISTICS, C UNIVERSITY OF PAISLEY, C HIGH ST., C PAISLEY, C SCOTLAND C PA1 2BE C C (e-mail: macl_ms0@paisley.ac.uk ) C C C LATEST UPDATE: C 21 NOVEMBER, 1996 C INTEGER NTERM1,NTERM2,NTERM3 REAL 1 ARRFD1(0:12),ARRFD2(0:22),ARRFD3(0:55), 2 CHEVAL,CHV,EXPX,FIFTY,FORTY2, 3 GAM3P5,ONE,T,THREE,TWO,TWOE,X,XHIGH1, 4 XHIGH2,XMIN1,XMIN2,XSQ,XVALUE,ZERO DATA ARRFD1/1.9406 5492 1037 8650 E 0, 1 -0.2878 6747 5518 043 E -1, 2 0.8509 1579 5231 3 E -3, 3 -0.3327 8452 5669 E -4, 4 0.1517 1202 058 E -5, 5 -0.7622 0087 4 E -7, 6 0.4095 5489 E -8, 7 -0.2311 964 E -9, 8 0.1355 37 E -10, 9 -0.8187 E -12, X 0.507 E -13, 1 -0.32 E -14, 2 0.2 E -15/ DATA ARRFD2( 0)/ 3.5862 2516 1563 4306 E 0/ DATA ARRFD2( 1)/ 1.8518 2900 5626 5751 E 0/ DATA ARRFD2( 2)/ 0.4612 3491 0241 7150 E 0/ DATA ARRFD2( 3)/ 0.5793 0397 6126 881 E -1/ DATA ARRFD2( 4)/ 0.1704 3790 5548 75 E -2/ DATA ARRFD2( 5)/-0.3970 5201 2249 6 E -3/ DATA ARRFD2( 6)/-0.7070 2491 890 E -5/ DATA ARRFD2( 7)/ 0.7659 9748 792 E -5/ DATA ARRFD2( 8)/-0.1857 8113 33 E -6/ DATA ARRFD2( 9)/-0.1832 2379 56 E -6/ DATA ARRFD2(10)/ 0.1392 4949 5 E -7/ DATA ARRFD2(11)/ 0.4670 2027 E -8/ DATA ARRFD2(12)/-0.6671 984 E -9/ DATA ARRFD2(13)/-0.1161 292 E -9/ DATA ARRFD2(14)/ 0.2844 38 E -10/ DATA ARRFD2(15)/ 0.2490 6 E -11/ DATA ARRFD2(16)/-0.1143 1 E -11/ DATA ARRFD2(17)/-0.279 E -13/ DATA ARRFD2(18)/ 0.439 E -13/ DATA ARRFD2(19)/-0.14 E -14/ DATA ARRFD2(20)/-0.16 E -14/ DATA ARRFD2(21)/ 0.1 E -15/ DATA ARRFD2(22)/ 0.1 E -15/ DATA ARRFD3( 0)/12.1307 5817 3688 4627 E 0/ DATA ARRFD3( 1)/-0.1547 5011 1128 7255 E 0/ DATA ARRFD3( 2)/-0.7390 0738 8850 999 E -1/ DATA ARRFD3( 3)/-0.3072 3537 7959 258 E -1/ DATA ARRFD3( 4)/-0.1145 4857 9330 328 E -1/ DATA ARRFD3( 5)/-0.4056 7636 8095 39 E -2/ DATA ARRFD3( 6)/-0.1398 0158 3732 27 E -2/ DATA ARRFD3( 7)/-0.4454 9018 1015 3 E -3/ DATA ARRFD3( 8)/-0.1173 9461 1270 4 E -3/ DATA ARRFD3( 9)/-0.1484 0898 0093 E -4/ DATA ARRFD3(10)/ 0.1188 9515 4223 E -4/ DATA ARRFD3(11)/ 0.1464 7695 8178 E -4/ DATA ARRFD3(12)/ 0.1132 2874 1730 E -4/ DATA ARRFD3(13)/ 0.7576 2292 948 E -5/ DATA ARRFD3(14)/ 0.4712 0400 466 E -5/ DATA ARRFD3(15)/ 0.2813 2628 202 E -5/ DATA ARRFD3(16)/ 0.1637 0517 341 E -5/ DATA ARRFD3(17)/ 0.9351 0762 72 E -6/ DATA ARRFD3(18)/ 0.5278 6892 10 E -6/ DATA ARRFD3(19)/ 0.2951 0798 70 E -6/ DATA ARRFD3(20)/ 0.1638 6001 90 E -6/ DATA ARRFD3(21)/ 0.9052 0540 9 E -7/ DATA ARRFD3(22)/ 0.4977 5697 5 E -7/ DATA ARRFD3(23)/ 0.2729 5586 3 E -7/ DATA ARRFD3(24)/ 0.1492 1458 5 E -7/ DATA ARRFD3(25)/ 0.8142 0359 E -8/ DATA ARRFD3(26)/ 0.4434 9200 E -8/ DATA ARRFD3(27)/ 0.2411 6032 E -8/ DATA ARRFD3(28)/ 0.1310 5018 E -8/ DATA ARRFD3(29)/ 0.7109 736 E -9/ DATA ARRFD3(30)/ 0.3856 721 E -9/ DATA ARRFD3(31)/ 0.2089 529 E -9/ DATA ARRFD3(32)/ 0.1131 735 E -9/ DATA ARRFD3(33)/ 0.6127 85 E -10/ DATA ARRFD3(34)/ 0.3314 48 E -10/ DATA ARRFD3(35)/ 0.1794 19 E -10/ DATA ARRFD3(36)/ 0.9695 3 E -11/ DATA ARRFD3(37)/ 0.5246 3 E -11/ DATA ARRFD3(38)/ 0.2834 3 E -11/ DATA ARRFD3(39)/ 0.1532 3 E -11/ DATA ARRFD3(40)/ 0.8284 E -12/ DATA ARRFD3(41)/ 0.4472 E -12/ DATA ARRFD3(42)/ 0.2421 E -12/ DATA ARRFD3(43)/ 0.1304 E -12/ DATA ARRFD3(44)/ 0.707 E -13/ DATA ARRFD3(45)/ 0.381 E -13/ DATA ARRFD3(46)/ 0.206 E -13/ DATA ARRFD3(47)/ 0.111 E -13/ DATA ARRFD3(48)/ 0.60 E -14/ DATA ARRFD3(49)/ 0.33 E -14/ DATA ARRFD3(50)/ 0.17 E -14/ DATA ARRFD3(51)/ 0.11 E -14/ DATA ARRFD3(52)/ 0.5 E -15/ DATA ARRFD3(53)/ 0.3 E -15/ DATA ARRFD3(54)/ 0.1 E -15/ DATA ARRFD3(55)/ 0.1 E -15/ DATA ZERO,ONE,TWO/ 0.0 E 0 , 1.0 E 0 , 2.0 E 0/ DATA THREE,FORTY2,FIFTY/ 3.0 E 0 , 42.0 E 0 , 50.0 E 0/ DATA GAM3P5/0.3323 3509 7044 7843 E 1/ DATA TWOE/5.4365 6365 6918 0905 E 0/ C C Machine-dependent constants (suitable for IEEE machines) C DATA NTERM1,NTERM2,NTERM3/7,13,28/ DATA XMIN1,XMIN2/-14.9027E0,-87.33654E0/ DATA XHIGH1,XHIGH2/7193.411E0,4.136236E15/ C C Start calculation C X=XVALUE C C Test for error condition C IF ( X .GT. XHIGH2 ) THEN PRINT*,'** ERROR ** - X TOO LARGE FOR FDP1P5' FDP1P5 = ZERO RETURN ENDIF C C Code for x < -1 C IF ( X .LT. -ONE ) THEN IF ( X .GT. XMIN1 ) THEN EXPX = EXP(X) T = TWOE * EXPX - ONE FDP1P5 = EXPX * CHEVAL ( NTERM1 , ARRFD1 , T ) ELSE IF ( X .LT. XMIN2 ) THEN FDP1P5 = ZERO ELSE FDP1P5 = EXP(X) ENDIF ENDIF ELSE C C Code for -1 <= x <= 2 C IF ( X .LE. TWO ) THEN T = ( TWO * X - ONE ) / THREE FDP1P5 = CHEVAL ( NTERM2 , ARRFD2 , T ) ELSE C C Code for x > 2 C FDP1P5 = X * X * SQRT(X) / GAM3P5 IF ( X .LE. XHIGH1 ) THEN XSQ = X * X T = ( FIFTY - XSQ ) / ( FORTY2 + XSQ ) CHV = CHEVAL ( NTERM3 , ARRFD3 , T ) FDP1P5 = FDP1P5 * ( ONE + CHV / XSQ ) ENDIF ENDIF ENDIF RETURN END REAL FUNCTION FDP2P5(XVALUE) C C DESCRIPTION: C C This function computes the Fermi-Dirac function of C order 5/2, defined as C C Int{0 to inf} t**(5/2) / (1+exp(t-x)) dt C FDP2P5(x) = ----------------------------------------- C Gamma(7/2) C C The function uses Chebyshev expansions which are given to C 16 decimal places for x <= 2, but only 10 decimal places C for x > 2. C C C ERROR RETURNS: C C If XVALUE too large and positive, the function value C will overflow. An error message is printed and the function C returns the value 0.0. C C C MACHINE-DEPENDENT CONSTANTS: C C NTERMS1 - INTEGER - The number of terms used from the array C ARRFD1. The recommended value is such that C ABS(ARRFD1(NTERMS1)) < EPS/10 C subject to 1 <= NTERMS1 <= 11. C C NTERMS2 - INTEGER - The number of terms used from the array C ARRFD2. The recommended value is such that C ABS(ARRFD2(NTERMS2)) < EPS/10 C subject to 1 <= NTERMS1 <= 21. C C NTERMS3 - INTEGER - The number of terms used from the array C ARRFD3. The recommended value is such that C ABS(ARRFD3(NTERMS3)) < EPS/10 C subject to 1 <= NTERMS3 <= 39. C C XMIN1 - REAL - The value of x below which C FDP2P5(x) = exp(x) C to machine precision. The recommended value C is 3.5*LN(2) + LN(EPSNEG) C C XMIN2 - REAL - The value of x below which C FDP2P5(x) = 0.0 C to machine precision. The recommended value C is LN ( XMIN ) C C XHIGH1 - REAL - The value of x above which C FDP2P5(x) = x**(7/2)/GAMMA(9/2) C to machine precision. The recommended value C is pi * SQRT(35/(12*EPS)) C C XHIGH2 - REAL - The value of x above which FDP2P5 would C overflow. The reommended value is C (11.6317*XMAX)**(2/7) C C For values of EPS, EPSNEG, and XMIN the user should refer to the C paper by Cody in ACM. Trans. Math. Soft. Vol. 14 (1988) p303-311. C C This code is provided with single and double precision values C of the machine-dependent parameters, suitable for machines C which satisfy the IEEE floating-point standard. C C C C AUTHOR: C DR. ALLAN MACLEOD, C DEPT. OF MATHEMATICS AND STATISTICS, C UNIVERSITY OF PAISLEY, C HIGH ST., C PAISLEY, C SCOTLAND C PA1 2BE C C (e-mail: macl-ms0@paisley.ac.uk ) C C C LATEST UPDATE: C 21 NOVEMBER, 1996 C INTEGER NTERM1,NTERM2,NTERM3 REAL 1 ARRFD1(0:11),ARRFD2(0:21),ARRFD3(0:61), 2 CHEVAL,CHV,EXPX,FIFTY,FORTY2, 3 GAM4P5,ONE,T,THREE,TWO,TWOE,X,XHIGH1, 4 XHIGH2,XMIN1,XMIN2,XSQ,XVALUE,ZERO DATA ARRFD1/1.9694 4166 8589 6693 E 0, 1 -0.1496 9179 4643 492 E -1, 2 0.3006 9558 1662 7 E -3, 3 -0.8946 2485 950 E -5, 4 0.3298 0720 25 E -6, 5 -0.1392 3929 8 E -7, 6 0.6455 885 E -9, 7 -0.3206 23 E -10, 8 0.1678 3 E -11, 9 -0.916 E -13, X 0.52 E -14, 1 -0.3 E -15/ DATA ARRFD2( 0)/ 4.2642 8383 9865 5301 E 0/ DATA ARRFD2( 1)/ 2.3437 4268 8491 2867 E 0/ DATA ARRFD2( 2)/ 0.6727 1197 8005 2076 E 0/ DATA ARRFD2( 3)/ 0.1148 8263 2796 5569 E 0/ DATA ARRFD2( 4)/ 0.1093 6396 8046 758 E -1/ DATA ARRFD2( 5)/ 0.2567 1739 5701 5 E -3/ DATA ARRFD2( 6)/-0.5058 8998 3911 E -4/ DATA ARRFD2( 7)/-0.7376 2157 74 E -6/ DATA ARRFD2( 8)/ 0.7352 9987 58 E -6/ DATA ARRFD2( 9)/-0.1664 2173 6 E -7/ DATA ARRFD2(10)/-0.1409 2049 9 E -7/ DATA ARRFD2(11)/ 0.9949 192 E -9/ DATA ARRFD2(12)/ 0.2991 457 E -9/ DATA ARRFD2(13)/-0.4013 32 E -10/ DATA ARRFD2(14)/-0.6354 6 E -11/ DATA ARRFD2(15)/ 0.1479 3 E -11/ DATA ARRFD2(16)/ 0.1181 E -12/ DATA ARRFD2(17)/-0.524 E -13/ DATA ARRFD2(18)/-0.11 E -14/ DATA ARRFD2(19)/ 0.18 E -14/ DATA ARRFD2(20)/-0.1 E -15/ DATA ARRFD2(21)/-0.1 E -15/ DATA ARRFD3( 0)/30.2895 6768 5980 2579 E 0/ DATA ARRFD3( 1)/ 1.1678 9766 4206 0562 E 0/ DATA ARRFD3( 2)/ 0.6420 5918 0082 1472 E 0/ DATA ARRFD3( 3)/ 0.3461 7238 6840 7417 E 0/ DATA ARRFD3( 4)/ 0.1840 8167 9078 1889 E 0/ DATA ARRFD3( 5)/ 0.9730 9243 5354 509 E -1/ DATA ARRFD3( 6)/ 0.5139 7329 2675 393 E -1/ DATA ARRFD3( 7)/ 0.2717 0980 1041 757 E -1/ DATA ARRFD3( 8)/ 0.1438 3327 1401 165 E -1/ DATA ARRFD3( 9)/ 0.7626 4863 9521 55 E -2/ DATA ARRFD3(10)/ 0.4050 3695 7672 02 E -2/ DATA ARRFD3(11)/ 0.2154 3961 4641 49 E -2/ DATA ARRFD3(12)/ 0.1147 5689 9017 77 E -2/ DATA ARRFD3(13)/ 0.6120 6223 6928 2 E -3/ DATA ARRFD3(14)/ 0.3268 3403 3785 9 E -3/ DATA ARRFD3(15)/ 0.1747 1455 2274 2 E -3/ DATA ARRFD3(16)/ 0.9348 7845 7860 E -4/ DATA ARRFD3(17)/ 0.5006 9221 2553 E -4/ DATA ARRFD3(18)/ 0.2683 7382 1846 E -4/ DATA ARRFD3(19)/ 0.1439 5719 1251 E -4/ DATA ARRFD3(20)/ 0.7727 2440 700 E -5/ DATA ARRFD3(21)/ 0.4150 3820 336 E -5/ DATA ARRFD3(22)/ 0.2230 5118 261 E -5/ DATA ARRFD3(23)/ 0.1199 3697 093 E -5/ DATA ARRFD3(24)/ 0.6452 3443 69 E -6/ DATA ARRFD3(25)/ 0.3472 8228 81 E -6/ DATA ARRFD3(26)/ 0.1869 9642 15 E -6/ DATA ARRFD3(27)/ 0.1007 3002 72 E -6/ DATA ARRFD3(28)/ 0.5428 0756 1 E -7/ DATA ARRFD3(29)/ 0.2926 0782 9 E -7/ DATA ARRFD3(30)/ 0.1577 8591 8 E -7/ DATA ARRFD3(31)/ 0.8511 0768 E -8/ DATA ARRFD3(32)/ 0.4592 2760 E -8/ DATA ARRFD3(33)/ 0.2478 5001 E -8/ DATA ARRFD3(34)/ 0.1338 0255 E -8/ DATA ARRFD3(35)/ 0.7225 103 E -9/ DATA ARRFD3(36)/ 0.3902 350 E -9/ DATA ARRFD3(37)/ 0.2108 157 E -9/ DATA ARRFD3(38)/ 0.1139 122 E -9/ DATA ARRFD3(39)/ 0.6156 38 E -10/ DATA ARRFD3(40)/ 0.3327 81 E -10/ DATA ARRFD3(41)/ 0.1799 19 E -10/ DATA ARRFD3(42)/ 0.9728 8 E -11/ DATA ARRFD3(43)/ 0.5261 7 E -11/ DATA ARRFD3(44)/ 0.2846 1 E -11/ DATA ARRFD3(45)/ 0.1539 7 E -11/ DATA ARRFD3(46)/ 0.8331 E -12/ DATA ARRFD3(47)/ 0.4508 E -12/ DATA ARRFD3(48)/ 0.2440 E -12/ DATA ARRFD3(49)/ 0.1321 E -12/ DATA ARRFD3(50)/ 0.715 E -13/ DATA ARRFD3(51)/ 0.387 E -13/ DATA ARRFD3(52)/ 0.210 E -13/ DATA ARRFD3(53)/ 0.114 E -13/ DATA ARRFD3(54)/ 0.61 E -14/ DATA ARRFD3(55)/ 0.33 E -14/ DATA ARRFD3(56)/ 0.18 E -14/ DATA ARRFD3(57)/ 0.11 E -14/ DATA ARRFD3(58)/ 0.5 E -15/ DATA ARRFD3(59)/ 0.3 E -15/ DATA ARRFD3(60)/ 0.2 E -15/ DATA ARRFD3(61)/ 0.1 E -15/ DATA ZERO,ONE,TWO/ 0.0 E 0 , 1.0 E 0 , 2.0 E 0/ DATA THREE,FORTY2,FIFTY/ 3.0 E 0 , 42.0 E 0 , 50.0 E 0/ DATA GAM4P5/0.1163 1728 3965 6745 E 2/ DATA TWOE/5.4365 6365 6918 0905 E 0/ C C Machine-dependent constants (suitable for IEEE machines) C DATA NTERM1,NTERM2,NTERM3/6,12,34/ DATA XMIN1,XMIN2/-14.20951E0,-87.33654E0/ DATA XHIGH1,XHIGH2/10988.12E0,2.058573E11/ C C Start calculation C X=XVALUE C C Test for error condition C IF ( X .GT. XHIGH2 ) THEN PRINT*,'** ERROR ** - X TOO LARGE FOR FDP2P5' FDP2P5 = ZERO RETURN ENDIF C C Code for x < -1 C IF ( X .LT. -ONE ) THEN IF ( X .GT. XMIN1 ) THEN EXPX = EXP(X) T = TWOE * EXPX - ONE FDP2P5 = EXPX * CHEVAL ( NTERM1 , ARRFD1 , T ) ELSE IF ( X .LT. XMIN2 ) THEN FDP2P5 = ZERO ELSE FDP2P5 = EXP(X) ENDIF ENDIF ELSE C C Code for -1 <= x <= 2 C IF ( X .LE. TWO ) THEN T = ( TWO * X - ONE ) / THREE FDP2P5 = CHEVAL ( NTERM2 , ARRFD2 , T ) ELSE C C Code for x > 2 C FDP2P5 = X * X * X * SQRT(X) / GAM4P5 IF ( X .LE. XHIGH1 ) THEN XSQ = X * X T = ( FIFTY - XSQ ) / ( FORTY2 + XSQ ) CHV = CHEVAL ( NTERM3 , ARRFD3 , T ) FDP2P5 = FDP2P5 * ( ONE + CHV / XSQ ) ENDIF ENDIF ENDIF RETURN END REAL FUNCTION CHEVAL(N,A,T) C C DESCRIPTION: C C This function evaluates a Chebyshev series, using the C Clenshaw method with Reinsch modification, as analysed C in the paper by Oliver. C C C INPUT PARAMETERS C C N - INTEGER - The no. of terms in the sequence C C A - REAL ARRAY, dimension 0 to N - The coefficients of C the Chebyshev series C C T - REAL - The value at which the series is to be C evaluated C C C REFERENCES C C "An error analysis of the modified Clenshaw method for C evaluating Chebyshev and Fourier series" J. Oliver, C J.I.M.A., vol. 20, 1977, pp379-391 C C C MACHINE-DEPENDENT CONSTANTS: NONE C C C INTRINSIC FUNCTIONS USED; C C ABS C C C C AUTHOR: Dr. Allan J. MacLeod, C Dept. of Mathematics and Statistics, C University of Paisley , C High St., C PAISLEY, C SCOTLAND C ( e-mail: macl-ms0@paisley.ac.uk ) C C C LATEST MODIFICATION: C 21 September , 1995 C C INTEGER I,N REAL 1 A(0:N),D1,D2,HALF,T,TEST,TT,TWO,U0,U1,U2,ZERO DATA ZERO,HALF/ 0.0 E 0 , 0.5 E 0 / DATA TEST,TWO/ 0.6 E 0 , 2.0 E 0 / U1 = ZERO C C If ABS ( T ) < 0.6 use the standard Clenshaw method C IF ( ABS( T ) .LT. TEST ) THEN U0 = ZERO TT = T + T DO 100 I = N , 0 , -1 U2 = U1 U1 = U0 U0 = TT * U1 + A( I ) - U2 100 CONTINUE CHEVAL = ( U0 - U2 ) / TWO ELSE C C If ABS ( T ) > = 0.6 use the Reinsch modification C D1 = ZERO C C T > = 0.6 code C IF ( T .GT. ZERO ) THEN TT = ( T - HALF ) - HALF TT = TT + TT DO 200 I = N , 0 , -1 D2 = D1 U2 = U1 D1 = TT * U2 + A( I ) + D2 U1 = D1 + U2 200 CONTINUE CHEVAL = ( D1 + D2 ) / TWO ELSE C C T < = -0.6 code C TT = ( T + HALF ) + HALF TT = TT + TT DO 300 I = N , 0 , -1 D2 = D1 U2 = U1 D1 = TT * U2 + A( I ) - D2 U1 = D1 - U2 300 CONTINUE CHEVAL = ( D1 - D2 ) / TWO ENDIF ENDIF RETURN END SHAR_EOF fi # end of overwriting check cd .. cd .. # End of shell archive exit 0