From 1e1406ff2ad64fdb9490eb962073d2a78384cf23 Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Thu, 11 Jun 2015 16:48:56 +0000 Subject: [PATCH] calcei module replaced with eierf, calerf removed from grayl. --- Makefile | 4 +- src/calcei_mod.f90 | 610 ------------------------------ src/dispersion.f90 | 2 +- src/eierf.f90 | 906 +++++++++++++++++++++++++++++++++++++++++++++ src/grayl.f | 250 ------------- 5 files changed, 909 insertions(+), 863 deletions(-) delete mode 100644 src/calcei_mod.f90 create mode 100644 src/eierf.f90 diff --git a/Makefile b/Makefile index 25bd673..c76cb49 100644 --- a/Makefile +++ b/Makefile @@ -3,7 +3,7 @@ EXE=gray # Objects list MAINOBJ=gray.o -OTHOBJ=dispersion.o calcei_mod.o dqagmv.o grayl.o reflections.o green_func_p.o \ +OTHOBJ=dispersion.o eierf.o dqagmv.o grayl.o reflections.o green_func_p.o \ const_and_precisions.o graydata_flags.o graydata_par.o graydata_anequil.o \ magsurf_data.o interp_eqprof.o @@ -28,7 +28,7 @@ gray.o: dispersion.o dqagmv.o green_func_p.o reflections.o const_and_precisions. graydata_flags.o graydata_par.o graydata_anequil.o magsurf_data.o interp_eqprof.o green_func_p.o: const_and_precisions.o reflections.o: const_and_precisions.o -dispersion.o: const_and_precisions.o calcei_mod.o dqagmv.o +dispersion.o: const_and_precisions.o eierf.o dqagmv.o graydata_flags.o: const_and_precisions.o graydata_par.o: const_and_precisions.o graydata_anequil.o: const_and_precisions.o diff --git a/src/calcei_mod.f90 b/src/calcei_mod.f90 deleted file mode 100644 index e57f502..0000000 --- a/src/calcei_mod.f90 +++ /dev/null @@ -1,610 +0,0 @@ -module calcei_mod - - implicit none - -contains - -! ====================================================================== -! nist guide to available math software. -! fullsource for module ei from package specfun. -! retrieved from netlib on fri mar 26 05:52:39 1999. -! ====================================================================== - subroutine calcei(arg,result,int) -!---------------------------------------------------------------------- -! -! this fortran 77 packet computes the exponential integrals eint(x), -! e1(x), and exp(-x)*eint(x) for real arguments x where -! -! integral (from t=-infinity to t=x) (exp(t)/t), x > 0, -! eint(x) = -! -integral (from t=-x to t=infinity) (exp(t)/t), x < 0, -! -! and where the first integral is a principal value integral. -! the packet contains three function type subprograms: ei, eone, -! and expei; and one subroutine type subprogram: calcei. the -! calling statements for the primary entries are -! -! y = eint(x), where x .ne. 0, -! -! y = eone(x), where x .gt. 0, -! and -! y = expei(x), where x .ne. 0, -! -! and where the entry points correspond to the functions eint(x), -! e1(x), and exp(-x)*eint(x), respectively. the routine calcei -! is intended for internal packet use only, all computations within -! the packet being concentrated in this routine. the function -! subprograms invoke calcei with the fortran statement -! call calcei(arg,result,int) -! where the parameter usage is as follows -! -! function parameters for calcei -! call arg result int -! -! eint(x) x .ne. 0 eint(x) 1 -! eone(x) x .gt. 0 -eint(-x) 2 -! expei(x) x .ne. 0 exp(-x)*eint(x) 3 -!---------------------------------------------------------------------- - integer i,int - double precision& - & a,arg,b,c,d,exp40,e,ei,f,four,fourty,frac,half,one,p,& - & plg,px,p037,p1,p2,q,qlg,qx,q1,q2,r,result,s,six,sump,& - & sumq,t,three,twelve,two,two4,w,x,xbig,xinf,xmax,xmx0,& - & x0,x01,x02,x11,y,ysq,zero - dimension a(7),b(6),c(9),d(9),e(10),f(10),p(10),q(10),r(10),& - & s(9),p1(10),q1(9),p2(10),q2(9),plg(4),qlg(4),px(10),qx(10) -!---------------------------------------------------------------------- -! mathematical constants -! exp40 = exp(40) -! x0 = zero of ei -! x01/x11 + x02 = zero of ei to extra precision -!---------------------------------------------------------------------- - data zero,p037,half,one,two/0.0d0,0.037d0,0.5d0,1.0d0,2.0d0/,& - & three,four,six,twelve,two4/3.0d0,4.0d0,6.0d0,12.d0,24.0d0/,& - & fourty,exp40/40.0d0,2.3538526683701998541d17/,& - & x01,x11,x02/381.5d0,1024.0d0,-5.1182968633365538008d-5/,& - & x0/3.7250741078136663466d-1/ -!---------------------------------------------------------------------- -! machine-dependent constants -!---------------------------------------------------------------------- - data xinf/1.79d+308/,xmax/716.351d0/,xbig/701.84d0/ -!---------------------------------------------------------------------- -! coefficients for -1.0 <= x < 0.0 -!---------------------------------------------------------------------- - data a/1.1669552669734461083368d2, 2.1500672908092918123209d3,& - & 1.5924175980637303639884d4, 8.9904972007457256553251d4,& - & 1.5026059476436982420737d5,-1.4815102102575750838086d5,& - & 5.0196785185439843791020d0/ - data b/4.0205465640027706061433d1, 7.5043163907103936624165d2,& - & 8.1258035174768735759855d3, 5.2440529172056355429883d4,& - & 1.8434070063353677359298d5, 2.5666493484897117319268d5/ -!---------------------------------------------------------------------- -! coefficients for -4.0 <= x < -1.0 -!---------------------------------------------------------------------- - data c/3.828573121022477169108d-1, 1.107326627786831743809d+1,& - & 7.246689782858597021199d+1, 1.700632978311516129328d+2,& - & 1.698106763764238382705d+2, 7.633628843705946890896d+1,& - & 1.487967702840464066613d+1, 9.999989642347613068437d-1,& - & 1.737331760720576030932d-8/ - data d/8.258160008564488034698d-2, 4.344836335509282083360d+0,& - & 4.662179610356861756812d+1, 1.775728186717289799677d+2,& - & 2.953136335677908517423d+2, 2.342573504717625153053d+2,& - & 9.021658450529372642314d+1, 1.587964570758947927903d+1,& - & 1.000000000000000000000d+0/ -!---------------------------------------------------------------------- -! coefficients for x < -4.0 -!---------------------------------------------------------------------- - data e/1.3276881505637444622987d+2,3.5846198743996904308695d+4,& - & 1.7283375773777593926828d+5,2.6181454937205639647381d+5,& - & 1.7503273087497081314708d+5,5.9346841538837119172356d+4,& - & 1.0816852399095915622498d+4,1.0611777263550331766871d03,& - & 5.2199632588522572481039d+1,9.9999999999999999087819d-1/ - data f/3.9147856245556345627078d+4,2.5989762083608489777411d+5,& - & 5.5903756210022864003380d+5,5.4616842050691155735758d+5,& - & 2.7858134710520842139357d+5,7.9231787945279043698718d+4,& - & 1.2842808586627297365998d+4,1.1635769915320848035459d+3,& - & 5.4199632588522559414924d+1,1.0d0/ -!---------------------------------------------------------------------- -! coefficients for rational approximation to ln(x/a), |1-x/a| < .1 -!---------------------------------------------------------------------- - data plg/-2.4562334077563243311d+01,2.3642701335621505212d+02,& - & -5.4989956895857911039d+02,3.5687548468071500413d+02/ - data qlg/-3.5553900764052419184d+01,1.9400230218539473193d+02,& - & -3.3442903192607538956d+02,1.7843774234035750207d+02/ -!---------------------------------------------------------------------- -! coefficients for 0.0 < x < 6.0, -! ratio of chebyshev polynomials -!---------------------------------------------------------------------- - data p/-1.2963702602474830028590d01,-1.2831220659262000678155d03,& - & -1.4287072500197005777376d04,-1.4299841572091610380064d06,& - & -3.1398660864247265862050d05,-3.5377809694431133484800d08,& - & 3.1984354235237738511048d08,-2.5301823984599019348858d10,& - & 1.2177698136199594677580d10,-2.0829040666802497120940d11/ - data q/ 7.6886718750000000000000d01,-5.5648470543369082846819d03,& - & 1.9418469440759880361415d05,-4.2648434812177161405483d06,& - & 6.4698830956576428587653d07,-7.0108568774215954065376d08,& - & 5.4229617984472955011862d09,-2.8986272696554495342658d10,& - & 9.8900934262481749439886d10,-8.9673749185755048616855d10/ -!---------------------------------------------------------------------- -! j-fraction coefficients for 6.0 <= x < 12.0 -!---------------------------------------------------------------------- - data r/-2.645677793077147237806d00,-2.378372882815725244124d00,& - & -2.421106956980653511550d01, 1.052976392459015155422d01,& - & 1.945603779539281810439d01,-3.015761863840593359165d01,& - & 1.120011024227297451523d01,-3.988850730390541057912d00,& - & 9.565134591978630774217d00, 9.981193787537396413219d-1/ - data s/ 1.598517957704779356479d-4, 4.644185932583286942650d00,& - & 3.697412299772985940785d02,-8.791401054875438925029d00,& - & 7.608194509086645763123d02, 2.852397548119248700147d01,& - & 4.731097187816050252967d02,-2.369210235636181001661d02,& - & 1.249884822712447891440d00/ -!---------------------------------------------------------------------- -! j-fraction coefficients for 12.0 <= x < 24.0 -!---------------------------------------------------------------------- - data p1/-1.647721172463463140042d00,-1.860092121726437582253d01,& - & -1.000641913989284829961d01,-2.105740799548040450394d01,& - & -9.134835699998742552432d-1,-3.323612579343962284333d01,& - & 2.495487730402059440626d01, 2.652575818452799819855d01,& - & -1.845086232391278674524d00, 9.999933106160568739091d-1/ - data q1/ 9.792403599217290296840d01, 6.403800405352415551324d01,& - & 5.994932325667407355255d01, 2.538819315630708031713d02,& - & 4.429413178337928401161d01, 1.192832423968601006985d03,& - & 1.991004470817742470726d02,-1.093556195391091143924d01,& - & 1.001533852045342697818d00/ -!---------------------------------------------------------------------- -! j-fraction coefficients for x .ge. 24.0 -!---------------------------------------------------------------------- - data p2/ 1.75338801265465972390d02,-2.23127670777632409550d02,& - & -1.81949664929868906455d01,-2.79798528624305389340d01,& - & -7.63147701620253630855d00,-1.52856623636929636839d01,& - & -7.06810977895029358836d00,-5.00006640413131002475d00,& - & -3.00000000320981265753d00, 1.00000000000000485503d00/ - data q2/ 3.97845977167414720840d04, 3.97277109100414518365d00,& - & 1.37790390235747998793d02, 1.17179220502086455287d02,& - & 7.04831847180424675988d01,-1.20187763547154743238d01,& - & -7.99243595776339741065d00,-2.99999894040324959612d00,& - & 1.99999999999048104167d00/ -!---------------------------------------------------------------------- - x = arg - if (x .eq. zero) then - ei = -xinf - if (int .eq. 2) ei = -ei - else if ((x .lt. zero) .or. (int .eq. 2)) then -!---------------------------------------------------------------------- -! calculate ei for negative argument or for e1. -!---------------------------------------------------------------------- - y = abs(x) - if (y .le. one) then - sump = a(7) * y + a(1) - sumq = y + b(1) - do 110 i = 2, 6 - sump = sump * y + a(i) - sumq = sumq * y + b(i) - 110 continue - ei = log(y) - sump / sumq - if (int .eq. 3) ei = ei * exp(y) - else if (y .le. four) then - w = one / y - sump = c(1) - sumq = d(1) - do 130 i = 2, 9 - sump = sump * w + c(i) - sumq = sumq * w + d(i) - 130 continue - ei = - sump / sumq - if (int .ne. 3) ei = ei * exp(-y) - else - if ((y .gt. xbig) .and. (int .lt. 3)) then - ei = zero - else - w = one / y - sump = e(1) - sumq = f(1) - do 150 i = 2, 10 - sump = sump * w + e(i) - sumq = sumq * w + f(i) - 150 continue - ei = -w * (one - w * sump / sumq ) - if (int .ne. 3) ei = ei * exp(-y) - end if - end if - if (int .eq. 2) ei = -ei - else if (x .lt. six) then -!---------------------------------------------------------------------- -! to improve conditioning, rational approximations are expressed -! in terms of chebyshev polynomials for 0 <= x < 6, and in -! continued fraction form for larger x. -!---------------------------------------------------------------------- - t = x + x - t = t / three - two - px(1) = zero - qx(1) = zero - px(2) = p(1) - qx(2) = q(1) - do 210 i = 2, 9 - px(i+1) = t * px(i) - px(i-1) + p(i) - qx(i+1) = t * qx(i) - qx(i-1) + q(i) - 210 continue - sump = half * t * px(10) - px(9) + p(10) - sumq = half * t * qx(10) - qx(9) + q(10) - frac = sump / sumq - xmx0 = (x - x01/x11) - x02 - if (abs(xmx0) .ge. p037) then - ei = log(x/x0) + xmx0 * frac - if (int .eq. 3) ei = exp(-x) * ei - else -!---------------------------------------------------------------------- -! special approximation to ln(x/x0) for x close to x0 -!---------------------------------------------------------------------- - y = xmx0 / (x + x0) - ysq = y*y - sump = plg(1) - sumq = ysq + qlg(1) - do 220 i = 2, 4 - sump = sump*ysq + plg(i) - sumq = sumq*ysq + qlg(i) - 220 continue - ei = (sump / (sumq*(x+x0)) + frac) * xmx0 - if (int .eq. 3) ei = exp(-x) * ei - end if - else if (x .lt. twelve) then - frac = zero - do 230 i = 1, 9 - frac = s(i) / (r(i) + x + frac) - 230 continue - ei = (r(10) + frac) / x - if (int .ne. 3) ei = ei * exp(x) - else if (x .le. two4) then - frac = zero - do 240 i = 1, 9 - frac = q1(i) / (p1(i) + x + frac) - 240 continue - ei = (p1(10) + frac) / x - if (int .ne. 3) ei = ei * exp(x) - else - if ((x .ge. xmax) .and. (int .lt. 3)) then - ei = xinf - else - y = one / x - frac = zero - do 250 i = 1, 9 - frac = q2(i) / (p2(i) + x + frac) - 250 continue - frac = p2(10) + frac - ei = y + y * y * frac - if (int .ne. 3) then - if (x .le. xmax-two4) then - ei = ei * exp(x) - else -!---------------------------------------------------------------------- -! calculation reformulated to avoid premature overflow -!---------------------------------------------------------------------- - ei = (ei * exp(x-fourty)) * exp40 - end if - end if - end if - end if - result = ei - return -!---------- last line of calcei ---------- - end - function eint(x) -!-------------------------------------------------------------------- -! -! this function program computes approximate values for the -! exponential integral eint(x), where x is real. -! -! author: w. j. cody -! -! latest modification: january 12, 1988 -! -!-------------------------------------------------------------------- - integer int - double precision eint, x, result -!-------------------------------------------------------------------- - int = 1 - call calcei(x,result,int) - eint = result - return -!---------- last line of ei ---------- - end - function expei(x) -!-------------------------------------------------------------------- -! -! this function program computes approximate values for the -! function exp(-x) * eint(x), where eint(x) is the exponential -! integral, and x is real. -! -! author: w. j. cody -! -! latest modification: january 12, 1988 -! -!-------------------------------------------------------------------- - integer int - double precision expei, x, result -!-------------------------------------------------------------------- - int = 3 - call calcei(x,result,int) - expei = result - return -!---------- last line of expei ---------- - end - function eone(x) -!-------------------------------------------------------------------- -! -! this function program computes approximate values for the -! exponential integral e1(x), where x is real. -! -! author: w. j. cody -! -! latest modification: january 12, 1988 -! -!-------------------------------------------------------------------- - integer int - double precision eone, x, result -!-------------------------------------------------------------------- - int = 2 - call calcei(x,result,int) - eone = result - return -!---------- last line of eone ---------- - end -! -! calcei3 = calcei for int=3 -! -! ====================================================================== - subroutine calcei3(arg,result) -!---------------------------------------------------------------------- -! -! this fortran 77 packet computes the exponential integrals eint(x), -! e1(x), and exp(-x)*eint(x) for real arguments x where -! -! integral (from t=-infinity to t=x) (exp(t)/t), x > 0, -! eint(x) = -! -integral (from t=-x to t=infinity) (exp(t)/t), x < 0, -! -! and where the first integral is a principal value integral. -! the packet contains three function type subprograms: ei, eone, -! and expei; and one subroutine type subprogram: calcei. the -! calling statements for the primary entries are -! -! y = eint(x), where x .ne. 0, -! -! y = eone(x), where x .gt. 0, -! and -! y = expei(x), where x .ne. 0, -! -! and where the entry points correspond to the functions eint(x), -! e1(x), and exp(-x)*eint(x), respectively. the routine calcei -! is intended for internal packet use only, all computations within -! the packet being concentrated in this routine. the function -! subprograms invoke calcei with the fortran statement -! call calcei(arg,result,int) -! where the parameter usage is as follows -! -! function parameters for calcei -! call arg result int -! -! eint(x) x .ne. 0 eint(x) 1 -! eone(x) x .gt. 0 -eint(-x) 2 -! expei(x) x .ne. 0 exp(-x)*eint(x) 3 -!---------------------------------------------------------------------- - integer i,int - double precision& - & a,arg,b,c,d,exp40,e,ei,f,four,fourty,frac,half,one,p,& - & plg,px,p037,p1,p2,q,qlg,qx,q1,q2,r,result,s,six,sump,& - & sumq,t,three,twelve,two,two4,w,x,xbig,xinf,xmax,xmx0,& - & x0,x01,x02,x11,y,ysq,zero - dimension a(7),b(6),c(9),d(9),e(10),f(10),p(10),q(10),r(10),& - & s(9),p1(10),q1(9),p2(10),q2(9),plg(4),qlg(4),px(10),qx(10) -!---------------------------------------------------------------------- -! mathematical constants -! exp40 = exp(40) -! x0 = zero of ei -! x01/x11 + x02 = zero of ei to extra precision -!---------------------------------------------------------------------- - data zero,p037,half,one,two/0.0d0,0.037d0,0.5d0,1.0d0,2.0d0/,& - & three,four,six,twelve,two4/3.0d0,4.0d0,6.0d0,12.d0,24.0d0/,& - & fourty,exp40/40.0d0,2.3538526683701998541d17/,& - & x01,x11,x02/381.5d0,1024.0d0,-5.1182968633365538008d-5/,& - & x0/3.7250741078136663466d-1/ -!---------------------------------------------------------------------- -! machine-dependent constants -!---------------------------------------------------------------------- - data xinf/1.79d+308/,xmax/716.351d0/,xbig/701.84d0/ -!---------------------------------------------------------------------- -! coefficients for -1.0 <= x < 0.0 -!---------------------------------------------------------------------- - data a/1.1669552669734461083368d2, 2.1500672908092918123209d3,& - & 1.5924175980637303639884d4, 8.9904972007457256553251d4,& - & 1.5026059476436982420737d5,-1.4815102102575750838086d5,& - & 5.0196785185439843791020d0/ - data b/4.0205465640027706061433d1, 7.5043163907103936624165d2,& - & 8.1258035174768735759855d3, 5.2440529172056355429883d4,& - & 1.8434070063353677359298d5, 2.5666493484897117319268d5/ -!---------------------------------------------------------------------- -! coefficients for -4.0 <= x < -1.0 -!---------------------------------------------------------------------- - data c/3.828573121022477169108d-1, 1.107326627786831743809d+1,& - & 7.246689782858597021199d+1, 1.700632978311516129328d+2,& - & 1.698106763764238382705d+2, 7.633628843705946890896d+1,& - & 1.487967702840464066613d+1, 9.999989642347613068437d-1,& - & 1.737331760720576030932d-8/ - data d/8.258160008564488034698d-2, 4.344836335509282083360d+0,& - & 4.662179610356861756812d+1, 1.775728186717289799677d+2,& - & 2.953136335677908517423d+2, 2.342573504717625153053d+2,& - & 9.021658450529372642314d+1, 1.587964570758947927903d+1,& - & 1.000000000000000000000d+0/ -!---------------------------------------------------------------------- -! coefficients for x < -4.0 -!---------------------------------------------------------------------- - data e/1.3276881505637444622987d+2,3.5846198743996904308695d+4,& - & 1.7283375773777593926828d+5,2.6181454937205639647381d+5,& - & 1.7503273087497081314708d+5,5.9346841538837119172356d+4,& - & 1.0816852399095915622498d+4,1.0611777263550331766871d03,& - & 5.2199632588522572481039d+1,9.9999999999999999087819d-1/ - data f/3.9147856245556345627078d+4,2.5989762083608489777411d+5,& - & 5.5903756210022864003380d+5,5.4616842050691155735758d+5,& - & 2.7858134710520842139357d+5,7.9231787945279043698718d+4,& - & 1.2842808586627297365998d+4,1.1635769915320848035459d+3,& - & 5.4199632588522559414924d+1,1.0d0/ -!---------------------------------------------------------------------- -! coefficients for rational approximation to ln(x/a), |1-x/a| < .1 -!---------------------------------------------------------------------- - data plg/-2.4562334077563243311d+01,2.3642701335621505212d+02,& - & -5.4989956895857911039d+02,3.5687548468071500413d+02/ - data qlg/-3.5553900764052419184d+01,1.9400230218539473193d+02,& - & -3.3442903192607538956d+02,1.7843774234035750207d+02/ -!---------------------------------------------------------------------- -! coefficients for 0.0 < x < 6.0, -! ratio of chebyshev polynomials -!---------------------------------------------------------------------- - data p/-1.2963702602474830028590d01,-1.2831220659262000678155d03,& - & -1.4287072500197005777376d04,-1.4299841572091610380064d06,& - & -3.1398660864247265862050d05,-3.5377809694431133484800d08,& - & 3.1984354235237738511048d08,-2.5301823984599019348858d10,& - & 1.2177698136199594677580d10,-2.0829040666802497120940d11/ - data q/ 7.6886718750000000000000d01,-5.5648470543369082846819d03,& - & 1.9418469440759880361415d05,-4.2648434812177161405483d06,& - & 6.4698830956576428587653d07,-7.0108568774215954065376d08,& - & 5.4229617984472955011862d09,-2.8986272696554495342658d10,& - & 9.8900934262481749439886d10,-8.9673749185755048616855d10/ -!---------------------------------------------------------------------- -! j-fraction coefficients for 6.0 <= x < 12.0 -!---------------------------------------------------------------------- - data r/-2.645677793077147237806d00,-2.378372882815725244124d00,& - & -2.421106956980653511550d01, 1.052976392459015155422d01,& - & 1.945603779539281810439d01,-3.015761863840593359165d01,& - & 1.120011024227297451523d01,-3.988850730390541057912d00,& - & 9.565134591978630774217d00, 9.981193787537396413219d-1/ - data s/ 1.598517957704779356479d-4, 4.644185932583286942650d00,& - & 3.697412299772985940785d02,-8.791401054875438925029d00,& - & 7.608194509086645763123d02, 2.852397548119248700147d01,& - & 4.731097187816050252967d02,-2.369210235636181001661d02,& - & 1.249884822712447891440d00/ -!---------------------------------------------------------------------- -! j-fraction coefficients for 12.0 <= x < 24.0 -!---------------------------------------------------------------------- - data p1/-1.647721172463463140042d00,-1.860092121726437582253d01,& - & -1.000641913989284829961d01,-2.105740799548040450394d01,& - & -9.134835699998742552432d-1,-3.323612579343962284333d01,& - & 2.495487730402059440626d01, 2.652575818452799819855d01,& - & -1.845086232391278674524d00, 9.999933106160568739091d-1/ - data q1/ 9.792403599217290296840d01, 6.403800405352415551324d01,& - & 5.994932325667407355255d01, 2.538819315630708031713d02,& - & 4.429413178337928401161d01, 1.192832423968601006985d03,& - & 1.991004470817742470726d02,-1.093556195391091143924d01,& - & 1.001533852045342697818d00/ -!---------------------------------------------------------------------- -! j-fraction coefficients for x .ge. 24.0 -!---------------------------------------------------------------------- - data p2/ 1.75338801265465972390d02,-2.23127670777632409550d02,& - & -1.81949664929868906455d01,-2.79798528624305389340d01,& - & -7.63147701620253630855d00,-1.52856623636929636839d01,& - & -7.06810977895029358836d00,-5.00006640413131002475d00,& - & -3.00000000320981265753d00, 1.00000000000000485503d00/ - data q2/ 3.97845977167414720840d04, 3.97277109100414518365d00,& - & 1.37790390235747998793d02, 1.17179220502086455287d02,& - & 7.04831847180424675988d01,-1.20187763547154743238d01,& - & -7.99243595776339741065d00,-2.99999894040324959612d00,& - & 1.99999999999048104167d00/ -!---------------------------------------------------------------------- - data int/ 3/ - x = arg - if (x .eq. zero) then - ei = -xinf - else if ((x .lt. zero)) then -!---------------------------------------------------------------------- -! calculate ei for negative argument or for e1. -!---------------------------------------------------------------------- - y = abs(x) - if (y .le. one) then - sump = a(7) * y + a(1) - sumq = y + b(1) - do 110 i = 2, 6 - sump = sump * y + a(i) - sumq = sumq * y + b(i) - 110 continue - ei = (log(y) - sump / sumq ) * exp(y) - else if (y .le. four) then - w = one / y - sump = c(1) - sumq = d(1) - do 130 i = 2, 9 - sump = sump * w + c(i) - sumq = sumq * w + d(i) - 130 continue - ei = - sump / sumq - else - w = one / y - sump = e(1) - sumq = f(1) - do 150 i = 2, 10 - sump = sump * w + e(i) - sumq = sumq * w + f(i) - 150 continue - ei = -w * (one - w * sump / sumq ) - end if - else if (x .lt. six) then -!---------------------------------------------------------------------- -! to improve conditioning, rational approximations are expressed -! in terms of chebyshev polynomials for 0 <= x < 6, and in -! continued fraction form for larger x. -!---------------------------------------------------------------------- - t = x + x - t = t / three - two - px(1) = zero - qx(1) = zero - px(2) = p(1) - qx(2) = q(1) - do 210 i = 2, 9 - px(i+1) = t * px(i) - px(i-1) + p(i) - qx(i+1) = t * qx(i) - qx(i-1) + q(i) - 210 continue - sump = half * t * px(10) - px(9) + p(10) - sumq = half * t * qx(10) - qx(9) + q(10) - frac = sump / sumq - xmx0 = (x - x01/x11) - x02 - if (abs(xmx0) .ge. p037) then - ei = exp(-x) * ( log(x/x0) + xmx0 * frac ) - else -!---------------------------------------------------------------------- -! special approximation to ln(x/x0) for x close to x0 -!---------------------------------------------------------------------- - y = xmx0 / (x + x0) - ysq = y*y - sump = plg(1) - sumq = ysq + qlg(1) - do 220 i = 2, 4 - sump = sump*ysq + plg(i) - sumq = sumq*ysq + qlg(i) - 220 continue - ei = exp(-x) * (sump / (sumq*(x+x0)) + frac) * xmx0 - end if - else if (x .lt. twelve) then - frac = zero - do 230 i = 1, 9 - frac = s(i) / (r(i) + x + frac) - 230 continue - ei = (r(10) + frac) / x - else if (x .le. two4) then - frac = zero - do 240 i = 1, 9 - frac = q1(i) / (p1(i) + x + frac) - 240 continue - ei = (p1(10) + frac) / x - else - y = one / x - frac = zero - do 250 i = 1, 9 - frac = q2(i) / (p2(i) + x + frac) - 250 continue - frac = p2(10) + frac - ei = y + y * y * frac - end if - result = ei - return -!---------- last line of calcei ---------- - end - -end module calcei_mod \ No newline at end of file diff --git a/src/dispersion.f90 b/src/dispersion.f90 index 34c8218..2dd5392 100644 --- a/src/dispersion.f90 +++ b/src/dispersion.f90 @@ -1,7 +1,7 @@ module dispersion ! use const_and_precisions, only : wp_,zero,one,im,czero,cunit,pi,sqrt_pi - use calcei_mod, only : calcei3 + use eierf, only : calcei3 implicit none ! local constants integer, parameter :: npts=500 diff --git a/src/eierf.f90 b/src/eierf.f90 new file mode 100644 index 0000000..ef3dd8f --- /dev/null +++ b/src/eierf.f90 @@ -0,0 +1,906 @@ +module eierf + + use const_and_precisions, only : wp_, zero, one + implicit none + real(wp_), parameter, private :: half=0.5_wp_, two=2.0_wp_, three=3.0_wp_, & + four=4.0_wp_, six=6.0_wp_, twelve=12._wp_, sixten=16.0_wp_, & + two4=24.0_wp_, fourty=40.0_wp_ + +contains + +! ====================================================================== +! nist guide to available math software. +! fullsource for module ei from package specfun. +! retrieved from netlib on fri mar 26 05:52:39 1999. +! ====================================================================== + subroutine calcei(arg,result,intt) +!---------------------------------------------------------------------- +! +! this fortran 77 packet computes the exponential integrals ei(x), +! e1(x), and exp(-x)*ei(x) for real arguments x where +! +! integral (from t=-infinity to t=x) (exp(t)/t), x > 0, +! ei(x) = +! -integral (from t=-x to t=infinity) (exp(t)/t), x < 0, +! +! and where the first integral is a principal value integral. +! the packet contains three function type subprograms: ei, eone, +! and expei; and one subroutine type subprogram: calcei. the +! calling statements for the primary entries are +! +! y = ei(x), where x /= 0, +! +! y = eone(x), where x > 0, +! and +! y = expei(x), where x /= 0, +! +! and where the entry points correspond to the functions ei(x), +! e1(x), and exp(-x)*ei(x), respectively. the routine calcei +! is intended for internal packet use only, all computations within +! the packet being concentrated in this routine. the function +! subprograms invoke calcei with the fortran statement +! call calcei(arg,result,intt) +! where the parameter usage is as follows +! +! function parameters for calcei +! call arg result intt +! +! ei(x) x /= 0 ei(x) 1 +! eone(x) x > 0 -ei(-x) 2 +! expei(x) x /= 0 exp(-x)*ei(x) 3 +!---------------------------------------------------------------------- + implicit none + integer, intent(in) :: intt + real(wp_), intent(in) :: arg + real(wp_), intent(out) :: result + integer :: i + real(wp_) :: ei,frac,sump,sumq,t,w,x,xmx0,y,ysq + real(wp_), dimension(10) :: px,qx +!---------------------------------------------------------------------- +! mathematical constants +! exp40 = exp(40) +! x0 = zero of ei +! x01/x11 + x02 = zero of ei to extra precision +!---------------------------------------------------------------------- + real(wp_), parameter :: p037=0.037_wp_, & + exp40=2.3538526683701998541e17_wp_, x01=381.5_wp_, x11=1024.0_wp_, & + x02=-5.1182968633365538008e-5_wp_, x0=3.7250741078136663466e-1_wp_ +!---------------------------------------------------------------------- +! machine-dependent constants +!---------------------------------------------------------------------- + real(wp_), parameter :: xinf=1.79e+308_wp_,xmax=716.351_wp_,xbig=701.84_wp_ +!---------------------------------------------------------------------- +! coefficients for -1.0 <= x < 0.0 +!---------------------------------------------------------------------- + real(wp_), dimension(7), parameter :: & + a=(/1.1669552669734461083368e2_wp_, 2.1500672908092918123209e3_wp_, & + 1.5924175980637303639884e4_wp_, 8.9904972007457256553251e4_wp_, & + 1.5026059476436982420737e5_wp_,-1.4815102102575750838086e5_wp_, & + 5.0196785185439843791020_wp_/) + real(wp_), dimension(6), parameter :: & + b=(/4.0205465640027706061433e1_wp_, 7.5043163907103936624165e2_wp_, & + 8.1258035174768735759855e3_wp_, 5.2440529172056355429883e4_wp_, & + 1.8434070063353677359298e5_wp_, 2.5666493484897117319268e5_wp_/) +!---------------------------------------------------------------------- +! coefficients for -4.0 <= x < -1.0 +!---------------------------------------------------------------------- + real(wp_), dimension(9), parameter :: & + c=(/3.828573121022477169108e-1_wp_, 1.107326627786831743809e+1_wp_, & + 7.246689782858597021199e+1_wp_, 1.700632978311516129328e+2_wp_, & + 1.698106763764238382705e+2_wp_, 7.633628843705946890896e+1_wp_, & + 1.487967702840464066613e+1_wp_, 9.999989642347613068437e-1_wp_, & + 1.737331760720576030932e-8_wp_/), & + d=(/8.258160008564488034698e-2_wp_, 4.344836335509282083360e+0_wp_, & + 4.662179610356861756812e+1_wp_, 1.775728186717289799677e+2_wp_, & + 2.953136335677908517423e+2_wp_, 2.342573504717625153053e+2_wp_, & + 9.021658450529372642314e+1_wp_, 1.587964570758947927903e+1_wp_, & + 1.000000000000000000000e+0_wp_/) +!---------------------------------------------------------------------- +! coefficients for x < -4.0 +!---------------------------------------------------------------------- + real(wp_), dimension(10), parameter :: & + e=(/1.3276881505637444622987e+2_wp_,3.5846198743996904308695e+4_wp_, & + 1.7283375773777593926828e+5_wp_,2.6181454937205639647381e+5_wp_, & + 1.7503273087497081314708e+5_wp_,5.9346841538837119172356e+4_wp_, & + 1.0816852399095915622498e+4_wp_,1.0611777263550331766871e03_wp_, & + 5.2199632588522572481039e+1_wp_,9.9999999999999999087819e-1_wp_/),& + f=(/3.9147856245556345627078e+4_wp_,2.5989762083608489777411e+5_wp_, & + 5.5903756210022864003380e+5_wp_,5.4616842050691155735758e+5_wp_, & + 2.7858134710520842139357e+5_wp_,7.9231787945279043698718e+4_wp_, & + 1.2842808586627297365998e+4_wp_,1.1635769915320848035459e+3_wp_, & + 5.4199632588522559414924e+1_wp_,1.0_wp_/) +!---------------------------------------------------------------------- +! coefficients for rational approximation to ln(x/a), |1-x/a| < .1 +!---------------------------------------------------------------------- + real(wp_), dimension(4), parameter :: & + plg=(/-2.4562334077563243311e+01_wp_,2.3642701335621505212e+02_wp_, & + -5.4989956895857911039e+02_wp_,3.5687548468071500413e+02_wp_/), & + qlg=(/-3.5553900764052419184e+01_wp_,1.9400230218539473193e+02_wp_, & + -3.3442903192607538956e+02_wp_,1.7843774234035750207e+02_wp_/) +!---------------------------------------------------------------------- +! coefficients for 0.0 < x < 6.0, +! ratio of chebyshev polynomials +!---------------------------------------------------------------------- + real(wp_), dimension(10), parameter :: & + p=(/-1.2963702602474830028590e01_wp_,-1.2831220659262000678155e03_wp_, & + -1.4287072500197005777376e04_wp_,-1.4299841572091610380064e06_wp_, & + -3.1398660864247265862050e05_wp_,-3.5377809694431133484800e08_wp_, & + 3.1984354235237738511048e08_wp_,-2.5301823984599019348858e10_wp_, & + 1.2177698136199594677580e10_wp_,-2.0829040666802497120940e11_wp_/),& + q=(/ 7.6886718750000000000000e01_wp_,-5.5648470543369082846819e03_wp_, & + 1.9418469440759880361415e05_wp_,-4.2648434812177161405483e06_wp_, & + 6.4698830956576428587653e07_wp_,-7.0108568774215954065376e08_wp_, & + 5.4229617984472955011862e09_wp_,-2.8986272696554495342658e10_wp_, & + 9.8900934262481749439886e10_wp_,-8.9673749185755048616855e10_wp_/) +!---------------------------------------------------------------------- +! j-fraction coefficients for 6.0 <= x < 12.0 +!---------------------------------------------------------------------- + real(wp_), dimension(10), parameter :: & + r=(/-2.645677793077147237806_wp_,-2.378372882815725244124_wp_, & + -2.421106956980653511550e01_wp_, 1.052976392459015155422e01_wp_, & + 1.945603779539281810439e01_wp_,-3.015761863840593359165e01_wp_, & + 1.120011024227297451523e01_wp_,-3.988850730390541057912_wp_, & + 9.565134591978630774217_wp_, 9.981193787537396413219e-1_wp_/) + real(wp_), dimension(9), parameter :: & + s=(/ 1.598517957704779356479e-4_wp_, 4.644185932583286942650_wp_, & + 3.697412299772985940785e02_wp_,-8.791401054875438925029_wp_, & + 7.608194509086645763123e02_wp_, 2.852397548119248700147e01_wp_, & + 4.731097187816050252967e02_wp_,-2.369210235636181001661e02_wp_, & + 1.249884822712447891440_wp_/) +!---------------------------------------------------------------------- +! j-fraction coefficients for 12.0 <= x < 24.0 +!---------------------------------------------------------------------- + real(wp_), dimension(10), parameter :: & + p1=(/-1.647721172463463140042_wp_,-1.860092121726437582253e01_wp_, & + -1.000641913989284829961e01_wp_,-2.105740799548040450394e01_wp_, & + -9.134835699998742552432e-1_wp_,-3.323612579343962284333e01_wp_, & + 2.495487730402059440626e01_wp_, 2.652575818452799819855e01_wp_, & + -1.845086232391278674524_wp_, 9.999933106160568739091e-1_wp_/) + real(wp_), dimension(9), parameter :: & + q1=(/ 9.792403599217290296840e01_wp_, 6.403800405352415551324e01_wp_, & + 5.994932325667407355255e01_wp_, 2.538819315630708031713e02_wp_, & + 4.429413178337928401161e01_wp_, 1.192832423968601006985e03_wp_, & + 1.991004470817742470726e02_wp_,-1.093556195391091143924e01_wp_, & + 1.001533852045342697818_wp_/) +!---------------------------------------------------------------------- +! j-fraction coefficients for x >= 24.0 +!---------------------------------------------------------------------- + real(wp_), dimension(10), parameter :: & + p2=(/ 1.75338801265465972390e02_wp_,-2.23127670777632409550e02_wp_, & + -1.81949664929868906455e01_wp_,-2.79798528624305389340e01_wp_, & + -7.63147701620253630855_wp_,-1.52856623636929636839e01_wp_, & + -7.06810977895029358836_wp_,-5.00006640413131002475_wp_, & + -3.00000000320981265753_wp_, 1.00000000000000485503_wp_/) + real(wp_), dimension(9), parameter :: & + q2=(/ 3.97845977167414720840e04_wp_, 3.97277109100414518365_wp_, & + 1.37790390235747998793e02_wp_, 1.17179220502086455287e02_wp_, & + 7.04831847180424675988e01_wp_,-1.20187763547154743238e01_wp_, & + -7.99243595776339741065_wp_,-2.99999894040324959612_wp_, & + 1.99999999999048104167_wp_/) +!---------------------------------------------------------------------- + x = arg + if (x == zero) then + ei = -xinf + if (intt == 2) ei = -ei + else if ((x < zero) .or. (intt == 2)) then +!---------------------------------------------------------------------- +! calculate ei for negative argument or for e1. +!---------------------------------------------------------------------- + y = abs(x) + if (y <= one) then + sump = a(7) * y + a(1) + sumq = y + b(1) + do i = 2, 6 + sump = sump * y + a(i) + sumq = sumq * y + b(i) + end do + ei = log(y) - sump / sumq + if (intt == 3) ei = ei * exp(y) + else if (y <= four) then + w = one / y + sump = c(1) + sumq = d(1) + do i = 2, 9 + sump = sump * w + c(i) + sumq = sumq * w + d(i) + end do + ei = - sump / sumq + if (intt /= 3) ei = ei * exp(-y) + else + if ((y > xbig) .and. (intt < 3)) then + ei = zero + else + w = one / y + sump = e(1) + sumq = f(1) + do i = 2, 10 + sump = sump * w + e(i) + sumq = sumq * w + f(i) + end do + ei = -w * (one - w * sump / sumq ) + if (intt /= 3) ei = ei * exp(-y) + end if + end if + if (intt == 2) ei = -ei + else if (x < six) then +!---------------------------------------------------------------------- +! to improve conditioning, rational approximations are expressed +! in terms of chebyshev polynomials for 0 <= x < 6, and in +! continued fraction form for larger x. +!---------------------------------------------------------------------- + t = x + x + t = t / three - two + px(1) = zero + qx(1) = zero + px(2) = p(1) + qx(2) = q(1) + do i = 2, 9 + px(i+1) = t * px(i) - px(i-1) + p(i) + qx(i+1) = t * qx(i) - qx(i-1) + q(i) + end do + sump = half * t * px(10) - px(9) + p(10) + sumq = half * t * qx(10) - qx(9) + q(10) + frac = sump / sumq + xmx0 = (x - x01/x11) - x02 + if (abs(xmx0) >= p037) then + ei = log(x/x0) + xmx0 * frac + if (intt == 3) ei = exp(-x) * ei + else +!---------------------------------------------------------------------- +! special approximation to ln(x/x0) for x close to x0 +!---------------------------------------------------------------------- + y = xmx0 / (x + x0) + ysq = y*y + sump = plg(1) + sumq = ysq + qlg(1) + do i = 2, 4 + sump = sump*ysq + plg(i) + sumq = sumq*ysq + qlg(i) + end do + ei = (sump / (sumq*(x+x0)) + frac) * xmx0 + if (intt == 3) ei = exp(-x) * ei + end if + else if (x < twelve) then + frac = zero + do i = 1, 9 + frac = s(i) / (r(i) + x + frac) + end do + ei = (r(10) + frac) / x + if (intt /= 3) ei = ei * exp(x) + else if (x <= two4) then + frac = zero + do i = 1, 9 + frac = q1(i) / (p1(i) + x + frac) + end do + ei = (p1(10) + frac) / x + if (intt /= 3) ei = ei * exp(x) + else + if ((x >= xmax) .and. (intt < 3)) then + ei = xinf + else + y = one / x + frac = zero + do i = 1, 9 + frac = q2(i) / (p2(i) + x + frac) + end do + frac = p2(10) + frac + ei = y + y * y * frac + if (intt /= 3) then + if (x <= xmax-two4) then + ei = ei * exp(x) + else +!---------------------------------------------------------------------- +! calculation reformulated to avoid premature overflow +!---------------------------------------------------------------------- + ei = (ei * exp(x-fourty)) * exp40 + end if + end if + end if + end if + result = ei + end subroutine calcei + + function ei(x) +!-------------------------------------------------------------------- +! +! this function program computes approximate values for the +! exponential integral ei(x), where x is real. +! +! author: w. j. cody +! +! latest modification: january 12, 1988 +! +!-------------------------------------------------------------------- + implicit none + integer :: intt + real(wp_) :: ei + real(wp_), intent(in) :: x + real(wp_) :: result +!-------------------------------------------------------------------- + intt = 1 + call calcei(x,result,intt) + ei = result + end function ei + + function expei(x) +!-------------------------------------------------------------------- +! +! this function program computes approximate values for the +! function exp(-x) * ei(x), where ei(x) is the exponential +! integral, and x is real. +! +! author: w. j. cody +! +! latest modification: january 12, 1988 +! +!-------------------------------------------------------------------- + implicit none + integer :: intt + real(wp_) :: expei + real(wp_), intent(in) :: x + real(wp_) :: result +!-------------------------------------------------------------------- + intt = 3 + call calcei(x,result,intt) + expei = result + end function expei + + function eone(x) +!-------------------------------------------------------------------- +! +! this function program computes approximate values for the +! exponential integral e1(x), where x is real. +! +! author: w. j. cody +! +! latest modification: january 12, 1988 +! +!-------------------------------------------------------------------- + implicit none + integer :: intt + real(wp_) :: eone + real(wp_), intent(in) :: x + real(wp_) :: result +!-------------------------------------------------------------------- + intt = 2 + call calcei(x,result,intt) + eone = result + end function eone + +! ====================================================================== +! calcei3 = calcei for int=3 +! ====================================================================== + subroutine calcei3(arg,result) +!---------------------------------------------------------------------- +! +! this fortran 77 packet computes the exponential integrals ei(x), +! e1(x), and exp(-x)*ei(x) for real arguments x where +! +! integral (from t=-infinity to t=x) (exp(t)/t), x > 0, +! ei(x) = +! -integral (from t=-x to t=infinity) (exp(t)/t), x < 0, +! +! and where the first integral is a principal value integral. +! the packet contains three function type subprograms: ei, eone, +! and expei; and one subroutine type subprogram: calcei. the +! calling statements for the primary entries are +! +! y = ei(x), where x /= 0, +! +! y = eone(x), where x > 0, +! and +! y = expei(x), where x /= 0, +! +! and where the entry points correspond to the functions ei(x), +! e1(x), and exp(-x)*ei(x), respectively. the routine calcei +! is intended for internal packet use only, all computations within +! the packet being concentrated in this routine. the function +! subprograms invoke calcei with the fortran statement +! call calcei(arg,result,int) +! where the parameter usage is as follows +! +! function parameters for calcei +! call arg result int +! +! ei(x) x /= 0 ei(x) 1 +! eone(x) x > 0 -ei(-x) 2 +! expei(x) x /= 0 exp(-x)*ei(x) 3 +!---------------------------------------------------------------------- + implicit none + real(wp_), intent(in) :: arg + real(wp_), intent(out) :: result + integer :: i + real(wp_) :: ei,frac,sump,sumq,t,w,x,xmx0,y,ysq + real(wp_), dimension(10) :: px,qx +!---------------------------------------------------------------------- +! mathematical constants +! exp40 = exp(40) +! x0 = zero of ei +! x01/x11 + x02 = zero of ei to extra precision +!---------------------------------------------------------------------- + real(wp_), parameter :: p037=0.037_wp_, & + x01=381.5_wp_, x11=1024.0_wp_, x02=-5.1182968633365538008e-5_wp_, & + x0=3.7250741078136663466e-1_wp_ +!---------------------------------------------------------------------- +! machine-dependent constants +!---------------------------------------------------------------------- + real(wp_), parameter :: xinf=1.79e+308_wp_ +!---------------------------------------------------------------------- +! coefficients for -1.0 <= x < 0.0 +!---------------------------------------------------------------------- + real(wp_), dimension(7), parameter :: & + a=(/1.1669552669734461083368e2_wp_, 2.1500672908092918123209e3_wp_, & + 1.5924175980637303639884e4_wp_, 8.9904972007457256553251e4_wp_, & + 1.5026059476436982420737e5_wp_,-1.4815102102575750838086e5_wp_, & + 5.0196785185439843791020_wp_/) + real(wp_), dimension(6), parameter :: & + b=(/4.0205465640027706061433e1_wp_, 7.5043163907103936624165e2_wp_, & + 8.1258035174768735759855e3_wp_, 5.2440529172056355429883e4_wp_, & + 1.8434070063353677359298e5_wp_, 2.5666493484897117319268e5_wp_/) +!---------------------------------------------------------------------- +! coefficients for -4.0 <= x < -1.0 +!---------------------------------------------------------------------- + real(wp_), dimension(9), parameter :: & + c=(/3.828573121022477169108e-1_wp_, 1.107326627786831743809e+1_wp_, & + 7.246689782858597021199e+1_wp_, 1.700632978311516129328e+2_wp_, & + 1.698106763764238382705e+2_wp_, 7.633628843705946890896e+1_wp_, & + 1.487967702840464066613e+1_wp_, 9.999989642347613068437e-1_wp_, & + 1.737331760720576030932e-8_wp_/), & + d=(/8.258160008564488034698e-2_wp_, 4.344836335509282083360e+0_wp_, & + 4.662179610356861756812e+1_wp_, 1.775728186717289799677e+2_wp_, & + 2.953136335677908517423e+2_wp_, 2.342573504717625153053e+2_wp_, & + 9.021658450529372642314e+1_wp_, 1.587964570758947927903e+1_wp_, & + 1.000000000000000000000e+0_wp_/) +!---------------------------------------------------------------------- +! coefficients for x < -4.0 +!---------------------------------------------------------------------- + real(wp_), dimension(10), parameter :: & + e=(/1.3276881505637444622987e+2_wp_,3.5846198743996904308695e+4_wp_, & + 1.7283375773777593926828e+5_wp_,2.6181454937205639647381e+5_wp_, & + 1.7503273087497081314708e+5_wp_,5.9346841538837119172356e+4_wp_, & + 1.0816852399095915622498e+4_wp_,1.0611777263550331766871e03_wp_, & + 5.2199632588522572481039e+1_wp_,9.9999999999999999087819e-1_wp_/), & + f=(/3.9147856245556345627078e+4_wp_,2.5989762083608489777411e+5_wp_, & + 5.5903756210022864003380e+5_wp_,5.4616842050691155735758e+5_wp_, & + 2.7858134710520842139357e+5_wp_,7.9231787945279043698718e+4_wp_, & + 1.2842808586627297365998e+4_wp_,1.1635769915320848035459e+3_wp_, & + 5.4199632588522559414924e+1_wp_,1.0_wp_/) +!---------------------------------------------------------------------- +! coefficients for rational approximation to ln(x/a), |1-x/a| < .1 +!---------------------------------------------------------------------- + real(wp_), dimension(4), parameter :: & + plg=(/-2.4562334077563243311e+01_wp_,2.3642701335621505212e+02_wp_, & + -5.4989956895857911039e+02_wp_,3.5687548468071500413e+02_wp_/), & + qlg=(/-3.5553900764052419184e+01_wp_,1.9400230218539473193e+02_wp_, & + -3.3442903192607538956e+02_wp_,1.7843774234035750207e+02_wp_/) +!---------------------------------------------------------------------- +! coefficients for 0.0 < x < 6.0, +! ratio of chebyshev polynomials +!---------------------------------------------------------------------- + real(wp_), dimension(10), parameter :: & + p=(/-1.2963702602474830028590e01_wp_,-1.2831220659262000678155e03_wp_, & + -1.4287072500197005777376e04_wp_,-1.4299841572091610380064e06_wp_, & + -3.1398660864247265862050e05_wp_,-3.5377809694431133484800e08_wp_, & + 3.1984354235237738511048e08_wp_,-2.5301823984599019348858e10_wp_, & + 1.2177698136199594677580e10_wp_,-2.0829040666802497120940e11_wp_/),& + q=(/ 7.6886718750000000000000e01_wp_,-5.5648470543369082846819e03_wp_, & + 1.9418469440759880361415e05_wp_,-4.2648434812177161405483e06_wp_, & + 6.4698830956576428587653e07_wp_,-7.0108568774215954065376e08_wp_, & + 5.4229617984472955011862e09_wp_,-2.8986272696554495342658e10_wp_, & + 9.8900934262481749439886e10_wp_,-8.9673749185755048616855e10_wp_/) +!---------------------------------------------------------------------- +! j-fraction coefficients for 6.0 <= x < 12.0 +!---------------------------------------------------------------------- + real(wp_), dimension(10), parameter :: & + r=(/-2.645677793077147237806_wp_,-2.378372882815725244124_wp_, & + -2.421106956980653511550e01_wp_, 1.052976392459015155422e01_wp_, & + 1.945603779539281810439e01_wp_,-3.015761863840593359165e01_wp_, & + 1.120011024227297451523e01_wp_,-3.988850730390541057912_wp_, & + 9.565134591978630774217_wp_, 9.981193787537396413219e-1_wp_/) + real(wp_), dimension(9), parameter :: & + s=(/ 1.598517957704779356479e-4_wp_, 4.644185932583286942650_wp_, & + 3.697412299772985940785e02_wp_,-8.791401054875438925029_wp_, & + 7.608194509086645763123e02_wp_, 2.852397548119248700147e01_wp_, & + 4.731097187816050252967e02_wp_,-2.369210235636181001661e02_wp_, & + 1.249884822712447891440_wp_/) +!---------------------------------------------------------------------- +! j-fraction coefficients for 12.0 <= x < 24.0 +!---------------------------------------------------------------------- + real(wp_), dimension(10), parameter :: & + p1=(/-1.647721172463463140042_wp_,-1.860092121726437582253e01_wp_, & + -1.000641913989284829961e01_wp_,-2.105740799548040450394e01_wp_, & + -9.134835699998742552432e-1_wp_,-3.323612579343962284333e01_wp_, & + 2.495487730402059440626e01_wp_, 2.652575818452799819855e01_wp_, & + -1.845086232391278674524_wp_, 9.999933106160568739091e-1_wp_/) + real(wp_), dimension(9), parameter :: & + q1=(/ 9.792403599217290296840e01_wp_, 6.403800405352415551324e01_wp_, & + 5.994932325667407355255e01_wp_, 2.538819315630708031713e02_wp_, & + 4.429413178337928401161e01_wp_, 1.192832423968601006985e03_wp_, & + 1.991004470817742470726e02_wp_,-1.093556195391091143924e01_wp_, & + 1.001533852045342697818_wp_/) +!---------------------------------------------------------------------- +! j-fraction coefficients for x >= 24.0 +!---------------------------------------------------------------------- + real(wp_), dimension(10), parameter :: & + p2=(/ 1.75338801265465972390e02_wp_,-2.23127670777632409550e02_wp_, & + -1.81949664929868906455e01_wp_,-2.79798528624305389340e01_wp_, & + -7.63147701620253630855_wp_,-1.52856623636929636839e01_wp_, & + -7.06810977895029358836_wp_,-5.00006640413131002475_wp_, & + -3.00000000320981265753_wp_, 1.00000000000000485503_wp_/) + real(wp_), dimension(9), parameter :: & + q2=(/ 3.97845977167414720840e04_wp_, 3.97277109100414518365_wp_, & + 1.37790390235747998793e02_wp_, 1.17179220502086455287e02_wp_, & + 7.04831847180424675988e01_wp_,-1.20187763547154743238e01_wp_, & + -7.99243595776339741065_wp_,-2.99999894040324959612_wp_, & + 1.99999999999048104167_wp_/) +!---------------------------------------------------------------------- + x = arg + if (x == zero) then + ei = -xinf + else if ((x < zero)) then +!---------------------------------------------------------------------- +! calculate ei for negative argument or for e1. +!---------------------------------------------------------------------- + y = abs(x) + if (y <= one) then + sump = a(7) * y + a(1) + sumq = y + b(1) + do i = 2, 6 + sump = sump * y + a(i) + sumq = sumq * y + b(i) + end do + ei = (log(y) - sump / sumq ) * exp(y) + else if (y <= four) then + w = one / y + sump = c(1) + sumq = d(1) + do i = 2, 9 + sump = sump * w + c(i) + sumq = sumq * w + d(i) + end do + ei = - sump / sumq + else + w = one / y + sump = e(1) + sumq = f(1) + do i = 2, 10 + sump = sump * w + e(i) + sumq = sumq * w + f(i) + end do + ei = -w * (one - w * sump / sumq ) + end if + else if (x < six) then +!---------------------------------------------------------------------- +! to improve conditioning, rational approximations are expressed +! in terms of chebyshev polynomials for 0 <= x < 6, and in +! continued fraction form for larger x. +!---------------------------------------------------------------------- + t = x + x + t = t / three - two + px(1) = zero + qx(1) = zero + px(2) = p(1) + qx(2) = q(1) + do i = 2, 9 + px(i+1) = t * px(i) - px(i-1) + p(i) + qx(i+1) = t * qx(i) - qx(i-1) + q(i) + end do + sump = half * t * px(10) - px(9) + p(10) + sumq = half * t * qx(10) - qx(9) + q(10) + frac = sump / sumq + xmx0 = (x - x01/x11) - x02 + if (abs(xmx0) >= p037) then + ei = exp(-x) * ( log(x/x0) + xmx0 * frac ) + else +!---------------------------------------------------------------------- +! special approximation to ln(x/x0) for x close to x0 +!---------------------------------------------------------------------- + y = xmx0 / (x + x0) + ysq = y*y + sump = plg(1) + sumq = ysq + qlg(1) + do i = 2, 4 + sump = sump*ysq + plg(i) + sumq = sumq*ysq + qlg(i) + end do + ei = exp(-x) * (sump / (sumq*(x+x0)) + frac) * xmx0 + end if + else if (x < twelve) then + frac = zero + do i = 1, 9 + frac = s(i) / (r(i) + x + frac) + end do + ei = (r(10) + frac) / x + else if (x <= two4) then + frac = zero + do i = 1, 9 + frac = q1(i) / (p1(i) + x + frac) + end do + ei = (p1(10) + frac) / x + else + y = one / x + frac = zero + do i = 1, 9 + frac = q2(i) / (p2(i) + x + frac) + end do + frac = p2(10) + frac + ei = y + y * y * frac + end if + result = ei + end subroutine calcei3 + + subroutine calerf(arg,result,jintt) +!------------------------------------------------------------------ +! +! this packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x) +! for a real argument x. it contains three function type +! subprograms: erf, erfc, and erfcx (or derf, derfc, and derfcx), +! and one subroutine type subprogram, calerf. the calling +! statements for the primary entries are: +! +! y=erf(x) (or y=derf(x)), +! +! y=erfc(x) (or y=derfc(x)), +! and +! y=erfcx(x) (or y=derfcx(x)). +! +! the routine calerf is intended for internal packet use only, +! all computations within the packet being concentrated in this +! routine. the function subprograms invoke calerf with the +! statement +! +! call calerf(arg,result,jintt) +! +! where the parameter usage is as follows +! +! function parameters for calerf +! call arg result jintt +! +! erf(arg) any real argument erf(arg) 0 +! erfc(arg) abs(arg) < xbig erfc(arg) 1 +! erfcx(arg) xneg < arg < xmax erfcx(arg) 2 +! +!******************************************************************* +!******************************************************************* +! +! Explanation of machine-dependent constants +! +! XMIN = the smallest positive floating-point number. +! XINF = the largest positive finite floating-point number. +! XNEG = the largest negative argument acceptable to ERFCX; +! the negative of the solution to the equation +! 2*exp(x*x) = XINF. +! XSMALL = argument below which erf(x) may be represented by +! 2*x/sqrt(pi) and above which x*x will not underflow. +! A conservative value is the largest machine number X +! such that 1.0 + X = 1.0 to machine precision. +! XBIG = largest argument acceptable to ERFC; solution to +! the equation: W(x) * (1-0.5/x**2) = XMIN, where +! W(x) = exp(-x*x)/[x*sqrt(pi)]. +! XHUGE = argument above which 1.0 - 1/(2*x*x) = 1.0 to +! machine precision. A conservative value is +! 1/[2*sqrt(XSMALL)] +! XMAX = largest acceptable argument to ERFCX; the minimum +! of XINF and 1/[sqrt(pi)*XMIN]. +! +!******************************************************************* +!******************************************************************* +! +! error returns +! +! the program returns erfc = 0 for arg >= xbig; +! +! erfcx = xinf for arg < xneg; +! and +! erfcx = 0 for arg >= xmax. +! +! +! intrinsic functions required are: +! +! abs, aint, exp +! +! +! author: w. j. cody +! mathematics and computer science division +! argonne national laboratory +! argonne, il 60439 +! +! latest modification: march 19, 1990 +! +!------------------------------------------------------------------ + implicit none + real(wp_), intent(in) :: arg + real(wp_), intent(out) :: result + integer, intent(in) :: jintt + integer :: i + real(wp_) :: del,x,xden,xnum,y,ysq +!------------------------------------------------------------------ +! mathematical constants +!------------------------------------------------------------------ + real(wp_), parameter :: sqrpi=5.6418958354775628695e-1_wp_, & + thresh=0.46875_wp_ +!------------------------------------------------------------------ +! machine-dependent constants +!------------------------------------------------------------------ + real(wp_), parameter :: xinf=1.79e308_wp_, & ! ~huge + xneg=-26.628_wp_, & ! ? + xsmall=1.11e-16_wp_, & ! ~epsilon/2 + xbig=26.543_wp_, & ! ? + xhuge=6.71e7_wp_, & ! ~1/sqrt(epsilon) + xmax=2.53e307_wp_ ! ? +!------------------------------------------------------------------ +! coefficients for approximation to erf in first interval +!------------------------------------------------------------------ + real(wp_), dimension(5), parameter :: & + a=(/3.16112374387056560_wp_,1.13864154151050156e02_wp_, & + 3.77485237685302021e02_wp_,3.20937758913846947e03_wp_, & + 1.85777706184603153e-1_wp_/) + real(wp_), dimension(4), parameter :: & + b=(/2.36012909523441209e01_wp_,2.44024637934444173e02_wp_, & + 1.28261652607737228e03_wp_,2.84423683343917062e03_wp_/) +!------------------------------------------------------------------ +! coefficients for approximation to erfc in second interval +!------------------------------------------------------------------ + real(wp_), dimension(9), parameter :: & + c=(/5.64188496988670089e-1_wp_,8.88314979438837594_wp_, & + 6.61191906371416295e01_wp_,2.98635138197400131e02_wp_, & + 8.81952221241769090e02_wp_,1.71204761263407058e03_wp_, & + 2.05107837782607147e03_wp_,1.23033935479799725e03_wp_, & + 2.15311535474403846e-8_wp_/) + real(wp_), dimension(8), parameter :: & + d=(/1.57449261107098347e01_wp_,1.17693950891312499e02_wp_, & + 5.37181101862009858e02_wp_,1.62138957456669019e03_wp_, & + 3.29079923573345963e03_wp_,4.36261909014324716e03_wp_, & + 3.43936767414372164e03_wp_,1.23033935480374942e03_wp_/) +!------------------------------------------------------------------ +! coefficients for approximation to erfc in third interval +!------------------------------------------------------------------ + real(wp_), dimension(6), parameter :: & + p=(/3.05326634961232344e-1_wp_,3.60344899949804439e-1_wp_, & + 1.25781726111229246e-1_wp_,1.60837851487422766e-2_wp_, & + 6.58749161529837803e-4_wp_,1.63153871373020978e-2_wp_/) + real(wp_), dimension(5), parameter :: & + q=(/2.56852019228982242_wp_,1.87295284992346047_wp_, & + 5.27905102951428412e-1_wp_,6.05183413124413191e-2_wp_, & + 2.33520497626869185e-3_wp_/) +!------------------------------------------------------------------ + x = arg + y = abs(x) + if (y <= thresh) then +!------------------------------------------------------------------ +! evaluate erf for |x| <= 0.46875 +!------------------------------------------------------------------ + ysq = zero + if (y > xsmall) ysq = y * y + xnum = a(5)*ysq + xden = ysq + do i = 1, 3 + xnum = (xnum + a(i)) * ysq + xden = (xden + b(i)) * ysq + end do + result = x * (xnum + a(4)) / (xden + b(4)) + if (jintt /= 0) result = one - result + if (jintt == 2) result = exp(ysq) * result + return +!------------------------------------------------------------------ +! evaluate erfc for 0.46875 <= |x| <= 4.0 +!------------------------------------------------------------------ + else if (y <= four) then + xnum = c(9)*y + xden = y + do i = 1, 7 + xnum = (xnum + c(i)) * y + xden = (xden + d(i)) * y + end do + result = (xnum + c(8)) / (xden + d(8)) + if (jintt /= 2) then + ysq = aint(y*sixten)/sixten + del = (y-ysq)*(y+ysq) + result = exp(-ysq*ysq) * exp(-del) * result + end if +!------------------------------------------------------------------ +! evaluate erfc for |x| > 4.0 +!------------------------------------------------------------------ + else if (y < xbig .or. (y < xmax .and. jintt == 2)) then + ysq = one / (y * y) + xnum = p(6)*ysq + xden = ysq + do i = 1, 4 + xnum = (xnum + p(i)) * ysq + xden = (xden + q(i)) * ysq + end do + result = ysq *(xnum + p(5)) / (xden + q(5)) + result = (sqrpi - result) / y + if (jintt /= 2) then + ysq = aint(y*sixten)/sixten + del = (y-ysq)*(y+ysq) + result = exp(-ysq*ysq) * exp(-del) * result + end if + else if (y >= xhuge) then + result = sqrpi / y + else + result = zero + end if +!------------------------------------------------------------------ +! fix up for negative argument, erf, etc. +!------------------------------------------------------------------ + if (jintt == 0) then + result = (half - result) + half + if (x < zero) result = -result + else if (jintt == 1) then + if (x < zero) result = two - result + else + if (x < zero) then + if (x < xneg) then + result = xinf + else + ysq = aint(x*sixten)/sixten + del = (x-ysq)*(x+ysq) + y = exp(ysq*ysq) * exp(del) + result = (y+y) - result + end if + end if + end if + end subroutine calerf + + function derf(x) +!-------------------------------------------------------------------- +! +! this subprogram computes approximate values for erf(x). +! (see comments heading calerf). +! +! author/date: w. j. cody, january 8, 1985 +! +!-------------------------------------------------------------------- + implicit none + real(wp_) :: derf + real(wp_), intent(in) :: x + integer :: jintt + real(wp_) :: result +!------------------------------------------------------------------ + jintt = 0 + call calerf(x,result,jintt) + derf = result + end function derf + + function derfc(x) +!-------------------------------------------------------------------- +! +! this subprogram computes approximate values for erfc(x). +! (see comments heading calerf). +! +! author/date: w. j. cody, january 8, 1985 +! +!-------------------------------------------------------------------- + implicit none + real(wp_) :: derfc + real(wp_), intent(in) :: x + integer :: jintt + real(wp_) :: result +!------------------------------------------------------------------ + jintt = 1 + call calerf(x,result,jintt) + derfc = result + end function derfc + + function derfcx(x) +!------------------------------------------------------------------ +! +! this subprogram computes approximate values for exp(x*x) * erfc(x). +! (see comments heading calerf). +! +! author/date: w. j. cody, march 30, 1987 +! +!------------------------------------------------------------------ + implicit none + real(wp_) :: derfcx + real(wp_), intent(in) :: x + integer :: jintt + real(wp_) :: result +!------------------------------------------------------------------ + jintt = 2 + call calerf(x,result,jintt) + derfcx = result + end function derfcx + +end module eierf \ No newline at end of file diff --git a/src/grayl.f b/src/grayl.f index b85530a..1184f52 100644 --- a/src/grayl.f +++ b/src/grayl.f @@ -1651,256 +1651,6 @@ c end c c -c -c -c -c -c - subroutine calerf(arg,result,jint) -c------------------------------------------------------------------ -c -c this packet evaluates erf(x), erfc(x), and exp(x*x)*erfc(x) -c for a real argument x. it contains three function type -c subprograms: erf, erfc, and erfcx (or derf, derfc, and derfcx), -c and one subroutine type subprogram, calerf. the calling -c statements for the primary entries are: -c -c y=erf(x) (or y=derf(x)), -c -c y=erfc(x) (or y=derfc(x)), -c and -c y=erfcx(x) (or y=derfcx(x)). -c -c the routine calerf is intended for internal packet use only, -c all computations within the packet being concentrated in this -c routine. the function subprograms invoke calerf with the -c statement -c -c call calerf(arg,result,jint) -c -c where the parameter usage is as follows -c -c function parameters for calerf -c call arg result jint -c -c erf(arg) any real argument erf(arg) 0 -c erfc(arg) abs(arg) .lt. xbig erfc(arg) 1 -c erfcx(arg) xneg .lt. arg .lt. xmax erfcx(arg) 2 -c -c******************************************************************* -c -c error returns -c -c the program returns erfc = 0 for arg .ge. xbig; -c -c erfcx = xinf for arg .lt. xneg; -c and -c erfcx = 0 for arg .ge. xmax. -c -c -c intrinsic functions required are: -c -c abs, aint, exp -c -c -c author: w. j. cody -c mathematics and computer science division -c argonne national laboratory -c argonne, il 60439 -c -c latest modification: march 19, 1990 -c -c------------------------------------------------------------------ - integer i,jint - double precision - 1 a,arg,b,c,d,del,four,half,p,one,q,result,sixten,sqrpi, - 2 two,thresh,x,xbig,xden,xhuge,xinf,xmax,xneg,xnum,xsmall, - 3 y,ysq,zero - dimension a(5),b(4),c(9),d(8),p(6),q(5) -c------------------------------------------------------------------ -c mathematical constants -c------------------------------------------------------------------ - data four,one,half,two,zero/4.0d0,1.0d0,0.5d0,2.0d0,0.0d0/, - 1 sqrpi/5.6418958354775628695d-1/,thresh/0.46875d0/, - 2 sixten/16.0d0/ -c------------------------------------------------------------------ -c machine-dependent constants -c------------------------------------------------------------------ - data xinf,xneg,xsmall/1.79d308,-26.628d0,1.11d-16/, - 1 xbig,xhuge,xmax/26.543d0,6.71d7,2.53d307/ -c------------------------------------------------------------------ -c coefficients for approximation to erf in first interval -c------------------------------------------------------------------ - data a/3.16112374387056560d00,1.13864154151050156d02, - 1 3.77485237685302021d02,3.20937758913846947d03, - 2 1.85777706184603153d-1/ - data b/2.36012909523441209d01,2.44024637934444173d02, - 1 1.28261652607737228d03,2.84423683343917062d03/ -c------------------------------------------------------------------ -c coefficients for approximation to erfc in second interval -c------------------------------------------------------------------ - data c/5.64188496988670089d-1,8.88314979438837594d0, - 1 6.61191906371416295d01,2.98635138197400131d02, - 2 8.81952221241769090d02,1.71204761263407058d03, - 3 2.05107837782607147d03,1.23033935479799725d03, - 4 2.15311535474403846d-8/ - data d/1.57449261107098347d01,1.17693950891312499d02, - 1 5.37181101862009858d02,1.62138957456669019d03, - 2 3.29079923573345963d03,4.36261909014324716d03, - 3 3.43936767414372164d03,1.23033935480374942d03/ -c------------------------------------------------------------------ -c coefficients for approximation to erfc in third interval -c------------------------------------------------------------------ - data p/3.05326634961232344d-1,3.60344899949804439d-1, - 1 1.25781726111229246d-1,1.60837851487422766d-2, - 2 6.58749161529837803d-4,1.63153871373020978d-2/ - data q/2.56852019228982242d00,1.87295284992346047d00, - 1 5.27905102951428412d-1,6.05183413124413191d-2, - 2 2.33520497626869185d-3/ -c------------------------------------------------------------------ - x = arg - y = abs(x) - if (y .le. thresh) then -c------------------------------------------------------------------ -c evaluate erf for |x| <= 0.46875 -c------------------------------------------------------------------ - ysq = zero - if (y .gt. xsmall) ysq = y * y - xnum = a(5)*ysq - xden = ysq - do 20 i = 1, 3 - xnum = (xnum + a(i)) * ysq - xden = (xden + b(i)) * ysq - 20 continue - result = x * (xnum + a(4)) / (xden + b(4)) - if (jint .ne. 0) result = one - result - if (jint .eq. 2) result = exp(ysq) * result - go to 800 -c------------------------------------------------------------------ -c evaluate erfc for 0.46875 <= |x| <= 4.0 -c------------------------------------------------------------------ - else if (y .le. four) then - xnum = c(9)*y - xden = y - do 120 i = 1, 7 - xnum = (xnum + c(i)) * y - xden = (xden + d(i)) * y - 120 continue - result = (xnum + c(8)) / (xden + d(8)) - if (jint .ne. 2) then - ysq = aint(y*sixten)/sixten - del = (y-ysq)*(y+ysq) - result = exp(-ysq*ysq) * exp(-del) * result - end if -c------------------------------------------------------------------ -c evaluate erfc for |x| > 4.0 -c------------------------------------------------------------------ - else - result = zero - if (y .ge. xbig) then - if ((jint .ne. 2) .or. (y .ge. xmax)) go to 300 - if (y .ge. xhuge) then - result = sqrpi / y - go to 300 - end if - end if - ysq = one / (y * y) - xnum = p(6)*ysq - xden = ysq - do 240 i = 1, 4 - xnum = (xnum + p(i)) * ysq - xden = (xden + q(i)) * ysq - 240 continue - result = ysq *(xnum + p(5)) / (xden + q(5)) - result = (sqrpi - result) / y - if (jint .ne. 2) then - ysq = aint(y*sixten)/sixten - del = (y-ysq)*(y+ysq) - result = exp(-ysq*ysq) * exp(-del) * result - end if - end if -c------------------------------------------------------------------ -c fix up for negative argument, erf, etc. -c------------------------------------------------------------------ - 300 if (jint .eq. 0) then - result = (half - result) + half - if (x .lt. zero) result = -result - else if (jint .eq. 1) then - if (x .lt. zero) result = two - result - else - if (x .lt. zero) then - if (x .lt. xneg) then - result = xinf - else - ysq = aint(x*sixten)/sixten - del = (x-ysq)*(x+ysq) - y = exp(ysq*ysq) * exp(del) - result = (y+y) - result - end if - end if - end if - 800 return -c---------- last card of calerf ---------- - end -c - double precision function derf(x) -c-------------------------------------------------------------------- -c -c this subprogram computes approximate values for erf(x). -c (see comments heading calerf). -c -c author/date: w. j. cody, january 8, 1985 -c -c-------------------------------------------------------------------- - integer jint - double precision x, result -c------------------------------------------------------------------ - jint = 0 - call calerf(x,result,jint) - derf = result - return -c---------- last card of derf ---------- - end -c - double precision function derfc(x) -c-------------------------------------------------------------------- -c -c this subprogram computes approximate values for erfc(x). -c (see comments heading calerf). -c -c author/date: w. j. cody, january 8, 1985 -c -c-------------------------------------------------------------------- - integer jint - double precision x, result -c------------------------------------------------------------------ - jint = 1 - call calerf(x,result,jint) - derfc = result - return -c---------- last card of derfc ---------- - end -c - double precision function derfcx(x) -c------------------------------------------------------------------ -c -c this subprogram computes approximate values for exp(x*x) * erfc(x). -c (see comments heading calerf). -c -c author/date: w. j. cody, march 30, 1987 -c -c------------------------------------------------------------------ - double precision x, result - integer jint -c------------------------------------------------------------------ - jint = 2 - call calerf(x,result,jint) - derfcx = result - return -c---------- last card of derfcx ---------- - end -c -c c Integration routine dqags.f from quadpack and dependencies: BEGIN c c