istp/gray
2
0
Fork 0
Code Issues Pull Requests Packages Projects Releases 1 Wiki Activity

dispersion, calcei module added

Browse Source
This commit is contained in:
Daniele Micheletti 2015-05-22 09:05:37 +00:00
parent 81b03f9df3
commit 0736f9793a
7 changed files with 1627 additions and 1419 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

3
Makefile
View File

@ -97,10 +97,11 @@ clean:
# Dependencies
# ------------
gray_main.o: const_and_precisions.o
gray-externals.o: green_func_p.o reflections.o beamdata.o const_and_precisions.o
gray-externals.o: green_func_p.o reflections.o beamdata.o const_and_precisions.o dispersion.o
green_func_p.o: const_and_precisions.o
scatterspl.o: const_and_precisions.o
beamdata.o: const_and_precisions.o
dispersion.o: calcei.o
# Special rule to preprocess source file and include svn revision
# ---------------------------------------------------------------

6
Makefile.single
View File

@ -3,7 +3,8 @@ EXE=gray-single
# Objects list
OBJ= gray-single.o gray_main.o gray-externals.o grayl.o beamdata.o \
const_and_precisions.o green_func_p.o reflections.o scatterspl.o
const_and_precisions.o green_func_p.o reflections.o scatterspl.o \
dispersion.o calcei.o
# Alternative search paths
vpath %.f90 src
@ -63,10 +64,11 @@ clean:
# ------------
gray-single.o: gray_main.o grayl.o
gray_main.o: const_and_precisions.o
gray-externals.o: green_func_p.o reflections.o beamdata.o const_and_precisions.o
gray-externals.o: green_func_p.o reflections.o beamdata.o const_and_precisions.o dispersion.o
green_func_p.o: const_and_precisions.o
scatterspl.o: const_and_precisions.o
beamdata.o: const_and_precisions.o
dispersion.o: calcei.o
## library name
## ------------

5
Makefile.standalone
View File

@ -4,7 +4,7 @@ LIBFILE=lib$(EXE).a
# Objects list
OBJ=gray_main.o gray-externals.o grayl.o reflections.o scatterspl.o \
beamdata.o green_func_p.o const_and_precisions.o
beamdata.o green_func_p.o const_and_precisions.o dispersion.o calcei.o
# Alternative search paths
vpath %.f90 src
@ -30,10 +30,11 @@ $(LIBFILE): $(OBJ)
# Dependencies on modules
main.o: const_and_precisions.o
gray_main.o: const_and_precisions.o
gray-externals.o: green_func_p.o reflections.o beamdata.o
gray-externals.o: green_func_p.o reflections.o beamdata.o dispersion.o
green_func_p.o: const_and_precisions.o
scatterspl.o: const_and_precisions.o
beamdata.o: const_and_precisions.o
dispersion.o: calcei.o
# General object compilation command
%.o: %.f90

608
src/calcei.f90 Normal file
View File

@ -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

958
src/dispersion.f90 Normal file
View File

@ -0,0 +1,958 @@
module dispersion
!
use calcei_mod
implicit none
integer, parameter :: dp=kind(1.0d0)
integer, parameter :: tmax=5,npts=500
real(dp), save :: ttv(npts+1),extv(npts+1)
real(dp), parameter :: 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,nprf2,n_f2,dnl,del
!
npl2=npl*npl
dnl=1.0d0-npl2
yg2=yg*yg
del=sqrt(dnl*dnl+4.0d0*npl2*(1.0d0-xg)/yg2)
n_f2=1.0d0-xg-xg*yg2*(1.0d0+npl2+sox*del)/(1.0d0-xg-yg2)/2.0d0
!
nprf2=n_f2-npl2
nprf=sqrt(nprf2)
!
end subroutine colddisp
!
!
!
subroutine larmornumber(yg,mu,npl,larm,nharm,nhf,lrm,fast)
!
! computation of local harmonic number
!
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 :: nharm,nhf ! n=number of the armonic
real(dp) :: ygn ! n*Y
real(dp) :: npl ! parallel N
real(dp) :: gg ! gamma=sqrt(1+u^2)
real(dp) :: dnl,rdu2
real(dp) :: argexp
real(dp) :: mu
integer :: lrm,larm,fast
real(dp) :: u1,u2,uu2,g1,g2
!----------------------------------------
integer :: nn
!----------------------------------------
!
dnl=1.0d0-npl**2
!
if(fast.eq.1) then
ygc=1.0d0-npl**2/2.0d0
else
ygc=sqrt(1.0d0-npl**2)
end if
nharm=int(ygc/yg-eps)+1
if(nharm.eq.0) return
!
nhf=0
!----------------------------------------
! do
! ygn=nharm*yg
!----------------------------------------
do nn=nharm,nharm+4
ygn=nn*yg
!----------------------------------------
if(fast.eq.1) then
rdu2=2.0d0*(ygn-ygc)
u1=npl+sqrt(rdu2)
u2=npl-sqrt(rdu2)
uu2=min(u1*u1,u2*u2)
argexp=mu*uu2/2.0d0
else
rdu2=ygn**2-ygc**2
g1=(ygn+npl*sqrt(rdu2))/dnl
g2=(ygn-npl*sqrt(rdu2))/dnl
gg=min(g1,g2)
argexp=mu*(gg-1.0d0)
u1=(ygn*npl+sqrt(rdu2))/dnl
u2=(ygn*npl-sqrt(rdu2))/dnl
end if
!----------------------------------------
if(argexp.le.expcr) nhf=nn
!----------------------------------------
! if(argexp.gt.expcr) then
! exit
! else
! nhf=nharm
! end if
! nharm=nharm+1
!----------------------------------------
end do
!----------------------------------------
lrm=larm
if(lrm.lt.nhf) lrm=nhf
!----------------------------------------
! lrm=max(larm,nharm)
!----------------------------------------
end subroutine larmornumber
!
!
subroutine warmdisp(xg,yg,mu,npl,nprf,sox,lrm,err,nprr,npri,fast,imx,ex,ey,ez)
implicit none
!
! integer, parameter :: imx=-20 ! number of iterations
real(dp), parameter :: one=1.0d0
complex(dp) :: cc0,cc2,cc4,discr
complex(dp) :: npra2,npra
complex(dp) :: npr,npr2
complex(dp) :: epsl(3,3,lrm),sepsl(3,3),e330
complex(dp) :: e11,e22,e12,e13,e23
complex(dp) :: a13,a31,a23,a32,a33
real(dp) :: errnpr
real(dp) :: nprf,xg,yg,npl,npl2,s,nprr,npri,mu,cr,ci,sox
integer :: err,lrm,ilrm,fast
integer :: i,j,k,imx,imxx
real(dp) :: ez2,ey2,ex2,enx2
complex(dp) :: ex,ey,ez,den
!
!
err=0
errnpr=1.0d0
npra2=nprf**2
npl2=npl*npl
imxx=abs(imx)
!
if (fast.eq.1) then
! call diel_tens_wr(e330,epsl,lrm)
call diel_tens_wr(xg,yg,mu,npl,e330,epsl,lrm)
else
call diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast)
end if
!
do
do i=1,imxx
!
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)*npra2**(ilrm-1)
end do
end do
end do
!
npra=sqrt(npra2)
!
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
! e33=e330+npra2*a33
e13=npra*a13
e23=npra*a23
! e21=-e12
! e31=e13
! e32=-e23
!
! if(i.gt.2.and.errnpr.lt.1.0d-3) exit
!
cc4=(e11-npl2)*(1.0d0-a33)+(a13+npl)*(a31+npl)
cc2=-e12*e12*(1.0d0-a33) &
-a32*e12*(a13+npl)+a23*e12*(a31+npl) &
-(a23*a32+e330+(e22-npl2)*(1.0d0-a33))*(e11-npl2) &
-(a13+npl)*(a31+npl)*(e22-npl2)
cc0=e330*((e11-npl2)*(e22-npl2)+e12*e12)
!
discr=cc2*cc2-4.0d0*cc0*cc4
!
if(yg.gt.1.0d0) then
s=sox
if(dimag(discr).LE.0.0d0) s=-s
else
s=-sox
if(dble(discr).le.0.d0.and.dimag(discr).ge.0.d0) s=-s
end if
!
npr2=(-cc2+s*sqrt(discr))/(2.0d0*cc4)
!
errnpr=abs(1.0d0-abs(npr2)/abs(npra2))
if(i.gt.1.and.errnpr.lt.1.0d-5) exit
npra2=npr2
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(npr2)),npl
imxx=1
else
write(*,"(' X =',f7.4,' Y =',f10.7,' Nperp =',f7.4,': convergence failed.',e12.5)") xg,yg,sqrt(abs(npr2)),npl
exit
end if
else
exit
end if
print*,yg,imx,imxx
end do
!
! if(i.gt.imx) print*,' i>imx ',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)
!
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)
!
dt=2.0d0*tmax/dble(npts)
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*dt
!
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 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

864
src/gray-externals.f
View File

@ -3120,16 +3120,20 @@ c if(ip.eq.nintp) ip=nintp-1
end
subroutine vectinit
c
use const_and_precisions, only : r8
use beamdata, only : jmx,kmx,nmx,psjki,tauv,alphav,pdjki,ppabs,
. currj,didst,ccci,iiv,iop,iow,tau1v,
. nrayr,nrayth,nstep,alloc_beam
implicit real*8 (a-h,o-z)
parameter(tmax=5,npts=500)
real*8 ttv(npts+1),extv(npts+1)
common/ttex/ttv,extv
. currj,didst,ccci,iiv,iop,iow,tau1v,nrayr,nrayth,nstep,alloc_beam
use dispersion, only : expinit
c
implicit none
c internal
integer i,j,k
real(r8) dt
c common
integer ierr,iwarm,ilarm
c
common/warm/iwarm,ilarm
c
common/ierr/ierr
c
if(nstep.gt.nmx) nstep=nmx
@ -3157,17 +3161,12 @@ 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
subroutine vectfree
use beamdata, only : dealloc_beam
implicit none
@ -4996,31 +4995,44 @@ c
c
c
subroutine ecrh_cd
implicit real*8 (a-h,o-z)
real*8 mc2,mesi
parameter(taucr=12.0d0,xxcr=16.0d0,eps=1.d-8)
parameter(qesi=1.602176487d-19,mesi=9.10938215d-31)
parameter(vcsi=2.99792458d+8)
parameter(mc2 = mesi*vcsi*vcsi/qesi*1d-3)
c
use const_and_precisions, only : r8
use dispersion, only : larmornumber,warmdisp
c
implicit none
c parameters
real(r8) taucr,qesi,mesi,vcsi,mc2,xxcr,eps
parameter(taucr=12.0d0,qesi=1.602176487d-19)
parameter(mesi=9.10938215d-31,vcsi=2.99792458d+8)
parameter(mc2 = mesi*vcsi*vcsi/qesi*1d-3)
parameter(xxcr=16.0d0,eps=1.d-8)
c internal
integer lrm,ierr
real(r8) ratiovgr,fzeff,temperature
c common
integer ithn,nharm,nhf,err,iwarm,ilarm,ieccd,imx
real(r8) yg,anpl,anprf,vgm,derdnm,anprre,anprim,alpha,effjcd,
. akim,tau,psinv,tekev,amu,xg,zeff,ak0,akinv,fhz,sox
complex*16 :: ex,ey,ez
C
common/ithn/ithn
common/nharm/nharm,nhf
common/warm/iwarm,ilarm
common/ieccd/ieccd
c
common/ygyg/yg
common/nplr/anpl,anpr
common/nplr/anpl,anprf
common/vgv/vgm,derdnm
common/parwv/ak0,akinv,fhz
c
common/nprw/anprre,anprim
c
common/absor/alpha,effjcd,akim,tau
c
common/psinv/psinv
common/tete/tekev
common/amut/amu
common/xgxg/xg
common/epolar/ex,ey,ez
common/warm/iwarm,ilarm
common/ieccd/ieccd
common/zz/Zeff
common/parwv/ak0,akinv,fhz
common/mode/sox
common/imx/imx
c
c absorption computation
c
@ -5030,50 +5042,19 @@ c
c
tekev=temperature(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
call larmornumber(yg,amu,anpl,ilarm,nharm,nhf,lrm,iwarm)
if(nhf.eq.0.or.nharm.eq.0) return
nharm=nharm-1
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)
call warmdisp(xg,yg,amu,anpl,anprf,sox,
. lrm,err,anprre,anprim,iwarm,imx,ex,ey,ez)
c
akim=ak0*anprim
ratiovgr=2.0d0*anpr/derdnm
c ratiovgr=2.0d0*anpr/derdnm*vgm
ratiovgr=2.0d0*anprf/derdnm
c ratiovgr=2.0d0*anprf/derdnm*vgm
alpha=2.0d0*akim*ratiovgr
if(alpha.lt.0.0d0) then
ierr=94
@ -5081,420 +5062,15 @@ c ratiovgr=2.0d0*anpr/derdnm*vgm
c
ithn=1
if(lrm.gt.nharm) ithn=2
c
zeff=fzeff(psinv)
if(ieccd.gt.0) call eccd(effjcd)
return
999 format(12(1x,e12.5))
end
c
c
c
subroutine dispersion(lrm)
implicit real*8(a-h,o-z)
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/amut/amu
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
99 format(20(1x,e12.5))
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)
implicit real*8(a-h,o-z)
complex*16 ui
complex*16 e330,epsl(3,3,lrm)
complex*16 ca11,ca12,ca22,ca13,ca23,ca33
complex*16 cq0p,cq0m,cq1p,cq1m,cq2p
parameter(pi=3.14159265358979d0,rpi=1.77245385090552d0)
parameter(ui=(0.0d0,1.0d0))
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
@ -5641,6 +5217,9 @@ c
function fhermit(t)
use calcei_mod, only : calcei3
implicit real*8 (a-h,o-z)
c
common/ygyg/yg
@ -5682,345 +5261,6 @@ c
c
c
c
subroutine antihermitian(ri,lrm)
implicit none
integer lmx,nmx,lrm,n,k,m,mm
real*8 rpi
parameter(rpi=1.77245385090552d0)
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)
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)
parameter(cui=(0.0d0,1.0d0))
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)
implicit real*8 (a-b,d-h,o-z)
implicit complex*16 (c)
parameter(apsicr=0.7d0)
parameter(cui=(0.0d0,1.0d0))
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

602
src/grayl.f
View File

@ -1652,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 eint(x),
c e1(x), and exp(-x)*eint(x) for real arguments x where
c
c integral (from t=-infinity to t=x) (exp(t)/t), x > 0,
c eint(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 = eint(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 eint(x),
c e1(x), and exp(-x)*eint(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 eint(x) x .ne. 0 eint(x) 1
c eone(x) x .gt. 0 -eint(-x) 2
c expei(x) x .ne. 0 exp(-x)*eint(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 eint(x)
c--------------------------------------------------------------------
c
c this function program computes approximate values for the
c exponential integral eint(x), where x is real.
c
c author: w. j. cody
c
c latest modification: january 12, 1988
c
c--------------------------------------------------------------------
integer int
double precision eint, x, result
c--------------------------------------------------------------------
int = 1
call calcei(x,result,int)
eint = 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) * eint(x), where eint(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 eint(x),
c e1(x), and exp(-x)*eint(x) for real arguments x where
c
c integral (from t=-infinity to t=x) (exp(t)/t), x > 0,
c eint(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 = eint(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 eint(x),
c e1(x), and exp(-x)*eint(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 eint(x) x .ne. 0 eint(x) 1
c eone(x) x .gt. 0 -eint(-x) 2
c expei(x) x .ne. 0 exp(-x)*eint(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