From 9a56975e75296d157ba2de9d286c727ac1408df3 Mon Sep 17 00:00:00 2001 From: Daniele Micheletti Date: Fri, 22 May 2015 10:14:36 +0000 Subject: [PATCH] hermitian_2 and fhermit moved to dispersion module, dqagmv module added --- Makefile | 2 +- Makefile.single | 4 +- Makefile.standalone | 5 +- src/dispersion.f90 | 217 +++++- src/dqagmv.f | 1695 ++++++++++++++++++++++++++++++++++++++++++ src/gray-externals.f | 193 ----- 6 files changed, 1917 insertions(+), 199 deletions(-) create mode 100644 src/dqagmv.f diff --git a/Makefile b/Makefile index 0157a6c..f33d9c0 100644 --- a/Makefile +++ b/Makefile @@ -101,7 +101,7 @@ gray-externals.o: green_func_p.o reflections.o beamdata.o const_and_precisions.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 +dispersion.o: calcei.o dqagmv.o # Special rule to preprocess source file and include svn revision # --------------------------------------------------------------- diff --git a/Makefile.single b/Makefile.single index 2131e4b..4b4893b 100644 --- a/Makefile.single +++ b/Makefile.single @@ -4,7 +4,7 @@ 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 \ - dispersion.o calcei.o + dispersion.o calcei.o dqagmv.o # Alternative search paths vpath %.f90 src @@ -68,7 +68,7 @@ gray-externals.o: green_func_p.o reflections.o beamdata.o const_and_precisions.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 +dispersion.o: calcei.o dqagmv.o ## library name ## ------------ diff --git a/Makefile.standalone b/Makefile.standalone index a4cfa81..3eece41 100644 --- a/Makefile.standalone +++ b/Makefile.standalone @@ -4,7 +4,8 @@ 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 dispersion.o calcei.o + beamdata.o green_func_p.o const_and_precisions.o dispersion.o \ + calcei.o dqagmv.o # Alternative search paths vpath %.f90 src @@ -34,7 +35,7 @@ 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 +dispersion.o: calcei.o dqagmv.o # General object compilation command %.o: %.f90 diff --git a/src/dispersion.f90 b/src/dispersion.f90 index 7656c2a..ef073d1 100644 --- a/src/dispersion.f90 +++ b/src/dispersion.f90 @@ -306,7 +306,19 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast) end do end do ! - call hermitian(rr,yg,mu,npl,cr,ci,fast,lrm) +! 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(*,*) 'fast dispersion flag unknown value:', fast + + end select ! call antihermitian(ri,yg,mu,npl,cr,ci,lrm) ! @@ -572,8 +584,211 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast) end subroutine hermitian ! ! + subroutine hermitian_2(rr,yg,mu,npl,cr,ci,fast,lrm) +! + implicit none +! parameters + integer :: lw,liw,npar + parameter(lw=5000,liw=lw/4,npar=7) + real(dp) :: epsa,epsr + 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 + 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 + apar(5) = real(n,dp) + apar(6) = real(m,dp) + apar(7) = real(ih,dp) +! + do n=n1,llm + nn=abs(n) + do m=nn,llm + if(n.eq.0.and.m.eq.0) then + ih=2 +! 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,2,m) = resfh + else + do ih=0,2 +! 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 if + 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 +! + return + end ! ! + function fhermit(t,apar,npar) +! + implicit none +! arguments + integer npar + real(dp) t,fhermit + real(dp), dimension(npar) :: 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 + 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) + 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 + return + end ! ! subroutine antihermitian(ri,yg,mu,npl,cr,ci,lrm) diff --git a/src/dqagmv.f b/src/dqagmv.f new file mode 100644 index 0000000..eda8b61 --- /dev/null +++ b/src/dqagmv.f @@ -0,0 +1,1695 @@ +c +c +c Integration routine dqags.f from quadpack and dependencies: BEGIN +c +c + subroutine dqagsmv(f,a,b,apar,np,epsabs,epsrel,result,abserr, + * neval,ier,limit,lenw,last,iwork,work) +c***begin prologue dqagsmv +c***date written 800101 (yymmdd) +c***revision date 830518 (yymmdd) +c***category no. h2a1a1 +c***keywords automatic integrator, general-purpose, +c (end-point) singularities, extrapolation, +c globally adaptive +c***author piessens,robert,appl. math. & progr. div. - k.u.leuven +c de doncker,elise,appl. math. & prog. div. - k.u.leuven +c***purpose the routine calculates an approximation result to a given +c definite integral i = integral of f over (a,b), +c hopefully satisfying following claim for accuracy +c abs(i-result).le.max(epsabs,epsrel*abs(i)). +c***description +c +c computation of a definite integral +c standard fortran subroutine +c double precision version +c +c +c parameters +c on entry +c f - double precision +c function subprogram defining the integrand +c function f(x,apar,np). the actual name for f needs +c to be declared e x t e r n a l in the driver +c program. +c +c a - double precision +c lower limit of integration +c +c b - double precision +c upper limit of integration +c +c apar - array of parameters of the integrand function f +c +c np - number of parameters. size of apar +c +c epsabs - double precision +c absolute accuracy requested +c epsrel - double precision +c relative accuracy requested +c if epsabs.le.0 +c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), +c the routine will end with ier = 6. +c +c on return +c result - double precision +c approximation to the integral +c +c abserr - double precision +c estimate of the modulus of the absolute error, +c which should equal or exceed abs(i-result) +c +c neval - integer +c number of integrand evaluations +c +c ier - integer +c ier = 0 normal and reliable termination of the +c routine. it is assumed that the requested +c accuracy has been achieved. +c ier.gt.0 abnormal termination of the routine +c the estimates for integral and error are +c less reliable. it is assumed that the +c requested accuracy has not been achieved. +c error messages +c ier = 1 maximum number of subdivisions allowed +c has been achieved. one can allow more sub- +c divisions by increasing the value of limit +c (and taking the according dimension +c adjustments into account. however, if +c this yields no improvement it is advised +c to analyze the integrand in order to +c determine the integration difficulties. if +c the position of a local difficulty can be +c determined (e.g. singularity, +c discontinuity within the interval) one +c will probably gain from splitting up the +c interval at this point and calling the +c integrator on the subranges. if possible, +c an appropriate special-purpose integrator +c should be used, which is designed for +c handling the type of difficulty involved. +c = 2 the occurrence of roundoff error is detec- +c ted, which prevents the requested +c tolerance from being achieved. +c the error may be under-estimated. +c = 3 extremely bad integrand behaviour +c occurs at some points of the integration +c interval. +c = 4 the algorithm does not converge. +c roundoff error is detected in the +c extrapolation table. it is presumed that +c the requested tolerance cannot be +c achieved, and that the returned result is +c the best which can be obtained. +c = 5 the integral is probably divergent, or +c slowly convergent. it must be noted that +c divergence can occur with any other value +c of ier. +c = 6 the input is invalid, because +c (epsabs.le.0 and +c epsrel.lt.max(50*rel.mach.acc.,0.5d-28) +c or limit.lt.1 or lenw.lt.limit*4. +c result, abserr, neval, last are set to +c zero.except when limit or lenw is invalid, +c iwork(1), work(limit*2+1) and +c work(limit*3+1) are set to zero, work(1) +c is set to a and work(limit+1) to b. +c +c dimensioning parameters +c limit - integer +c dimensioning parameter for iwork +c limit determines the maximum number of subintervals +c in the partition of the given integration interval +c (a,b), limit.ge.1. +c if limit.lt.1, the routine will end with ier = 6. +c +c lenw - integer +c dimensioning parameter for work +c lenw must be at least limit*4. +c if lenw.lt.limit*4, the routine will end +c with ier = 6. +c +c last - integer +c on return, last equals the number of subintervals +c produced in the subdivision process, detemines the +c number of significant elements actually in the work +c arrays. +c +c work arrays +c iwork - integer +c vector of dimension at least limit, the first k +c elements of which contain pointers +c to the error estimates over the subintervals +c such that work(limit*3+iwork(1)),... , +c work(limit*3+iwork(k)) form a decreasing +c sequence, with k = last if last.le.(limit/2+2), +c and k = limit+1-last otherwise +c +c work - double precision +c vector of dimension at least lenw +c on return +c work(1), ..., work(last) contain the left +c end-points of the subintervals in the +c partition of (a,b), +c work(limit+1), ..., work(limit+last) contain +c the right end-points, +c work(limit*2+1), ..., work(limit*2+last) contain +c the integral approximations over the subintervals, +c work(limit*3+1), ..., work(limit*3+last) +c contain the error estimates. +c +c***references (none) +c***routines called dqagsemv,xerror +c***end prologue dqagsmv +c +c + double precision a,abserr,b,epsabs,epsrel,f,result,work,apar + integer ier,iwork,last,lenw,limit,lvl,l1,l2,l3,neval,np +c + dimension iwork(limit),work(lenw),apar(np) +c + external f +c +c check validity of limit and lenw. +c +c***first executable statement dqags + ier = 6 + neval = 0 + last = 0 + result = 0.0d+00 + abserr = 0.0d+00 + if(limit.lt.1.or.lenw.lt.limit*4) go to 10 +c +c prepare call for dqagse. +c + l1 = limit+1 + l2 = limit+l1 + l3 = limit+l2 +c + call dqagsemv(f,a,b,apar,np,epsabs,epsrel,limit,result,abserr, + * neval,ier,work(1),work(l1),work(l2),work(l3),iwork,last) +c +c call error handler if necessary. +c + lvl = 0 +10 if(ier.eq.6) lvl = 1 + if(ier.ne.0) print*,'habnormal return from dqags',ier,lvl + return + end +c + subroutine dqagsemv(f,a,b,apar,np,epsabs,epsrel,limit,result, + * abserr,neval,ier,alist,blist,rlist,elist,iord,last) +c***begin prologue dqagsemv +c***date written 800101 (yymmdd) +c***revision date 830518 (yymmdd) +c***category no. h2a1a1 +c***keywords automatic integrator, general-purpose, +c (end point) singularities, extrapolation, +c globally adaptive +c***author piessens,robert,appl. math. & progr. div. - k.u.leuven +c de doncker,elise,appl. math. & progr. div. - k.u.leuven +c***purpose the routine calculates an approximation result to a given +c definite integral i = integral of f over (a,b), +c hopefully satisfying following claim for accuracy +c abs(i-result).le.max(epsabs,epsrel*abs(i)). +c***description +c +c computation of a definite integral +c standard fortran subroutine +c double precision version +c +c parameters +c on entry +c f - double precision +c function subprogram defining the integrand +c function f(x,apar,np). the actual name for f needs +c to be declared e x t e r n a l in the driver +c program. +c +c a - double precision +c lower limit of integration +c +c b - double precision +c upper limit of integration +c +c apar - array of parameters of the integrand function f +c +c np - number of parameters. size of apar +c +c epsabs - double precision +c absolute accuracy requested +c epsrel - double precision +c relative accuracy requested +c if epsabs.le.0 +c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), +c the routine will end with ier = 6. +c +c limit - integer +c gives an upperbound on the number of subintervals +c in the partition of (a,b) +c +c on return +c result - double precision +c approximation to the integral +c +c abserr - double precision +c estimate of the modulus of the absolute error, +c which should equal or exceed abs(i-result) +c +c neval - integer +c number of integrand evaluations +c +c ier - integer +c ier = 0 normal and reliable termination of the +c routine. it is assumed that the requested +c accuracy has been achieved. +c ier.gt.0 abnormal termination of the routine +c the estimates for integral and error are +c less reliable. it is assumed that the +c requested accuracy has not been achieved. +c error messages +c = 1 maximum number of subdivisions allowed +c has been achieved. one can allow more sub- +c divisions by increasing the value of limit +c (and taking the according dimension +c adjustments into account). however, if +c this yields no improvement it is advised +c to analyze the integrand in order to +c determine the integration difficulties. if +c the position of a local difficulty can be +c determined (e.g. singularity, +c discontinuity within the interval) one +c will probably gain from splitting up the +c interval at this point and calling the +c integrator on the subranges. if possible, +c an appropriate special-purpose integrator +c should be used, which is designed for +c handling the type of difficulty involved. +c = 2 the occurrence of roundoff error is detec- +c ted, which prevents the requested +c tolerance from being achieved. +c the error may be under-estimated. +c = 3 extremely bad integrand behaviour +c occurs at some points of the integration +c interval. +c = 4 the algorithm does not converge. +c roundoff error is detected in the +c extrapolation table. +c it is presumed that the requested +c tolerance cannot be achieved, and that the +c returned result is the best which can be +c obtained. +c = 5 the integral is probably divergent, or +c slowly convergent. it must be noted that +c divergence can occur with any other value +c of ier. +c = 6 the input is invalid, because +c epsabs.le.0 and +c epsrel.lt.max(50*rel.mach.acc.,0.5d-28). +c result, abserr, neval, last, rlist(1), +c iord(1) and elist(1) are set to zero. +c alist(1) and blist(1) are set to a and b +c respectively. +c +c alist - double precision +c vector of dimension at least limit, the first +c last elements of which are the left end points +c of the subintervals in the partition of the +c given integration range (a,b) +c +c blist - double precision +c vector of dimension at least limit, the first +c last elements of which are the right end points +c of the subintervals in the partition of the given +c integration range (a,b) +c +c rlist - double precision +c vector of dimension at least limit, the first +c last elements of which are the integral +c approximations on the subintervals +c +c elist - double precision +c vector of dimension at least limit, the first +c last elements of which are the moduli of the +c absolute error estimates on the subintervals +c +c iord - integer +c vector of dimension at least limit, the first k +c elements of which are pointers to the +c error estimates over the subintervals, +c such that elist(iord(1)), ..., elist(iord(k)) +c form a decreasing sequence, with k = last +c if last.le.(limit/2+2), and k = limit+1-last +c otherwise +c +c last - integer +c number of subintervals actually produced in the +c subdivision process +c +c***references (none) +c***routines called d1mach,dqelg,dqk21mv,dqpsrt +c***end prologue dqagsemv +c + double precision a,abseps,abserr,alist,area,area1,area12,area2,a1, + * a2,b,blist,b1,b2,correc,dabs,defabs,defab1,defab2,d1mach,dmax1, + * dres,elist,epmach,epsabs,epsrel,erlarg,erlast,errbnd,errmax, + * error1,error2,erro12,errsum,ertest,f,oflow,resabs,reseps,result, + * res3la,rlist,rlist2,small,uflow,apar + integer id,ier,ierro,iord,iroff1,iroff2,iroff3,jupbnd,k,ksgn, + * ktmin,last,limit,maxerr,neval,nres,nrmax,numrl2,np + logical extrap,noext +c + dimension alist(limit),blist(limit),elist(limit),iord(limit), + * res3la(3),rlist(limit),rlist2(52),apar(np) +c + external f +c +c the dimension of rlist2 is determined by the value of +c limexp in subroutine dqelg (rlist2 should be of dimension +c (limexp+2) at least). +c +c list of major variables +c ----------------------- +c +c alist - list of left end points of all subintervals +c considered up to now +c blist - list of right end points of all subintervals +c considered up to now +c rlist(i) - approximation to the integral over +c (alist(i),blist(i)) +c rlist2 - array of dimension at least limexp+2 containing +c the part of the epsilon table which is still +c needed for further computations +c elist(i) - error estimate applying to rlist(i) +c maxerr - pointer to the interval with largest error +c estimate +c errmax - elist(maxerr) +c erlast - error on the interval currently subdivided +c (before that subdivision has taken place) +c area - sum of the integrals over the subintervals +c errsum - sum of the errors over the subintervals +c errbnd - requested accuracy max(epsabs,epsrel* +c abs(result)) +c *****1 - variable for the left interval +c *****2 - variable for the right interval +c last - index for subdivision +c nres - number of calls to the extrapolation routine +c numrl2 - number of elements currently in rlist2. if an +c appropriate approximation to the compounded +c integral has been obtained it is put in +c rlist2(numrl2) after numrl2 has been increased +c by one. +c small - length of the smallest interval considered up +c to now, multiplied by 1.5 +c erlarg - sum of the errors over the intervals larger +c than the smallest interval considered up to now +c extrap - logical variable denoting that the routine is +c attempting to perform extrapolation i.e. before +c subdividing the smallest interval we try to +c decrease the value of erlarg. +c noext - logical variable denoting that extrapolation +c is no longer allowed (true value) +c +c machine dependent constants +c --------------------------- +c +c epmach is the largest relative spacing. +c uflow is the smallest positive magnitude. +c oflow is the largest positive magnitude. +c +c***first executable statement dqagse + epmach = d1mach(4) +c +c test on validity of parameters +c ------------------------------ + ier = 0 + neval = 0 + last = 0 + result = 0.0d+00 + abserr = 0.0d+00 + alist(1) = a + blist(1) = b + rlist(1) = 0.0d+00 + elist(1) = 0.0d+00 + if(epsabs.le.0.0d+00.and.epsrel.lt.dmax1(0.5d+02*epmach,0.5d-28)) + * ier = 6 + if(ier.eq.6) go to 999 +c +c first approximation to the integral +c ----------------------------------- +c + uflow = d1mach(1) + oflow = d1mach(2) + ierro = 0 + call dqk21mv(f,a,b,apar,np,result,abserr,defabs,resabs) +c +c test on accuracy. +c + dres = dabs(result) + errbnd = dmax1(epsabs,epsrel*dres) + last = 1 + rlist(1) = result + elist(1) = abserr + iord(1) = 1 + if(abserr.le.1.0d+02*epmach*defabs.and.abserr.gt.errbnd) ier = 2 + if(limit.eq.1) ier = 1 + if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs).or. + * abserr.eq.0.0d+00) go to 140 +c +c initialization +c -------------- +c + rlist2(1) = result + errmax = abserr + maxerr = 1 + area = result + errsum = abserr + abserr = oflow + nrmax = 1 + nres = 0 + numrl2 = 2 + ktmin = 0 + extrap = .false. + noext = .false. + iroff1 = 0 + iroff2 = 0 + iroff3 = 0 + ksgn = -1 + if(dres.ge.(0.1d+01-0.5d+02*epmach)*defabs) ksgn = 1 +c +c main do-loop +c ------------ +c + do 90 last = 2,limit +c +c bisect the subinterval with the nrmax-th largest error +c estimate. +c + a1 = alist(maxerr) + b1 = 0.5d+00*(alist(maxerr)+blist(maxerr)) + a2 = b1 + b2 = blist(maxerr) + erlast = errmax + call dqk21mv(f,a1,b1,apar,np,area1,error1,resabs,defab1) + call dqk21mv(f,a2,b2,apar,np,area2,error2,resabs,defab2) +c +c improve previous approximations to integral +c and error and test for accuracy. +c + area12 = area1+area2 + erro12 = error1+error2 + errsum = errsum+erro12-errmax + area = area+area12-rlist(maxerr) + if(defab1.eq.error1.or.defab2.eq.error2) go to 15 + if(dabs(rlist(maxerr)-area12).gt.0.1d-04*dabs(area12) + * .or.erro12.lt.0.99d+00*errmax) go to 10 + if(extrap) iroff2 = iroff2+1 + if(.not.extrap) iroff1 = iroff1+1 + 10 if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1 + 15 rlist(maxerr) = area1 + rlist(last) = area2 + errbnd = dmax1(epsabs,epsrel*dabs(area)) +c +c test for roundoff error and eventually set error flag. +c + if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2 + if(iroff2.ge.5) ierro = 3 +c +c set error flag in the case that the number of subintervals +c equals limit. +c + if(last.eq.limit) ier = 1 +c +c set error flag in the case of bad integrand behaviour +c at a point of the integration range. +c + if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03*epmach)* + * (dabs(a2)+0.1d+04*uflow)) ier = 4 +c +c append the newly-created intervals to the list. +c + if(error2.gt.error1) go to 20 + alist(last) = a2 + blist(maxerr) = b1 + blist(last) = b2 + elist(maxerr) = error1 + elist(last) = error2 + go to 30 + 20 alist(maxerr) = a2 + alist(last) = a1 + blist(last) = b1 + rlist(maxerr) = area2 + rlist(last) = area1 + elist(maxerr) = error2 + elist(last) = error1 +c +c call subroutine dqpsrt to maintain the descending ordering +c in the list of error estimates and select the subinterval +c with nrmax-th largest error estimate (to be bisected next). +c + 30 call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) +c ***jump out of do-loop + if(errsum.le.errbnd) go to 115 +c ***jump out of do-loop + if(ier.ne.0) go to 100 + if(last.eq.2) go to 80 + if(noext) go to 90 + erlarg = erlarg-erlast + if(dabs(b1-a1).gt.small) erlarg = erlarg+erro12 + if(extrap) go to 40 +c +c test whether the interval to be bisected next is the +c smallest interval. +c + if(dabs(blist(maxerr)-alist(maxerr)).gt.small) go to 90 + extrap = .true. + nrmax = 2 + 40 if(ierro.eq.3.or.erlarg.le.ertest) go to 60 +c +c the smallest interval has the largest error. +c before bisecting decrease the sum of the errors over the +c larger intervals (erlarg) and perform extrapolation. +c + id = nrmax + jupbnd = last + if(last.gt.(2+limit/2)) jupbnd = limit+3-last + do 50 k = id,jupbnd + maxerr = iord(nrmax) + errmax = elist(maxerr) +c ***jump out of do-loop + if(dabs(blist(maxerr)-alist(maxerr)).gt.small) go to 90 + nrmax = nrmax+1 + 50 continue +c +c perform extrapolation. +c + 60 numrl2 = numrl2+1 + rlist2(numrl2) = area + call dqelg(numrl2,rlist2,reseps,abseps,res3la,nres) + ktmin = ktmin+1 + if(ktmin.gt.5.and.abserr.lt.0.1d-02*errsum) ier = 5 + if(abseps.ge.abserr) go to 70 + ktmin = 0 + abserr = abseps + result = reseps + correc = erlarg + ertest = dmax1(epsabs,epsrel*dabs(reseps)) +c ***jump out of do-loop + if(abserr.le.ertest) go to 100 +c +c prepare bisection of the smallest interval. +c + 70 if(numrl2.eq.1) noext = .true. + if(ier.eq.5) go to 100 + maxerr = iord(1) + errmax = elist(maxerr) + nrmax = 1 + extrap = .false. + small = small*0.5d+00 + erlarg = errsum + go to 90 + 80 small = dabs(b-a)*0.375d+00 + erlarg = errsum + ertest = errbnd + rlist2(2) = area + 90 continue +c +c set final result and error estimate. +c ------------------------------------ +c + 100 if(abserr.eq.oflow) go to 115 + if(ier+ierro.eq.0) go to 110 + if(ierro.eq.3) abserr = abserr+correc + if(ier.eq.0) ier = 3 + if(result.ne.0.0d+00.and.area.ne.0.0d+00) go to 105 + if(abserr.gt.errsum) go to 115 + if(area.eq.0.0d+00) go to 130 + go to 110 + 105 if(abserr/dabs(result).gt.errsum/dabs(area)) go to 115 +c +c test on divergence. +c + 110 if(ksgn.eq.(-1).and.dmax1(dabs(result),dabs(area)).le. + * defabs*0.1d-01) go to 130 + if(0.1d-01.gt.(result/area).or.(result/area).gt.0.1d+03 + * .or.errsum.gt.dabs(area)) ier = 6 + go to 130 +c +c compute global integral sum. +c + 115 result = 0.0d+00 + do 120 k = 1,last + result = result+rlist(k) + 120 continue + abserr = errsum + 130 if(ier.gt.2) ier = ier-1 + 140 neval = 42*last-21 + 999 return + end +c + subroutine dqk21mv(f,a,b,apar,np,result,abserr,resabs,resasc) +c***begin prologue dqk21mv +c***date written 800101 (yymmdd) +c***revision date 830518 (yymmdd) +c***category no. h2a1a2 +c***keywords 21-point gauss-kronrod rules +c***author piessens,robert,appl. math. & progr. div. - k.u.leuven +c de doncker,elise,appl. math. & progr. div. - k.u.leuven +c***purpose to compute i = integral of f over (a,b), with error +c estimate +c j = integral of abs(f) over (a,b) +c***description +c +c integration rules +c standard fortran subroutine +c double precision version +c +c parameters +c on entry +c f - double precision +c function subprogram defining the integrand +c function f(x,apar,np). the actual name for f +c needs to be declared e x t e r n a l in the +c driver program. +c +c a - double precision +c lower limit of integration +c +c b - double precision +c upper limit of integration +c +c apar - array of parameters of the integrand function f +c +c np - number of parameters. size of apar +c +c on return +c result - double precision +c approximation to the integral i +c result is computed by applying the 21-point +c kronrod rule (resk) obtained by optimal addition +c of abscissae to the 10-point gauss rule (resg). +c +c abserr - double precision +c estimate of the modulus of the absolute error, +c which should not exceed abs(i-result) +c +c resabs - double precision +c approximation to the integral j +c +c resasc - double precision +c approximation to the integral of abs(f-i/(b-a)) +c over (a,b) +c +c***references (none) +c***routines called d1mach +c***end prologue dqk21mv +c + double precision a,absc,abserr,b,centr,dabs,dhlgth,dmax1,dmin1, + * d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth,resabs,resasc, + * resg,resk,reskh,result,uflow,wg,wgk,xgk,apar + integer j,jtw,jtwm1,np + external f +c + dimension fv1(10),fv2(10),wg(5),wgk(11),xgk(11),apar(np) +c +c the abscissae and weights are given for the interval (-1,1). +c because of symmetry only the positive abscissae and their +c corresponding weights are given. +c +c xgk - abscissae of the 21-point kronrod rule +c xgk(2), xgk(4), ... abscissae of the 10-point +c gauss rule +c xgk(1), xgk(3), ... abscissae which are optimally +c added to the 10-point gauss rule +c +c wgk - weights of the 21-point kronrod rule +c +c wg - weights of the 10-point gauss rule +c +c +c gauss quadrature weights and kronron quadrature abscissae and weights +c as evaluated with 80 decimal digit arithmetic by l. w. fullerton, +c bell labs, nov. 1981. +c + data wg ( 1) / 0.0666713443 0868813759 3568809893 332 d0 / + data wg ( 2) / 0.1494513491 5058059314 5776339657 697 d0 / + data wg ( 3) / 0.2190863625 1598204399 5534934228 163 d0 / + data wg ( 4) / 0.2692667193 0999635509 1226921569 469 d0 / + data wg ( 5) / 0.2955242247 1475287017 3892994651 338 d0 / +c + data xgk ( 1) / 0.9956571630 2580808073 5527280689 003 d0 / + data xgk ( 2) / 0.9739065285 1717172007 7964012084 452 d0 / + data xgk ( 3) / 0.9301574913 5570822600 1207180059 508 d0 / + data xgk ( 4) / 0.8650633666 8898451073 2096688423 493 d0 / + data xgk ( 5) / 0.7808177265 8641689706 3717578345 042 d0 / + data xgk ( 6) / 0.6794095682 9902440623 4327365114 874 d0 / + data xgk ( 7) / 0.5627571346 6860468333 9000099272 694 d0 / + data xgk ( 8) / 0.4333953941 2924719079 9265943165 784 d0 / + data xgk ( 9) / 0.2943928627 0146019813 1126603103 866 d0 / + data xgk ( 10) / 0.1488743389 8163121088 4826001129 720 d0 / + data xgk ( 11) / 0.0000000000 0000000000 0000000000 000 d0 / +c + data wgk ( 1) / 0.0116946388 6737187427 8064396062 192 d0 / + data wgk ( 2) / 0.0325581623 0796472747 8818972459 390 d0 / + data wgk ( 3) / 0.0547558965 7435199603 1381300244 580 d0 / + data wgk ( 4) / 0.0750396748 1091995276 7043140916 190 d0 / + data wgk ( 5) / 0.0931254545 8369760553 5065465083 366 d0 / + data wgk ( 6) / 0.1093871588 0229764189 9210590325 805 d0 / + data wgk ( 7) / 0.1234919762 6206585107 7958109831 074 d0 / + data wgk ( 8) / 0.1347092173 1147332592 8054001771 707 d0 / + data wgk ( 9) / 0.1427759385 7706008079 7094273138 717 d0 / + data wgk ( 10) / 0.1477391049 0133849137 4841515972 068 d0 / + data wgk ( 11) / 0.1494455540 0291690566 4936468389 821 d0 / +c +c +c list of major variables +c ----------------------- +c +c centr - mid point of the interval +c hlgth - half-length of the interval +c absc - abscissa +c fval* - function value +c resg - result of the 10-point gauss formula +c resk - result of the 21-point kronrod formula +c reskh - approximation to the mean value of f over (a,b), +c i.e. to i/(b-a) +c +c +c machine dependent constants +c --------------------------- +c +c epmach is the largest relative spacing. +c uflow is the smallest positive magnitude. +c +c***first executable statement dqk21 + epmach = d1mach(4) + uflow = d1mach(1) +c + centr = 0.5d+00*(a+b) + hlgth = 0.5d+00*(b-a) + dhlgth = dabs(hlgth) +c +c compute the 21-point kronrod approximation to +c the integral, and estimate the absolute error. +c + resg = 0.0d+00 + fc = f(centr,apar,np) + resk = wgk(11)*fc + resabs = dabs(resk) + do 10 j=1,5 + jtw = 2*j + absc = hlgth*xgk(jtw) + fval1 = f(centr-absc,apar,np) + fval2 = f(centr+absc,apar,np) + fv1(jtw) = fval1 + fv2(jtw) = fval2 + fsum = fval1+fval2 + resg = resg+wg(j)*fsum + resk = resk+wgk(jtw)*fsum + resabs = resabs+wgk(jtw)*(dabs(fval1)+dabs(fval2)) + 10 continue + do 15 j = 1,5 + jtwm1 = 2*j-1 + absc = hlgth*xgk(jtwm1) + fval1 = f(centr-absc,apar,np) + fval2 = f(centr+absc,apar,np) + fv1(jtwm1) = fval1 + fv2(jtwm1) = fval2 + fsum = fval1+fval2 + resk = resk+wgk(jtwm1)*fsum + resabs = resabs+wgk(jtwm1)*(dabs(fval1)+dabs(fval2)) + 15 continue + reskh = resk*0.5d+00 + resasc = wgk(11)*dabs(fc-reskh) + do 20 j=1,10 + resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh)) + 20 continue + result = resk*hlgth + resabs = resabs*dhlgth + resasc = resasc*dhlgth + abserr = dabs((resk-resg)*hlgth) + if(resasc.ne.0.0d+00.and.abserr.ne.0.0d+00) + * abserr = resasc*dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00) + if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1 + * ((epmach*0.5d+02)*resabs,abserr) + return + end +c + subroutine dqagimv(f,bound,inf,apar,np,epsabs,epsrel,result, + * abserr,neval,ier,limit,lenw,last,iwork,work) +c***begin prologue dqagimv +c***date written 800101 (yymmdd) +c***revision date 830518 (yymmdd) +c***category no. h2a3a1,h2a4a1 +c***keywords automatic integrator, infinite intervals, +c general-purpose, transformation, extrapolation, +c globally adaptive +c***author piessens,robert,appl. math. & progr. div. - k.u.leuven +c de doncker,elise,appl. math. & progr. div. -k.u.leuven +c***purpose the routine calculates an approximation result to a given +c integral i = integral of f over (bound,+infinity) +c or i = integral of f over (-infinity,bound) +c or i = integral of f over (-infinity,+infinity) +c hopefully satisfying following claim for accuracy +c abs(i-result).le.max(epsabs,epsrel*abs(i)). +c***description +c +c integration over infinite intervals +c standard fortran subroutine +c +c parameters +c on entry +c f - double precision +c function subprogram defining the integrand +c function f(x,apar,np). the actual name for f needs +c to be declared e x t e r n a l in the driver +c program. +c +c bound - double precision +c finite bound of integration range +c (has no meaning if interval is doubly-infinite) +c +c inf - integer +c indicating the kind of integration range involved +c inf = 1 corresponds to (bound,+infinity), +c inf = -1 to (-infinity,bound), +c inf = 2 to (-infinity,+infinity). +c +c apar - array of parameters of the integrand function f +c +c np - number of parameters. size of apar +c +c epsabs - double precision +c absolute accuracy requested +c epsrel - double precision +c relative accuracy requested +c if epsabs.le.0 +c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), +c the routine will end with ier = 6. +c +c +c on return +c result - double precision +c approximation to the integral +c +c abserr - double precision +c estimate of the modulus of the absolute error, +c which should equal or exceed abs(i-result) +c +c neval - integer +c number of integrand evaluations +c +c ier - integer +c ier = 0 normal and reliable termination of the +c routine. it is assumed that the requested +c accuracy has been achieved. +c - ier.gt.0 abnormal termination of the routine. the +c estimates for result and error are less +c reliable. it is assumed that the requested +c accuracy has not been achieved. +c error messages +c ier = 1 maximum number of subdivisions allowed +c has been achieved. one can allow more +c subdivisions by increasing the value of +c limit (and taking the according dimension +c adjustments into account). however, if +c this yields no improvement it is advised +c to analyze the integrand in order to +c determine the integration difficulties. if +c the position of a local difficulty can be +c determined (e.g. singularity, +c discontinuity within the interval) one +c will probably gain from splitting up the +c interval at this point and calling the +c integrator on the subranges. if possible, +c an appropriate special-purpose integrator +c should be used, which is designed for +c handling the type of difficulty involved. +c = 2 the occurrence of roundoff error is +c detected, which prevents the requested +c tolerance from being achieved. +c the error may be under-estimated. +c = 3 extremely bad integrand behaviour occurs +c at some points of the integration +c interval. +c = 4 the algorithm does not converge. +c roundoff error is detected in the +c extrapolation table. +c it is assumed that the requested tolerance +c cannot be achieved, and that the returned +c result is the best which can be obtained. +c = 5 the integral is probably divergent, or +c slowly convergent. it must be noted that +c divergence can occur with any other value +c of ier. +c = 6 the input is invalid, because +c (epsabs.le.0 and +c epsrel.lt.max(50*rel.mach.acc.,0.5d-28)) +c or limit.lt.1 or leniw.lt.limit*4. +c result, abserr, neval, last are set to +c zero. exept when limit or leniw is +c invalid, iwork(1), work(limit*2+1) and +c work(limit*3+1) are set to zero, work(1) +c is set to a and work(limit+1) to b. +c +c dimensioning parameters +c limit - integer +c dimensioning parameter for iwork +c limit determines the maximum number of subintervals +c in the partition of the given integration interval +c (a,b), limit.ge.1. +c if limit.lt.1, the routine will end with ier = 6. +c +c lenw - integer +c dimensioning parameter for work +c lenw must be at least limit*4. +c if lenw.lt.limit*4, the routine will end +c with ier = 6. +c +c last - integer +c on return, last equals the number of subintervals +c produced in the subdivision process, which +c determines the number of significant elements +c actually in the work arrays. +c +c work arrays +c iwork - integer +c vector of dimension at least limit, the first +c k elements of which contain pointers +c to the error estimates over the subintervals, +c such that work(limit*3+iwork(1)),... , +c work(limit*3+iwork(k)) form a decreasing +c sequence, with k = last if last.le.(limit/2+2), and +c k = limit+1-last otherwise +c +c work - double precision +c vector of dimension at least lenw +c on return +c work(1), ..., work(last) contain the left +c end points of the subintervals in the +c partition of (a,b), +c work(limit+1), ..., work(limit+last) contain +c the right end points, +c work(limit*2+1), ...,work(limit*2+last) contain the +c integral approximations over the subintervals, +c work(limit*3+1), ..., work(limit*3) +c contain the error estimates. +c***references (none) +c***routines called dqagiemv,xerror +c***end prologue dqagimv +c + double precision abserr,bound,epsabs,epsrel,f,result,work,apar + integer ier,inf,iwork,last,lenw,limit,lvl,l1,l2,l3,neval,np +c + dimension iwork(limit),work(lenw),apar(np) +c + external f +c +c check validity of limit and lenw. +c +c***first executable statement dqagi + ier = 6 + neval = 0 + last = 0 + result = 0.0d+00 + abserr = 0.0d+00 + if(limit.lt.1.or.lenw.lt.limit*4) go to 10 +c +c prepare call for dqagie. +c + l1 = limit+1 + l2 = limit+l1 + l3 = limit+l2 +c + call dqagiemv(f,bound,inf,apar,np,epsabs,epsrel,limit,result, + * abserr,neval,ier,work(1),work(l1),work(l2),work(l3),iwork,last) +c +c call error handler if necessary. +c + lvl = 0 +10 if(ier.eq.6) lvl = 1 + if(ier.ne.0) print*,'habnormal return from dqagi' + return + end +c + subroutine dqagiemv(f,bound,inf,apar,np,epsabs,epsrel,limit, + * result,abserr,neval,ier,alist,blist,rlist,elist,iord,last) +c***begin prologue dqagiemv +c***date written 800101 (yymmdd) +c***revision date 830518 (yymmdd) +c***category no. h2a3a1,h2a4a1 +c***keywords automatic integrator, infinite intervals, +c general-purpose, transformation, extrapolation, +c globally adaptive +c***author piessens,robert,appl. math & progr. div - k.u.leuven +c de doncker,elise,appl. math & progr. div - k.u.leuven +c***purpose the routine calculates an approximation result to a given +c integral i = integral of f over (bound,+infinity) +c or i = integral of f over (-infinity,bound) +c or i = integral of f over (-infinity,+infinity), +c hopefully satisfying following claim for accuracy +c abs(i-result).le.max(epsabs,epsrel*abs(i)) +c***description +c +c integration over infinite intervals +c standard fortran subroutine +c +c f - double precision +c function subprogram defining the integrand +c function f(x,apar,np). the actual name for f needs +c to be declared e x t e r n a l in the driver +c program. +c +c bound - double precision +c finite bound of integration range +c (has no meaning if interval is doubly-infinite) +c +c inf - double precision +c indicating the kind of integration range involved +c inf = 1 corresponds to (bound,+infinity), +c inf = -1 to (-infinity,bound), +c inf = 2 to (-infinity,+infinity). +c +c apar - array of parameters of the integrand function f +c +c np - number of parameters. size of apar +c +c epsabs - double precision +c absolute accuracy requested +c epsrel - double precision +c relative accuracy requested +c if epsabs.le.0 +c and epsrel.lt.max(50*rel.mach.acc.,0.5d-28), +c the routine will end with ier = 6. +c +c limit - integer +c gives an upper bound on the number of subintervals +c in the partition of (a,b), limit.ge.1 +c +c on return +c result - double precision +c approximation to the integral +c +c abserr - double precision +c estimate of the modulus of the absolute error, +c which should equal or exceed abs(i-result) +c +c neval - integer +c number of integrand evaluations +c +c ier - integer +c ier = 0 normal and reliable termination of the +c routine. it is assumed that the requested +c accuracy has been achieved. +c - ier.gt.0 abnormal termination of the routine. the +c estimates for result and error are less +c reliable. it is assumed that the requested +c accuracy has not been achieved. +c error messages +c ier = 1 maximum number of subdivisions allowed +c has been achieved. one can allow more +c subdivisions by increasing the value of +c limit (and taking the according dimension +c adjustments into account). however,if +c this yields no improvement it is advised +c to analyze the integrand in order to +c determine the integration difficulties. +c if the position of a local difficulty can +c be determined (e.g. singularity, +c discontinuity within the interval) one +c will probably gain from splitting up the +c interval at this point and calling the +c integrator on the subranges. if possible, +c an appropriate special-purpose integrator +c should be used, which is designed for +c handling the type of difficulty involved. +c = 2 the occurrence of roundoff error is +c detected, which prevents the requested +c tolerance from being achieved. +c the error may be under-estimated. +c = 3 extremely bad integrand behaviour occurs +c at some points of the integration +c interval. +c = 4 the algorithm does not converge. +c roundoff error is detected in the +c extrapolation table. +c it is assumed that the requested tolerance +c cannot be achieved, and that the returned +c result is the best which can be obtained. +c = 5 the integral is probably divergent, or +c slowly convergent. it must be noted that +c divergence can occur with any other value +c of ier. +c = 6 the input is invalid, because +c (epsabs.le.0 and +c epsrel.lt.max(50*rel.mach.acc.,0.5d-28), +c result, abserr, neval, last, rlist(1), +c elist(1) and iord(1) are set to zero. +c alist(1) and blist(1) are set to 0 +c and 1 respectively. +c +c alist - double precision +c vector of dimension at least limit, the first +c last elements of which are the left +c end points of the subintervals in the partition +c of the transformed integration range (0,1). +c +c blist - double precision +c vector of dimension at least limit, the first +c last elements of which are the right +c end points of the subintervals in the partition +c of the transformed integration range (0,1). +c +c rlist - double precision +c vector of dimension at least limit, the first +c last elements of which are the integral +c approximations on the subintervals +c +c elist - double precision +c vector of dimension at least limit, the first +c last elements of which are the moduli of the +c absolute error estimates on the subintervals +c +c iord - integer +c vector of dimension limit, the first k +c elements of which are pointers to the +c error estimates over the subintervals, +c such that elist(iord(1)), ..., elist(iord(k)) +c form a decreasing sequence, with k = last +c if last.le.(limit/2+2), and k = limit+1-last +c otherwise +c +c last - integer +c number of subintervals actually produced +c in the subdivision process +c +c***references (none) +c***routines called d1mach,dqelg,dqk15imv,dqpsrt +c***end prologue dqagiemv + double precision abseps,abserr,alist,area,area1,area12,area2,a1, + * a2,blist,boun,bound,b1,b2,correc,dabs,defabs,defab1,defab2, + * dmax1,dres,d1mach,elist,epmach,epsabs,epsrel,erlarg,erlast, + * errbnd,errmax,error1,error2,erro12,errsum,ertest,f,oflow,resabs, + * reseps,result,res3la,rlist,rlist2,small,uflow,apar + integer id,ier,ierro,inf,iord,iroff1,iroff2,iroff3,jupbnd,k,ksgn, + * ktmin,last,limit,maxerr,neval,nres,nrmax,numrl2,np + logical extrap,noext +c + dimension alist(limit),blist(limit),elist(limit),iord(limit), + * res3la(3),rlist(limit),rlist2(52),apar(np) +c + external f +c +c the dimension of rlist2 is determined by the value of +c limexp in subroutine dqelg. +c +c +c list of major variables +c ----------------------- +c +c alist - list of left end points of all subintervals +c considered up to now +c blist - list of right end points of all subintervals +c considered up to now +c rlist(i) - approximation to the integral over +c (alist(i),blist(i)) +c rlist2 - array of dimension at least (limexp+2), +c containing the part of the epsilon table +c wich is still needed for further computations +c elist(i) - error estimate applying to rlist(i) +c maxerr - pointer to the interval with largest error +c estimate +c errmax - elist(maxerr) +c erlast - error on the interval currently subdivided +c (before that subdivision has taken place) +c area - sum of the integrals over the subintervals +c errsum - sum of the errors over the subintervals +c errbnd - requested accuracy max(epsabs,epsrel* +c abs(result)) +c *****1 - variable for the left subinterval +c *****2 - variable for the right subinterval +c last - index for subdivision +c nres - number of calls to the extrapolation routine +c numrl2 - number of elements currently in rlist2. if an +c appropriate approximation to the compounded +c integral has been obtained, it is put in +c rlist2(numrl2) after numrl2 has been increased +c by one. +c small - length of the smallest interval considered up +c to now, multiplied by 1.5 +c erlarg - sum of the errors over the intervals larger +c than the smallest interval considered up to now +c extrap - logical variable denoting that the routine +c is attempting to perform extrapolation. i.e. +c before subdividing the smallest interval we +c try to decrease the value of erlarg. +c noext - logical variable denoting that extrapolation +c is no longer allowed (true-value) +c +c machine dependent constants +c --------------------------- +c +c epmach is the largest relative spacing. +c uflow is the smallest positive magnitude. +c oflow is the largest positive magnitude. +c +c***first executable statement dqagie + epmach = d1mach(4) +c +c test on validity of parameters +c ----------------------------- +c + ier = 0 + neval = 0 + last = 0 + result = 0.0d+00 + abserr = 0.0d+00 + alist(1) = 0.0d+00 + blist(1) = 0.1d+01 + rlist(1) = 0.0d+00 + elist(1) = 0.0d+00 + iord(1) = 0 + if(epsabs.le.0.0d+00.and.epsrel.lt.dmax1(0.5d+02*epmach,0.5d-28)) + * ier = 6 + if(ier.eq.6) go to 999 +c +c +c first approximation to the integral +c ----------------------------------- +c +c determine the interval to be mapped onto (0,1). +c if inf = 2 the integral is computed as i = i1+i2, where +c i1 = integral of f over (-infinity,0), +c i2 = integral of f over (0,+infinity). +c + boun = bound + if(inf.eq.2) boun = 0.0d+00 + call dqk15imv(f,boun,inf,0.0d+00,0.1d+01,apar,np,result,abserr, + * defabs,resabs) +c +c test on accuracy +c + last = 1 + rlist(1) = result + elist(1) = abserr + iord(1) = 1 + dres = dabs(result) + errbnd = dmax1(epsabs,epsrel*dres) + if(abserr.le.1.0d+02*epmach*defabs.and.abserr.gt.errbnd) ier = 2 + if(limit.eq.1) ier = 1 + if(ier.ne.0.or.(abserr.le.errbnd.and.abserr.ne.resabs).or. + * abserr.eq.0.0d+00) go to 130 +c +c initialization +c -------------- +c + uflow = d1mach(1) + oflow = d1mach(2) + rlist2(1) = result + errmax = abserr + maxerr = 1 + area = result + errsum = abserr + abserr = oflow + nrmax = 1 + nres = 0 + ktmin = 0 + numrl2 = 2 + extrap = .false. + noext = .false. + ierro = 0 + iroff1 = 0 + iroff2 = 0 + iroff3 = 0 + ksgn = -1 + if(dres.ge.(0.1d+01-0.5d+02*epmach)*defabs) ksgn = 1 +c +c main do-loop +c ------------ +c + do 90 last = 2,limit +c +c bisect the subinterval with nrmax-th largest error estimate. +c + a1 = alist(maxerr) + b1 = 0.5d+00*(alist(maxerr)+blist(maxerr)) + a2 = b1 + b2 = blist(maxerr) + erlast = errmax + call dqk15imv(f,boun,inf,a1,b1,apar,np,area1,error1,resabs, + * defab1) + call dqk15imv(f,boun,inf,a2,b2,apar,np,area2,error2,resabs, + * defab2) +c +c improve previous approximations to integral +c and error and test for accuracy. +c + area12 = area1+area2 + erro12 = error1+error2 + errsum = errsum+erro12-errmax + area = area+area12-rlist(maxerr) + if(defab1.eq.error1.or.defab2.eq.error2)go to 15 + if(dabs(rlist(maxerr)-area12).gt.0.1d-04*dabs(area12) + * .or.erro12.lt.0.99d+00*errmax) go to 10 + if(extrap) iroff2 = iroff2+1 + if(.not.extrap) iroff1 = iroff1+1 + 10 if(last.gt.10.and.erro12.gt.errmax) iroff3 = iroff3+1 + 15 rlist(maxerr) = area1 + rlist(last) = area2 + errbnd = dmax1(epsabs,epsrel*dabs(area)) +c +c test for roundoff error and eventually set error flag. +c + if(iroff1+iroff2.ge.10.or.iroff3.ge.20) ier = 2 + if(iroff2.ge.5) ierro = 3 +c +c set error flag in the case that the number of +c subintervals equals limit. +c + if(last.eq.limit) ier = 1 +c +c set error flag in the case of bad integrand behaviour +c at some points of the integration range. +c + if(dmax1(dabs(a1),dabs(b2)).le.(0.1d+01+0.1d+03*epmach)* + * (dabs(a2)+0.1d+04*uflow)) ier = 4 +c +c append the newly-created intervals to the list. +c + if(error2.gt.error1) go to 20 + alist(last) = a2 + blist(maxerr) = b1 + blist(last) = b2 + elist(maxerr) = error1 + elist(last) = error2 + go to 30 + 20 alist(maxerr) = a2 + alist(last) = a1 + blist(last) = b1 + rlist(maxerr) = area2 + rlist(last) = area1 + elist(maxerr) = error2 + elist(last) = error1 +c +c call subroutine dqpsrt to maintain the descending ordering +c in the list of error estimates and select the subinterval +c with nrmax-th largest error estimate (to be bisected next). +c + 30 call dqpsrt(limit,last,maxerr,errmax,elist,iord,nrmax) + if(errsum.le.errbnd) go to 115 + if(ier.ne.0) go to 100 + if(last.eq.2) go to 80 + if(noext) go to 90 + erlarg = erlarg-erlast + if(dabs(b1-a1).gt.small) erlarg = erlarg+erro12 + if(extrap) go to 40 +c +c test whether the interval to be bisected next is the +c smallest interval. +c + if(dabs(blist(maxerr)-alist(maxerr)).gt.small) go to 90 + extrap = .true. + nrmax = 2 + 40 if(ierro.eq.3.or.erlarg.le.ertest) go to 60 +c +c the smallest interval has the largest error. +c before bisecting decrease the sum of the errors over the +c larger intervals (erlarg) and perform extrapolation. +c + id = nrmax + jupbnd = last + if(last.gt.(2+limit/2)) jupbnd = limit+3-last + do 50 k = id,jupbnd + maxerr = iord(nrmax) + errmax = elist(maxerr) + if(dabs(blist(maxerr)-alist(maxerr)).gt.small) go to 90 + nrmax = nrmax+1 + 50 continue +c +c perform extrapolation. +c + 60 numrl2 = numrl2+1 + rlist2(numrl2) = area + call dqelg(numrl2,rlist2,reseps,abseps,res3la,nres) + ktmin = ktmin+1 + if(ktmin.gt.5.and.abserr.lt.0.1d-02*errsum) ier = 5 + if(abseps.ge.abserr) go to 70 + ktmin = 0 + abserr = abseps + result = reseps + correc = erlarg + ertest = dmax1(epsabs,epsrel*dabs(reseps)) + if(abserr.le.ertest) go to 100 +c +c prepare bisection of the smallest interval. +c + 70 if(numrl2.eq.1) noext = .true. + if(ier.eq.5) go to 100 + maxerr = iord(1) + errmax = elist(maxerr) + nrmax = 1 + extrap = .false. + small = small*0.5d+00 + erlarg = errsum + go to 90 + 80 small = 0.375d+00 + erlarg = errsum + ertest = errbnd + rlist2(2) = area + 90 continue +c +c set final result and error estimate. +c ------------------------------------ +c + 100 if(abserr.eq.oflow) go to 115 + if((ier+ierro).eq.0) go to 110 + if(ierro.eq.3) abserr = abserr+correc + if(ier.eq.0) ier = 3 + if(result.ne.0.0d+00.and.area.ne.0.0d+00)go to 105 + if(abserr.gt.errsum)go to 115 + if(area.eq.0.0d+00) go to 130 + go to 110 + 105 if(abserr/dabs(result).gt.errsum/dabs(area))go to 115 +c +c test on divergence +c + 110 if(ksgn.eq.(-1).and.dmax1(dabs(result),dabs(area)).le. + * defabs*0.1d-01) go to 130 + if(0.1d-01.gt.(result/area).or.(result/area).gt.0.1d+03. + *or.errsum.gt.dabs(area)) ier = 6 + go to 130 +c +c compute global integral sum. +c + 115 result = 0.0d+00 + do 120 k = 1,last + result = result+rlist(k) + 120 continue + abserr = errsum + 130 neval = 30*last-15 + if(inf.eq.2) neval = 2*neval + if(ier.gt.2) ier=ier-1 + 999 return + end +c + subroutine dqk15imv(f,boun,inf,a,b,apar,np,result,abserr,resabs, + * resasc) +c***begin prologue dqk15imv +c***date written 800101 (yymmdd) +c***revision date 830518 (yymmdd) +c***category no. h2a3a2,h2a4a2 +c***keywords 15-point transformed gauss-kronrod rules +c***author piessens,robert,appl. math. & progr. div. - k.u.leuven +c de doncker,elise,appl. math. & progr. div. - k.u.leuven +c***purpose the original (infinite integration range is mapped +c onto the interval (0,1) and (a,b) is a part of (0,1). +c it is the purpose to compute +c i = integral of transformed integrand over (a,b), +c j = integral of abs(transformed integrand) over (a,b). +c***description +c +c integration rule +c standard fortran subroutine +c double precision version +c +c parameters +c on entry +c f - double precision +c fuction subprogram defining the integrand +c function f(x,apar,np). the actual name for f +c needs to be declared e x t e r n a l in the +c calling program. +c +c boun - double precision +c finite bound of original integration +c range (set to zero if inf = +2) +c +c inf - integer +c if inf = -1, the original interval is +c (-infinity,bound), +c if inf = +1, the original interval is +c (bound,+infinity), +c if inf = +2, the original interval is +c (-infinity,+infinity) and +c the integral is computed as the sum of two +c integrals, one over (-infinity,0) and one over +c (0,+infinity). +c +c a - double precision +c lower limit for integration over subrange +c of (0,1) +c +c b - double precision +c upper limit for integration over subrange +c of (0,1) +c +c apar - array of parameters of the integrand function f +c +c np - number of parameters. size of apar +c +c on return +c result - double precision +c approximation to the integral i +c result is computed by applying the 15-point +c kronrod rule(resk) obtained by optimal addition +c of abscissae to the 7-point gauss rule(resg). +c +c abserr - double precision +c estimate of the modulus of the absolute error, +c which should equal or exceed abs(i-result) +c +c resabs - double precision +c approximation to the integral j +c +c resasc - double precision +c approximation to the integral of +c abs((transformed integrand)-i/(b-a)) over (a,b) +c +c***references (none) +c***routines called d1mach +c***end prologue dqk15imv +c + double precision a,absc,absc1,absc2,abserr,b,boun,centr,dabs,dinf, + * dmax1,dmin1,d1mach,epmach,f,fc,fsum,fval1,fval2,fv1,fv2,hlgth, + * resabs,resasc,resg,resk,reskh,result,tabsc1,tabsc2,uflow,wg,wgk, + * xgk,apar + integer inf,j,np + external f +c + dimension fv1(7),fv2(7),xgk(8),wgk(8),wg(8),apar(np) +c +c the abscissae and weights are supplied for the interval +c (-1,1). because of symmetry only the positive abscissae and +c their corresponding weights are given. +c +c xgk - abscissae of the 15-point kronrod rule +c xgk(2), xgk(4), ... abscissae of the 7-point +c gauss rule +c xgk(1), xgk(3), ... abscissae which are optimally +c added to the 7-point gauss rule +c +c wgk - weights of the 15-point kronrod rule +c +c wg - weights of the 7-point gauss rule, corresponding +c to the abscissae xgk(2), xgk(4), ... +c wg(1), wg(3), ... are set to zero. +c + data wg(1) / 0.0d0 / + data wg(2) / 0.1294849661 6886969327 0611432679 082d0 / + data wg(3) / 0.0d0 / + data wg(4) / 0.2797053914 8927666790 1467771423 780d0 / + data wg(5) / 0.0d0 / + data wg(6) / 0.3818300505 0511894495 0369775488 975d0 / + data wg(7) / 0.0d0 / + data wg(8) / 0.4179591836 7346938775 5102040816 327d0 / +c + data xgk(1) / 0.9914553711 2081263920 6854697526 329d0 / + data xgk(2) / 0.9491079123 4275852452 6189684047 851d0 / + data xgk(3) / 0.8648644233 5976907278 9712788640 926d0 / + data xgk(4) / 0.7415311855 9939443986 3864773280 788d0 / + data xgk(5) / 0.5860872354 6769113029 4144838258 730d0 / + data xgk(6) / 0.4058451513 7739716690 6606412076 961d0 / + data xgk(7) / 0.2077849550 0789846760 0689403773 245d0 / + data xgk(8) / 0.0000000000 0000000000 0000000000 000d0 / +c + data wgk(1) / 0.0229353220 1052922496 3732008058 970d0 / + data wgk(2) / 0.0630920926 2997855329 0700663189 204d0 / + data wgk(3) / 0.1047900103 2225018383 9876322541 518d0 / + data wgk(4) / 0.1406532597 1552591874 5189590510 238d0 / + data wgk(5) / 0.1690047266 3926790282 6583426598 550d0 / + data wgk(6) / 0.1903505780 6478540991 3256402421 014d0 / + data wgk(7) / 0.2044329400 7529889241 4161999234 649d0 / + data wgk(8) / 0.2094821410 8472782801 2999174891 714d0 / +c +c +c list of major variables +c ----------------------- +c +c centr - mid point of the interval +c hlgth - half-length of the interval +c absc* - abscissa +c tabsc* - transformed abscissa +c fval* - function value +c resg - result of the 7-point gauss formula +c resk - result of the 15-point kronrod formula +c reskh - approximation to the mean value of the transformed +c integrand over (a,b), i.e. to i/(b-a) +c +c machine dependent constants +c --------------------------- +c +c epmach is the largest relative spacing. +c uflow is the smallest positive magnitude. +c +c***first executable statement dqk15imv + epmach = d1mach(4) + uflow = d1mach(1) + dinf = min0(1,inf) +c + centr = 0.5d+00*(a+b) + hlgth = 0.5d+00*(b-a) + tabsc1 = boun+dinf*(0.1d+01-centr)/centr + fval1 = f(tabsc1,apar,np) + if(inf.eq.2) fval1 = fval1+f(-tabsc1,apar,np) + fc = (fval1/centr)/centr +c +c compute the 15-point kronrod approximation to +c the integral, and estimate the error. +c + resg = wg(8)*fc + resk = wgk(8)*fc + resabs = dabs(resk) + do 10 j=1,7 + absc = hlgth*xgk(j) + absc1 = centr-absc + absc2 = centr+absc + tabsc1 = boun+dinf*(0.1d+01-absc1)/absc1 + tabsc2 = boun+dinf*(0.1d+01-absc2)/absc2 + fval1 = f(tabsc1,apar,np) + fval2 = f(tabsc2,apar,np) + if(inf.eq.2) fval1 = fval1+f(-tabsc1,apar,np) + if(inf.eq.2) fval2 = fval2+f(-tabsc2,apar,np) + fval1 = (fval1/absc1)/absc1 + fval2 = (fval2/absc2)/absc2 + fv1(j) = fval1 + fv2(j) = fval2 + fsum = fval1+fval2 + resg = resg+wg(j)*fsum + resk = resk+wgk(j)*fsum + resabs = resabs+wgk(j)*(dabs(fval1)+dabs(fval2)) + 10 continue + reskh = resk*0.5d+00 + resasc = wgk(8)*dabs(fc-reskh) + do 20 j=1,7 + resasc = resasc+wgk(j)*(dabs(fv1(j)-reskh)+dabs(fv2(j)-reskh)) + 20 continue + result = resk*hlgth + resasc = resasc*hlgth + resabs = resabs*hlgth + abserr = dabs((resk-resg)*hlgth) + if(resasc.ne.0.0d+00.and.abserr.ne.0.d0) abserr = resasc* + * dmin1(0.1d+01,(0.2d+03*abserr/resasc)**1.5d+00) + if(resabs.gt.uflow/(0.5d+02*epmach)) abserr = dmax1 + * ((epmach*0.5d+02)*resabs,abserr) + return + end diff --git a/src/gray-externals.f b/src/gray-externals.f index 3bb21e8..7ea0ca9 100644 --- a/src/gray-externals.f +++ b/src/gray-externals.f @@ -5071,199 +5071,6 @@ c c c c - - - - subroutine hermitian_2(rr,lrm) - implicit real*8(a-h,o-z) - parameter(tmax=5,npts=500) - 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/cri/cr,ci - 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 - bth=sqrt(bth2) - amu2=amu*amu - amu4=amu2*amu2 - amu6=amu4*amu2 -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) - - use calcei_mod, only : calcei3 - - 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 - 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) - 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 eccd(effjcd) use green_func_p implicit none