From df2faa8213b151af45cdf203be9fa66f0a32c5ca Mon Sep 17 00:00:00 2001 From: Lorenzo Figini Date: Fri, 22 May 2015 10:48:24 +0000 Subject: [PATCH] fixed parameters passing to fhermit in hermitian_2 --- src/dispersion.f90 | 191 +++++++++++++++++++++------------------------ 1 file changed, 90 insertions(+), 101 deletions(-) diff --git a/src/dispersion.f90 b/src/dispersion.f90 index ef073d1..85defcc 100644 --- a/src/dispersion.f90 +++ b/src/dispersion.f90 @@ -312,11 +312,11 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast) case(2:3) call hermitian(rr,yg,mu,npl,cr,ci,fast,lrm) - case(4) + case(4:) call hermitian_2(rr,yg,mu,npl,cr,ci,fast,lrm) case default - write(*,*) 'fast dispersion flag unknown value:', fast + write(*,*) "unexpected value for flag 'fast' in dispersion:", fast end select ! @@ -585,19 +585,17 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast) ! ! 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) + integer,parameter :: lw=5000,liw=lw/4 + integer, parameter :: npar=7 + real(dp), parameter :: epsa=0.0d0,epsr=1.0d-4 ! arguments integer :: lrm,fast real(dp) :: yg,mu,npl,cr,ci real(dp), dimension(-lrm:lrm,0:2,0:lrm) :: rr ! internal - integer :: n,m,ih,i,k,n1,nn,llm,neval,ier,last + integer :: n,m,ih,i,k,n1,nn,llm,neval,ier,last,ihmin integer, dimension(liw) :: iw real(dp) :: mu2,mu4,mu6,npl2,dt,bth,bth2,bth4,bth6 real(dp) :: t,x,s,sy1,sy2,sy3,resfh,epp @@ -620,36 +618,29 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast) mu4=mu2*mu2 mu6=mu4*mu2 ! - n1=1 - if(fast.gt.10) n1=-llm + 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) + apar(1) = yg + apar(2) = mu + apar(3) = npl + apar(4) = cr ! - 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 + do n=n1,llm + nn=abs(n) + do m=nn,llm + apar(5) = real(n,dp) + apar(6) = real(m,dp) + if(n.eq.0.and.m.eq.0) ihmin=2 + do ih=ihmin,2 + apar(7) = real(ih,dp) +! call dqagi(fhermit,bound,2,epsa,epsr,resfh, + call dqagsmv(fhermit,-tmax,tmax,apar,npar,epsa,epsr,resfh, & + epp,neval,ier,liw,lw,last,iw,w) + rr(n,ih,m) = resfh end do end do + end do if(fast.gt.10) return ! @@ -662,88 +653,88 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast) ! 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,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)) + 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)) + 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 + 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(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(-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)) + rr(-3,0,3) = -12.0d0*bth4/sy3* & + (1.0d0+bth2*(3.0d0-2.25d0/sy3+0.5d0*npl2/sy3**2)+ & + bth4*(3.75d0-8.71875d0/sy3+(6.1875d0+2.625d0*npl2) & + /sy3**2-4.125d0*npl2/sy3**3+0.75d0*npl2*npl2/sy3**4)) + rr(-3,1,3) = -6.0d0*npl*bth6/sy3**2* & + (1.0d0+bth2*(5.25d0-5.5d0/sy3+1.5d0*npl2/sy3**2)) + rr(-3,2,3) = -6.0d0*bth6/sy3* & + (1.0d0+bth2*(5.25d0-2.75d0/sy3+1.5d0*npl2/sy3**2)) ! - end if + end if end if ! - return - end + end subroutine hermitian_2 ! ! function fhermit(t,apar,npar) ! implicit none ! arguments - integer npar - real(dp) t,fhermit - real(dp), dimension(npar) :: apar + integer, intent(in) :: npar + real(dp), intent(in) :: t + real(dp) :: fhermit + real(dp), dimension(npar), intent(in) :: apar ! internal integer :: n,m,ih integer nn @@ -771,7 +762,6 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast) 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) @@ -787,8 +777,7 @@ subroutine diel_tens_fr(xg,yg,mu,npl,e330,cr,ci,epsl,lrm,fast) 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 + end function fhermit ! ! subroutine antihermitian(ri,yg,mu,npl,cr,ci,lrm)