diff --git a/Makefile b/Makefile index a4080fe..0ffe52b 100644 --- a/Makefile +++ b/Makefile @@ -3,7 +3,7 @@ EXE=gray # Objects list MAINOBJ=gray.o -OTHOBJ=dqagmv.o grayl.o reflections.o green_func_p.o \ +OTHOBJ=dispersion.o calcei_mod.o dqagmv.o grayl.o reflections.o green_func_p.o \ const_and_precisions.o # Alternative search paths @@ -23,9 +23,10 @@ $(EXE): $(MAINOBJ) $(OTHOBJ) $(FC) $(FFLAGS) -o $@ $^ # Dependencies on modules -gray.o: dqagmv.o green_func_p.o reflections.o const_and_precisions.o +gray.o: dispersion.o dqagmv.o green_func_p.o reflections.o const_and_precisions.o green_func_p.o: const_and_precisions.o reflections.o: const_and_precisions.o +dispersion.o: calcei_mod.o dqagmv.o # General object compilation command %.o: %.f90 diff --git a/src/calcei_mod.f90 b/src/calcei_mod.f90 new file mode 100644 index 0000000..4413918 --- /dev/null +++ b/src/calcei_mod.f90 @@ -0,0 +1,608 @@ +module calcei_mod + + 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 new file mode 100644 index 0000000..b41e087 --- /dev/null +++ b/src/dispersion.f90 @@ -0,0 +1,1146 @@ +module dispersion +! + use calcei_mod + implicit none + integer, parameter :: dp=kind(1.0d0) + integer, parameter :: npts=500 + real(dp), save :: ttv(npts+1),extv(npts+1) + real(dp), parameter :: tmax=5.0d0,dtex=2.0d0*tmax/dble(npts) +! +contains +! +subroutine colddisp(xg,yg,npl,nprf,sox) +! +! solution cold dispersion relation +! + implicit none +! + real(dp) :: xg ! X=omegap^2/omega^2 + real(dp) :: yg ! Y=omegac/omega + real(dp) :: npl ! N parallel to B + real(dp) :: nprf ! N perpendicular to B (cold) + real(dp) :: sox ! sox=-1 ==> O mode, sox=1 ==> X mode +! + real(dp) :: yg2,npl2,dnl,del,dxg +! + npl2=npl*npl + dnl=1.0d0-npl2 + dxg=1.0d0-xg + yg2=yg*yg + del=sqrt(dnl*dnl+4.0d0*npl2*dxg/yg2) + nprf=sqrt(dxg-npl2-xg*yg2*(1.0d0+npl2+sox*del)/(dxg-yg2)/2.0d0) +! +end subroutine colddisp +! +! +! +subroutine harmnumber(yg,mu,npl,nhmin,nhmax,iwr) +! +! computation of minimum and maximum harmonic +! + implicit none +! + real(dp), parameter :: expcr=16.d0 ! maximum value for mu*(gamma-1) above + ! which the distribution function is + ! considered to be 0 + real(dp), parameter :: eps=1.d-8 ! small number to have a correct rounding + ! when ygnc/yg is an integer + + real(dp) :: yg,ygc ! Y=omegac/omega + integer :: nhmin,nhmax ! n=number of the armonic + real(dp) :: ygn ! n*Y + real(dp) :: npl,npl2 ! parallel N + real(dp) :: gg ! gamma=sqrt(1+u^2) + real(dp) :: dnl,rdu2 + real(dp) :: argexp + real(dp) :: mu ! mu=mc^2/Te + integer :: iwr ! weakly (iwr=1) or fully relativistic approximation + real(dp) :: u1,u2,uu2,g1,g2 +!---------------------------------------- + integer :: nh,nhc +!---------------------------------------- +! + nhmin=0 + nhmax=0 + npl2=npl**2 + dnl=1.0d0-npl2 +! + if(iwr.eq.1) then + ygc=max(1.0d0-0.5d0*npl2,0.0d0) + else + ygc=sqrt(max(dnl,0.0d0)) + end if + nhc=int(ygc/yg) + if (nhc*ygimx ',yg,errnpr,i +! + if(sqrt(dble(npr2)).lt.0.0d0.or.npr2.ne.npr2.or.abs(npr2).eq.huge(one).or.abs(npr2).le.tiny(one)) then + write(*,"(' X =',f7.4,' Y =',f7.4,' Nperp =',f7.4,'!')") xg,yg,sqrt(abs(npr2)) + npr2=(0.0d0,0.0d0) + end if +! if(dble(npr2).lt.0.0d0) then +! npr2=0.0d0 +! print*,' Y =',yg,' npr2 < 0' +! err=99 +! end if +! +! write(11,99) yg,dble(npr2),dimag(npr2),nprf**2,dble(i) +! + npr=sqrt(npr2) + nprr=dble(npr) + npri=dimag(npr) +! + ex=dcmplx(0.0d0,0.0d0) + ey=dcmplx(0.0d0,0.0d0) + ez=dcmplx(0.0d0,0.0d0) +! + if (abs(npl).gt.1.0d-6) then + den=e12*e23-(e13+npr*npl)*(e22-npr2-npl2) + ey=-(e12*(e13+npr*npl)+(e11-npl2)*e23)/den + ez=(e12*e12+(e22-npr2-npl2)*(e11-npl2))/den + ez2=abs(ez)**2 + ey2=abs(ey)**2 + enx2=1.0d0/(1.0d0+ey2+ez2) + ex=dcmplx(sqrt(enx2),0.0d0) + ez2=ez2*enx2 + ey2=ey2*enx2 + ez=ez*ex + ey=ey*ex + else + if(sox.lt.0.0d0) then + ez=dcmplx(1.0d0,0.0d0) + ez2=abs(ez)**2 + else + ex2=1.0d0/(1.0d0+abs(-e11/e12)**2) + ex=sqrt(ex2) + ey=-ex*e11/e12 + ey2=abs(ey)**2 + ez2=0.0d0 + end if + end if +! +end subroutine warmdisp +! +! +! Fully relativistic case +! computation of dielectric tensor elements +! up to third order in Larmor radius for hermitian part +subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast) + implicit none + integer :: i,j,l,is,k,lm,fast + real(dp) :: xg,yg,mu,npl,npl2,dnl,cr,ci,w + real(dp) :: asl,bsl,cmxw,fal + integer :: lrm + complex(dp) :: e330,epsl(3,3,lrm) + complex(dp) :: ca11,ca12,ca22,ca13,ca23,ca33 + complex(dp) :: cq0p,cq0m,cq1p,cq1m,cq2p + real(dp), parameter :: pi=3.14159265358979d0,rpi=1.77245385090552d0 + complex(dp), parameter :: ui=(0.0d0,1.0d0) + real(dp), dimension(-lrm:lrm,0:2,0:lrm) :: rr + real(dp), dimension(lrm,0:2,lrm) :: ri +! + npl2=npl**2 + dnl=1.0d0-npl2 +! + cmxw=1.0d0+15.0d0/(8.0d0*mu)+105.0d0/(128.0d0*mu**2) + cr=-mu*mu/(rpi*cmxw) + ci=sqrt(2.0d0*pi*mu)*mu**2/cmxw +! + do l=1,lrm + do j=1,3 + do i=1,3 + epsl(i,j,l)=dcmplx(0.0d0,0.0d0) + end do + end do + end do +! +! call hermitian(rr,yg,mu,npl,cr,ci,fast,lrm) + select case(fast) + + case(2:3) + call hermitian(rr,yg,mu,npl,cr,ci,fast,lrm) + + case(4:) + call hermitian_2(rr,yg,mu,npl,cr,ci,fast,lrm) + + case default + write(*,*) "unexpected value for flag 'fast' in dispersion:", fast + + end select +! + call antihermitian(ri,yg,mu,npl,cr,ci,lrm) +! + do l=1,lrm + lm=l-1 + fal=-0.25d0**l*fact(2*l)/(fact(l)**2*yg**(2*lm)) + ca11=(0.0d0,0.0d0) + ca12=(0.0d0,0.0d0) + ca13=(0.0d0,0.0d0) + ca22=(0.0d0,0.0d0) + ca23=(0.0d0,0.0d0) + ca33=(0.0d0,0.0d0) + do is=0,l + k=l-is + w=(-1.0d0)**k +! + asl=w/(fact(is+l)*fact(l-is)) + bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0d0*l-1.0d0)) +! + if(is.gt.0) then + cq0p=rr(is,0,l)+rr(-is,0,l)+ui*ri(is,0,l) + cq0m=rr(is,0,l)-rr(-is,0,l)+ui*ri(is,0,l) + cq1p=rr(is,1,l)+rr(-is,1,l)+ui*ri(is,1,l) + cq1m=rr(is,1,l)-rr(-is,1,l)+ui*ri(is,1,l) + cq2p=rr(is,2,l)+rr(-is,2,l)+ui*ri(is,2,l) + else + cq0p=rr(is,0,l) + cq0m=rr(is,0,l) + cq1p=rr(is,1,l) + cq1m=rr(is,1,l) + cq2p=rr(is,2,l) + end if +! + ca11=ca11+is**2*asl*cq0p + ca12=ca12+is*l*asl*cq0m + ca22=ca22+bsl*cq0p + ca13=ca13+is*asl*cq1m/yg + ca23=ca23+l*asl*cq1p/yg + ca33=ca33+asl*cq2p/yg**2 + end do + epsl(1,1,l) = - xg*ca11*fal + epsl(1,2,l) = + ui*xg*ca12*fal + epsl(2,2,l) = - xg*ca22*fal + epsl(1,3,l) = - xg*ca13*fal + epsl(2,3,l) = - ui*xg*ca23*fal + epsl(3,3,l) = - xg*ca33*fal + end do +! + cq2p=rr(0,2,0) + e330=1.0d0+xg*cq2p +! + epsl(1,1,1) = 1.d0 + epsl(1,1,1) + epsl(2,2,1) = 1.d0 + epsl(2,2,1) +! + do l=1,lrm + epsl(2,1,l) = - epsl(1,2,l) + epsl(3,1,l) = epsl(1,3,l) + epsl(3,2,l) = - epsl(2,3,l) + end do +! +! + return + end subroutine diel_tens_fr +! +! +! + subroutine hermitian(rr,yg,mu,npl,cr,ci,fast,lrm) + implicit none + real(dp) :: yg,mu,mu2,mu4,mu6,npl,npl2,npl4,cr,ci + real(dp), dimension(-lrm:lrm,0:2,0:lrm) :: rr + real(dp) :: dt,bth,bth2,bth4,bth6,t,rxt,upl2,upl,gx,exdx + real(dp) :: x,gr,s,zm,zm2,zm3,fe0m,ffe,sy1,sy2,sy3 + integer :: i,k,n,n1,nn,m,lrm,llm,fast +! + do n=-lrm,lrm + do k=0,2 + do m=0,lrm + rr(n,k,m)=0.0d0 + end do + end do + end do +! + llm=min(3,lrm) +! + bth2=2.0d0/mu + bth=sqrt(bth2) + mu2=mu*mu + mu4=mu2*mu2 + mu6=mu4*mu2 +! + do i = 1, npts+1 + t = ttv(i) + rxt=sqrt(1.0d0+t*t/(2.0d0*mu)) + x = t*rxt + upl2=bth2*x**2 + upl=bth*x + gx=1.0d0+t*t/mu + exdx=cr*extv(i)*gx/rxt*dtex +! + n1=1 + if(fast.gt.2) n1=-llm +! + do n=n1,llm + nn=abs(n) + gr=npl*upl+n*yg + zm=-mu*(gx-gr) + s=mu*(gx+gr) + zm2=zm*zm + zm3=zm2*zm + call calcei3(zm,fe0m) +! + do m=nn,llm + if(n.eq.0.and.m.eq.0) then + rr(0,2,0) = rr(0,2,0) - exdx*fe0m*upl2 + else + if (m.eq.1) then + ffe=(1.0d0+s*(1.0d0-zm*fe0m))/mu2 + else if (m.eq.2) then + ffe=(6.0d0-2.0d0*zm+4.0d0*s+s*s*(1.0d0+zm-zm2*fe0m))/mu4 + else + ffe=(18.0d0*s*(s+4.0d0-zm)+6.0d0* & + (20.0d0-8.0d0*zm+zm2)+s**3*(2.0d0+zm+zm2-zm3*fe0m))/mu6 + end if +! + rr(n,0,m) = rr(n,0,m) + exdx*ffe + rr(n,1,m) = rr(n,1,m) + exdx*ffe*upl + rr(n,2,m) = rr(n,2,m) + exdx*ffe*upl2 +! + end if +! + end do + end do + end do +! + if(fast.gt.2) return +! + sy1=1.0d0+yg + sy2=1.0d0+yg*2.0d0 + sy3=1.0d0+yg*3.0d0 +! + bth4=bth2*bth2 + bth6=bth4*bth2 +! + npl2=npl*npl + npl4=npl2*npl2 +! + rr(0,2,0) = -(1.0d0+bth2*(-1.25d0+1.5d0*npl2) & + +bth4*(1.71875d0-6.0d0*npl2+3.75d0*npl2*npl2) & + +bth6*3.0d0*(-65.0d0+456.0d0*npl2-660.d0*npl4 & + + 280.0d0*npl2*npl4)/64.0d0 & + + bth6*bth2*15.0d0*(252.853d3-2850.816d3*npl2 & + +6942.720d3*npl4-6422.528d3*npl4*npl2 & + +2064.384d3*npl4*npl4)/524.288d3) + rr(0,1,1) = -npl*bth2*(1.0d0+bth2*(-2.25d0+1.5d0*npl2) & + +bth4*9.375d-2*(6.1d1-9.6d1*npl2+4.d1*npl4 & + +bth2*(-184.5d0+4.92d2*npl2-4.5d2*npl4 & + +1.4d2*npl2*npl4))) + rr(0,2,1) = -bth2*(1.0d0+bth2*(-0.5d0+1.5d0*npl2) & + +0.375d0*bth4*(3.d0-15.d0*npl2+10.d0*npl4) + & + 3.d0*bth6*(-61.d0+471.d0*npl2-680*npl4+280.d0*npl2*npl4)/64.d0) + rr(-1,0,1) = -2.0d0/sy1*(1.0d0+bth2/sy1*(-1.25d0+0.5d0*npl2/sy1) & + +bth4/sy1*(-0.46875d0+(2.1875d0+0.625d0*npl2)/sy1 & + -2.625d0*npl2/sy1**2+0.75d0*npl4/sy1**3) + bth6/sy1* & + (0.234375d0+(1.640625d0+0.234375d0*npl2)/sy1 + & + (-4.921875d0-4.921875d0*npl2)/sy1**2 + & + 2.25d0*npl2*(5.25d0 + npl2)/sy1**3 - 8.4375d0*npl4/sy1**4 + & + 1.875d0*npl2*npl4/sy1**5)+bth6*bth2/sy1*(0.019826889038*sy1 & + -0.06591796875d0+(-0.7177734375d0 - 0.1171875d0*npl2)/sy1+ & + (-5.537109375d0 - 2.4609375d0*npl2)/sy1**2 + & + (13.53515625d0 + 29.53125d0*npl2 + 2.8125d0*npl4)/sy1**3 + & + (-54.140625d0*npl2 - 32.6953125d0*npl4)/sy1**4 + & + (69.609375d0*npl4 + 9.84375d0*npl2*npl4)/sy1**5 - & + 36.09375d0*npl2*npl4/sy1**6 + 6.5625d0*npl4**2/sy1**7)) + rr(-1,1,1) = -npl*bth2/sy1**2*(1.0d0+bth2*(1.25d0-3.5d0/sy1+ & + 1.5d0*npl2/sy1**2)+bth4*9.375d-2*((5.0d0-7.d1/sy1+ & + (126.d0+48.d0*npl2)/sy1**2-144.d0*npl2/sy1**3+4.d1*npl4/sy1**4)+ & + bth2*(-2.5d0-3.5d1/sy1+(3.15d2+6.d1*npl2)/sy1**2+ & + (-4.62d2-5.58d2*npl2)/sy1**3+(9.9d2*npl2+2.1d2*npl4)/sy1**4 & + -6.6d2*npl4/sy1**5+1.4d2*npl4*npl2/sy1**6))) + rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(1.25d0-1.75d0/sy1+1.5d0*npl2/sy1**2)+ & + bth4*3.d0/32.d0*(5.d0-35.d0/sy1+(42.d0+48.d0*npl2)/sy1**2- & + 108.d0*npl2/sy1**3+40.0d0*npl4/sy1**4 + & + 0.5d0*bth2*(-5.d0-35.d0/sy1+(21.d1+12.d1*npl2)/sy1**2- & + (231.d0+837.d0*npl2)/sy1**3+12.d0*npl2*(99.d0+35.d0*npl2)/sy1**4 - & + 1100.d0*npl4/sy1**5 + 280.d0*npl2*npl4/sy1**6))) +! + if(llm.gt.1) then +! + rr(0,0,2) = -4.0d0*bth2*(1.0d0+bth2*(-0.5d0+0.5d0*npl2)+ & + bth4*(1.125d0-1.875d0*npl2+0.75d0*npl4)+bth6* & + 3.0d0*(-61.d0+157.d0*npl2-136.d0*npl4+40.d0*npl2*npl4)/64.d0) + rr(0,1,2) = -2.0d0*npl*bth4*(1.0d0+bth2*(-1.5d0+1.5d0*npl2)+ & + bth4*(39.d0-69.d0*npl2+30.0d0*npl4)/8.0d0) + rr(0,2,2) = -2.0d0*bth4*(1.0d0+bth2*(0.75d0+1.5d0*npl2)+ bth4* & + (13.d0-48.d0*npl2 +40.0d0*npl4)*3.0d0/32.0d0) + rr(-1,0,2) = -4.0d0*bth2/sy1*(1.0d0+bth2* & + (1.25d0-1.75d0/sy1+0.5d0*npl2/sy1**2)+bth4* & + (0.46875d0-3.28125d0/sy1+(3.9375d0+1.5d0*npl2)/sy1**2 & + -3.375d0*npl2/sy1**3+0.75d0*npl4/sy1**4)+bth4*bth2*3.0d0/64.d0* & + (-5.0d0-35.0d0/sy1+(210.d0+40.d0*npl2)/sy1**2-3.0d0*(77.d0+93.d0*npl2)/sy1**3+ & + (396.0*npl2+84.d0*npl4)/sy1**4-220.d0*npl4/sy1**5+40.d0*npl4*npl2/sy1**6)) + rr(-1,1,2) = -2.0d0*bth4*npl/sy1**2*(1.0d0+bth2* & + (3.0d0-4.5d0/sy1+1.5d0*npl2/sy1**2) + bth4* & + (20.d0-93.d0/sy1+(99.d0+42.d0*npl2)/sy1**2- & + 88.d0*npl2/sy1**3+20.d0*npl4/sy1**4)*3.0d0/16.0d0) + rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+bth2* & + (3.0d0-2.25d0/sy1+1.5d0*npl2/sy1**2)+ bth4* & + (40.d0*npl4/sy1**4-132.0d0*npl2/sy1**3+ & + (66.d0+84.d0*npl2)/sy1**2-93.d0/sy1+40.0d0)*3.0d0/32.0d0) + rr(-2,0,2) = -4.0d0*bth2/sy2*(1.0d0+bth2* & + (1.25d0-1.75d0/sy2+0.5d0*npl2/sy2**2)+bth4* & + (0.46875d0-3.28125d0/sy2+(3.9375d0+1.5d0*npl2)/sy2**2 & + -3.375d0*npl2/sy2**3+0.75d0*npl4/sy2**4)+bth4*bth2*3.0d0/64.d0* & + (-5.0d0-35.0d0/sy2+(210.d0+40.d0*npl2)/sy2**2-3.0d0*(77.d0+93.d0*npl2)/sy2**3+ & + (396.0*npl2+84.d0*npl4)/sy2**4-220.d0*npl4/sy2**5+40.d0*npl4*npl2/sy2**6)) + rr(-2,1,2) =-2.0d0*bth4*npl/sy2**2*(1.0d0+bth2* & + (3.0d0-4.5d0/sy2+1.5d0*npl2/sy2**2) + bth4* & + (20.d0-93.d0/sy2+(99.d0+42.d0*npl2)/sy2**2- & + 88.d0*npl2/sy2**3+20.d0*npl4/sy2**4)*3.0d0/16.0d0) + rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+bth2* & + (3.0d0-2.25d0/sy2+1.5d0*npl2/sy2**2) + bth4* & + (40.d0*npl4/sy2**4-132.0d0*npl2/sy2**3+ & + (66.d0+84.d0*npl2)/sy2**2-93.d0/sy2+40.0d0)*3.0d0/32.0d0) +! + if(llm.gt.2) then +! + rr(0,0,3) = -12.0d0*bth4*(1.0d0+bth2*(0.75d0+0.5d0*npl2)+bth4* & + (1.21875d0-1.5d0*npl2+0.75d0*npl2*npl2)) + rr(0,1,3) = -6.0d0*npl*bth6*(1+bth2*(-0.25d0+1.5d0*npl2)) + rr(0,2,3) = -6.0d0*bth6*(1.0d0+bth2*(2.5d0+1.5d0*npl2)) + rr(-1,0,3) = -12.0d0*bth4/sy1* & + (1.0d0+bth2*(3.0d0-2.25d0/sy1+0.5d0*npl2/sy1**2)+ & + bth4*(3.75d0-8.71875d0/sy1+(6.1875d0+2.625d0*npl2) & + /sy1**2-4.125d0*npl2/sy1**3+0.75*npl2*npl2/sy1**4)) + rr(-1,1,3) = -6.0d0*npl*bth6/sy1**2* & + (1.0d0+bth2*(5.25d0-5.5d0/sy1+1.5d0*npl2/sy1**2)) + rr(-1,2,3) = -6.0d0*bth6/sy1* & + (1.0d0+bth2*(5.25d0-2.75d0/sy1+1.5d0*npl2/sy1**2)) +! + rr(-2,0,3) = -12.0d0*bth4/sy2* & + (1.0d0+bth2*(3.0d0-2.25d0/sy2+0.5d0*npl2/sy2**2)+ & + bth4*(3.75d0-8.71875d0/sy2+(6.1875d0+2.625d0*npl2) & + /sy2**2-4.125d0*npl2/sy2**3+0.75*npl2*npl2/sy2**4)) + rr(-2,1,3) = -6.0d0*npl*bth6/sy2**2* & + (1.0d0+bth2*(5.25d0-5.5d0/sy2+1.5d0*npl2/sy2**2)) + rr(-2,2,3) = -6.0d0*bth6/sy2* & + (1.0d0+bth2*(5.25d0-2.75d0/sy2+1.5d0*npl2/sy2**2)) +! + rr(-3,0,3) = -12.0d0*bth4/sy3* & + (1.0d0+bth2*(3.0d0-2.25d0/sy3+0.5d0*npl2/sy3**2)+ & + bth4*(3.75d0-8.71875d0/sy3+(6.1875d0+2.625d0*npl2) & + /sy3**2-4.125d0*npl2/sy3**3+0.75*npl2*npl2/sy3**4)) + rr(-3,1,3) = -6.0d0*npl*bth6/sy3**2* & + (1.0d0+bth2*(5.25d0-5.5d0/sy3+1.5d0*npl2/sy3**2)) + rr(-3,2,3) = -6.0d0*bth6/sy3* & + (1.0d0+bth2*(5.25d0-2.75d0/sy3+1.5d0*npl2/sy3**2)) + +! + end if + end if +! + end subroutine hermitian +! +! + subroutine hermitian_2(rr,yg,mu,npl,cr,ci,fast,lrm) + implicit none +! parameters + integer,parameter :: lw=5000,liw=lw/4 + integer, parameter :: npar=7 + real(dp), parameter :: epsa=0.0d0,epsr=1.0d-4 +! arguments + integer :: lrm,fast + real(dp) :: yg,mu,npl,cr,ci + real(dp), dimension(-lrm:lrm,0:2,0:lrm) :: rr +! internal + integer :: n,m,ih,i,k,n1,nn,llm,neval,ier,last,ihmin + integer, dimension(liw) :: iw + real(dp) :: mu2,mu4,mu6,npl2,dt,bth,bth2,bth4,bth6 + real(dp) :: t,x,s,sy1,sy2,sy3,resfh,epp + real(dp), dimension(lw) :: w + real(dp), dimension(npar) :: apar +! + do n=-lrm,lrm + do k=0,2 + do m=0,lrm + rr(n,k,m)=0.0d0 + end do + end do + end do +! + llm=min(3,lrm) +! + bth2=2.0d0/mu + bth=sqrt(bth2) + mu2=mu*mu + mu4=mu2*mu2 + mu6=mu4*mu2 +! + n1=1 + if(fast.gt.10) n1=-llm +! + apar(1) = yg + apar(2) = mu + apar(3) = npl + apar(4) = cr +! + do n=n1,llm + nn=abs(n) + apar(5) = real(n,dp) + do m=nn,llm + apar(6) = real(m,dp) + ihmin=0 + if(n.eq.0.and.m.eq.0) ihmin=2 + do ih=ihmin,2 + apar(7) = real(ih,dp) +! call dqagi(fhermit,bound,2,epsa,epsr,resfh, + call dqagsmv(fhermit,-tmax,tmax,apar,npar,epsa,epsr,resfh, & + epp,neval,ier,liw,lw,last,iw,w) + rr(n,ih,m) = resfh + end do + end do + end do + + if(fast.gt.10) return +! + sy1=1.0d0+yg + sy2=1.0d0+yg*2.0d0 + sy3=1.0d0+yg*3.0d0 +! + bth4=bth2*bth2 + bth6=bth4*bth2 +! + npl2=npl*npl +! + rr(0,2,0) = -(1.0d0+bth2*(-1.25d0+1.5d0*npl2) & + +bth4*(1.71875d0-6.0d0*npl2+3.75d0*npl2*npl2)) + rr(0,1,1) = -npl*bth2*(1.0d0+bth2*(-2.25d0+1.5d0*npl2)) + rr(0,2,1) = -bth2*(1.0d0+bth2*(-0.5d0+1.5d0*npl2)) + rr(-1,0,1) = -2.0d0/sy1*(1.0d0+bth2/sy1*(-1.25d0+0.5d0*npl2 & + /sy1)+bth4/sy1*(-0.46875d0+(2.1875d0+0.625d0*npl2)/ & + sy1-2.625d0*npl2/sy1**2+0.75d0*npl2*npl2/sy1**3)) + rr(-1,1,1) = -npl*bth2/sy1**2*(1.0d0+bth2*(1.25d0-3.5d0/sy1+ & + 1.5d0*npl2/sy1**2)) + rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(1.25d0-1.75d0/sy1+1.5d0* & + npl2/sy1**2)) +! + if(llm.gt.1) then +! + rr(0,0,2) = -4.0d0*bth2*(1.0d0+bth2*(-0.5d0+0.5d0*npl2)+ & + bth4*(1.125d0-1.875d0*npl2+0.75d0*npl2*npl2)) + rr(0,1,2) = -2.0d0*npl*bth4*(1.0d0+bth2*(-1.5d0+1.5d0*npl2)) + rr(0,2,2) = -2.0d0*bth4*(1.0d0+bth2*(0.75d0+1.5d0*npl2)) + rr(-1,0,2) = -4.0d0*bth2/sy1*(1.0d0+bth2* & + (1.25d0-1.75d0/sy1+0.5d0*npl2/sy1**2)+bth4* & + (0.46875d0-3.28125d0/sy1+(3.9375d0+1.5d0*npl2)/sy1**2 & + -3.375d0*npl2/sy1**3+0.75d0*npl2*npl2/sy1**4)) + rr(-1,1,2) = -2.0d0*bth4*npl/sy1**2*(1.0d0+bth2* & + (3.0d0-4.5d0/sy1+1.5d0*npl2/sy1**2)) + rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+bth2* & + (3.0d0-2.25d0/sy1+1.5d0*npl2/sy1**2)) + rr(-2,0,2) = -4.0d0*bth2/sy2*(1.0d0+bth2* & + (1.25d0-1.75d0/sy2+0.5d0*npl2/sy2**2)+bth4* & + (0.46875d0-3.28125d0/sy2+(3.9375d0+1.5d0*npl2)/sy2**2 & + -3.375d0*npl2/sy2**3+0.75d0*npl2*npl2/sy2**4)) + rr(-2,1,2) = -2.0d0*bth4*npl/sy2**2*(1.0d0+bth2* & + (3.0d0-4.5d0/sy2+1.5d0*npl2/sy2**2)) + rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+bth2* & + (3.0d0-2.25d0/sy2+1.5d0*npl2/sy2**2)) +! + if(llm.gt.2) then +! + rr(0,0,3) = -12.0d0*bth4*(1+bth2*(0.75d0+0.5d0*npl2)+bth4* & + (1.21875d0-1.5d0*npl2+0.75d0*npl2*npl2)) + rr(0,1,3) = -6.0d0*npl*bth6*(1+bth2*(-0.25d0+1.5d0*npl2)) + rr(0,2,3) = -6.0d0*bth6*(1+bth2*(2.5d0+1.5d0*npl2)) + rr(-1,0,3) = -12.0d0*bth4/sy1*& + (1.0d0+bth2*(3.0d0-2.25d0/sy1+0.5d0*npl2/sy1**2)+ & + bth4*(3.75d0-8.71875d0/sy1+(6.1875d0+2.625d0*npl2) & + /sy1**2-4.125d0*npl2/sy1**3+0.75d0*npl2*npl2/sy1**4)) + rr(-1,1,3) = -6.0d0*npl*bth6/sy1**2*& + (1.0d0+bth2*(5.25d0-5.5d0/sy1+1.5d0*npl2/sy1**2)) + rr(-1,2,3) = -6.0d0*bth6/sy1*& + (1.0d0+bth2*(5.25d0-2.75d0/sy1+1.5d0*npl2/sy1**2)) +! + rr(-2,0,3) = -12.0d0*bth4/sy2*& + (1.0d0+bth2*(3.0d0-2.25d0/sy2+0.5d0*npl2/sy2**2)+ & + bth4*(3.75d0-8.71875d0/sy2+(6.1875d0+2.625d0*npl2) & + /sy2**2-4.125d0*npl2/sy2**3+0.75d0*npl2*npl2/sy2**4)) + rr(-2,1,3) = -6.0d0*npl*bth6/sy2**2* & + (1.0d0+bth2*(5.25d0-5.5d0/sy2+1.5d0*npl2/sy2**2)) + rr(-2,2,3) = -6.0d0*bth6/sy2* & + (1.0d0+bth2*(5.25d0-2.75d0/sy2+1.5d0*npl2/sy2**2)) +! + rr(-3,0,3) = -12.0d0*bth4/sy3* & + (1.0d0+bth2*(3.0d0-2.25d0/sy3+0.5d0*npl2/sy3**2)+ & + bth4*(3.75d0-8.71875d0/sy3+(6.1875d0+2.625d0*npl2) & + /sy3**2-4.125d0*npl2/sy3**3+0.75d0*npl2*npl2/sy3**4)) + rr(-3,1,3) = -6.0d0*npl*bth6/sy3**2* & + (1.0d0+bth2*(5.25d0-5.5d0/sy3+1.5d0*npl2/sy3**2)) + rr(-3,2,3) = -6.0d0*bth6/sy3* & + (1.0d0+bth2*(5.25d0-2.75d0/sy3+1.5d0*npl2/sy3**2)) +! + end if + end if +! + end subroutine hermitian_2 +! +! + function fhermit(t,apar,npar) +! + implicit none +! arguments + integer, intent(in) :: npar + real(dp), intent(in) :: t + real(dp) :: fhermit + real(dp), dimension(npar), intent(in) :: apar +! internal + integer :: n,m,ih + integer nn + real(dp) :: yg,mu,npl,cr + real(dp) :: mu2,mu4,mu6,bth,bth2,rxt,x,upl2,upl,gx + real(dp) :: exdxdt,gr,zm,s,zm2,zm3,fe0m,ffe,uplh +! + yg = apar(1) + mu = apar(2) + npl = apar(3) + cr = apar(4) + n = int(apar(5)) + m = int(apar(6)) + ih = int(apar(7)) +! + bth2=2.0d0/mu + bth=sqrt(bth2) + mu2=mu*mu + mu4=mu2*mu2 + mu6=mu4*mu2 +! + rxt=sqrt(1.0d0+t*t/(2.0d0*mu)) + x = t*rxt + upl2=bth2*x**2 + upl=bth*x + gx=1.0d0+t*t/mu + exdxdt=cr*exp(-t*t)*gx/rxt + gr=npl*upl+n*yg + zm=-mu*(gx-gr) + s=mu*(gx+gr) + zm2=zm*zm + zm3=zm2*zm + call calcei3(zm,fe0m) + ffe=0.0d0 + uplh=upl**ih + if(n.eq.0.and.m.eq.0) ffe=exdxdt*fe0m*upl2 + if(m.eq.1) ffe=(1.0d0+s*(1.0d0-zm*fe0m))*uplh/mu2 + if(m.eq.2) ffe=(6.0d0-2.0d0*zm+4.0d0*s & + & +s*s*(1.0d0+zm-zm2*fe0m))*uplh/mu4 + if(m.eq.3) ffe=(18.0d0*s*(s+4.0d0-zm)+6.0d0*(20.0d0-8.0d0*zm+zm2) & + & +s**3*(2.0d0+zm+zm2-zm3*fe0m))*uplh/mu6 + fhermit= exdxdt*ffe + end function fhermit +! +! + subroutine antihermitian(ri,yg,mu,npl,cr,ci,lrm) + implicit none + integer ::lrm,n,k,m,mm + real(dp), parameter :: rpi=1.77245385090552d0 + integer, parameter :: lmx=20,nmx=lmx+2 + real(dp) :: fsbi(nmx) + real(dp) :: ri(lrm,0:2,lrm) + real(dp) :: yg,mu,npl,cmu,npl2,dnl,cr,ci,ygn,rdu2,rdu + real(dp) :: du,ub,aa,up,um,gp,gm,xp,xm,eep,eem,ee,cm,cim + real(dp) :: fi0p0,fi1p0,fi2p0,fi0m0,fi1m0,fi2m0 + real(dp) :: fi0p1,fi1p1,fi2p1,fi0m1,fi1m1,fi2m1,fi0m,fi1m,fi2m +! + do n=1,lrm + do k=0,2 + do m=1,lrm + ri(n,k,m)=0.0d0 + end do + end do + end do +! + npl2=npl*npl + dnl=1.0d0-npl2 + cmu=npl*mu +! + do n=1,lrm + ygn=n*yg + rdu2=ygn**2-dnl + if(rdu2.gt.0.0d0) then + rdu=sqrt(rdu2) + du=rdu/dnl + ub=npl*ygn/dnl + aa=mu*npl*du + if (abs(aa).gt.5.0d0) then + up=ub+du + um=ub-du + gp=npl*up+ygn + gm=npl*um+ygn + xp=up+1.0d0/cmu + xm=um+1.0d0/cmu + eem=exp(-mu*(gm-1.0d0)) + eep=exp(-mu*(gp-1.0d0)) + fi0p0=-1.0d0/cmu + fi1p0=-xp/cmu + fi2p0=-(1.0d0/cmu**2+xp*xp)/cmu + fi0m0=-1.0d0/cmu + fi1m0=-xm/cmu + fi2m0=-(1.0d0/cmu**2+xm*xm)/cmu + do m=1,lrm + fi0p1=-2.0d0*m*(fi1p0-ub*fi0p0)/cmu + fi0m1=-2.0d0*m*(fi1m0-ub*fi0m0)/cmu + fi1p1=-((1.0d0+2.0d0*m)*fi2p0-2.0d0*(m+1.0d0)*ub*fi1p0 & + +up*um*fi0p0)/cmu + fi1m1=-((1.0d0+2.0d0*m)*fi2m0-2.0d0*(m+1.0d0)*ub*fi1m0 & + +up*um*fi0m0)/cmu + fi2p1=(2.0d0*(1.0d0+m)*fi1p1-2.0d0*m* & + (ub*fi2p0-up*um*fi1p0))/cmu + fi2m1=(2.0d0*(1.0d0+m)*fi1m1-2.0d0*m* & + (ub*fi2m0-up*um*fi1m0))/cmu + if(m.ge.n) then + ri(n,0,m)=0.5d0*ci*dnl**m*(fi0p1*eep-fi0m1*eem) + ri(n,1,m)=0.5d0*ci*dnl**m*(fi1p1*eep-fi1m1*eem) + ri(n,2,m)=0.5d0*ci*dnl**m*(fi2p1*eep-fi2m1*eem) + end if + fi0p0=fi0p1 + fi1p0=fi1p1 + fi2p0=fi2p1 + fi0m0=fi0m1 + fi1m0=fi1m1 + fi2m0=fi2m1 + end do + else + ee=exp(-mu*(ygn-1.0d0+npl*ub)) + call ssbi(aa,n,lrm,fsbi) + do m=n,lrm + cm=rpi*fact(m)*du**(2*m+1) + cim=0.5d0*ci*dnl**m + mm=m-n+1 + fi0m=cm*fsbi(mm) + fi1m=-0.5d0*aa*cm*fsbi(mm+1) + fi2m=0.5d0*cm*(fsbi(mm+1)+0.5d0*aa*aa*fsbi(mm+2)) + ri(n,0,m)=cim*ee*fi0m + ri(n,1,m)=cim*ee*(du*fi1m+ub*fi0m) + ri(n,2,m)=cim*ee*(du*du*fi2m+2.0d0*du*ub*fi1m+ub*ub*fi0m) + end do + end if + end if + end do +! + return +! + end subroutine antihermitian +! + subroutine ssbi(zz,n,l,fsbi) + implicit none + integer :: n,l,k,m,mm + real(dp) :: c0,c1,sbi,zz + real(dp) :: gamm + real(dp), parameter :: eps=1.0d-10 + integer, parameter :: lmx=20,nmx=lmx+2 + real(dp) :: fsbi(nmx) + do m=n,l+2 + c0=1.0d0/gamm(dble(m)+1.5d0) + sbi=c0 + do k=1,50 + c1=c0*0.25d0*zz**2/(dble(m+k)+0.5d0)/dble(k) + sbi=sbi+c1 + if(c1/sbi.lt.eps) exit + c0=c1 + end do + mm=m-n+1 + fsbi(mm)=sbi + end do + return +! + end subroutine ssbi +! +! +! +! +! +! + subroutine expinit + implicit none + integer :: i +! +! + do i = 1, npts+1 + ttv(i) = -tmax+dble(i-1)*dtex + extv(i)=exp(-ttv(i)*ttv(i)) + end do +! + end subroutine expinit +! +! +! + function fact(k) + integer :: i,k + real(dp) :: fact + fact=0.0d0 + if(k.lt.0) return + fact=1.0d0 + if(k.eq.0) return + do i=1,k + fact=fact*i + end do + return + end function fact +! +! +! +! ====================================================================== +! Weakly relativistic dielectric tensor +! computation of dielectric tensor elements +! Krivenki and Orefice, JPP 30,125 (1983) +! + subroutine diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm) +! + implicit none +! parameters + complex(dp), parameter :: ui=(0.0d0,1.0d0) +! arguments + integer :: lrm + real(dp) :: xg,yg,npl,mu + complex(dp) :: e330,epsl(3,3,lrm) +! internal + integer :: l,lm,is,k + real(dp) :: npl2,fcl,w,asl,bsl,fact + complex(dp) :: ca11,ca12,ca13,ca22,ca23,ca33 + complex(dp) :: cq0p,cq0m,cq1p,cq1m,cq2p + complex(dp), dimension(0:lrm,0:2) :: cefp,cefm +! + npl2=npl*npl +! + call fsup(cefp,cefm,lrm,yg,npl,mu) +! + do l=1,lrm + lm=l-1 + fcl=0.5d0**l*((1.0d0/yg)**2/mu)**lm*fact(2*l)/fact(l) + ca11=(0.d0,0.d0) + ca12=(0.d0,0.d0) + ca13=(0.d0,0.d0) + ca22=(0.d0,0.d0) + ca23=(0.d0,0.d0) + ca33=(0.d0,0.d0) + do is=0,l + k=l-is + w=(-1.0d0)**k +! + asl=w/(fact(is+l)*fact(l-is)) + bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0d0*l-1.0d0)) +! + cq0p=mu*cefp(is,0) + cq0m=mu*cefm(is,0) + cq1p=mu*npl*(cefp(is,0)-cefp(is,1)) + cq1m=mu*npl*(cefm(is,0)-cefm(is,1)) + cq2p=cefp(is,1)+mu*npl2*(cefp(is,2)+cefp(is,0)-2.0d0*cefp(is,1)) +! + ca11=ca11+is**2*asl*cq0p + ca12=ca12+is*l*asl*cq0m + ca22=ca22+bsl*cq0p + ca13=ca13+is*asl*cq1m/yg + ca23=ca23+l*asl*cq1p/yg + ca33=ca33+asl*cq2p/yg**2 + end do + epsl(1,1,l) = - xg*ca11*fcl + epsl(1,2,l) = + ui*xg*ca12*fcl + epsl(2,2,l) = - xg*ca22*fcl + epsl(1,3,l) = - xg*ca13*fcl + epsl(2,3,l) = - ui*xg*ca23*fcl + epsl(3,3,l) = - xg*ca33*fcl + end do +! + cq2p=cefp(0,1)+mu*npl2*(cefp(0,2)+cefp(0,0)-2.0d0*cefp(0,1)) + e330=1.0d0-xg*mu*cq2p +! + epsl(1,1,1) = 1.d0 + epsl(1,1,1) + epsl(2,2,1) = 1.d0 + epsl(2,2,1) +! + do l=1,lrm + epsl(2,1,l) = - epsl(1,2,l) + epsl(3,1,l) = epsl(1,3,l) + epsl(3,2,l) = - epsl(2,3,l) + end do +! + return + end +! + subroutine fsup(cefp,cefm,lrm,yg,npl,mu) +! + implicit none +! parameters + real(dp), parameter :: apsicr=0.7d0 + complex(dp), parameter :: ui=(0.0d0,1.0d0) +! arguments + integer :: lrm + real(dp) :: yg,npl,mu + complex(dp), dimension(0:lrm,0:2) :: cefp,cefm +! internal + integer :: is,l,iq,ir,iflag + real(dp) :: psi,apsi,alpha,phi2,phim + real(dp) :: xp,yp,xm,ym,x0,y0 + real(dp) :: zrp,zip,zrm,zim,zr0,zi0 + complex(dp) :: czp,czm,cf12,cf32,cphi + complex(dp) :: cz0,cdz0,cf0,cf1,cf2 +! + psi=sqrt(0.5d0*mu)*npl + apsi=abs(psi) +! + do is=0,lrm + alpha=npl*npl/2.0d0+is*yg-1.0d0 + phi2=mu*alpha + phim=sqrt(abs(phi2)) + if (alpha.ge.0) then + xp=psi-phim + yp=0.0d0 + xm=-psi-phim + ym=0.0d0 + x0=-phim + y0=0.0d0 + else + xp=psi + yp=phim + xm=-psi + ym=phim + x0=0.0d0 + y0=phim + end if + call zetac (xp,yp,zrp,zip,iflag) + call zetac (xm,ym,zrm,zim,iflag) +! + czp=dcmplx(zrp,zip) + czm=dcmplx(zrm,zim) + cf12=(0.0d0,0.0d0) + if (alpha.ge.0) then + if (alpha.ne.0) cf12=-(czp+czm)/(2.0d0*phim) + else + cf12=-ui*(czp+czm)/(2.0d0*phim) + end if +! + if(apsi.gt.apsicr) then + cf32=-(czp-czm)/(2.0d0*psi) + else + cphi=phim + if(alpha.lt.0) cphi=-ui*phim + call zetac (x0,y0,zr0,zi0,iflag) + cz0=dcmplx(zr0,zi0) + cdz0=2.0d0*(1.0d0-cphi*cz0) + cf32=cdz0 + end if +! + cf0=cf12 + cf1=cf32 + cefp(is,0)=cf32 + cefm(is,0)=cf32 + do l=1,is+2 + iq=l-1 + if(apsi.gt.apsicr) then + cf2=(1.0d0+phi2*cf0-(iq+0.5d0)*cf1)/psi**2 + else + cf2=(1.0d0+phi2*cf1)/dble(iq+1.5d0) + end if + ir=l-is + if(ir.ge.0) then + cefp(is,ir)=cf2 + cefm(is,ir)=cf2 + end if + cf0=cf1 + cf1=cf2 + end do +! + if(is.ne.0) then +! + alpha=npl*npl/2.0d0-is*yg-1.0d0 + phi2=mu*alpha + phim=sqrt(abs(phi2)) + if (alpha.ge.0.0d0) then + xp=psi-phim + yp=0.0d0 + xm=-psi-phim + ym=0.0d0 + x0=-phim + y0=0.0d0 + else + xp=psi + yp=phim + xm=-psi + ym=phim + x0=0.0d0 + y0=phim + end if + call zetac (xp,yp,zrp,zip,iflag) + call zetac (xm,ym,zrm,zim,iflag) +! + czp=dcmplx(zrp,zip) + czm=dcmplx(zrm,zim) +! + cf12=(0.0d0,0.0d0) + if (alpha.ge.0) then + if (alpha.ne.0.0d0) cf12=-(czp+czm)/(2.0d0*phim) + else + cf12=-ui*(czp+czm)/(2.0d0*phim) + end if + if(apsi.gt.apsicr) then + cf32=-(czp-czm)/(2.0d0*psi) + else + cphi=phim + if(alpha.lt.0) cphi=-ui*phim + call zetac (x0,y0,zr0,zi0,iflag) + cz0=dcmplx(zr0,zi0) + cdz0=2.0d0*(1.0d0-cphi*cz0) + cf32=cdz0 + end if +! + cf0=cf12 + cf1=cf32 + do l=1,is+2 + iq=l-1 + if(apsi.gt.apsicr) then + cf2=(1.0d0+phi2*cf0-(iq+0.5d0)*cf1)/psi**2 + else + cf2=(1.0d0+phi2*cf1)/dble(iq+1.5d0) + end if + ir=l-is + if(ir.ge.0) then + cefp(is,ir)=cefp(is,ir)+cf2 + cefm(is,ir)=cefm(is,ir)-cf2 + end if + cf0=cf1 + cf1=cf2 + end do +! + end if +! + end do +! + return + end +! ====================================================================== + +end module dispersion \ No newline at end of file diff --git a/src/gray.f b/src/gray.f index 3f89652..84853d9 100644 --- a/src/gray.f +++ b/src/gray.f @@ -682,7 +682,7 @@ c common/parban/b0,rr0m,zr0m,rpam common/parqq/q0,qa,alq common/parqte/te0,dte0,alt1,alt2 - common/zz/Zeff + common/zz/Zeffan common/bound/zbmin,zbmax c common/parbres/bres @@ -3265,6 +3265,7 @@ c spline interpolation of rhotor versus rhopol subroutine vectinit + use dispersion, only: expinit implicit real*8 (a-h,o-z) parameter(jmx=31,kmx=36,nmx=8000) dimension psjki(jmx,kmx,nmx) @@ -3273,10 +3274,6 @@ c spline interpolation of rhotor versus rhopol dimension currj(jmx,kmx,nmx),didst(jmx,kmx,nmx),ccci(jmx,kmx,nmx) dimension iiv(jmx,kmx),iop(jmx,kmx),iow(jmx,kmx),ihcd(jmx,kmx) dimension istore(jmx,kmx),anwcl(3),xwcl(3) - parameter(tmax=5,npts=500) - real*8 ttv(npts+1),extv(npts+1) - - common/ttex/ttv,extv c common/warm/iwarm,ilarm common/iiv/iiv @@ -3319,13 +3316,7 @@ c end do end do c - if(iwarm.gt.1) then - dt=2.0d0*tmax/dble(npts) - do i = 1, npts+1 - ttv(i) = -tmax+dble(i-1)*dt - extv(i)=exp(-ttv(i)*ttv(i)) - end do - end if + if(iwarm.gt.1) call expinit c return end @@ -4551,13 +4542,13 @@ c common/crad/psrad,derad,terad,zfc common/coffz/cz common/eqnn/nr,nz,npp,nintp - common/zz/Zeff + common/zz/Zeffan c fzeff=1 ps=arg if(ps.gt.1.0d0.and.ps.lt.0.0d0) return if(iprof.eq.0) then - fzeff=zeff + fzeff=zeffan else call locate(psrad,npp,ps,k) k=max(1,min(k,npp-1)) @@ -5366,20 +5357,23 @@ c c c subroutine ecrh_cd - use const_and_precisions, only : eqsi=>e_,mesi=>me_,vcsi=>c_, - * mc2=>mc2_ + use const_and_precisions, only : mc2=>mc2_ + use dispersion, only : harmnumber, warmdisp implicit real*8 (a-h,o-z) parameter(taucr=12.0d0,xxcr=16.0d0,eps=1.d-8) complex*16 ex,ey,ez c - common/nharm/nharm,nhf + common/nharm/nhmin,nhmax common/warm/iwarm,ilarm + common/imx/imx common/ieccd/ieccd c + common/xgxg/xg common/ygyg/yg common/nplr/anpl,anpr common/vgv/vgm,derdnm common/parwv/ak0,akinv,fhz + common/mode/sox c common/nprw/anprre,anprim common/epolar/ex,ey,ez @@ -5389,12 +5383,9 @@ c common/ierr/ierr common/iokh/iokhawa c - common/psival/psinv common/tete/tekev common/dens/dens,ddens - common/amut/amu - common/zz/Zeff common/btot/btot common/bmxmn/bmax,bmin common/fc/fc @@ -5406,47 +5397,16 @@ c effjcd=0.0d0 c tekev=temp(psinv) - amu=mc2/tekev - zeff=fzeff(psinv) c if(tekev.le.0.0d0.or.tau.gt.taucr) return c - dnl=1.0d0-anpl*anpl - if(iwarm.eq.1) then - ygc=1.0d0-anpl**2/2.0d0 - else - ygc=sqrt(1.0d0-anpl**2) - end if + amu=mc2/tekev + call harmnumber(yg,amu,anpl,nhmin,nhmax,iwarm) + if(nhmin.eq.0) return c - nharm=int(ygc/yg-eps)+1 -c - if(nharm.eq.0) return -c - nhf=0 - do nn=nharm,nharm+4 - ygn=nn*yg - if(iwarm.eq.1) then - rdu2=2.0d0*(ygn-ygc) - u1=anpl+sqrt(rdu2) - u2=anpl-sqrt(rdu2) - uu2=min(u1*u1,u2*u2) - argex=amu*uu2/2.0d0 - else - rdu2=ygn**2-ygc**2 - g1=(ygn+anpl*sqrt(rdu2))/dnl - g2=(ygn-anpl*sqrt(rdu2))/dnl - gg=min(g1,g2) - argex=amu*(gg-1.0d0) - u1=(ygn*anpl+sqrt(rdu2))/dnl - u2=(ygn*anpl-sqrt(rdu2))/dnl - end if - if(argex.le.xxcr) nhf=nn - end do - if(nhf.eq.0) return -c - lrm=ilarm - if(lrm.lt.nhf) lrm=nhf - call dispersion(lrm) + lrm=max(ilarm,nhmax) + call warmdisp(xg,yg,amu,anpl,anpr,sox,lrm,ierr,anprre,anprim, + * iwarm,imx,ex,ey,ez) c akim=ak0*anprim ratiovgr=2.0d0*anpr/derdnm @@ -5457,13 +5417,14 @@ c ratiovgr=2.0d0*anpr/derdnm*vgm print*,' IERR = ', ierr,' alpha negative' end if c - ithn=1 - if(lrm.gt.nharm) ithn=2 if(ieccd.gt.0) then + ithn=1 + if(lrm.gt.nhmin) ithn=2 + zeff=fzeff(psinv) rbn=btot/bmin rbx=btot/bmax call eccd(yg,anpl,anprre,amu,zeff,rbn,rbx,fc,ex,ey,ez, - * dens,psinv,ieccd,nharm,nhf,ithn,effjcd,iokhawa,ierr) + * dens,psinv,ieccd,nhmin,nhmax,ithn,effjcd,iokhawa,ierr) end if return @@ -5471,931 +5432,6 @@ c c c c - subroutine dispersion(lrm) - implicit real*8(a-h,o-z) - parameter(one=1.0d0) - complex*16 cc0,cc2,cc4,rr - complex*16 anpr2a,anpra - complex*16 anpr,anpr2,ex,ey,ez,den - complex*16 epsl(3,3,lrm),sepsl(3,3),e330 - complex*16 e11,e22,e12,e13,e23 -c complex*16 e33,e21,e31,e32 - complex*16 a13,a31,a23,a32,a33 -c - common/ygyg/yg - common/xgxg/xg - common/nplr/anpl,anprf - common/mode/sox - common/warm/iwarm,ilarm - common/nprw/anprr,anpri - common/epolar/ex,ey,ez - common/imx/imx - - errnpr=1.0d0 - anpr2a=anprf**2 - anpl2=anpl*anpl -c - if (iwarm.eq.1) then - call diel_tens_wr(e330,epsl,lrm) - else - call diel_tens_fr(e330,epsl,lrm) - end if -c - imxx=abs(imx) -c loop to disable convergence if failure detected - do - do i=1,imxx -c - do j=1,3 - do k=1,3 - sepsl(k,j)=dcmplx(0.0d0,0.0d0) - do ilrm=1,lrm - sepsl(k,j)=sepsl(k,j)+epsl(k,j,ilrm)*anpr2a**(ilrm-1) - end do - end do - end do -c - anpra=sqrt(anpr2a) -c - e11=sepsl(1,1) - e22=sepsl(2,2) - e12=sepsl(1,2) - a33=sepsl(3,3) - a13=sepsl(1,3) - a23=sepsl(2,3) - a31=a13 - a32=-a23 -c e33=e330+anpr2a*a33 - e13=anpra*a13 - e23=anpra*a23 -c e21=-e12 -c e31=e13 -c e32=-e23 -c - cc4=(e11-anpl2)*(1.0d0-a33)+(a13+anpl)*(a31+anpl) - cc2=-e12*e12*(1.0d0-a33) - . -a32*e12*(a13+anpl)+a23*e12*(a31+anpl) - . -(a23*a32+e330+(e22-anpl2)*(1.0d0-a33))*(e11-anpl2) - . -(a13+anpl)*(a31+anpl)*(e22-anpl2) - cc0=e330*((e11-anpl2)*(e22-anpl2)+e12*e12) -c - rr=cc2*cc2-4.0d0*cc0*cc4 -c - if(yg.gt.1.0d0) then - s=sox - if(dimag(rr).LE.0.0d0) s=-s - else - s=-sox - if(dble(rr).le.0.d0.and.dimag(rr).ge.0.d0) s=-s - end if -c - anpr2=(-cc2+s*sqrt(rr))/(2.0d0*cc4) -c - errnpr=abs(1.0d0-abs(anpr2)/abs(anpr2a)) - if(i.gt.1.and.errnpr.lt.1.0d-5) exit - anpr2a=anpr2 - end do - if(i.gt.imxx .and. imxx.gt.1) then - if (imx.lt.0) then - write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,"// - ."': convergence disabled.',e12.5)") xg,yg,sqrt(abs(anpr2)),anpl - imxx=1 - else - write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,"// - ."': convergence failed.',e12.5)") xg,yg,sqrt(abs(anpr2)),anpl - exit - end if - else - exit - end if - print*,yg,imx,imxx - end do -c end loop to disable convergence - if(sqrt(dble(anpr2)).lt.0.0d0 .or. anpr2.ne.anpr2 - . .or. abs(anpr2).eq.huge(one) .or. abs(anpr2).le.tiny(one)) then - write(*,"(' X =',f7.4,' Y =',f7.4,"// - . "' Nperp =',f7.4,'!')") xg,yg,sqrt(abs(anpr2)) - ierr=99 - anpr2=(0.0d0,0.0d0) - end if -c - anpr=sqrt(anpr2) - anprr=dble(anpr) - anpri=dimag(anpr) -c - ex=dcmplx(0.0d0,0.0d0) - ey=dcmplx(0.0d0,0.0d0) - ez=dcmplx(0.0d0,0.0d0) -c - if (abs(anpl).gt.1.0d-6) then - den=e12*e23-(e13+anpr*anpl)*(e22-anpr2-anpl2) - ey=-(e12*(e13+anpr*anpl)+(e11-anpl2)*e23)/den - ez=(e12*e12+(e22-anpr2-anpl2)*(e11-anpl2))/den - ez2=abs(ez)**2 - ey2=abs(ey)**2 - enx2=1.0d0/(1.0d0+ey2+ez2) - ex=dcmplx(sqrt(enx2),0.0d0) - ez2=ez2*enx2 - ey2=ey2*enx2 - ez=ez*ex - ey=ey*ex - else - if(sox.lt.0.0d0) then - ez=dcmplx(1.0d0,0.0d0) - ez2=abs(ez)**2 - else - ex2=1.0d0/(1.0d0+abs(-e11/e12)**2) - ex=sqrt(ex2) - ey=-ex*e11/e12 - ey2=abs(ey)**2 - ez2=0.0d0 - end if - end if -c - return - end -c -c Fully relativistic case -c computation of dielectric tensor elements -c up to third order in Larmor radius for hermitian part -c - subroutine diel_tens_fr(e330,epsl,lrm) - use const_and_precisions, only : pi,rpi=>sqrt_pi,ui=>im - implicit real*8(a-h,o-z) - complex*16 e330,epsl(3,3,lrm) - complex*16 ca11,ca12,ca22,ca13,ca23,ca33 - complex*16 cq0p,cq0m,cq1p,cq1m,cq2p - dimension rr(-lrm:lrm,0:2,0:lrm),ri(lrm,0:2,lrm) -c - common/xgxg/xg - common/ygyg/yg - common/amut/amu - common/nplr/anpl,anpr - common/resah/anpl2,dnl -c - common/cri/cr,ci - common/warm/iwarm,ilarm -c - anpl2=anpl**2 - dnl=1.0d0-anpl2 -c - cmxw=1.0d0+15.0d0/(8.0d0*amu)+105.0d0/(128.0d0*amu**2) - cr=-amu*amu/(rpi*cmxw) - ci=sqrt(2.0d0*pi*amu)*amu**2/cmxw -c - do l=1,lrm - do j=1,3 - do i=1,3 - epsl(i,j,l)=dcmplx(0.0d0,0.0d0) - end do - end do - end do -c - if(iwarm.eq.2) call hermitian(rr,lrm) - if(iwarm.eq.4) call hermitian_2(rr,lrm) -c - call antihermitian(ri,lrm) -c - do l=1,lrm - lm=l-1 - fal=-0.25d0**l*fact(2*l)/(fact(l)**2*yg**(2*lm)) - ca11=(0.0d0,0.0d0) - ca12=(0.0d0,0.0d0) - ca13=(0.0d0,0.0d0) - ca22=(0.0d0,0.0d0) - ca23=(0.0d0,0.0d0) - ca33=(0.0d0,0.0d0) - do is=0,l - k=l-is - w=(-1.0d0)**k -c - asl=w/(fact(is+l)*fact(l-is)) - bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0d0*l-1.0d0)) -c - if(is.gt.0) then - cq0p=rr(is,0,l)+rr(-is,0,l)+ui*ri(is,0,l) - cq0m=rr(is,0,l)-rr(-is,0,l)+ui*ri(is,0,l) - cq1p=rr(is,1,l)+rr(-is,1,l)+ui*ri(is,1,l) - cq1m=rr(is,1,l)-rr(-is,1,l)+ui*ri(is,1,l) - cq2p=rr(is,2,l)+rr(-is,2,l)+ui*ri(is,2,l) - else - cq0p=rr(is,0,l) - cq0m=rr(is,0,l) - cq1p=rr(is,1,l) - cq1m=rr(is,1,l) - cq2p=rr(is,2,l) - end if -c - ca11=ca11+is**2*asl*cq0p - ca12=ca12+is*l*asl*cq0m - ca22=ca22+bsl*cq0p - ca13=ca13+is*asl*cq1m/yg - ca23=ca23+l*asl*cq1p/yg - ca33=ca33+asl*cq2p/yg**2 - end do - epsl(1,1,l) = - xg*ca11*fal - epsl(1,2,l) = + ui*xg*ca12*fal - epsl(2,2,l) = - xg*ca22*fal - epsl(1,3,l) = - xg*ca13*fal - epsl(2,3,l) = - ui*xg*ca23*fal - epsl(3,3,l) = - xg*ca33*fal - end do -c - cq2p=rr(0,2,0) - e330=1.0d0+xg*cq2p -c - epsl(1,1,1) = 1.d0 + epsl(1,1,1) - epsl(2,2,1) = 1.d0 + epsl(2,2,1) -c - do l=1,lrm - epsl(2,1,l) = - epsl(1,2,l) - epsl(3,1,l) = epsl(1,3,l) - epsl(3,2,l) = - epsl(2,3,l) - end do -c - return - end -c -c -c - subroutine hermitian(rr,lrm) - implicit real*8(a-h,o-z) - parameter(tmax=5,npts=500) - dimension rr(-lrm:lrm,0:2,0:lrm) - real*8 ttv(npts+1),extv(npts+1) -c - common/ttex/ttv,extv -c - common/ygyg/yg - common/amut/amu - common/nplr/anpl,anpr - common/warm/iwarm,ilarm - common/cri/cr,ci -c - do n=-lrm,lrm - do k=0,2 - do m=0,lrm - rr(n,k,m)=0.0d0 - end do - end do - end do -c - llm=min(3,lrm) -c - dt=2.0d0*tmax/dble(npts) - bth2=2.0d0/amu - bth=sqrt(bth2) - amu2=amu*amu - amu4=amu2*amu2 - amu6=amu4*amu2 -c - do i = 1, npts+1 - t = ttv(i) - rxt=sqrt(1.0d0+t*t/(2.0d0*amu)) - x = t*rxt - upl2=bth2*x**2 - upl=bth*x - gx=1.0d0+t*t/amu - exdx=cr*extv(i)*gx/rxt*dt -c - n1=1 - if(iwarm.gt.2) n1=-llm -c - do n=n1,llm - nn=abs(n) - gr=anpl*upl+n*yg - zm=-amu*(gx-gr) - s=amu*(gx+gr) - zm2=zm*zm - zm3=zm2*zm - call calcei3(zm,fe0m) -c - do m=nn,llm - if(n.eq.0.and.m.eq.0) then - rr(0,2,0) = rr(0,2,0) - exdx*fe0m*upl2 - else - ffe=0.0d0 - if(m.eq.1) ffe=(1.0d0+s*(1.0d0-zm*fe0m))/amu2 - if(m.eq.2) ffe=(6.0d0-2.0d0*zm+4.0d0*s - . +s*s*(1.0d0+zm-zm2*fe0m))/amu4 - if(m.eq.3) ffe=(18.0d0*s*(s+4.0d0-zm)+6.0d0* - . (20.0d0-8.0d0*zm+zm2)+s**3*(2.0d0+zm+zm2-zm3*fe0m))/amu6 -c - rr(n,0,m) = rr(n,0,m) + exdx*ffe - rr(n,1,m) = rr(n,1,m) + exdx*ffe*upl - rr(n,2,m) = rr(n,2,m) + exdx*ffe*upl2 - end if -c - end do - end do - end do -c - if(iwarm.gt.2) return -c - sy1=1.0d0+yg - sy2=1.0d0+yg*2.0d0 - sy3=1.0d0+yg*3.0d0 -c - bth4=bth2*bth2 - bth6=bth4*bth2 -c - anpl2=anpl*anpl -c - rr(0,2,0) = -(1.0d0+bth2*(-1.25d0+1.5d0*anpl2) - . +bth4*(1.71875d0-6.0d0*anpl2+3.75d0*anpl2*anpl2)) - rr(0,1,1) = -anpl*bth2*(1.0d0+bth2*(-2.25d0+1.5d0*anpl2)) - rr(0,2,1) = -bth2*(1.0d0+bth2*(-0.5d0+1.5d0*anpl2)) - rr(-1,0,1) = -2.0d0/sy1*(1.0d0+bth2/sy1*(-1.25d0+0.5d0*anpl2 - . /sy1)+bth4/sy1*(-0.46875d0+(2.1875d0+0.625d0*anpl2)/ - . sy1-2.625d0*anpl2/sy1**2+0.75d0*anpl2*anpl2/sy1**3)) - rr(-1,1,1) = -anpl*bth2/sy1**2*(1.0d0+bth2*(1.25d0-3.5d0/sy1+ - . 1.5d0*anpl2/sy1**2)) - rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(1.25d0-1.75d0/sy1+1.5d0* - . anpl2/sy1**2)) -c - if(llm.gt.1) then -c - rr(0,0,2) = -4.0d0*bth2*(1.0d0+bth2*(-0.5d0+0.5d0*anpl2)+ - . bth4*(1.125d0-1.875d0*anpl2+0.75d0*anpl2*anpl2)) - rr(0,1,2) = -2.0d0*anpl*bth4*(1.0d0+bth2*(-1.5d0+1.5d0*anpl2)) - rr(0,2,2) = -2.0d0*bth4*(1.0d0+bth2*(0.75d0+1.5d0*anpl2)) - rr(-1,0,2) = -4.0d0*bth2/sy1*(1.0d0+bth2* - . (1.25d0-1.75d0/sy1+0.5d0*anpl2/sy1**2)+bth4* - . (0.46875d0-3.28125d0/sy1+(3.9375d0+1.5d0*anpl2)/sy1**2 - . -3.375d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4)) - rr(-1,1,2) = -2.0d0*bth4*anpl/sy1**2*(1.0d0+bth2* - . (3.0d0-4.5d0/sy1+1.5d0*anpl2/sy1**2)) - rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+bth2* - . (3.0d0-2.25d0/sy1+1.5d0*anpl2/sy1**2)) - rr(-2,0,2) = -4.0d0*bth2/sy2*(1.0d0+bth2* - . (1.25d0-1.75d0/sy2+0.5d0*anpl2/sy2**2)+bth4* - . (0.46875d0-3.28125d0/sy2+(3.9375d0+1.5d0*anpl2)/sy2**2 - . -3.375d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4)) - rr(-2,1,2) =-2.0d0*bth4*anpl/sy2**2*(1.0d0+bth2* - . (3.0d0-4.5d0/sy2+1.5d0*anpl2/sy2**2)) - rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+bth2* - . (3.0d0-2.25d0/sy2+1.5d0*anpl2/sy2**2)) -c - if(llm.gt.2) then -c - rr(0,0,3) = -12.0d0*bth4*(1+bth2*(0.75d0+0.5d0*anpl2)+bth4* - . (1.21875d0-1.5d0*anpl2+0.75d0*anpl2*anpl2)) - rr(0,1,3) = -6.0d0*anpl*bth6*(1+bth2*(-0.25d0+1.5d0*anpl2)) - rr(0,2,3) = -6.0d0*bth6*(1+bth2*(2.5d0+1.5d0*anpl2)) - rr(-1,0,3) = -12.0d0*bth4/sy1* - . (1.0d0+bth2*(3.0d0-2.25d0/sy1+0.5d0*anpl2/sy1**2)+ - . bth4*(3.75d0-8.71875d0/sy1+(6.1875d0+2.625d0*anpl2) - . /sy1**2-4.125d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4)) - rr(-1,1,3) = -6.0d0*anpl*bth6/sy1**2* - . (1.0d0+bth2*(5.25d0-5.5d0/sy1+1.5d0*anpl2/sy1**2)) - rr(-1,2,3) = -6.0d0*bth6/sy1* - . (1.0d0+bth2*(5.25d0-2.75d0/sy1+1.5d0*anpl2/sy1**2)) -c - rr(-2,0,3) = -12.0d0*bth4/sy2* - . (1.0d0+bth2*(3.0d0-2.25d0/sy2+0.5d0*anpl2/sy2**2)+ - . bth4*(3.75d0-8.71875d0/sy2+(6.1875d0+2.625d0*anpl2) - . /sy2**2-4.125d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4)) - rr(-2,1,3) = -6.0d0*anpl*bth6/sy2**2* - . (1.0d0+bth2*(5.25d0-5.5d0/sy2+1.5d0*anpl2/sy2**2)) - rr(-2,2,3) = -6.0d0*bth6/sy2* - . (1.0d0+bth2*(5.25d0-2.75d0/sy2+1.5d0*anpl2/sy2**2)) -c - rr(-3,0,3) = -12.0d0*bth4/sy3* - . (1.0d0+bth2*(3.0d0-2.25d0/sy3+0.5d0*anpl2/sy3**2)+ - . bth4*(3.75d0-8.71875d0/sy3+(6.1875d0+2.625d0*anpl2) - . /sy3**2-4.125d0*anpl2/sy3**3+0.75d0*anpl2*anpl2/sy3**4)) - rr(-3,1,3) = -6.0d0*anpl*bth6/sy3**2* - . (1.0d0+bth2*(5.25d0-5.5d0/sy3+1.5d0*anpl2/sy3**2)) - rr(-3,2,3) = -6.0d0*bth6/sy3* - . (1.0d0+bth2*(5.25d0-2.75d0/sy3+1.5d0*anpl2/sy3**2)) -c - end if - end if -c - return - end - - - - subroutine hermitian_2(rr,lrm) - implicit real*8(a-h,o-z) - parameter(tmax=5) - dimension rr(-lrm:lrm,0:2,0:lrm) - parameter(epsa=0.0d0,epsr=1.0d-4) - parameter (lw=5000,liw=lw/4) - dimension w(lw),iw(liw) - external fhermit -c - common/ygyg/yg - common/amut/amu - common/nplr/anpl,anpr - common/warm/iwarm,ilarm - common/nmhermit/n,m,ih -c - do n=-lrm,lrm - do k=0,2 - do m=0,lrm - rr(n,k,m)=0.0d0 - end do - end do - end do -c - llm=min(3,lrm) -c - bth2=2.0d0/amu -c - n1=1 - if(iwarm.gt.10) n1=-llm -c - do n=n1,llm - nn=abs(n) - do m=nn,llm - if(n.eq.0.and.m.eq.0) then - ih=2 -c call dqagi(fhermit,bound,2,epsa,epsr,resfh, - call dqags(fhermit,-tmax,tmax,epsa,epsr,resfh, - . epp,neval,ier,liw,lw,last,iw,w) - rr(n,2,m) = resfh - else - do ih=0,2 -c call dqagi(fhermit,bound,2,epsa,epsr,resfh, - call dqags(fhermit,-tmax,tmax,epsa,epsr,resfh, - . epp,neval,ier,liw,lw,last,iw,w) - rr(n,ih,m) = resfh - end do - end if - end do - end do - - if(iwarm.gt.10) return -c - sy1=1.0d0+yg - sy2=1.0d0+yg*2.0d0 - sy3=1.0d0+yg*3.0d0 -c - bth4=bth2*bth2 - bth6=bth4*bth2 -c - anpl2=anpl*anpl -c - rr(0,2,0) = -(1.0d0+bth2*(-1.25d0+1.5d0*anpl2) - . +bth4*(1.71875d0-6.0d0*anpl2+3.75d0*anpl2*anpl2)) - rr(0,1,1) = -anpl*bth2*(1.0d0+bth2*(-2.25d0+1.5d0*anpl2)) - rr(0,2,1) = -bth2*(1.0d0+bth2*(-0.5d0+1.5d0*anpl2)) - rr(-1,0,1) = -2.0d0/sy1*(1.0d0+bth2/sy1*(-1.25d0+0.5d0*anpl2 - . /sy1)+bth4/sy1*(-0.46875d0+(2.1875d0+0.625d0*anpl2)/ - . sy1-2.625d0*anpl2/sy1**2+0.75d0*anpl2*anpl2/sy1**3)) - rr(-1,1,1) = -anpl*bth2/sy1**2*(1.0d0+bth2*(1.25d0-3.5d0/sy1+ - . 1.5d0*anpl2/sy1**2)) - rr(-1,2,1) = -bth2/sy1*(1.0d0+bth2*(1.25d0-1.75d0/sy1+1.5d0* - . anpl2/sy1**2)) -c - if(llm.gt.1) then -c - rr(0,0,2) = -4.0d0*bth2*(1.0d0+bth2*(-0.5d0+0.5d0*anpl2)+ - . bth4*(1.125d0-1.875d0*anpl2+0.75d0*anpl2*anpl2)) - rr(0,1,2) = -2.0d0*anpl*bth4*(1.0d0+bth2*(-1.5d0+1.5d0*anpl2)) - rr(0,2,2) = -2.0d0*bth4*(1.0d0+bth2*(0.75d0+1.5d0*anpl2)) - rr(-1,0,2) = -4.0d0*bth2/sy1*(1.0d0+bth2* - . (1.25d0-1.75d0/sy1+0.5d0*anpl2/sy1**2)+bth4* - . (0.46875d0-3.28125d0/sy1+(3.9375d0+1.5d0*anpl2)/sy1**2 - . -3.375d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4)) - rr(-1,1,2) = -2.0d0*bth4*anpl/sy1**2*(1.0d0+bth2* - . (3.0d0-4.5d0/sy1+1.5d0*anpl2/sy1**2)) - rr(-1,2,2) = -2.0d0*bth4/sy1*(1.0d0+bth2* - . (3.0d0-2.25d0/sy1+1.5d0*anpl2/sy1**2)) - rr(-2,0,2) = -4.0d0*bth2/sy2*(1.0d0+bth2* - . (1.25d0-1.75d0/sy2+0.5d0*anpl2/sy2**2)+bth4* - . (0.46875d0-3.28125d0/sy2+(3.9375d0+1.5d0*anpl2)/sy2**2 - . -3.375d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4)) - rr(-2,1,2) =-2.0d0*bth4*anpl/sy2**2*(1.0d0+bth2* - . (3.0d0-4.5d0/sy2+1.5d0*anpl2/sy2**2)) - rr(-2,2,2) = -2.0d0*bth4/sy2*(1.0d0+bth2* - . (3.0d0-2.25d0/sy2+1.5d0*anpl2/sy2**2)) -c - if(llm.gt.2) then -c - rr(0,0,3) = -12.0d0*bth4*(1+bth2*(0.75d0+0.5d0*anpl2)+bth4* - . (1.21875d0-1.5d0*anpl2+0.75d0*anpl2*anpl2)) - rr(0,1,3) = -6.0d0*anpl*bth6*(1+bth2*(-0.25d0+1.5d0*anpl2)) - rr(0,2,3) = -6.0d0*bth6*(1+bth2*(2.5d0+1.5d0*anpl2)) - rr(-1,0,3) = -12.0d0*bth4/sy1* - . (1.0d0+bth2*(3.0d0-2.25d0/sy1+0.5d0*anpl2/sy1**2)+ - . bth4*(3.75d0-8.71875d0/sy1+(6.1875d0+2.625d0*anpl2) - . /sy1**2-4.125d0*anpl2/sy1**3+0.75d0*anpl2*anpl2/sy1**4)) - rr(-1,1,3) = -6.0d0*anpl*bth6/sy1**2* - . (1.0d0+bth2*(5.25d0-5.5d0/sy1+1.5d0*anpl2/sy1**2)) - rr(-1,2,3) = -6.0d0*bth6/sy1* - . (1.0d0+bth2*(5.25d0-2.75d0/sy1+1.5d0*anpl2/sy1**2)) -c - rr(-2,0,3) = -12.0d0*bth4/sy2* - . (1.0d0+bth2*(3.0d0-2.25d0/sy2+0.5d0*anpl2/sy2**2)+ - . bth4*(3.75d0-8.71875d0/sy2+(6.1875d0+2.625d0*anpl2) - . /sy2**2-4.125d0*anpl2/sy2**3+0.75d0*anpl2*anpl2/sy2**4)) - rr(-2,1,3) = -6.0d0*anpl*bth6/sy2**2* - . (1.0d0+bth2*(5.25d0-5.5d0/sy2+1.5d0*anpl2/sy2**2)) - rr(-2,2,3) = -6.0d0*bth6/sy2* - . (1.0d0+bth2*(5.25d0-2.75d0/sy2+1.5d0*anpl2/sy2**2)) -c - rr(-3,0,3) = -12.0d0*bth4/sy3* - . (1.0d0+bth2*(3.0d0-2.25d0/sy3+0.5d0*anpl2/sy3**2)+ - . bth4*(3.75d0-8.71875d0/sy3+(6.1875d0+2.625d0*anpl2) - . /sy3**2-4.125d0*anpl2/sy3**3+0.75d0*anpl2*anpl2/sy3**4)) - rr(-3,1,3) = -6.0d0*anpl*bth6/sy3**2* - . (1.0d0+bth2*(5.25d0-5.5d0/sy3+1.5d0*anpl2/sy3**2)) - rr(-3,2,3) = -6.0d0*bth6/sy3* - . (1.0d0+bth2*(5.25d0-2.75d0/sy3+1.5d0*anpl2/sy3**2)) -c - end if - end if -c - return - end - - - - function fhermit(t) - implicit real*8 (a-h,o-z) -c - common/ygyg/yg - common/amut/amu - common/nplr/anpl,anpr - common/cri/cr,ci - common/nmhermit/n,m,ih - - bth2=2.0d0/amu - bth=sqrt(bth2) - amu2=amu*amu - amu4=amu2*amu2 - amu6=amu4*amu2 - - rxt=sqrt(1.0d0+t*t/(2.0d0*amu)) - x = t*rxt - upl2=bth2*x**2 - upl=bth*x - gx=1.0d0+t*t/amu - exdxdt=cr*exp(-t*t)*gx/rxt - gr=anpl*upl+n*yg - zm=-amu*(gx-gr) - s=amu*(gx+gr) - zm2=zm*zm - zm3=zm2*zm - call calcei3(zm,fe0m) - ffe=0.0d0 - uplh=upl**ih - if(n.eq.0.and.m.eq.0) ffe=exdxdt*fe0m*upl2 - if(m.eq.1) ffe=(1.0d0+s*(1.0d0-zm*fe0m))*uplh/amu2 - if(m.eq.2) ffe=(6.0d0-2.0d0*zm+4.0d0*s - . +s*s*(1.0d0+zm-zm2*fe0m))*uplh/amu4 - if(m.eq.3) ffe=(18.0d0*s*(s+4.0d0-zm)+6.0d0*(20.0d0-8.0d0*zm+zm2) - . +s**3*(2.0d0+zm+zm2-zm3*fe0m))*uplh/amu6 - fhermit= exdxdt*ffe - return - end -c -c -c - subroutine antihermitian(ri,lrm) - use const_and_precisions, only : rpi=>sqrt_pi - implicit none - integer lmx,nmx,lrm,n,k,m,mm - parameter(lmx=20,nmx=lmx+2) - real*8 fsbi(nmx) - real*8 ri(lrm,0:2,lrm) - real*8 yg,amu,anpl,cmu,anpl2,dnl,cr,ci,ygn,rdu2,rdu,anpr - real*8 du,ub,aa,up,um,gp,gm,xp,xm,eep,eem,ee,cm,cim - real*8 fi0p0,fi1p0,fi2p0,fi0m0,fi1m0,fi2m0 - real*8 fi0p1,fi1p1,fi2p1,fi0m1,fi1m1,fi2m1,fi0m,fi1m,fi2m - real*8 fact -c - common/ygyg/yg - common/amut/amu - common/nplr/anpl,anpr - common/uu/ub,cmu - common/resah/anpl2,dnl -c - common/cri/cr,ci -c - do n=1,lrm - do k=0,2 - do m=1,lrm - ri(n,k,m)=0.0d0 - end do - end do - end do -c - dnl=1.0d0-anpl2 - cmu=anpl*amu -c - do n=1,lrm - ygn=n*yg - rdu2=ygn**2-dnl - if(rdu2.gt.0.0d0) then - rdu=sqrt(rdu2) - du=rdu/dnl - ub=anpl*ygn/dnl - aa=amu*anpl*du - if (abs(aa).gt.5.0d0) then - up=ub+du - um=ub-du - gp=anpl*up+ygn - gm=anpl*um+ygn - xp=up+1.0d0/cmu - xm=um+1.0d0/cmu - eem=exp(-amu*(gm-1.0d0)) - eep=exp(-amu*(gp-1.0d0)) - fi0p0=-1.0d0/cmu - fi1p0=-xp/cmu - fi2p0=-(1.0d0/cmu**2+xp*xp)/cmu - fi0m0=-1.0d0/cmu - fi1m0=-xm/cmu - fi2m0=-(1.0d0/cmu**2+xm*xm)/cmu - do m=1,lrm - fi0p1=-2.0d0*m*(fi1p0-ub*fi0p0)/cmu - fi0m1=-2.0d0*m*(fi1m0-ub*fi0m0)/cmu - fi1p1=-((1.0d0+2.0d0*m)*fi2p0-2.0d0*(m+1.0d0)*ub*fi1p0 - . +up*um*fi0p0)/cmu - fi1m1=-((1.0d0+2.0d0*m)*fi2m0-2.0d0*(m+1.0d0)*ub*fi1m0 - . +up*um*fi0m0)/cmu - fi2p1=(2.0d0*(1.0d0+m)*fi1p1-2.0d0*m* - . (ub*fi2p0-up*um*fi1p0))/cmu - fi2m1=(2.0d0*(1.0d0+m)*fi1m1-2.0d0*m* - . (ub*fi2m0-up*um*fi1m0))/cmu - if(m.ge.n) then - ri(n,0,m)=0.5d0*ci*dnl**m*(fi0p1*eep-fi0m1*eem) - ri(n,1,m)=0.5d0*ci*dnl**m*(fi1p1*eep-fi1m1*eem) - ri(n,2,m)=0.5d0*ci*dnl**m*(fi2p1*eep-fi2m1*eem) - end if - fi0p0=fi0p1 - fi1p0=fi1p1 - fi2p0=fi2p1 - fi0m0=fi0m1 - fi1m0=fi1m1 - fi2m0=fi2m1 - end do - else - ee=exp(-amu*(ygn-1.0d0+anpl*ub)) - call ssbi(aa,n,lrm,fsbi) - do m=n,lrm - cm=rpi*fact(m)*du**(2*m+1) - cim=0.5d0*ci*dnl**m - mm=m-n+1 - fi0m=cm*fsbi(mm) - fi1m=-0.5d0*aa*cm*fsbi(mm+1) - fi2m=0.5d0*cm*(fsbi(mm+1)+0.5d0*aa*aa*fsbi(mm+2)) - ri(n,0,m)=cim*ee*fi0m - ri(n,1,m)=cim*ee*(du*fi1m+ub*fi0m) - ri(n,2,m)=cim*ee*(du*du*fi2m+2.0d0*du*ub*fi1m+ub*ub*fi0m) - end do - end if - end if - end do -c - return - end -c - subroutine ssbi(zz,n,l,fsbi) - implicit none - integer n,l,nmx,lmx,k,m,mm - real*8 eps,c0,c1,sbi,zz - real*8 gamm - parameter(eps=1.0d-10,lmx=20,nmx=lmx+2) - real*8 fsbi(nmx) - do m=n,l+2 - c0=1.0d0/gamm(dble(m)+1.5d0) - sbi=c0 - do k=1,50 - c1=c0*0.25d0*zz**2/(dble(m+k)+0.5d0)/dble(k) - sbi=sbi+c1 - if(c1/sbi.lt.eps) exit - c0=c1 - end do - mm=m-n+1 - fsbi(mm)=sbi - end do - return - end -c -c Weakly relativistic dielectric tensor -c computation of dielectric tensor elements -c Krivenki and Orefice, JPP 30,125 (1983) -c - subroutine diel_tens_wr(ce330,cepsl,lrm) - use const_and_precisions, only : cui=>im - implicit real*8 (a-b,d-h,o-z) - implicit complex*16 (c) - dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2),cepsl(3,3,lrm) -c - common/xgxg/xg - common/ygyg/yg - common/nplr/anpl,anprf - common/amut/amu -c - anpl2=anpl*anpl -c - call fsup(cefp,cefm,lrm) -c - do l=1,lrm - lm=l-1 - fcl=0.5d0**l*((1.0d0/yg)**2/amu)**lm*fact(2*l)/fact(l) - ca11=(0.d0,0.d0) - ca12=(0.d0,0.d0) - ca13=(0.d0,0.d0) - ca22=(0.d0,0.d0) - ca23=(0.d0,0.d0) - ca33=(0.d0,0.d0) - do is=0,l - k=l-is - w=(-1.0d0)**k -c - asl=w/(fact(is+l)*fact(l-is)) - bsl=asl*(is*is+dble(2*k*lm*(l+is))/(2.0d0*l-1.0d0)) -c - cq0p=amu*cefp(is,0) - cq0m=amu*cefm(is,0) - cq1p=amu*anpl*(cefp(is,0)-cefp(is,1)) - cq1m=amu*anpl*(cefm(is,0)-cefm(is,1)) - cq2p=cefp(is,1)+amu*anpl2*(cefp(is,2) - . +cefp(is,0)-2.0d0*cefp(is,1)) -c - ca11=ca11+is**2*asl*cq0p - ca12=ca12+is*l*asl*cq0m - ca22=ca22+bsl*cq0p - ca13=ca13+is*asl*cq1m/yg - ca23=ca23+l*asl*cq1p/yg - ca33=ca33+asl*cq2p/yg**2 - end do - cepsl(1,1,l) = - xg*ca11*fcl - cepsl(1,2,l) = + cui*xg*ca12*fcl - cepsl(2,2,l) = - xg*ca22*fcl - cepsl(1,3,l) = - xg*ca13*fcl - cepsl(2,3,l) = - cui*xg*ca23*fcl - cepsl(3,3,l) = - xg*ca33*fcl - end do -c - cq2p=cefp(0,1)+amu*anpl2*(cefp(0,2)+cefp(0,0)-2.0d0*cefp(0,1)) - ce330=1.0d0-xg*amu*cq2p -c - cepsl(1,1,1) = 1.d0 + cepsl(1,1,1) - cepsl(2,2,1) = 1.d0 + cepsl(2,2,1) -c - do l=1,lrm - cepsl(2,1,l) = - cepsl(1,2,l) - cepsl(3,1,l) = cepsl(1,3,l) - cepsl(3,2,l) = - cepsl(2,3,l) - end do -c - return - end -c -c -c - subroutine fsup(cefp,cefm,lrm) - use const_and_precisions, only : cui=>im - implicit real*8 (a-b,d-h,o-z) - implicit complex*16 (c) - parameter(apsicr=0.7d0) - dimension cefp(0:lrm,0:2),cefm(0:lrm,0:2) -c - common/ygyg/yg - common/nplr/anpl,anprf - common/amut/amu -c - psi=sqrt(0.5d0*amu)*anpl - apsi=abs(psi) -c - do is=0,lrm - alpha=anpl*anpl/2.0d0+is*yg-1.0d0 - phi2=amu*alpha - phim=sqrt(abs(phi2)) - if (alpha.ge.0) then - xp=psi-phim - yp=0.0d0 - xm=-psi-phim - ym=0.0d0 - x0=-phim - y0=0.0d0 - else - xp=psi - yp=phim - xm=-psi - ym=phim - x0=0.0d0 - y0=phim - end if - call zetac (xp,yp,zrp,zip,iflag) - call zetac (xm,ym,zrm,zim,iflag) -c - czp=dcmplx(zrp,zip) - czm=dcmplx(zrm,zim) - cf12=(0.0d0,0.0d0) - if (alpha.ge.0) then - if (alpha.ne.0) cf12=-(czp+czm)/(2.0d0*phim) - else - cf12=-cui*(czp+czm)/(2.0d0*phim) - end if -c - if(apsi.gt.apsicr) then - cf32=-(czp-czm)/(2.0d0*psi) - else - cphi=phim - if(alpha.lt.0) cphi=-cui*phim - call zetac (x0,y0,zr0,zi0,iflag) - cz0=dcmplx(zr0,zi0) - cdz0=2.0d0*(1.0d0-cphi*cz0) - cf32=cdz0 - end if -c - cf0=cf12 - cf1=cf32 - cefp(is,0)=cf32 - cefm(is,0)=cf32 - do l=1,is+2 - iq=l-1 - if(apsi.gt.apsicr) then - cf2=(1.0d0+phi2*cf0-(iq+0.5d0)*cf1)/psi**2 - else - cf2=(1.0d0+phi2*cf1)/dble(iq+1.5d0) - end if - ir=l-is - if(ir.ge.0) then - cefp(is,ir)=cf2 - cefm(is,ir)=cf2 - end if - cf0=cf1 - cf1=cf2 - end do -c - if(is.ne.0) then -c - alpha=anpl*anpl/2.0d0-is*yg-1.0d0 - phi2=amu*alpha - phim=sqrt(abs(phi2)) - if (alpha.ge.0.0d0) then - xp=psi-phim - yp=0.0d0 - xm=-psi-phim - ym=0.0d0 - x0=-phim - y0=0.0d0 - else - xp=psi - yp=phim - xm=-psi - ym=phim - x0=0.0d0 - y0=phim - end if - call zetac (xp,yp,zrp,zip,iflag) - call zetac (xm,ym,zrm,zim,iflag) -c - czp=dcmplx(zrp,zip) - czm=dcmplx(zrm,zim) -c - cf12=(0.0d0,0.0d0) - if (alpha.ge.0) then - if (alpha.ne.0.0d0) cf12=-(czp+czm)/(2.0d0*phim) - else - cf12=-cui*(czp+czm)/(2.0d0*phim) - end if - if(apsi.gt.apsicr) then - cf32=-(czp-czm)/(2.0d0*psi) - else - cphi=phim - if(alpha.lt.0) cphi=-cui*phim - call zetac (x0,y0,zr0,zi0,iflag) - cz0=dcmplx(zr0,zi0) - cdz0=2.0d0*(1.0d0-cphi*cz0) - cf32=cdz0 - end if -c - cf0=cf12 - cf1=cf32 - do l=1,is+2 - iq=l-1 - if(apsi.gt.apsicr) then - cf2=(1.0d0+phi2*cf0-(iq+0.5d0)*cf1)/psi**2 - else - cf2=(1.0d0+phi2*cf1)/dble(iq+1.5d0) - end if - ir=l-is - if(ir.ge.0) then - cefp(is,ir)=cefp(is,ir)+cf2 - cefm(is,ir)=cefm(is,ir)-cf2 - end if - cf0=cf1 - cf1=cf2 - end do -c - end if -c - end do -c - return - end - - - subroutine eccd(yg,anpl,anprre,amu,zeff,rbn,rbx,fc,ex,ey,ez, * dens,psinv,ieccd,nhmn,nhmx,ithn,effjcd,iokhawa,ierr) use green_func_p, only: SpitzFuncCoeff diff --git a/src/grayl.f b/src/grayl.f index 0f0a9f6..b85530a 100644 --- a/src/grayl.f +++ b/src/grayl.f @@ -257,6 +257,25 @@ c end subroutine simpson c c +c + subroutine trapezoid(n,xi,fi,s) +c subroutine for integration with the trapezoidal rule. +c fi: integrand f(x); xi: abscissa x; +c s: integral Int_{xi(1)}^{xi(n)} f(x)dx + implicit none + integer n + real*8 xi(n),fi(n) + real*8 s + integer i +c + s = 0.0d0 + do i = 1, n-1 + s = s+(xi(i+1)-xi(i))*(fi(i+1)-fi(i)) + end do + s = 0.5d0*s + end subroutine trapezoid +c +c c spline routines: begin c c @@ -1633,608 +1652,6 @@ c c c c -* ====================================================================== -* 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) -c---------------------------------------------------------------------- -c -c this fortran 77 packet computes the exponential integrals ei(x), -c e1(x), and exp(-x)*ei(x) for real arguments x where -c -c integral (from t=-infinity to t=x) (exp(t)/t), x > 0, -c ei(x) = -c -integral (from t=-x to t=infinity) (exp(t)/t), x < 0, -c -c and where the first integral is a principal value integral. -c the packet contains three function type subprograms: ei, eone, -c and expei; and one subroutine type subprogram: calcei. the -c calling statements for the primary entries are -c -c y = ei(x), where x .ne. 0, -c -c y = eone(x), where x .gt. 0, -c and -c y = expei(x), where x .ne. 0, -c -c and where the entry points correspond to the functions ei(x), -c e1(x), and exp(-x)*ei(x), respectively. the routine calcei -c is intended for internal packet use only, all computations within -c the packet being concentrated in this routine. the function -c subprograms invoke calcei with the fortran statement -c call calcei(arg,result,int) -c where the parameter usage is as follows -c -c function parameters for calcei -c call arg result int -c -c ei(x) x .ne. 0 ei(x) 1 -c eone(x) x .gt. 0 -ei(-x) 2 -c expei(x) x .ne. 0 exp(-x)*ei(x) 3 -c---------------------------------------------------------------------- - integer i,int - double precision - 1 a,arg,b,c,d,exp40,e,ei,f,four,fourty,frac,half,one,p, - 2 plg,px,p037,p1,p2,q,qlg,qx,q1,q2,r,result,s,six,sump, - 3 sumq,t,three,twelve,two,two4,w,x,xbig,xinf,xmax,xmx0, - 4 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), - 1 s(9),p1(10),q1(9),p2(10),q2(9),plg(4),qlg(4),px(10),qx(10) -c---------------------------------------------------------------------- -c mathematical constants -c exp40 = exp(40) -c x0 = zero of ei -c x01/x11 + x02 = zero of ei to extra precision -c---------------------------------------------------------------------- - data zero,p037,half,one,two/0.0d0,0.037d0,0.5d0,1.0d0,2.0d0/, - 1 three,four,six,twelve,two4/3.0d0,4.0d0,6.0d0,12.d0,24.0d0/, - 2 fourty,exp40/40.0d0,2.3538526683701998541d17/, - 3 x01,x11,x02/381.5d0,1024.0d0,-5.1182968633365538008d-5/, - 4 x0/3.7250741078136663466d-1/ -c---------------------------------------------------------------------- -c machine-dependent constants -c---------------------------------------------------------------------- - data xinf/1.79d+308/,xmax/716.351d0/,xbig/701.84d0/ -c---------------------------------------------------------------------- -c coefficients for -1.0 <= x < 0.0 -c---------------------------------------------------------------------- - data a/1.1669552669734461083368d2, 2.1500672908092918123209d3, - 1 1.5924175980637303639884d4, 8.9904972007457256553251d4, - 2 1.5026059476436982420737d5,-1.4815102102575750838086d5, - 3 5.0196785185439843791020d0/ - data b/4.0205465640027706061433d1, 7.5043163907103936624165d2, - 1 8.1258035174768735759855d3, 5.2440529172056355429883d4, - 2 1.8434070063353677359298d5, 2.5666493484897117319268d5/ -c---------------------------------------------------------------------- -c coefficients for -4.0 <= x < -1.0 -c---------------------------------------------------------------------- - data c/3.828573121022477169108d-1, 1.107326627786831743809d+1, - 1 7.246689782858597021199d+1, 1.700632978311516129328d+2, - 2 1.698106763764238382705d+2, 7.633628843705946890896d+1, - 3 1.487967702840464066613d+1, 9.999989642347613068437d-1, - 4 1.737331760720576030932d-8/ - data d/8.258160008564488034698d-2, 4.344836335509282083360d+0, - 1 4.662179610356861756812d+1, 1.775728186717289799677d+2, - 2 2.953136335677908517423d+2, 2.342573504717625153053d+2, - 3 9.021658450529372642314d+1, 1.587964570758947927903d+1, - 4 1.000000000000000000000d+0/ -c---------------------------------------------------------------------- -c coefficients for x < -4.0 -c---------------------------------------------------------------------- - data e/1.3276881505637444622987d+2,3.5846198743996904308695d+4, - 1 1.7283375773777593926828d+5,2.6181454937205639647381d+5, - 2 1.7503273087497081314708d+5,5.9346841538837119172356d+4, - 3 1.0816852399095915622498d+4,1.0611777263550331766871d03, - 4 5.2199632588522572481039d+1,9.9999999999999999087819d-1/ - data f/3.9147856245556345627078d+4,2.5989762083608489777411d+5, - 1 5.5903756210022864003380d+5,5.4616842050691155735758d+5, - 2 2.7858134710520842139357d+5,7.9231787945279043698718d+4, - 3 1.2842808586627297365998d+4,1.1635769915320848035459d+3, - 4 5.4199632588522559414924d+1,1.0d0/ -c---------------------------------------------------------------------- -c coefficients for rational approximation to ln(x/a), |1-x/a| < .1 -c---------------------------------------------------------------------- - data plg/-2.4562334077563243311d+01,2.3642701335621505212d+02, - 1 -5.4989956895857911039d+02,3.5687548468071500413d+02/ - data qlg/-3.5553900764052419184d+01,1.9400230218539473193d+02, - 1 -3.3442903192607538956d+02,1.7843774234035750207d+02/ -c---------------------------------------------------------------------- -c coefficients for 0.0 < x < 6.0, -c ratio of chebyshev polynomials -c---------------------------------------------------------------------- - data p/-1.2963702602474830028590d01,-1.2831220659262000678155d03, - 1 -1.4287072500197005777376d04,-1.4299841572091610380064d06, - 2 -3.1398660864247265862050d05,-3.5377809694431133484800d08, - 3 3.1984354235237738511048d08,-2.5301823984599019348858d10, - 4 1.2177698136199594677580d10,-2.0829040666802497120940d11/ - data q/ 7.6886718750000000000000d01,-5.5648470543369082846819d03, - 1 1.9418469440759880361415d05,-4.2648434812177161405483d06, - 2 6.4698830956576428587653d07,-7.0108568774215954065376d08, - 3 5.4229617984472955011862d09,-2.8986272696554495342658d10, - 4 9.8900934262481749439886d10,-8.9673749185755048616855d10/ -c---------------------------------------------------------------------- -c j-fraction coefficients for 6.0 <= x < 12.0 -c---------------------------------------------------------------------- - data r/-2.645677793077147237806d00,-2.378372882815725244124d00, - 1 -2.421106956980653511550d01, 1.052976392459015155422d01, - 2 1.945603779539281810439d01,-3.015761863840593359165d01, - 3 1.120011024227297451523d01,-3.988850730390541057912d00, - 4 9.565134591978630774217d00, 9.981193787537396413219d-1/ - data s/ 1.598517957704779356479d-4, 4.644185932583286942650d00, - 1 3.697412299772985940785d02,-8.791401054875438925029d00, - 2 7.608194509086645763123d02, 2.852397548119248700147d01, - 3 4.731097187816050252967d02,-2.369210235636181001661d02, - 4 1.249884822712447891440d00/ -c---------------------------------------------------------------------- -c j-fraction coefficients for 12.0 <= x < 24.0 -c---------------------------------------------------------------------- - data p1/-1.647721172463463140042d00,-1.860092121726437582253d01, - 1 -1.000641913989284829961d01,-2.105740799548040450394d01, - 2 -9.134835699998742552432d-1,-3.323612579343962284333d01, - 3 2.495487730402059440626d01, 2.652575818452799819855d01, - 4 -1.845086232391278674524d00, 9.999933106160568739091d-1/ - data q1/ 9.792403599217290296840d01, 6.403800405352415551324d01, - 1 5.994932325667407355255d01, 2.538819315630708031713d02, - 2 4.429413178337928401161d01, 1.192832423968601006985d03, - 3 1.991004470817742470726d02,-1.093556195391091143924d01, - 4 1.001533852045342697818d00/ -c---------------------------------------------------------------------- -c j-fraction coefficients for x .ge. 24.0 -c---------------------------------------------------------------------- - data p2/ 1.75338801265465972390d02,-2.23127670777632409550d02, - 1 -1.81949664929868906455d01,-2.79798528624305389340d01, - 2 -7.63147701620253630855d00,-1.52856623636929636839d01, - 3 -7.06810977895029358836d00,-5.00006640413131002475d00, - 4 -3.00000000320981265753d00, 1.00000000000000485503d00/ - data q2/ 3.97845977167414720840d04, 3.97277109100414518365d00, - 1 1.37790390235747998793d02, 1.17179220502086455287d02, - 2 7.04831847180424675988d01,-1.20187763547154743238d01, - 3 -7.99243595776339741065d00,-2.99999894040324959612d00, - 4 1.99999999999048104167d00/ -c---------------------------------------------------------------------- - 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 -c---------------------------------------------------------------------- -c calculate ei for negative argument or for e1. -c---------------------------------------------------------------------- - 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 -c---------------------------------------------------------------------- -c to improve conditioning, rational approximations are expressed -c in terms of chebyshev polynomials for 0 <= x < 6, and in -c continued fraction form for larger x. -c---------------------------------------------------------------------- - 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 -c---------------------------------------------------------------------- -c special approximation to ln(x/x0) for x close to x0 -c---------------------------------------------------------------------- - 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 -c---------------------------------------------------------------------- -c calculation reformulated to avoid premature overflow -c---------------------------------------------------------------------- - ei = (ei * exp(x-fourty)) * exp40 - end if - end if - end if - end if - result = ei - return -c---------- last line of calcei ---------- - end - function ei(x) -c-------------------------------------------------------------------- -c -c this function program computes approximate values for the -c exponential integral ei(x), where x is real. -c -c author: w. j. cody -c -c latest modification: january 12, 1988 -c -c-------------------------------------------------------------------- - integer int - double precision ei, x, result -c-------------------------------------------------------------------- - int = 1 - call calcei(x,result,int) - ei = result - return -c---------- last line of ei ---------- - end - function expei(x) -c-------------------------------------------------------------------- -c -c this function program computes approximate values for the -c function exp(-x) * ei(x), where ei(x) is the exponential -c integral, and x is real. -c -c author: w. j. cody -c -c latest modification: january 12, 1988 -c -c-------------------------------------------------------------------- - integer int - double precision expei, x, result -c-------------------------------------------------------------------- - int = 3 - call calcei(x,result,int) - expei = result - return -c---------- last line of expei ---------- - end - function eone(x) -c-------------------------------------------------------------------- -c -c this function program computes approximate values for the -c exponential integral e1(x), where x is real. -c -c author: w. j. cody -c -c latest modification: january 12, 1988 -c -c-------------------------------------------------------------------- - integer int - double precision eone, x, result -c-------------------------------------------------------------------- - int = 2 - call calcei(x,result,int) - eone = result - return -c---------- last line of eone ---------- - end -c -c calcei3 = calcei for int=3 -c -* ====================================================================== - subroutine calcei3(arg,result) -c---------------------------------------------------------------------- -c -c this fortran 77 packet computes the exponential integrals ei(x), -c e1(x), and exp(-x)*ei(x) for real arguments x where -c -c integral (from t=-infinity to t=x) (exp(t)/t), x > 0, -c ei(x) = -c -integral (from t=-x to t=infinity) (exp(t)/t), x < 0, -c -c and where the first integral is a principal value integral. -c the packet contains three function type subprograms: ei, eone, -c and expei; and one subroutine type subprogram: calcei. the -c calling statements for the primary entries are -c -c y = ei(x), where x .ne. 0, -c -c y = eone(x), where x .gt. 0, -c and -c y = expei(x), where x .ne. 0, -c -c and where the entry points correspond to the functions ei(x), -c e1(x), and exp(-x)*ei(x), respectively. the routine calcei -c is intended for internal packet use only, all computations within -c the packet being concentrated in this routine. the function -c subprograms invoke calcei with the fortran statement -c call calcei(arg,result,int) -c where the parameter usage is as follows -c -c function parameters for calcei -c call arg result int -c -c ei(x) x .ne. 0 ei(x) 1 -c eone(x) x .gt. 0 -ei(-x) 2 -c expei(x) x .ne. 0 exp(-x)*ei(x) 3 -c---------------------------------------------------------------------- - integer i,int - double precision - 1 a,arg,b,c,d,exp40,e,ei,f,four,fourty,frac,half,one,p, - 2 plg,px,p037,p1,p2,q,qlg,qx,q1,q2,r,result,s,six,sump, - 3 sumq,t,three,twelve,two,two4,w,x,xbig,xinf,xmax,xmx0, - 4 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), - 1 s(9),p1(10),q1(9),p2(10),q2(9),plg(4),qlg(4),px(10),qx(10) -c---------------------------------------------------------------------- -c mathematical constants -c exp40 = exp(40) -c x0 = zero of ei -c x01/x11 + x02 = zero of ei to extra precision -c---------------------------------------------------------------------- - data zero,p037,half,one,two/0.0d0,0.037d0,0.5d0,1.0d0,2.0d0/, - 1 three,four,six,twelve,two4/3.0d0,4.0d0,6.0d0,12.d0,24.0d0/, - 2 fourty,exp40/40.0d0,2.3538526683701998541d17/, - 3 x01,x11,x02/381.5d0,1024.0d0,-5.1182968633365538008d-5/, - 4 x0/3.7250741078136663466d-1/ -c---------------------------------------------------------------------- -c machine-dependent constants -c---------------------------------------------------------------------- - data xinf/1.79d+308/,xmax/716.351d0/,xbig/701.84d0/ -c---------------------------------------------------------------------- -c coefficients for -1.0 <= x < 0.0 -c---------------------------------------------------------------------- - data a/1.1669552669734461083368d2, 2.1500672908092918123209d3, - 1 1.5924175980637303639884d4, 8.9904972007457256553251d4, - 2 1.5026059476436982420737d5,-1.4815102102575750838086d5, - 3 5.0196785185439843791020d0/ - data b/4.0205465640027706061433d1, 7.5043163907103936624165d2, - 1 8.1258035174768735759855d3, 5.2440529172056355429883d4, - 2 1.8434070063353677359298d5, 2.5666493484897117319268d5/ -c---------------------------------------------------------------------- -c coefficients for -4.0 <= x < -1.0 -c---------------------------------------------------------------------- - data c/3.828573121022477169108d-1, 1.107326627786831743809d+1, - 1 7.246689782858597021199d+1, 1.700632978311516129328d+2, - 2 1.698106763764238382705d+2, 7.633628843705946890896d+1, - 3 1.487967702840464066613d+1, 9.999989642347613068437d-1, - 4 1.737331760720576030932d-8/ - data d/8.258160008564488034698d-2, 4.344836335509282083360d+0, - 1 4.662179610356861756812d+1, 1.775728186717289799677d+2, - 2 2.953136335677908517423d+2, 2.342573504717625153053d+2, - 3 9.021658450529372642314d+1, 1.587964570758947927903d+1, - 4 1.000000000000000000000d+0/ -c---------------------------------------------------------------------- -c coefficients for x < -4.0 -c---------------------------------------------------------------------- - data e/1.3276881505637444622987d+2,3.5846198743996904308695d+4, - 1 1.7283375773777593926828d+5,2.6181454937205639647381d+5, - 2 1.7503273087497081314708d+5,5.9346841538837119172356d+4, - 3 1.0816852399095915622498d+4,1.0611777263550331766871d03, - 4 5.2199632588522572481039d+1,9.9999999999999999087819d-1/ - data f/3.9147856245556345627078d+4,2.5989762083608489777411d+5, - 1 5.5903756210022864003380d+5,5.4616842050691155735758d+5, - 2 2.7858134710520842139357d+5,7.9231787945279043698718d+4, - 3 1.2842808586627297365998d+4,1.1635769915320848035459d+3, - 4 5.4199632588522559414924d+1,1.0d0/ -c---------------------------------------------------------------------- -c coefficients for rational approximation to ln(x/a), |1-x/a| < .1 -c---------------------------------------------------------------------- - data plg/-2.4562334077563243311d+01,2.3642701335621505212d+02, - 1 -5.4989956895857911039d+02,3.5687548468071500413d+02/ - data qlg/-3.5553900764052419184d+01,1.9400230218539473193d+02, - 1 -3.3442903192607538956d+02,1.7843774234035750207d+02/ -c---------------------------------------------------------------------- -c coefficients for 0.0 < x < 6.0, -c ratio of chebyshev polynomials -c---------------------------------------------------------------------- - data p/-1.2963702602474830028590d01,-1.2831220659262000678155d03, - 1 -1.4287072500197005777376d04,-1.4299841572091610380064d06, - 2 -3.1398660864247265862050d05,-3.5377809694431133484800d08, - 3 3.1984354235237738511048d08,-2.5301823984599019348858d10, - 4 1.2177698136199594677580d10,-2.0829040666802497120940d11/ - data q/ 7.6886718750000000000000d01,-5.5648470543369082846819d03, - 1 1.9418469440759880361415d05,-4.2648434812177161405483d06, - 2 6.4698830956576428587653d07,-7.0108568774215954065376d08, - 3 5.4229617984472955011862d09,-2.8986272696554495342658d10, - 4 9.8900934262481749439886d10,-8.9673749185755048616855d10/ -c---------------------------------------------------------------------- -c j-fraction coefficients for 6.0 <= x < 12.0 -c---------------------------------------------------------------------- - data r/-2.645677793077147237806d00,-2.378372882815725244124d00, - 1 -2.421106956980653511550d01, 1.052976392459015155422d01, - 2 1.945603779539281810439d01,-3.015761863840593359165d01, - 3 1.120011024227297451523d01,-3.988850730390541057912d00, - 4 9.565134591978630774217d00, 9.981193787537396413219d-1/ - data s/ 1.598517957704779356479d-4, 4.644185932583286942650d00, - 1 3.697412299772985940785d02,-8.791401054875438925029d00, - 2 7.608194509086645763123d02, 2.852397548119248700147d01, - 3 4.731097187816050252967d02,-2.369210235636181001661d02, - 4 1.249884822712447891440d00/ -c---------------------------------------------------------------------- -c j-fraction coefficients for 12.0 <= x < 24.0 -c---------------------------------------------------------------------- - data p1/-1.647721172463463140042d00,-1.860092121726437582253d01, - 1 -1.000641913989284829961d01,-2.105740799548040450394d01, - 2 -9.134835699998742552432d-1,-3.323612579343962284333d01, - 3 2.495487730402059440626d01, 2.652575818452799819855d01, - 4 -1.845086232391278674524d00, 9.999933106160568739091d-1/ - data q1/ 9.792403599217290296840d01, 6.403800405352415551324d01, - 1 5.994932325667407355255d01, 2.538819315630708031713d02, - 2 4.429413178337928401161d01, 1.192832423968601006985d03, - 3 1.991004470817742470726d02,-1.093556195391091143924d01, - 4 1.001533852045342697818d00/ -c---------------------------------------------------------------------- -c j-fraction coefficients for x .ge. 24.0 -c---------------------------------------------------------------------- - data p2/ 1.75338801265465972390d02,-2.23127670777632409550d02, - 1 -1.81949664929868906455d01,-2.79798528624305389340d01, - 2 -7.63147701620253630855d00,-1.52856623636929636839d01, - 3 -7.06810977895029358836d00,-5.00006640413131002475d00, - 4 -3.00000000320981265753d00, 1.00000000000000485503d00/ - data q2/ 3.97845977167414720840d04, 3.97277109100414518365d00, - 1 1.37790390235747998793d02, 1.17179220502086455287d02, - 2 7.04831847180424675988d01,-1.20187763547154743238d01, - 3 -7.99243595776339741065d00,-2.99999894040324959612d00, - 4 1.99999999999048104167d00/ -c---------------------------------------------------------------------- - data int/ 3/ - x = arg - if (x .eq. zero) then - ei = -xinf - else if ((x .lt. zero)) then -c---------------------------------------------------------------------- -c calculate ei for negative argument or for e1. -c---------------------------------------------------------------------- - 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 -c---------------------------------------------------------------------- -c to improve conditioning, rational approximations are expressed -c in terms of chebyshev polynomials for 0 <= x < 6, and in -c continued fraction form for larger x. -c---------------------------------------------------------------------- - 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 -c---------------------------------------------------------------------- -c special approximation to ln(x/x0) for x close to x0 -c---------------------------------------------------------------------- - 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 -c---------- last line of calcei ---------- - end c c c