! ! numerical constants ! ! ! last revised: 15 January 2011 ! module numconstants implicit none integer :: print_intermediate_results integer, allocatable :: monen(:) integer, private :: nmax=0 real(8) :: pi real(8), allocatable :: bcof(:,:),fnr(:),vwh_coef(:,:,:,:) real(8), allocatable :: vcc_const(:,:,:),fnm1_const(:,:),fn_const(:,:),fnp1_const(:,:) data pi/3.141592653589793/ contains subroutine init(notd) implicit none integer :: notd,l,n,ierr,nbc,m,mm1,mp1,np1,nm1,nn1,mn real(8) :: fnorm1,fnorm2 ! ! bcof(n,l)=((n+l)!/(n!l!))^(1/2) ! if(notd.le.nmax) return nmax=max(nmax,notd) nbc=6*notd+6 if(allocated(fnr)) deallocate(monen,fnr,bcof) allocate (monen(0:2*notd),bcof(0:nbc,0:nbc),fnr(0:2*nbc),stat=ierr) ! write(*,'('' nmax, bcof status:'',2i5)') nmax,ierr do n=0,2*notd monen(n)=(-1)**n enddo fnr(0)=0.d0 do n=1,2*nbc fnr(n)=dsqrt(dble(n)) enddo bcof(0,0)=1.d0 do n=0,nbc-1 do l=n+1,nbc bcof(n,l)=fnr(n+l)*bcof(n,l-1)/fnr(l) bcof(l,n)=bcof(n,l) enddo bcof(n+1,n+1)=fnr(n+n+2)*fnr(n+n+1)*bcof(n,n)/fnr(n+1)/fnr(n+1) enddo if(allocated(vwh_coef)) deallocate(vwh_coef) allocate(vwh_coef(-notd:notd,1:notd,-1:1,-1:1)) ! ! constants used for calculation of svwf functions. ! do n=1,notd nn1=n*(n+1) np1=n+1 nm1=n-1 fnorm1=-.5d0/fnr(n+n+1)/fnr(n)/fnr(n+1) fnorm2=-.5d0*fnr(n+n+1)/fnr(n)/fnr(n+1) m=-n mp1=m+1 mm1=m-1 vwh_coef(m,n, 1, 1)=-fnorm1*n*fnr(np1+m)*fnr(np1+mp1) vwh_coef(m,n, 1,-1)=fnorm1*np1*fnr(n-m)*fnr(nm1-m) vwh_coef(m,n,-1, 1)=fnorm1*n*fnr(np1-m)*fnr(np1-mm1) vwh_coef(m,n,-1,-1)=0.d0 vwh_coef(m,n, 0, 1)=fnorm1*n*fnr(np1+m)*fnr(np1-m) vwh_coef(m,n, 0,-1)=0.d0 vwh_coef(m,n, 1, 0)=-fnorm2*fnr(n-m)*fnr(np1+m) vwh_coef(m,n,-1, 0)=-0.d0 vwh_coef(m,n, 0, 0)=-fnorm2*m do m=-n+1,-1 mp1=m+1 mm1=m-1 vwh_coef(m,n, 1, 1)=-fnorm1*n*fnr(np1+m)*fnr(np1+mp1) vwh_coef(m,n, 1,-1)=fnorm1*np1*fnr(n-m)*fnr(nm1-m) vwh_coef(m,n,-1, 1)=fnorm1*n*fnr(np1-m)*fnr(np1-mm1) vwh_coef(m,n,-1,-1)=-fnorm1*np1*fnr(n+m)*fnr(nm1+m) vwh_coef(m,n, 0, 1)=fnorm1*n*fnr(np1+m)*fnr(np1-m) vwh_coef(m,n, 0,-1)=fnorm1*np1*fnr(n+m)*fnr(n-m) vwh_coef(m,n, 1, 0)=-fnorm2*fnr(n-m)*fnr(np1+m) vwh_coef(m,n,-1, 0)=-fnorm2*fnr(n+m)*fnr(np1-m) vwh_coef(m,n, 0, 0)=-fnorm2*m enddo do m=0,n-1 mp1=m+1 mm1=m-1 vwh_coef(m,n, 1, 1)=-fnorm1*n*fnr(np1+m)*fnr(np1+mp1) vwh_coef(m,n, 1,-1)=fnorm1*np1*fnr(n-m)*fnr(nm1-m) vwh_coef(m,n,-1, 1)=fnorm1*n*fnr(np1-m)*fnr(np1-mm1) vwh_coef(m,n,-1,-1)=-fnorm1*np1*fnr(n+m)*fnr(nm1+m) vwh_coef(m,n, 0, 1)=fnorm1*n*fnr(np1+m)*fnr(np1-m) vwh_coef(m,n, 0,-1)=fnorm1*np1*fnr(n+m)*fnr(n-m) vwh_coef(m,n, 1, 0)=-fnorm2*fnr(n-m)*fnr(np1+m) vwh_coef(m,n,-1, 0)=-fnorm2*fnr(n+m)*fnr(np1-m) vwh_coef(m,n, 0, 0)=-fnorm2*m enddo m=n mp1=m+1 mm1=m-1 vwh_coef(m,n, 1, 1)=-fnorm1*n*fnr(np1+m)*fnr(np1+mp1) vwh_coef(m,n, 1,-1)=0.d0 vwh_coef(m,n,-1, 1)=fnorm1*n*fnr(np1-m)*fnr(np1-mm1) vwh_coef(m,n,-1,-1)=-fnorm1*np1*fnr(n+m)*fnr(nm1+m) vwh_coef(m,n, 0, 1)=fnorm1*n*fnr(np1+m)*fnr(np1-m) vwh_coef(m,n, 0,-1)=0.d0 vwh_coef(m,n, 1, 0)=-0.d0 vwh_coef(m,n,-1, 0)=-fnorm2*fnr(n+m)*fnr(np1-m) vwh_coef(m,n, 0, 0)=-fnorm2*m enddo end subroutine init end module numconstants ! ! special function for the multiple sphere problem ! module specialfuncs implicit none contains subroutine timewrite(iunit,char1,time) use intrinsics implicit none integer :: iunit real(8) :: time,time2 character(*) :: char1 if(time.gt.3600.d0) then time2=time/3600.d0 write(iunit,'(a,f9.3,'' hours'')') char1,time2 elseif(time.gt.60.d0) then time2=time/60.d0 write(iunit,'(a,f9.2,'' min'')') char1,time2 else write(iunit,'(a,f9.2,'' sec'')') char1,time endif call flush(iunit) end subroutine timewrite ! ! ricatti-bessel function psi(n), real argument ! subroutine ricbessel(n,ds,eps,nmax,psi) implicit none integer :: n,nmax,ns,i real(8) :: ds,dns,sn,psi(0:n),psit,ds2,sum,eps,err if(int(ds).lt.n) then ns=nint(ds+4.*(ds**.3333d0)+17) ns=max(n+10,ns) dns=0.d0 do i=ns-1,n,-1 sn=dble(i+1)/ds dns=sn-1.d0/(dns+sn) enddo psi(n)=dns psi(n-1)=dble(n)/ds-1.d0/(dns+dble(n)/ds) do i=n-2,1,-1 sn=dble(i+1)/ds psi(i)=sn-1.d0/(psi(i+1)+sn) enddo psit=dsin(ds) psi(0)=psit ds2=ds*ds sum=psit*psit/ds2 do i=1,n psit=psit/(dble(i)/ds+psi(i)) sum=sum+dble(i+i+1)*psit*psit/ds2 err=dabs(1.d0-sum) psi(i)=psit if(err.lt.eps) then nmax=i return endif enddo nmax=n else psi(0)=dsin(ds) psi(1)=psi(0)/ds-dcos(ds) do i=1,n-1 sn=dble(i+i+1)/ds psi(i+1)=sn*psi(i)-psi(i-1) enddo nmax=n endif end subroutine ricbessel ! ! ricatti-hankel function xi(n), real argument ! ! ! last revised: 15 January 2011 ! subroutine richankel(n,ds,xi) implicit none integer :: n,i,ns real(8) :: ds,dns,sn,chi0,chi1,chi2,psi,psi0,psi1 complex(8) :: xi(0:n) if(int(ds).lt.n) then ns=nint(ds+4.*(ds**.3333)+17) ns=max(n+10,ns) dns=0.d0 do i=ns-1,n,-1 sn=dble(i+1)/ds dns=sn-1.d0/(dns+sn) enddo xi(n)=dns xi(n-1)=dble(n)/ds-1.d0/(dns+dble(n)/ds) do i=n-2,1,-1 sn=dble(i+1)/ds xi(i)=sn-1.d0/(xi(i+1)+sn) enddo chi0=-dcos(ds) psi=dsin(ds) chi1=chi0/ds-psi xi(0)=dcmplx(psi,chi0) do i=1,n chi2=dble(i+i+1)/ds*chi1-chi0 psi=psi/(dble(i)/ds+xi(i)) xi(i)=dcmplx(psi,chi1) chi0=chi1 chi1=chi2 enddo return else chi0=-dcos(ds) psi0=dsin(ds) chi1=chi0/ds-psi0 psi1=psi0/ds+chi0 xi(0)=dcmplx(psi0,chi0) xi(1)=dcmplx(psi1,chi1) do i=1,n-1 sn=dble(i+i+1)/ds xi(i+1)=sn*xi(i)-xi(i-1) enddo return endif end subroutine richankel ! ! ricatti-bessel function psi(n), complex argument ! ! ! last revised: 15 January 2011 ! subroutine cricbessel(n,ds,psi) implicit none integer :: n,i complex(8) :: ds,psi(0:n),chi(0:n) call cspherebessel(n,ds,psi,chi) do i=0,n psi(i)=psi(i)*ds enddo return end subroutine cricbessel ! ! ricatti-hankel function psi(n), complex argument ! ! ! last revised: 15 January 2011 ! 7 october 2011: forces upwards recurrence for real argument ds ! subroutine crichankel(n,ds,xi) implicit none integer :: n,i,i1 complex(8) :: ds,psi(0:n),chi(0:n),xi(0:n),ci data ci/(0.d0,1.d0)/ xi(0)=-ci*cdexp(ci*ds) xi(1)=-cdexp(ci*ds)*(ci+ds)/ds if(dimag(ds).eq.0.d0) then do i=1,n-1 i1=i+1 xi(i1)=dble(i+i1)/ds*xi(i)-xi(i-1) enddo return endif if(cdabs(xi(0)).lt.1.d-10) then do i=1,n-1 i1=i+1 xi(i1)=dble(i+i1)/ds*xi(i)-xi(i-1) enddo return else call cspherebessel(n,ds,psi,chi) do i=1,n-1 i1=i+1 xi(i1)=(psi(i1)+ci*chi(i1))*ds enddo return endif end subroutine crichankel ! ! ========================================================== ! Purpose: Compute spherical Bessel functions jn(z) & yn(z) ! for a complex argument ! Input : z --- Complex argument ! n --- Order of jn(z) & yn(z) ( n = 0,1,2,... ) ! Output: CSJ(n) --- jn(z) ! CSY(n) --- yn(z) ! NM --- Highest order computed ! Routines called: ! MSTA1 and MSTA2 for computing the starting ! point for backward recurrence ! ========================================================== ! ! obtained from, and copywrited by, Jian-Ming Jin ! http://jin.ece.uiuc.edu/ ! ! ! last revised: 15 January 2011 ! subroutine cspherebessel(n,z,csj,csy) implicit none integer :: n,nm,k,m real(8) :: a0 complex(8) :: z,csj(0:n),csy(0:n),csa,csb,cs,cf0,cf1,cf a0=cdabs(z) nm=n if (a0.lt.1.0d-60) then csj=(0.d0,0.d0) csy=(-1.d300,0.d0) csy(0)=(1.d0,0.d0) return endif csj=(0.d0,0.d0) csj(0)=cdsin(z)/z csj(1)=(csj(0)-cdcos(z))/z if (n.ge.2) then csa=csj(0) csb=csj(1) m=msta1(a0,200) if (m.lt.n) then nm=m else m=msta2(a0,n,15) endif cf0=0.0d0 cf1=1.0d0-100 do k=m,0,-1 cf=(2.0d0*k+3.0d0)*cf1/z-cf0 if (k.le.nm) csj(k)=cf cf0=cf1 cf1=cf enddo if (cdabs(csa).gt.cdabs(csb)) cs=csa/cf if (cdabs(csa).le.cdabs(csb)) cs=csb/cf0 do k=0,min(nm,n) csj(k)=cs*csj(k) enddo endif csy=(1.d200,0.d0) csy(0)=-cdcos(z)/z csy(1)=(csy(0)-cdsin(z))/z do k=2,min(nm,n) if (cdabs(csj(k-1)).gt.cdabs(csj(k-2))) then csy(k)=(csj(k)*csy(k-1)-1.0d0/(z*z))/csj(k-1) else csy(k)=(csj(k)*csy(k-2)-(2.0d0*k-1.0d0)/z**3)/csj(k-2) endif enddo end subroutine cspherebessel ! ! =================================================== ! Purpose: Determine the starting point for backward ! recurrence such that the magnitude of ! Jn(x) at that point is about 10^(-MP) ! Input : x --- Argument of Jn(x) ! MP --- Value of magnitude ! Output: MSTA1 --- Starting point ! =================================================== ! ! ! last revised: 15 January 2011 ! integer function msta1(x,mp) implicit none integer :: mp,n0,n1,it,nn real(8) :: x, a0,f1,f,f0 a0=dabs(x) n0=int(1.1*a0)+1 f0=envj(n0,a0)-mp n1=n0+5 f1=envj(n1,a0)-mp do it=1,20 nn=n1-(n1-n0)/(1.0d0-f0/f1) f=envj(nn,a0)-mp if(abs(nn-n1).lt.1) exit n0=n1 f0=f1 n1=nn f1=f enddo msta1=nn end function msta1 ! ! =================================================== ! Purpose: Determine the starting point for backward ! recurrence such that all Jn(x) has MP ! significant digits ! Input : x --- Argument of Jn(x) ! n --- Order of Jn(x) ! MP --- Significant digit ! Output: MSTA2 --- Starting point ! =================================================== ! ! ! last revised: 15 January 2011 ! integer function msta2(x,n,mp) implicit none integer :: n,mp,n0,n1,it,nn real(8) :: x,a0,hmp,ejn,obj,f0,f1,f a0=dabs(x) hmp=0.5d0*dble(mp) ejn=envj(n,a0) if (ejn.le.hmp) then obj=mp n0=int(1.1*a0) else obj=hmp+ejn n0=n endif f0=envj(n0,a0)-obj n1=n0+5 f1=envj(n1,a0)-obj do it=1,20 nn=n1-(n1-n0)/(1.0d0-f0/f1) f=envj(nn,a0)-obj if (abs(nn-n1).lt.1) exit n0=n1 f0=f1 n1=nn f1=f enddo msta2=nn+10 end function msta2 real(8) function envj(n,x) implicit none integer :: n real(8) :: x n=max(1,abs(n)) envj=0.5d0*dlog10(6.28d0*n)-n*dlog10(1.36d0*x/n) end function envj ! ! vector coupling coefficients vc(w) = C(m,n|k,l|m+k,w), w = |n-l|,... n+l ! uses downwards and upwards recurrence ! ! ! last revised: 15 January 2011 ! subroutine vcfunc(m,n,k,l,vcn) use numconstants implicit none integer :: m,n,k,l,wmax,wmin,w,mk real(8) :: vcn(0:n+l),t1,t2,t3,vcmax,vctest,rat vcn=0.d0 wmax=n+l wmin=max(abs(n-l),abs(m+k)) vcn(wmax)=bcof(n+m,l+k)*bcof(n-m,l-k)/bcof(n+n,l+l) if(wmin.eq.wmax) return vcn(wmax-1)=vcn(wmax)*(l*m-k*n)*fnr(2*(l+n)-1)/fnr(l)/fnr(n)& & /fnr(n+l+m+k)/fnr(n+l-m-k) if(wmin.eq.wmax-1) return mk=m+k vcmax=abs(vcn(wmax))+abs(vcn(wmax-1)) ! ! a downwards recurrence is used initially ! do w=wmax,wmin+2,-1 t1=2*w*fnr(w+w+1)*fnr(w+w-1)/(fnr(w+mk)*fnr(w-mk)& & *fnr(n-l+w)*fnr(l-n+w)*fnr(n+l-w+1)*fnr(n+l+w+1)) t2=dble((m-k)*w*(w-1)-mk*n*(n+1)+mk*l*(l+1))& & /dble(2*w*(w-1)) t3=fnr(w-mk-1)*fnr(w+mk-1)*fnr(l-n+w-1)*fnr(n-l+w-1)& & *fnr(n+l-w+2)*fnr(n+l+w)/(dble(2*(w-1))*fnr(2*w-3)& & *fnr(2*w-1)) vcn(w-2)=(t2*vcn(w-1)-vcn(w)/t1)/t3 if(mod(wmax-w,2).eq.1) then vctest=abs(vcn(w-2))+abs(vcn(w-1)) vcmax=max(vcmax,vctest) rat=vctest/vcmax ! ! if/when the coefficients start to decrease in magnitude, an upwards recurrence takes over ! if(rat.lt.0.01d0) exit endif enddo if(w-2.gt.wmin) then wmax=w-3 call vcfuncuprec(m,n,k,l,wmax,vcn) endif end subroutine vcfunc ! ! upwards VC coefficient recurrence ! ! ! last revised: 15 January 2011 ! subroutine vcfuncuprec(m,n,k,l,wmax,vcn) use numconstants implicit none integer :: m,n,k,l,wmax,wmin,w,mk,nl,m1,n1,l1,k1,w1,w2 real(8) :: vcn(0:n+l),t1,t2,t3,vc1 mk=abs(m+k) nl=abs(n-l) if(nl.ge.mk) then w=nl if(n.ge.l) then m1=m n1=n l1=l k1=k else m1=k n1=l k1=m l1=n endif vc1=(-1)**(k1+l1)*bcof(l1+k1,w-m1-k1) & *bcof(l1-k1,w+m1+k1)/bcof(l1+l1,w+w+1) else w=mk if(m+k.ge.0) then vc1=(-1)**(n+m)*bcof(n-l+w,l-k)*bcof(l-n+w,n-m) & /bcof(w+w+1,n+l-w) else vc1=(-1)**(l+k)*bcof(n-l+w,l+k)*bcof(l-n+w,n+m) & /bcof(w+w+1,n+l-w) endif endif w1=w vcn(w)=vc1 w=w1+1 mk=m+k w2=min(wmax,n+l) if(w2.gt.w1) then t1=2*w*fnr(w+w+1)*fnr(w+w-1)/(fnr(w+mk)*fnr(w-mk) & *fnr(n-l+w)*fnr(l-n+w)*fnr(n+l-w+1)*fnr(n+l+w+1)) if(w1.eq.0) then t2=.5*dble(m-k) else t2=dble((m-k)*w*(w-1)-mk*n*(n+1)+mk*l*(l+1)) & /dble(2*w*(w-1)) endif vcn(w)=t1*t2*vcn(w1) endif do w=w1+2,w2 t1=2*w*fnr(w+w+1)*fnr(w+w-1)/(fnr(w+mk)*fnr(w-mk) & *fnr(n-l+w)*fnr(l-n+w)*fnr(n+l-w+1)*fnr(n+l+w+1)) t2=dble((m-k)*w*(w-1)-mk*n*(n+1)+mk*l*(l+1)) & /dble(2*w*(w-1)) t3=fnr(w-mk-1)*fnr(w+mk-1)*fnr(l-n+w-1)*fnr(n-l+w-1) & *fnr(n+l-w+2)*fnr(n+l+w)/(dble(2*(w-1))*fnr(2*w-3) & *fnr(2*w-1)) vcn(w)=t1*(t2*vcn(w-1)-t3*vcn(w-2)) enddo end subroutine vcfuncuprec ! ! Normalized associated legendre functions ! ! ! last revised: 15 January 2011 ! subroutine normalizedlegendre(cbe,mmax,nmax,dc) use numconstants implicit none integer :: nmax,mmax,m,n,np1,nm1,im real(8) :: dc(-mmax:mmax,0:nmax),cbe,sbe sbe=dsqrt((1.d0+cbe)*(1.d0-cbe)) dc=0.d0 do m=0,mmax dc(m,m)=(-1)**m*(0.5d0*sbe)**m*bcof(m,m) if(m.eq.nmax) exit dc(m,m+1)=fnr(m+m+1)*cbe*dc(m,m) do n=m+1,nmax-1 dc(m,n+1)=(-fnr(n-m)*fnr(n+m)*dc(m,n-1)+dble(n+n+1)*cbe*dc(m,n)) & /(fnr(n+1-m)*fnr(n+1+m)) enddo enddo do m=1,mmax im=(-1)**m do n=m,nmax dc(-m,n)=im*dc(m,n) enddo enddo end subroutine normalizedlegendre ! ! Generalized spherical functions ! ! dc(m,n*(n+1)+k)=(-1)^(m + k)((n - k)!(n + k)!/(n - m)!/(n + m)!)^(1/2) ! ((1 + x)/2)^((m + k)/2)((1 - x)/2)^((k - m)/2)JacobiP[n - k, k - m, k + m, x] ! ! for |m| <= kmax, n=0,1,...nmax, |k| <= n ! ! ! last revised: 15 January 2011 ! subroutine rotcoef(cbe,kmax,nmax,dc) use numconstants implicit none integer :: kmax,nmax,k,m,in,n,knmax,nn1,kn,im,m1 real(8) :: cbe,sbe,dc(-kmax:kmax,0:nmax*(nmax+2)),cbe2,sbe2,dk0(-nmax-1:nmax+1),& dk01(-nmax-1:nmax+1),sben,dkt,fmn,dkm0,dkm1,dkn1 sbe=dsqrt((1.d0+cbe)*(1.d0-cbe)) cbe2=.5d0*(1.d0+cbe) sbe2=.5d0*(1.d0-cbe) in=1 dk0(0)=1.d0 sben=1.d0 dc(0,0)=1.d0 dk01(0)=0. do n=1,nmax knmax=min(n,kmax) nn1=n*(n+1) in=-in sben=sben*sbe/2.d0 dk0(n)=in*sben*bcof(n,n) dk0(-n)=in*dk0(n) dk01(n)=0. dk01(-n)=0. dc(0,nn1+n)=dk0(n) dc(0,nn1-n)=dk0(-n) do k=-n+1,n-1 kn=nn1+k dkt=dk01(k) dk01(k)=dk0(k) dk0(k)=(cbe*dble(n+n-1)*dk01(k)-fnr(n-k-1)*fnr(n+k-1)*dkt)& /(fnr(n+k)*fnr(n-k)) dc(0,kn)=dk0(k) enddo im=1 do m=1,knmax im=-im fmn=1.d0/fnr(n-m+1)/fnr(n+m) m1=m-1 dkm0=0. do k=-n,n kn=nn1+k dkm1=dkm0 dkm0=dc(m1,kn) if(k.eq.n) then dkn1=0. else dkn1=dc(m1,kn+1) endif dc(m,kn)=(fnr(n+k)*fnr(n-k+1)*cbe2*dkm1 & -fnr(n-k)*fnr(n+k+1)*sbe2*dkn1 & -dble(k)*sbe*dc(m1,kn))*fmn dc(-m,nn1-k)=dc(m,kn)*(-1)**(k)*im enddo enddo enddo end subroutine rotcoef subroutine rotcoefvecarg(narg,cbe,kmax,nmax,dc) use numconstants implicit none integer :: kmax,nmax,k,m,in,n,knmax,nn1,kn,im,m1,narg real(8) :: cbe(narg),sbe(narg),dc(-kmax:kmax,0:nmax*(nmax+2),narg), & cbe2(narg),sbe2(narg),dk0(-nmax-1:nmax+1,narg),& dk01(-nmax-1:nmax+1,narg),sben(narg),dkt(narg), & fmn,dkm0(narg),dkm1(narg),dkn1(narg) sbe=sqrt((1.d0+cbe)*(1.d0-cbe)) cbe2=.5d0*(1.d0+cbe) sbe2=.5d0*(1.d0-cbe) in=1 dk0(0,:)=1.d0 sben=1.d0 dc(0,0,:)=1.d0 dk01(0,:)=0. do n=1,nmax knmax=min(n,kmax) nn1=n*(n+1) in=-in sben=sben*sbe/2.d0 dk0(n,:)=in*sben(:)*bcof(n,n) dk0(-n,:)=in*dk0(n,:) dk01(n,:)=0. dk01(-n,:)=0. dc(0,nn1+n,:)=dk0(n,:) dc(0,nn1-n,:)=dk0(-n,:) do k=-n+1,n-1 kn=nn1+k dkt(:)=dk01(k,:) dk01(k,:)=dk0(k,:) dk0(k,:)=(cbe(:)*dble(n+n-1)*dk01(k,:)-fnr(n-k-1)*fnr(n+k-1)*dkt(:)) & /(fnr(n+k)*fnr(n-k)) dc(0,kn,:)=dk0(k,:) enddo im=1 do m=1,knmax im=-im fmn=1.d0/fnr(n-m+1)/fnr(n+m) m1=m-1 dkm0=0. do k=-n,n kn=nn1+k dkm1=dkm0 dkm0(:)=dc(m1,kn,:) if(k.eq.n) then dkn1=0. else dkn1(:)=dc(m1,kn+1,:) endif dc(m,kn,:)=(fnr(n+k)*fnr(n-k+1)*cbe2(:)*dkm1(:) & -fnr(n-k)*fnr(n+k+1)*sbe2(:)*dkn1(:) & -dble(k)*sbe(:)*dc(m1,kn,:))*fmn dc(-m,nn1-k,:)=dc(m,kn,:)*(-1)**(k)*im enddo enddo enddo end subroutine rotcoefvecarg ! ! tau are the vector spherical harmonic functions, normalized ! ! ! last revised: 15 January 2011 ! subroutine taufunc(cb,nmax,tau) use numconstants implicit none integer :: nmax,n,m,p,nn1,mn real(8) :: drot(-1:1,0:nmax*(nmax+2)),tau(0:nmax+1,nmax,2),cb,fnm call rotcoef(cb,1,nmax,drot) do n=1,nmax nn1=n*(n+1) fnm=sqrt(dble(n+n+1)/2.d0)/4.d0 do m=-n,-1 mn=nn1+m tau(n+1,-m,1)=-fnm*(-drot(-1,mn)+drot(1,mn)) tau(n+1,-m,2)=-fnm*(drot(-1,mn)+drot(1,mn)) enddo do m=0,n mn=nn1+m tau(m,n,1)=-fnm*(-drot(-1,mn)+drot(1,mn)) tau(m,n,2)=-fnm*(drot(-1,mn)+drot(1,mn)) enddo enddo end subroutine taufunc ! ! vector spherical harmonic function ! november 2011 ! subroutine pifunc(cb,ephi,nmax,ndim,pivec) use numconstants implicit none integer :: nmax,n,m,p,nn1,mn,ndim real(8) :: drot(-1:1,0:nmax*(nmax+2)),tau(2),cb,fnm complex(8) :: pivec(0:ndim+1,ndim,2),ephi,ephim(-nmax:nmax),cin call rotcoef(cb,1,nmax,drot) ephim(0)=1.d0 do m=1,nmax ephim(m)=ephi*ephim(m-1) ephim(-m)=dconjg(ephim(m)) enddo do n=1,nmax cin=(0.d0,-1.d0)**(n+1) nn1=n*(n+1) fnm=sqrt(dble(n+n+1)/2.d0)/4.d0 do m=-n,-1 mn=nn1+m tau(1)=-fnm*(-drot(-1,mn)+drot(1,mn)) tau(2)=-fnm*(drot(-1,mn)+drot(1,mn)) pivec(n+1,-m,1)=cin*tau(1)*ephim(m) pivec(n+1,-m,2)=cin*tau(2)*ephim(m) enddo do m=0,n mn=nn1+m tau(1)=-fnm*(-drot(-1,mn)+drot(1,mn)) tau(2)=-fnm*(drot(-1,mn)+drot(1,mn)) pivec(m,n,1)=cin*tau(1)*ephim(m) pivec(m,n,2)=cin*tau(2)*ephim(m) enddo enddo end subroutine pifunc ! ! regular vswf expansion coefficients for a plane wave. ! alpha, beta: incident azimuth and polar angles. ! ! ! last revised: 15 January 2011 ! subroutine planewavecoef(alpha,beta,nodr,pmnp0) use numconstants implicit none integer :: nodr,m,n,p,k,ierr real(8) :: alpha,beta,cb,sb,ca,sa real(8), allocatable :: tau(:,:,:) complex(8) :: ealpha,ci,cin complex(8), allocatable :: ealpham(:) complex(8) :: pmnp0(0:nodr+1,nodr,2,2) data ci/(0.d0,1.d0)/ call init(nodr) allocate(ealpham(-nodr:nodr)) allocate(tau(0:nodr+1,nodr,2)) cb=cos(beta) sb=sqrt((1.d0-cb)*(1.d0+cb)) ca=cos(alpha) sa=sin(alpha) ealpha=dcmplx(ca,sa) call taufunc(cb,nodr,tau) call ephicoef(ealpha,nodr,ealpham) do n=1,nodr cin=4.d0*ci**(n+1) do p=1,2 do m=-n,-1 pmnp0(n+1,-m,p,1)=-cin*tau(n+1,-m,p)*ealpham(-m) pmnp0(n+1,-m,p,2)=ci*cin*tau(n+1,-m,3-p)*ealpham(-m) enddo do m=0,n pmnp0(m,n,p,1)=-cin*tau(m,n,p)*ealpham(-m) pmnp0(m,n,p,2)=ci*cin*tau(m,n,3-p)*ealpham(-m) enddo enddo enddo deallocate(ealpham,tau) end subroutine planewavecoef ! ! regular vswf expansion coefficients for a gaussian beam, localized approximation. ! cbeam = 1/(k omega) ! ! ! last revised: 15 January 2011 ! subroutine gaussianbeamcoef(alpha,beta,cbeam,nodr,pmnp0) use numconstants implicit none integer :: nodr,m,n,p,k,ierr real(8) :: alpha,beta,cbeam,gbn complex(8) :: pmnp0(0:nodr+1,nodr,2,2) call planewavecoef(alpha,beta,nodr,pmnp0) do n=1,nodr gbn=dexp(-((dble(n)+.5d0)*cbeam)**2.) do p=1,2 do k=1,2 do m=-n,-1 pmnp0(n+1,-m,p,k)=pmnp0(n+1,-m,p,k)*gbn enddo do m=0,n pmnp0(m,n,p,k)=pmnp0(m,n,p,k)*gbn enddo enddo enddo enddo end subroutine gaussianbeamcoef ! ! plane wave expansion coefficients at sphere origins. uses a phase shift. ! ! ! last revised: 15 January 2011 ! subroutine sphereplanewavecoef(nsphere,neqns,nodr,nodrmax,alpha,beta,rpos,pmnp) implicit none integer :: m,n,p,nsphere,i,l,nodr(nsphere),nblk,nboff,nodrmax,neqns,k real(8) :: alpha,beta,cb,sb,ca,sa,rpos(3,nsphere) complex(8) :: ci,phasefac, pmnp(neqns,2) complex(8) :: pmnp0(0:nodrmax+1,nodrmax,2,2) data ci/(0.d0,1.d0)/ call planewavecoef(alpha,beta,nodrmax,pmnp0) cb=cos(beta) sb=sqrt((1.d0-cb)*(1.d0+cb)) ca=cos(alpha) sa=sin(alpha) l=0 do i=1,nsphere phasefac=cdexp(ci*((ca*rpos(1,i)+sa*rpos(2,i))*sb+rpos(3,i)*cb)) do p=1,2 do n=1,nodr(i) do m=0,nodr(i)+1 l=l+1 do k=1,2 pmnp(l,k)=phasefac*pmnp0(m,n,p,k) enddo enddo enddo enddo enddo end subroutine sphereplanewavecoef ! ! this computes the normalized translation coefficients for an ! axial translation of positive distance r. For itype=1 or 3, the translation ! uses the spherical Bessel or Hankel functions as a basis function, ! respectively. They are related to the coefficients appearing in ! M&M JOSA 96 by ! ! J^{ij}_{mnp mlq} = (E_{ml}/E_{mn})^(1/2) ac(s,n,l*(l+1)+m) ! ! where ! ! E_{mn} = n(n+1)(n+m)!/((2n+1)(n-m)!) ! s=mod(p+q,2)+1 (i.e., s=1 for the A coefficient, =2 for the B ! coefficient) ! ! The calculation procedure is based on the derivation ! of the addition theorem for vector harmonics, appearing in ! Fuller and Mackowski, proc. Light Scattering by Nonspherical ! Particles, NASA/GISS Sept. 1998. ! ! revised: 10 october 2011: used F90 vector arithmetic and precalculation ! of various constants. ! subroutine axialtrancoef(itype,r,ri,nmax,lmax,ac) use numconstants implicit none integer :: itype,nmax,lmax,n,l,m,p,w,n21,ll1,nlmin,lblk,wmin,wmax,ml integer :: iadd,nlmax integer, save :: nlmax0 real(8) :: r complex(8) :: ri,ci,z,xi(0:nmax+lmax) complex(8) :: ac(nmax,lmax*(lmax+3)/2,2) data ci,nlmax0/(0.d0,1.d0),0/ nlmax=max(nmax,lmax) if(nlmax.gt.nlmax0) then nlmax0=nlmax call axialtrancoefinit(nlmax) endif if(r.eq.0.d0) then ac=(0.d0,0.d0) if(itype.ne.1) return do m=0,min(nmax,lmax) do n=max(1,m),min(nmax,lmax) iadd=atcadd(m,n,lmax) ac(n,iadd,l)=1. enddo enddo return endif z=r*ri if(itype.eq.1) then call cricbessel(nmax+lmax,z,xi) else call crichankel(nmax+lmax,z,xi) endif xi=xi/z do n=1,nmax do l=1,lmax wmin=abs(n-l) wmax=n+l do m=0,min(n,l) iadd=atcadd(m,l,lmax) ml=l*(l+1)/2+m ac(n,iadd,1)=sum(vcc_const(n,ml,wmin:wmax:2)*xi(wmin:wmax:2)) ac(n,iadd,2)=ci*sum(vcc_const(n,ml,wmin+1:wmax-1:2)*xi(wmin+1:wmax-1:2)) enddo enddo enddo end subroutine axialtrancoef ! ! axial translation coefficients calculated by the diamond recurrence formula ! new: 10 october 2011 ! subroutine axialtrancoefrecurrence(itype,r,ri,nmax,lmax,ac) use numconstants implicit none integer :: itype,nmax,lmax,n,l,m,p,q,w,n21,ll1,nlmin,lblk, & wmin,wmax,ml,m1,np1,nm1,iaddp1,iaddm1,lm1,lp1 integer :: iadd,nlmax,iadd0,iadd1 integer, save :: nlmax0 real(8) :: r,fnp1,fn,fnm1,flp1,fl,flm1 complex(8) :: ri,ci,z,xi(0:nmax+lmax) complex(8) :: ac(nmax,lmax*(lmax+3)/2,2) data ci,nlmax0/(0.d0,1.d0),0/ nlmax=max(nmax,lmax) nlmin=min(nmax,lmax) if(nlmax.gt.nlmax0) then nlmax0=nlmax call axialtrancoefinit(nlmax) endif if(r.eq.0.d0) then ac=(0.d0,0.d0) if(itype.ne.1) return do m=0,nlmin m1=max(1,m) do n=m1,nlmin iadd=atcadd(m,n,lmax) ac(n,iadd,l)=1. enddo enddo return endif z=r*ri if(itype.eq.1) then call cricbessel(nmax+lmax,z,xi) else call crichankel(nmax+lmax,z,xi) endif xi=xi/z lm1=lmax-1 do m=0,nlmin m1=max(1,abs(m)) lp1=m1+1 iadd0=atcadd(m,m1,lmax) iadd1=atcadd(m,lmax,lmax) iaddp1=iadd0+1 iaddm1=iadd1-1 iadd=iadd0-1 n=m1 do l=m1,lmax wmin=abs(n-l) wmax=n+l iadd=iadd+1 ml=l*(l+1)/2+m ac(n,iadd,1)=sum(vcc_const(n,ml,wmin:wmax:2)*xi(wmin:wmax:2)) ac(n,iadd,2)=ci*sum(vcc_const(n,ml,wmin+1:wmax-1:2)*xi(wmin+1:wmax-1:2)) enddo l=lmax iadd=iadd1 ml=l*(l+1)/2+m do n=m1+1,nmax wmin=abs(n-l) wmax=n+l ac(n,iadd,1)=sum(vcc_const(n,ml,wmin:wmax:2)*xi(wmin:wmax:2)) ac(n,iadd,2)=ci*sum(vcc_const(n,ml,wmin+1:wmax-1:2)*xi(wmin+1:wmax-1:2)) enddo if(m1.eq.nlmin) cycle do n=m1,nmax-1 np1=n+1 nm1=n-1 do p=1,2 q=3-p ac(np1,iadd0:iaddm1,p)= & - ac(n,iaddp1:iadd1,p)*fnp1_const(m,m1:lm1) & + (fn_const(m,m1:lm1)-fn_const(m,n))*ci*ac(n,iadd0:iaddm1,q) ac(np1,iaddp1:iaddm1,p)=ac(np1,iaddp1:iaddm1,p) & + ac(n,iadd0:iadd1-2,p)*fnm1_const(m,lp1:lm1) if(n.gt.m1) then ac(np1,iadd0:iaddm1,p)=ac(np1,iadd0:iaddm1,p) & + ac(nm1,iadd0:iaddm1,p)*fnm1_const(m,n) endif ac(np1,iadd0:iaddm1,p)=ac(np1,iadd0:iaddm1,p)/fnp1_const(m,n) enddo enddo enddo end subroutine axialtrancoefrecurrence ! ! constants for translation coefficient calculation ! subroutine axialtrancoefinit(nmax) use numconstants implicit none integer :: nmax,m,n,l,w,n21,ml,ll1,wmin,wmax,nlmin,lp1,lm1 real(8) :: c1,c2,vc1(0:2*nmax),vc2(0:2*nmax),alnw complex(8) :: ci,inlw data ci/(0.d0,1.d0)/ if(allocated(vcc_const)) deallocate(vcc_const,fnm1_const,fn_const,fnp1_const) allocate(vcc_const(nmax,nmax*(nmax+1)/2+nmax,0:2*nmax),fnm1_const(0:nmax,nmax), & fn_const(0:nmax,nmax),fnp1_const(0:nmax,nmax)) do n=1,nmax n21=n+n+1 do l=1,nmax c1=fnr(n21)*fnr(l+l+1) ll1=l*(l+1)/2 call vcfunc(-1,n,1,l,vc2) wmin=abs(n-l) wmax=n+l nlmin=min(l,n) do m=0,nlmin ml=ll1+m c2=-c1*(-1)**m call vcfunc(-m,n,m,l,vc1) do w=wmin,wmax inlw=ci**(n-l+w) vcc_const(n,ml,w)=c2*vc1(w)*vc2(w)*(dble(inlw)+dimag(inlw)) enddo enddo enddo enddo fnm1_const=0. fn_const=0. fnp1_const=0. do m=0,nmax do l=max(1,m),nmax lp1=l+1 lm1=l-1 fnm1_const(m,l)=fnr(lm1)*fnr(lp1)*fnr(l-m)*fnr(l+m)/fnr(lm1+l)/fnr(l+lp1)/dble(l) fn_const(m,l)=dble(m)/dble(l)/dble(lp1) fnp1_const(m,l)=fnr(l)*fnr(l+2)*fnr(lp1-m)*fnr(lp1+m)/fnr(l+lp1)/fnr(l+l+3)/dble(lp1) enddo enddo end subroutine axialtrancoefinit ! ! test to determine convergence of regular vswf addition theorem for max. order lmax ! and translation distance r w/ refractive index ri. ! ! ! last revised: 15 January 2011 ! subroutine tranordertest(r,ri,lmax,eps,nmax) use numconstants implicit none integer :: itype,nmax,lmax,n,l,m,p,w,n21,ll1,nlmin,lblk,wmin,wmax integer, parameter :: nlim=200 integer :: iadd real(8) :: r,alnw,sum,eps real(8) :: vc1(0:nlim+lmax) complex(8) :: ri,ci,z,a,b,c complex(8) :: xi(0:nlim+lmax) data ci/(0.d0,1.d0)/ if(r.eq.0.d0) then nmax=lmax return endif z=r*ri sum=0.d0 do n=1,nlim call init(n+lmax) call cricbessel(n+lmax,z,xi) do l=0,n+lmax xi(l)=xi(l)/z*ci**l enddo n21=n+n+1 l=lmax c=fnr(n21)*fnr(l+l+1)*ci**(n-l) call vcfunc(-1,n,1,l,vc1) wmin=abs(n-l) wmax=n+l m=1 a=0. b=0. do w=wmin,wmax alnw=vc1(w)*vc1(w) if(mod(n+l+w,2).eq.0) then a=a+alnw*xi(w) else b=b+alnw*xi(w) endif enddo a=c*a b=c*b sum=sum+a*conjg(a)+b*conjg(b) if(abs(1.d0-sum).lt.eps) exit enddo nmax=min(n,nlim) nmax=max(nmax,lmax) end subroutine tranordertest ! ! address for axial translation coefficient ! ! ! last revised: 15 January 2011 ! integer function atcadd(m,n,ntot) implicit none integer :: m,n,ntot atcadd=n-ntot+(max(1,m)*(1+2*ntot-max(1,m)))/2+ntot*min(1,m) end function atcadd ! ! gentrancoef: calculates the vwh translation coefficients for ! a general translation from one origin to another ! ! input: itype: integer, =1, regular, =3, outgoing type harmonics ! xptran: real, dim 3 vector: x,y,z components of translation, in units ! of 1/k ! ri: complex, refractive index of medium ! nrow0,nrow1,ncol0,ncol1: integer, starting and stopping row and column order ! iaddrow0,iaddcol0: address offset for row and column order (see below) ! output: ac(p,mn,kl): complex translation matrix. calculated for mode p=1,2 (A or B type), ! order n=nrow0,nrow1, degree m=-n,n ! order l=ncol0,ncol1, degree k=-n,n ! address is given by ! mn=m+n*(n+1)-(nrow0-1)*(nrow0+1)+iaddrow0 ! kl=k+l*(l+1)-(ncol0-1)*(ncol0+1)+iaddcol0 ! that is, if iaddrow0=0 the address is mn=1 for n=nrow0 and m=-n. ! ! ! last revised: 15 January 2011 ! subroutine gentrancoef(itype,xptran,ri,nrow0,nrow1,ncol0,ncol1, & iaddrow0,iaddcol0,ac) use numconstants implicit none integer :: itype,nrow0,nrow1,ncol0,ncol1,iaddrow0,iaddcol0,kmax integer :: ntot,nblkr0,nblkr1,nblkc0,nblkc1 integer :: v,vw,w,wmax,wmin,n,l,m,k,p,nn1,ll1,mn,kl,m1m real(8) :: vc1(0:nrow1+ncol1),vc2(0:nrow1+ncol1),& xptran(3),r,ct,ct0 real(8) :: drot(0:0,0:(nrow1+ncol1)*(nrow1+ncol1+2)) complex(8) :: ri,ci,ephi,ac(2,nrow1*(nrow1+2)-(nrow0-1)*(nrow0+1)-iaddrow0,& ncol1*(ncol1+2)-(ncol0-1)*(ncol0+1)-iaddcol0),& z,c,a,b complex(8) :: ephim(-(nrow1+ncol1):nrow1+ncol1),jnc(0:nrow1+ncol1) data ci/(0.d0,1.d0)/ call cartosphere(xptran,r,ct,ephi) ntot=nrow1+ncol1 nblkr0=(nrow0-1)*(nrow0+1) nblkr1=nrow1*(nrow1+2) nblkc0=(ncol0-1)*(ncol0+1) nblkc1=ncol1*(ncol1+2) if(r.eq.0.d0) then do n=nblkr0+1,nblkr1 mn=n-nblkr0+iaddrow0 do l=nblkc0+1,nblkc1 kl=l-nblkc0+iaddcol0 do p=1,2 ac(p,mn,kl)=0.d0 enddo enddo if(n.gt.nblkc0.and.n.le.nblkc1.and.itype.eq.1) then ac(1,mn,n-nblkc0+iaddcol0)=1.d0 endif enddo return endif kmax=0 ct0=ct call rotcoef(ct0,kmax,ntot,drot) call ephicoef(ephi,ntot,ephim) z=ri*r if(itype.eq.1) then call cricbessel(ntot,z,jnc) else call crichankel(ntot,z,jnc) endif do n=0,ntot c=ci**n jnc(n)=c*jnc(n)/z enddo do l=ncol0,ncol1 ll1=l*(l+1) do n=nrow0,nrow1 nn1=n*(n+1) wmax=n+l call vcfunc(-1,n,1,l,vc2) c=-ci**(n-l)*fnr(n+n+1)*fnr(l+l+1) do k=-l,l kl=ll1+k-nblkc0+iaddcol0 do m=-n,n m1m=(-1)**m mn=nn1+m-nblkr0+iaddrow0 v=k-m call vcfunc(-m,n,k,l,vc1) a=0. b=0. wmin=max(abs(v),abs(n-l)) do w=wmax,wmin,-1 vw=w*(w+1)+v if(mod(wmax-w,2).eq.0) then a=a+vc1(w)*vc2(w)*jnc(w)*drot(0,vw) else b=b+vc1(w)*vc2(w)*jnc(w)*drot(0,vw) endif enddo ac(1,mn,kl)=a*c*m1m*ephim(v) ac(2,mn,kl)=b*c*m1m*ephim(v) enddo enddo enddo enddo return end subroutine gentrancoef ! ! cartosphere takes the cartesian point (x,y,z) = xp(1), xp(2), xp(3) ! and converts to polar form: r: radius, ct: cos(theta), ep = exp(i phi) ! ! ! last revised: 15 January 2011 ! subroutine cartosphere(xp,r,ct,ep) implicit none real(8) :: xp(3),r,ct complex(8) :: ep r=xp(1)*xp(1)+xp(2)*xp(2)+xp(3)*xp(3) if(r.eq.0.d0) then ct=1.d0 ep=(1.d0,0.d0) return endif r=sqrt(r) ct=xp(3)/r if(xp(1).eq.0.d0.and.xp(2).eq.0.d0) then ep=(1.d0,0.d0) else ep=dcmplx(xp(1),xp(2))/sqrt(xp(1)*xp(1)+xp(2)*xp(2)) endif return end subroutine cartosphere ! ! ephicoef returns the complex array epm(m) = exp(i m phi) for ! m=-nodr,nodr. ep =exp(i phi), and epm is dimensioned epm(-nd:nd) ! ! ! last revised: 15 January 2011 ! subroutine ephicoef(ep,nodr,epm) implicit none integer :: nodr,m complex(8) :: ep,epm(-nodr:nodr) epm(0)=(1.d0,0.d0) do m=1,nodr epm(m)=ep*epm(m-1) epm(-m)=dconjg(epm(m)) enddo return end subroutine ephicoef ! ! test to determine max order of vswf expansion of a plane wave at distance r ! ! ! last revised: 15 January 2011 ! subroutine planewavetruncationorder(r,eps,nodr) implicit none integer :: nodr,n1,n real(8) :: r,eps,err real(8), allocatable :: jn(:) complex(8) :: sum, ci,eir data ci/(0.d0,1.d0)/ n1=max(10,int(3.*r+1)) allocate(jn(0:n1)) call ricbessel(n1,r,-1.d0,n1,jn) jn(0:n1)=jn(0:n1)/r eir=cdexp(-ci*r) sum=jn(0)*eir do n=1,n1 sum=sum+ci**n*dble(n+n+1)*jn(n)*eir err=cdabs(1.d0-sum) if(err.lt.eps) then nodr=n deallocate(jn) return endif enddo nodr=n1 deallocate(jn) end subroutine planewavetruncationorder ! ! calculates the cartesian components of the vswf at position rpos, in ref. index ri. ! ! ! original: 15 January 2011 ! revised: 23 February 2011: multiplied by root 2 ! subroutine vwhcalc(rpos,ri,nodr,itype,vwh) use numconstants implicit none integer :: nodr,itype,m,n,p,nodrp1,nodrm1,nn1,mn,np1,nm1,mp1,mm1,ndim, & nblkp integer, save :: nodrmax real(8) :: rpos(3),r,ct,fnorm1,fnorm2 real(8) pmn(0:0,0:(nodr+1)*(nodr+3)) complex(8) :: ci,vwh(3,2,1:*),ri,ephi,a,b,a1,b1,z1,a2,b2,z2 complex(8) :: a1vec(-nodr:nodr), & b1vec(-nodr:nodr),z1vec(-nodr:nodr),a2vec(-nodr:nodr), & b2vec(-nodr:nodr),z2vec(-nodr:nodr) complex(8) :: umn(-nodr-2:nodr+2,0:nodr+1), hn(0:nodr+1), ephim(-nodr-1:nodr+1) data ci,nodrmax/(0.d0,1.d0),0/ if(nodr.gt.nodrmax) then nodrmax=nodr call init(nodr+2) endif call cartosphere(rpos,r,ct,ephi) if(r.le.1.d-4) then vwh(:,:,1:nodr*(nodr+1))=(0.d0,0.d0) if(itype.eq.3) return vwh(1,1,1)=.5d0*fnr(2)/fnr(3) vwh(2,1,1)=-.5d0*ci*fnr(2)/fnr(3) vwh(3,1,2)=1.d0*fnr(2)/fnr(6) vwh(1,1,3)=-.5d0*fnr(2)/fnr(3) vwh(2,1,3)=-.5d0*ci*fnr(2)/fnr(3) return endif nodrp1=nodr+1 nodrm1=nodr-1 a=ri*r if(itype.eq.1) then call cricbessel(nodrp1,a,hn) else call crichankel(nodrp1,a,hn) endif hn(0:nodrp1)=hn(0:nodrp1)/a call rotcoef(ct,0,nodrp1,pmn) call ephicoef(ephi,nodrp1,ephim) umn=0.d0 umn(0,0)=hn(0)*fnr(2) do n=1,nodrp1 nn1=n*(n+1) umn(-n:n,n)=fnr(2)*pmn(0,nn1-n:nn1+n)*ephim(-n:n)*hn(n) umn(-n-1,n)=0.d0 umn(n+1,n)=0.d0 enddo do n=1,nodr nn1=n*(n+1) np1=n+1 nm1=n-1 a1vec(-n:n)=vwh_coef(-n:n,n,1,1)*umn(-nm1:np1,np1) & +vwh_coef(-n:n,n,1,-1)*umn(-nm1:np1,nm1) b1vec(-n:n)=vwh_coef(-n:n,n,-1,1)*umn(-np1:nm1,np1) & +vwh_coef(-n:n,n,-1,-1)*umn(-np1:nm1,nm1) z1vec(-n:n)=vwh_coef(-n:n,n,0,1)*umn(-n:n,np1) & +vwh_coef(-n:n,n,0,-1)*umn(-n:n,nm1) a2vec(-n:n)=vwh_coef(-n:n,n,1,0)*umn(-nm1:np1,n) b2vec(-n:n)=vwh_coef(-n:n,n,-1,0)*umn(-np1:nm1,n) z2vec(-n:n)=vwh_coef(-n:n,n,0,0)*umn(-n:n,n) vwh(1,1,nn1-n:nn1+n)=-0.5d0*(a1vec(-n:n)+b1vec(-n:n)) vwh(2,1,nn1-n:nn1+n)=-0.5d0*ci*(-a1vec(-n:n)+b1vec(-n:n)) vwh(3,1,nn1-n:nn1+n)=-z1vec(-n:n) vwh(1,2,nn1-n:nn1+n)=-0.5d0*ci*(a2vec(-n:n)+b2vec(-n:n)) vwh(2,2,nn1-n:nn1+n)=-0.5d0*(a2vec(-n:n)-b2vec(-n:n)) vwh(3,2,nn1-n:nn1+n)=-ci*z2vec(-n:n) enddo end subroutine vwhcalc ! ! svwf calculation for an axial translation ! ! ! original: 15 January 2011 ! revised: 23 February 2011: multiplied by root 2 ! subroutine vwhaxialcalc(rpos,ri,nodr,itype,vwh) use numconstants implicit none integer :: nodr,itype,m,n,p,nodrp1,nodrm1,nn1,mn,np1,nm1,mp1,mm1,ndim, & nblkp integer, save :: nodrmax real(8) :: rpos(3),r,ct real(8) pmn(-2:2,0:nodr+1) complex(8) :: ci,vwh(3,2,2,1:nodr),ri,ephi,a,b,a1,b1,z1,a2,b2,z2 complex(8) :: umn(-2:2,0:nodr+1), hn(0:nodr+1), ephim(-2:2) data ci,nodrmax/(0.d0,1.d0),0/ if(nodr.gt.nodrmax) then nodrmax=nodr call init(nodr+2) endif call cartosphere(rpos,r,ct,ephi) if(r.le.1.d-4) then vwh(:,:,:,1:nodr)=(0.d0,0.d0) if(itype.eq.3) return vwh(1,1,1,1)=.5d0*fnr(2)/fnr(3) vwh(2,1,1,1)=-.5d0*ci*fnr(2)/fnr(3) vwh(1,1,2,1)=-.5d0*fnr(2)/fnr(3) vwh(2,1,2,1)=-.5d0*ci*fnr(2)/fnr(3) return endif nodrp1=nodr+1 nodrm1=nodr-1 a=ri*r if(itype.eq.1) then call cricbessel(nodrp1,a,hn) else call crichankel(nodrp1,a,hn) endif hn(0:nodrp1)=hn(0:nodrp1)/a call normalizedlegendre(ct,2,nodrp1,pmn) call ephicoef(ephi,2,ephim) umn(-2:2,0:nodrp1)=0.d0 umn(0,0)=hn(0)*fnr(2) do n=1,nodrp1 p=min(n,2) do m=-p,p umn(m,n)=fnr(2)*pmn(m,n)*ephim(m)*hn(n) enddo enddo vwh(:,:,:,1:nodr)=0.d0 do n=1,nodr np1=n+1 nm1=n-1 m=-1 mp1=m+1 mm1=m-1 a1=vwh_coef(m,n,1,1)*umn(mp1,np1) & +vwh_coef(m,n,1,-1)*umn(mp1,nm1) b1=vwh_coef(m,n,-1,1)*umn(mm1,np1) & +vwh_coef(m,n,-1,-1)*umn(mm1,nm1) z1=vwh_coef(m,n,0,1)*umn(m,np1) & +vwh_coef(m,n,0,-1)*umn(m,nm1) a2=vwh_coef(m,n,1,0)*umn(mp1,n) b2=vwh_coef(m,n,-1,0)*umn(mm1,n) z2=vwh_coef(m,n,0,0)*umn(m,n) vwh(1,1,1,n)=-0.5d0*(a1+b1) vwh(2,1,1,n)=-0.5d0*ci*(-a1+b1) vwh(3,1,1,n)=-z1 vwh(1,2,1,n)=-0.5d0*ci*(a2+b2) vwh(2,2,1,n)=-0.5d0*(a2-b2) vwh(3,2,1,n)=-ci*z2 m=1 mp1=m+1 mm1=m-1 a1=vwh_coef(m,n,1,1)*umn(mp1,np1) & +vwh_coef(m,n,1,-1)*umn(mp1,nm1) b1=vwh_coef(m,n,-1,1)*umn(mm1,np1) & +vwh_coef(m,n,-1,-1)*umn(mm1,nm1) z1=vwh_coef(m,n,0,1)*umn(m,np1) & +vwh_coef(m,n,0,-1)*umn(m,nm1) a2=vwh_coef(m,n,1,0)*umn(mp1,n) b2=vwh_coef(m,n,-1,0)*umn(mm1,n) z2=vwh_coef(m,n,0,0)*umn(m,n) vwh(1,1,2,n)=-0.5d0*(a1+b1) vwh(2,1,2,n)=-0.5d0*ci*(-a1+b1) vwh(3,1,2,n)=-z1 vwh(1,2,2,n)=-0.5d0*ci*(a2+b2) vwh(2,2,2,n)=-0.5d0*(a2-b2) vwh(3,2,2,n)=-ci*z2 enddo return end subroutine vwhaxialcalc end module specialfuncs ! ! module mpidata ! ! ! last revised: 15 January 2011 ! module mpidata implicit none integer :: group_comm,root_group_comm,base_rank,group_rank,root_group_rank, & base_group,number_groups,proc_per_group,number_proc integer, allocatable :: mpi_sphere_index(:), mpi_sphere_number(:) contains ! ! allocates the processors into groups ! ! last revised: 15 January 2011: original ! 20 April 2011: fixedorran=0 now looks for 2 groups. ! 10 october 2011: option for not storing matrices. If fixorran=0, 2 groups, else ! nproc groups ! november 2011: near and far field translation differentiation ! subroutine mpisetup(nsphere,nodr,rpos,fixorran,maxmbperproc,istore, & nfdistance,fftranpresent,iunit) use mpidefs use intrinsics use specialfuncs implicit none integer :: nsphere,numprocs,ierr,i,iunit,nodr(nsphere),fixorran, & nodrmax,nodrmin,temp_comm,newgroup,j,rank,maxmbperproc, & istore,nfspheres,fftranpresent,ffspheres integer, allocatable :: grouplist1(:),grouplist2(:) real(8) :: memrow(nsphere),memtot,maxmemproc,memperproc real(8) :: fp,sum,rpos(3,*),nfdistance,rij(3),r,avenfspheres,rmax, & nfdistancei,aveffspheres maxmemproc=maxmbperproc*1.d6 call ms_mpi(mpi_command='size',mpi_size=numprocs) call ms_mpi(mpi_command='rank',mpi_rank=rank) call ms_mpi(mpi_command='group',mpi_group=base_group) base_rank=rank number_proc=numprocs memrow=0.d0 memtot=0.d0 ! ! compute the memory storage requirements ! avenfspheres=0. aveffspheres=0. rmax=0. if(1.eq.1) then do i=1,nsphere nfspheres=0 do j=1,nsphere rij(:)=rpos(:,i)-rpos(:,j) if(j.ne.i) then if(nfdistance.lt.0.) then nfdistancei=(.5*dble(nodr(i)+nodr(j)))**2. else nfdistancei=nfdistance endif r=sqrt(dot_product(rij,rij)) rmax=max(rmax,r) if(r.le.nfdistancei) then nfspheres=nfspheres+1 nodrmax=max(nodr(j),nodr(i)) nodrmin=min(nodr(j),nodr(i)) memrow(i)=memrow(i)+(2*nodrmin+1)*(1+nodrmax*(nodrmax+2))*8.d0 memrow(i)=memrow(i)+nodr(i)*nodr(j)*(nodr(j)+3)*16.d0 memrow(i)=memrow(i)+(2*nodrmax+1)*16.d0 endif if(r.gt.nfdistancei.and.istore.eq.2) then memrow(i)=memrow(i)+2*nodrmax*(nodrmax+2)*16.d0 endif endif enddo ffspheres=nsphere-1-nfspheres avenfspheres=avenfspheres+nfspheres aveffspheres=aveffspheres+ffspheres memtot=memtot+memrow(i) enddo if(aveffspheres.eq.0) then fftranpresent=0 else fftranpresent=1 endif proc_per_group=ceiling(memtot/maxmemproc) proc_per_group=min(proc_per_group,numprocs) proc_per_group=max(proc_per_group,1) do if(mod(numprocs,proc_per_group).eq.0) exit if(proc_per_group.eq.numprocs) exit proc_per_group=proc_per_group+1 enddo endif avenfspheres=avenfspheres/dble(nsphere) if(rank.eq.0) then write(iunit,'('' average near field translations per sphere:'', f10.1)') avenfspheres call flush(iunit) endif ! ! no-store option ! if(istore.eq.0) then if(fixorran.eq.0) then proc_per_group=max(floor(dble(numprocs)/2.),1) else proc_per_group=1 endif memrow=1.d0 memtot=dble(nsphere) else ! ! only one or two groups for fixed orientation ! if(fixorran.eq.0) proc_per_group=max(floor(dble(numprocs)/2.),proc_per_group) endif number_groups=numprocs/proc_per_group if(allocated(mpi_sphere_index)) deallocate(mpi_sphere_index) if(allocated(mpi_sphere_number)) deallocate(mpi_sphere_number) allocate(mpi_sphere_index(0:proc_per_group-1),mpi_sphere_number(0:proc_per_group-1), & grouplist1(proc_per_group),grouplist2(number_groups)) memperproc=memtot/dble(proc_per_group) ! ! associate the spheres with the processors in a group ! mpi_sphere_index(0)=0 do j=1,proc_per_group-1 memtot=0.d0 do i=1,nsphere memtot=memtot+memrow(i) if(memtot.gt.dble(j)*memperproc) then mpi_sphere_index(j)=i-1 exit endif enddo enddo do i=0,proc_per_group-2 mpi_sphere_number(i)=mpi_sphere_index(i+1)-mpi_sphere_index(i) enddo mpi_sphere_number(proc_per_group-1)=nsphere-mpi_sphere_index(proc_per_group-1) ! ! assign the sphere-based groups ! do i=0,number_groups-1 do j=0,proc_per_group-1 grouplist1(j+1)=i*proc_per_group+j enddo call ms_mpi(mpi_command='incl',mpi_group=base_group,mpi_size=proc_per_group, & mpi_new_group_list=grouplist1,mpi_new_group=newgroup) call ms_mpi(mpi_command='create',mpi_group=newgroup,mpi_new_comm=temp_comm) if(rank.ge.grouplist1(1).and.rank.le.grouplist1(proc_per_group)) then group_comm=temp_comm endif grouplist2(i+1)=i*proc_per_group enddo ! ! make a group associated with the rank 0 members of the sphere groups ! call ms_mpi(mpi_command='incl',mpi_group=base_group,mpi_size=number_groups, & mpi_new_group_list=grouplist2,mpi_new_group=newgroup) call ms_mpi(mpi_command='create',mpi_group=newgroup,mpi_new_comm=temp_comm) group_rank=mod(rank,proc_per_group) root_group_rank=floor(dble(rank)/dble(proc_per_group)) if(group_rank.eq.0) root_group_comm=temp_comm if(rank.eq.0) then if(istore.ge.1) then write(iunit,'('' number of processors, number groups, mb mem/processor:'',2i5,f9.3)') & numprocs,number_groups,memperproc*1.d-6 if(memperproc.gt.maxmemproc) then write(iunit,'('' warning: set maximum memory/processor is exceeded!!'')') endif else write(iunit,'('' number of processors, number groups:'',2i5)') & numprocs,number_groups endif call flush(iunit) endif deallocate(grouplist1,grouplist2) end subroutine mpisetup end module mpidata ! ! module spheredata: used to 1) input sphere data, 2) dimension sphere data ! arrays, and 3) provide common access to the data in other subroutines. ! ! ! last revised: 15 January 2011 ! ! 30 March 2011: added optical activity ! module spheredata use specialfuncs use mpidata implicit none integer, private :: numberspheres,numberiterations,fixedorrandom,numbertheta, & calcnf,nfplane,calctmatrix,runprintunit,calcamn,maxmemperproc, & trackiterations,nfoutdata,normalizesm,storetranmat,niterstep, & fftranpresent real(8), private :: lengthscalefactor,realriscalefactor,imriscalefactor,epsmie, & epstran,epssoln,phideg,thetamindeg,thetamaxdeg,alphadeg, & betadeg,epstcon,nfplanepos,nfplanevert(2,2),deltax,gammadeg,epspw, & cgaussbeam,gaussbeamfocus(3),realchiralfactor,imchiralfactor,nfdistance character(30), private :: positionfile,outputfile,nfoutputfile,tmatrixfile,printfile, & amnfile real(8), private :: xspmax,xvsp real(8), private, allocatable :: rpos(:,:),xsp(:) complex(8), private, allocatable :: ri(:,:) data numberiterations,fixedorrandom,numbertheta/2000,0,181/ data calcamn,trackiterations,niterstep/1,1,20/ data lengthscalefactor,realriscalefactor,imriscalefactor,epsmie, & epstran,epssoln,phideg,thetamindeg,thetamaxdeg,alphadeg, & betadeg,epstcon/1.d0,1.d0,1.d0,1.d-4,1.d-6,1.d-10,0.d0,0.d0, & 180.d0,0.d0,0.d0,1.d-6/ data realchiralfactor,imchiralfactor/0.d0,0.d0/ data normalizesm,storetranmat,nfdistance/0,1,-1.d0/ data maxmemperproc/1500/ data cgaussbeam/0.d0/ data gaussbeamfocus/0.d0,0.d0,0.d0/ data calcnf,calctmatrix,nfoutdata/0,1,1/ data runprintunit/6/ data positionfile,outputfile,tmatrixfile,printfile/'at_bottom','test.dat','tmatrix-temp.dat',' '/ data nfoutputfile/'nf-temp.dat'/ data amnfile/'amn-temp.dat'/ contains ! ! Find the number of data points in input unit iunit, and reposition the unit to the ! point after record containing parmid ! ! ! last revised: 15 January 2011 ! subroutine numberinrecord(iunit,parmid,numrec) implicit none integer :: numrec,iunit character*1 :: a character*35 :: parmid character*10 :: rec numrec=0 do read(iunit,"(a)",advance="no",err=100,eor=100) a if(a.ne.' '.and.a.ne.',') then ! ! start of a number ! numrec=numrec+1 ! ! look for the delimeter ! do read(iunit,"(a)",advance="no",err=100,eor=100) a if(a.eq.' '.or.a.eq.',') exit enddo endif enddo 100 if(parmid.eq.'rewind') then rewind(iunit) else backspace(iunit) backspace(iunit) backspace(iunit) do read(iunit,'(a10)') rec if(rec.eq.parmid(1:10)) exit enddo endif end subroutine numberinrecord ! ! inputdata: reads parameters from inputfile ! reads sphere data from position file ! ! ! original: 15 January 2011 ! revised: 21 February 2011: fix output file initialization. ! 30 March 2011: added optical activity ! ! subroutine inputdata(inputfile,printdata) integer :: imax,i,j,ierr,iunit,numrec,nsphere,printdata real(8) :: rmax,rtoi,rposmean(3),rireal,riimag,dtemp,betareal,betaimag, & rij,xij(3),rijmax real(8), allocatable :: sdat(:) complex(8) :: ribulk,beta character*35 :: parmid character*30 :: inputfile ! ! cycle through parameter input operations ! open(1,file=inputfile) do read(1,'(a)',end=10) parmid parmid=parmid(:index(parmid,' ')) if(parmid.eq.'number_spheres') then read(1,*) numberspheres cycle endif if(parmid.eq.'sphere_position_file') then read(1,'(a)') positionfile positionfile=positionfile(:index(positionfile,' ')) cycle endif if(parmid.eq.'output_file') then read(1,'(a)') outputfile outputfile=outputfile(:index(outputfile,' ')) cycle endif if(parmid.eq.'run_print_file') then read(1,'(a)') printfile printfile=printfile(:index(printfile,' ')) if(printdata.eq.1) then if((printfile.eq.' '.or.printfile.eq.'console')) then printfile=' ' runprintunit=6 else runprintunit=4 open(runprintunit,file=printfile) endif else runprintunit=6 endif cycle endif if(parmid.eq.'length_scale_factor') then read(1,*) lengthscalefactor cycle endif if(parmid.eq.'real_ref_index_scale_factor') then read(1,*) realriscalefactor cycle endif if(parmid.eq.'imag_ref_index_scale_factor') then read(1,*) imriscalefactor cycle endif if(parmid.eq.'real_chiral_factor') then read(1,*) realchiralfactor cycle endif if(parmid.eq.'imag_chiral_factor') then read(1,*) imchiralfactor cycle endif if(parmid.eq.'mie_epsilon') then read(1,*) epsmie cycle endif if(parmid.eq.'translation_epsilon') then read(1,*) epstran cycle endif if(parmid.eq.'solution_epsilon') then read(1,*) epssoln cycle endif if(parmid.eq.'max_number_iterations') then read(1,*) numberiterations cycle endif if(parmid.eq.'max_memory_per_processor') then read(1,*) maxmemperproc cycle endif if(parmid.eq.'store_translation_matrix') then read(1,*) storetranmat cycle endif if(parmid.eq.'near_field_distance') then read(1,*) nfdistance cycle endif if(parmid.eq.'iterations_per_correction') then read(1,*) niterstep cycle endif if(parmid.eq.'fixed_or_random_orientation') then read(1,*) fixedorrandom cycle endif if(parmid.eq.'scattering_plane_angle_deg') then read(1,*) phideg cycle endif if(parmid.eq.'min_scattering_angle_deg') then read(1,*) thetamindeg cycle endif if(parmid.eq.'max_scattering_angle_deg') then read(1,*) thetamaxdeg cycle endif if(parmid.eq.'number_scattering_angles') then read(1,*) numbertheta cycle endif if(parmid.eq.'normalize_scattering_matrix') then read(1,*) normalizesm cycle endif if(parmid.eq.'incident_azimuth_angle_deg') then read(1,*) alphadeg cycle endif if(parmid.eq.'incident_polar_angle_deg') then read(1,*) betadeg cycle endif if(parmid.eq.'calculate_scattering_coefficients') then read(1,*) calcamn cycle endif if(parmid.eq.'scattering_coefficient_file') then read(1,'(a)') amnfile if(amnfile.eq.' ') then amnfile='amn-temp.dat' else amnfile=amnfile(:index(amnfile,' ')) endif cycle endif if(parmid.eq.'track_iterations') then read(1,*) trackiterations cycle endif if(parmid.eq.'calculate_near_field') then read(1,*) calcnf cycle endif if(parmid.eq.'near_field_plane_coord') then read(1,*) nfplane cycle endif if(parmid.eq.'near_field_plane_position') then read(1,*) nfplanepos cycle endif if(parmid.eq.'near_field_plane_vertices') then read(1,*) nfplanevert cycle endif if(parmid.eq.'spacial_step_size') then read(1,*) deltax cycle endif if(parmid.eq.'polarization_angle_deg') then read(1,*) gammadeg cycle endif if(parmid.eq.'near_field_output_file') then read(1,'(a)') nfoutputfile if(nfoutputfile.eq.' ') then nfoutputfile='nf-temp.dat' else nfoutputfile=nfoutputfile(:index(nfoutputfile,' ')) endif cycle endif if(parmid.eq.'near_field_output_data') then read(1,*) nfoutdata cycle endif if(parmid.eq.'plane_wave_epsilon') then read(1,*) epspw cycle endif if(parmid.eq.'gaussian_beam_constant') then read(1,*) cgaussbeam cycle endif if(parmid.eq.'gaussian_beam_focal_point') then read(1,*) gaussbeamfocus cycle endif if(parmid.eq.'t_matrix_convergence_epsilon') then read(1,*) epstcon cycle endif if(parmid.eq.'calculate_t_matrix') then read(1,*) calctmatrix cycle endif if(parmid.eq.'t_matrix_file') then read(1,'(a)') tmatrixfile if(tmatrixfile.eq.' ') then tmatrixfile='tmatrix-temp.dat' else tmatrixfile=tmatrixfile(:index(tmatrixfile,' ')) endif cycle endif if(parmid.eq.'sphere_sizes_and_positions') exit if(parmid.eq.'end_of_options') exit write(*,'('' warning: unknown parameter ID:'',a35)') parmid enddo ! ! end of parameter input options. Input of sphere data follows ! 10 write(runprintunit,'('' input file is '',a30)') inputfile if(positionfile.ne.'at_bottom'.and.positionfile.ne.' ') then close(1) open(1,file=positionfile) parmid='rewind' endif ! ! find number of records in position file ! call numberinrecord(1,parmid,numrec) if(printdata.eq.1) write(runprintunit,'('' position data has '',i3,'' records'')') numrec nsphere=numberspheres iunit=1 allocate(sdat(numrec)) allocate(xsp(0:nsphere),rpos(3,0:nsphere),ri(2,0:nsphere),stat=ierr) xvsp=0.d0 do i=1,nsphere read(iunit,*,end=20) sdat xsp(i)=sdat(1)*lengthscalefactor rpos(1:3,i)=sdat(2:4)*lengthscalefactor if(numrec.gt.4) then rireal=sdat(5)*realriscalefactor riimag=sdat(6)*imriscalefactor else rireal=realriscalefactor riimag=imriscalefactor endif if(numrec.gt.6) then betareal=sdat(7)*realchiralfactor betaimag=sdat(8)*imchiralfactor else betareal=realchiralfactor betaimag=imchiralfactor endif ribulk=dcmplx(rireal,riimag) beta=dcmplx(betareal,betaimag) if(beta.eq.(0.d0,0.d0)) then ri(1,i)=ribulk ri(2,i)=ribulk else ri(1,i)=ribulk/(1.d0-beta*ribulk) ri(2,i)=ribulk/(1.d0+beta*ribulk) endif xvsp=xvsp+xsp(i)**3.d0 enddo 20 nsphere=min(nsphere,i-1) close(iunit) deallocate(sdat) if(nsphere.ne.numberspheres.and.printdata.eq.1) then write(runprintunit,'('' warning: insufficient position points in file.'')') write(runprintunit,'('' number of spheres truncated to:'',i5)') nsphere endif ! ! check for overlapping spheres, and find maximum translation ! rijmax=0. do i=1,nsphere do j=i+1,nsphere xij=rpos(:,i)-rpos(:,j) rij=sqrt(dot_product(xij,xij)) rijmax=max(rijmax,rij) if(rij/(xsp(i)+xsp(j)).lt..999d0) then write(runprintunit,'('' warning: spheres '',i4,'' and '',i4 '' overlap. '',& & '' scaled distance:'' f8.4)') i,j,rij/(xsp(i)+xsp(j)) endif enddo enddo if(rijmax.gt.nfdistance) then fftranpresent=1 else fftranpresent=0 endif numberspheres=nsphere xvsp=xvsp**(1.d0/3.d0) gaussbeamfocus=gaussbeamfocus*lengthscalefactor if(nsphere.eq.1) then rposmean=rpos(:,1) rpos(:,1)=0.d0 xspmax=xsp(1) else rposmean=0.d0 do i=1,nsphere rposmean=rposmean+rpos(:,i) enddo rposmean=rposmean/dble(nsphere) rmax=0.d0 ! ! the target origin is defined as the GB focal point. ! do i=1,nsphere ! rpos(1:3,i)=rpos(1:3,i)-rposmean(1:3) rpos(1:3,i)=rpos(1:3,i)-gaussbeamfocus(1:3) rtoi=dot_product(rpos(:,i),rpos(:,i)) if(rtoi.gt.rmax) then rmax=rtoi imax=i endif enddo xspmax=sqrt(rmax)+xsp(imax) endif ! ! xsp(0) is the circumscribing sphere size parameter ! xsp(0)=xspmax ri(1,0)=(1.d0,0.d0) ri(2,0)=(1.d0,0.d0) rpos(:,0)=0.d0 ! ! write run data to run file and output file ! if(printdata.eq.1) then call writerundata(runprintunit) call flush(runprintunit) open(1,file=outputfile,status='replace',action='write') call writerundata(1) close(1) endif end subroutine inputdata ! ! writes run data to output unit iunit ! ! ! last revised: 15 January 2011 ! 30 March 2011: added optical activity ! subroutine writerundata(iunit) implicit none integer :: iunit,i character*1 :: lf if(iunit.ne.1) then lf = ' ' else lf = '/' endif write(iunit,'('' number of spheres, volume size parameter:'' '//lf//',i5,e13.5)') & numberspheres,xvsp write(iunit,'('' position file:'' '//lf//',a)') positionfile write(iunit,'('' output file:'' '//lf//',a)') outputfile write(iunit,'('' length, ref. indx. scale factors:'' '//lf//',3f8.3)') lengthscalefactor, & realriscalefactor,imriscalefactor write(iunit,'('' chiral factors:'' '//lf//',2e13.5)') & realchiralfactor,imchiralfactor write(iunit,'('' thetamin, thetamax, num. theta:'' '//lf//',2f9.1,i5)') & thetamindeg,thetamaxdeg,numbertheta write(iunit,'('' epsmie, epssoln, max number iterations:'' '//lf//',2e12.4,i5)') epsmie, & epssoln, numberiterations if(fftranpresent.eq.1) then write(iunit,'('' far field kr, iterations/correction:'' '//lf//',e12.4,i5)') & nfdistance,niterstep else write(iunit,'('' all translations computed exactly'' '//lf//')') endif if(cgaussbeam.ne.0.d0) then write(iunit,'('' gaussian incident beam: 1/width:'' '//lf//',f9.4,)') cgaussbeam write(iunit,'('' beam focal point:'' '//lf//',3f9.3,)') gaussbeamfocus else write(iunit,'('' plane wave incidence'')') endif if(fixedorrandom.eq.0) then write(iunit,'('' fixed orientation calculations'')') write(iunit,'('' scattering plane, incident alpha, beta:'' '//lf//',3f9.2)') & phideg,alphadeg,betadeg write(iunit,'('' common expansion epsilon:'' '//lf//',e12.4)') epstran if(calcamn.eq.0) then write(iunit,'('' scattering coefficients read from file '' '//lf//',a)') amnfile else write(iunit,'('' scattering coefficients calculated, stored in file '' '//lf//',a)') amnfile endif if(calcnf.eq.1) then write(iunit,'('' near field calculated, stored in file '' '//lf//',a)') nfoutputfile write(iunit,'('' near field data output option: '' '//lf//',i4)') nfoutdata write(iunit,'('' near field plane, position: '' '//lf//', i4,f9.3)') nfplane, nfplanepos write(iunit,'('' near field plane vertices: '' '//lf//',4f9.3)') nfplanevert write(iunit,'('' spacial step size:'' '//lf//',f9.4)') deltax write(iunit,'('' polarization angle, deg.:'' '//lf//',f9.2)') gammadeg write(iunit,'('' plane wave epsilon:'' '//lf//',e13.5)') epspw endif else write(iunit,'('' random orientation calculations'')') if(calctmatrix.eq.0) then write(iunit,'('' t matrix read from file '' '//lf//',a)') tmatrixfile elseif(calctmatrix.eq.1) then write(iunit,'('' t matrix calculated, stored in file '' '//lf//',a)') tmatrixfile write(iunit,'('' t matrix convergence epsilon:'' '//lf//',e12.4)') epstcon else write(iunit,'('' t matrix calculated from end of file '' '//lf//',a)') tmatrixfile write(iunit,'('' t matrix convergence epsilon:'' '//lf//',e12.4)') epstcon endif endif end subroutine writerundata ! ! getspheredata: retrieves sphere data ! ! ! last revised: 15 January 2011 ! 30 March 2011: added optical activity ! subroutine getspheredata(number_spheres, sphere_size_parameters, sphere_positions, & sphere_refractive_indices, volume_size_parameter) implicit none integer, optional :: number_spheres real(8), optional :: sphere_size_parameters(numberspheres), & sphere_positions(3,numberspheres), volume_size_parameter complex(8), optional :: sphere_refractive_indices(2,numberspheres) if (present(number_spheres)) number_spheres=numberspheres if (present(sphere_size_parameters)) sphere_size_parameters(1:numberspheres)=xsp(1:numberspheres) if (present(sphere_positions)) sphere_positions(:,1:numberspheres)=rpos(:,1:numberspheres) if (present(sphere_refractive_indices)) & sphere_refractive_indices(:,1:numberspheres)=ri(:,1:numberspheres) if (present(volume_size_parameter)) volume_size_parameter=xvsp end subroutine getspheredata subroutine getspheredataone(sphere,sphere_size_parameter, sphere_position, & sphere_refractive_index) implicit none integer :: sphere real(8), optional :: sphere_size_parameter,sphere_position(3) complex(8), optional :: sphere_refractive_index(2) if (present(sphere_size_parameter)) sphere_size_parameter=xsp(sphere) if (present(sphere_position)) sphere_position(:)=rpos(:,sphere) if (present(sphere_refractive_index)) & sphere_refractive_index(:)=ri(:,sphere) end subroutine getspheredataone ! ! setspheredata: sets sphere data ! subroutine setspheredata(number_spheres, sphere_size_parameters, sphere_positions, & sphere_refractive_indices, volume_size_parameter) implicit none integer :: i integer, optional :: number_spheres real(8), optional :: sphere_size_parameters(*), & sphere_positions(3,*), volume_size_parameter complex(8), optional :: sphere_refractive_indices(2,*) if (present(number_spheres)) then numberspheres=number_spheres if(allocated(xsp)) deallocate(xsp,rpos,ri) allocate(xsp(0:numberspheres),rpos(3,0:numberspheres),ri(2,0:numberspheres)) endif if (present(sphere_size_parameters)) xsp(1:numberspheres) =sphere_size_parameters(1:numberspheres) if (present(sphere_positions)) rpos(:,1:numberspheres) =sphere_positions(:,1:numberspheres) if (present(sphere_refractive_indices)) ri(:,1:numberspheres) =sphere_refractive_indices(:,1:numberspheres) if (present(volume_size_parameter)) xvsp =volume_size_parameter end subroutine setspheredata ! ! getrunparameters: retrieves run parameters read from input file ! ! ! last revised: 15 January 2011 ! 30 March 2011: added optical activity ! subroutine getrunparameters(number_spheres,sphere_position_file,output_file, & length_scale_factor,real_ref_index_scale_factor, & imag_ref_index_scale_factor,mie_epsilon,translation_epsilon,solution_epsilon, & max_number_iterations,fixed_or_random_orientation,scattering_plane_angle_deg, & min_scattering_angle_deg,max_scattering_angle_deg,number_scattering_angles, & incident_azimuth_angle_deg,incident_polar_angle_deg,calculate_near_field, & near_field_plane_coord,near_field_plane_position,near_field_plane_vertices, & spacial_step_size,polarization_angle_deg,near_field_output_file, & plane_wave_epsilon,t_matrix_convergence_epsilon,gaussian_beam_constant, & gaussian_beam_focal_point,calculate_t_matrix,t_matrix_file,run_print_file, & run_print_unit,calculate_scattering_coefficients,scattering_coefficient_file, & max_memory_per_processor,track_iterations,near_field_output_data, & real_chiral_factor,imag_chiral_factor,normalize_scattering_matrix, & store_translation_matrix,near_field_distance, & iterations_per_correction) implicit none integer, optional :: number_spheres,max_number_iterations,fixed_or_random_orientation, & number_scattering_angles,calculate_near_field,near_field_plane_coord, & calculate_t_matrix,run_print_unit,calculate_scattering_coefficients, & max_memory_per_processor,track_iterations,near_field_output_data, & normalize_scattering_matrix,store_translation_matrix, & iterations_per_correction real(8), optional :: length_scale_factor,real_ref_index_scale_factor, & imag_ref_index_scale_factor,mie_epsilon,translation_epsilon,solution_epsilon, & scattering_plane_angle_deg, & min_scattering_angle_deg,max_scattering_angle_deg, & incident_azimuth_angle_deg,incident_polar_angle_deg,t_matrix_convergence_epsilon, & near_field_plane_position,near_field_plane_vertices(2,2),spacial_step_size, & polarization_angle_deg,plane_wave_epsilon,gaussian_beam_constant, & gaussian_beam_focal_point(3),real_chiral_factor,imag_chiral_factor, & near_field_distance character*30, optional :: sphere_position_file,output_file,near_field_output_file, & t_matrix_file,run_print_file,scattering_coefficient_file if(present(number_spheres)) number_spheres =numberspheres if(present(sphere_position_file)) sphere_position_file =positionfile if(present(output_file)) output_file =outputfile if(present(length_scale_factor)) length_scale_factor =lengthscalefactor if(present(real_ref_index_scale_factor)) real_ref_index_scale_factor =realriscalefactor if(present(imag_ref_index_scale_factor)) imag_ref_index_scale_factor =imriscalefactor if(present(mie_epsilon)) mie_epsilon =epsmie if(present(translation_epsilon)) translation_epsilon =epstran if(present(solution_epsilon)) solution_epsilon =epssoln if(present(max_number_iterations)) max_number_iterations =numberiterations if(present(track_iterations)) track_iterations =trackiterations if(present(max_memory_per_processor)) max_memory_per_processor =maxmemperproc if(present(fixed_or_random_orientation)) fixed_or_random_orientation =fixedorrandom if(present(scattering_plane_angle_deg)) scattering_plane_angle_deg =phideg if(present(min_scattering_angle_deg)) min_scattering_angle_deg =thetamindeg if(present(max_scattering_angle_deg)) max_scattering_angle_deg =thetamaxdeg if(present(number_scattering_angles)) number_scattering_angles =numbertheta if(present(normalize_scattering_matrix)) normalize_scattering_matrix =normalizesm if(present(incident_azimuth_angle_deg)) incident_azimuth_angle_deg =alphadeg if(present(incident_polar_angle_deg)) incident_polar_angle_deg =betadeg if(present(t_matrix_convergence_epsilon)) t_matrix_convergence_epsilon =epstcon if(present(calculate_near_field)) calculate_near_field =calcnf if(present(near_field_plane_coord)) near_field_plane_coord =nfplane if(present(near_field_plane_position)) near_field_plane_position =nfplanepos if(present(near_field_plane_vertices)) near_field_plane_vertices =nfplanevert if(present(spacial_step_size)) spacial_step_size =deltax if(present(polarization_angle_deg)) polarization_angle_deg =gammadeg if(present(near_field_output_file)) near_field_output_file =nfoutputfile if(present(near_field_output_data)) near_field_output_data =nfoutdata if(present(plane_wave_epsilon)) plane_wave_epsilon =epspw if(present(gaussian_beam_constant)) gaussian_beam_constant =cgaussbeam if(present(gaussian_beam_focal_point)) gaussian_beam_focal_point =gaussbeamfocus if(present(t_matrix_file)) t_matrix_file =tmatrixfile if(present(calculate_t_matrix)) calculate_t_matrix =calctmatrix if(present(run_print_file)) run_print_file =printfile if(present(run_print_unit)) run_print_unit =runprintunit if(present(calculate_scattering_coefficients)) calculate_scattering_coefficients =calcamn if(present(scattering_coefficient_file)) scattering_coefficient_file =amnfile if(present(real_chiral_factor)) real_chiral_factor =realchiralfactor if(present(imag_chiral_factor)) imag_chiral_factor =imchiralfactor if(present(store_translation_matrix)) store_translation_matrix =storetranmat if(present(near_field_distance)) near_field_distance =nfdistance if(present(iterations_per_correction)) iterations_per_correction =niterstep end subroutine getrunparameters ! ! set run parameters: set run parameters ! ! ! last revised: 15 January 2011 ! 30 March 2011: added optical activity ! subroutine setrunparameters(number_spheres,sphere_position_file,output_file, & length_scale_factor,real_ref_index_scale_factor, & imag_ref_index_scale_factor,mie_epsilon,translation_epsilon,solution_epsilon, & max_number_iterations,fixed_or_random_orientation,scattering_plane_angle_deg, & min_scattering_angle_deg,max_scattering_angle_deg,number_scattering_angles, & incident_azimuth_angle_deg,incident_polar_angle_deg,calculate_near_field, & near_field_plane_coord,near_field_plane_position,near_field_plane_vertices, & spacial_step_size,polarization_angle_deg,near_field_output_file, & plane_wave_epsilon,t_matrix_convergence_epsilon,gaussian_beam_constant, & gaussian_beam_focal_point,calculate_t_matrix,t_matrix_file,run_print_file, & run_print_unit,calculate_scattering_coefficients,scattering_coefficient_file, & max_memory_per_processor,track_iterations,near_field_output_data, & real_chiral_factor,imag_chiral_factor,store_translation_matrix, & near_field_distance,iterations_per_correction) implicit none integer, optional :: number_spheres,max_number_iterations,fixed_or_random_orientation, & number_scattering_angles,calculate_near_field,near_field_plane_coord, & calculate_t_matrix,run_print_unit,calculate_scattering_coefficients, & max_memory_per_processor,track_iterations,near_field_output_data, & store_translation_matrix,iterations_per_correction real(8), optional :: length_scale_factor,real_ref_index_scale_factor, & imag_ref_index_scale_factor,mie_epsilon,translation_epsilon,solution_epsilon, & scattering_plane_angle_deg,near_field_distance,& min_scattering_angle_deg,max_scattering_angle_deg, & incident_azimuth_angle_deg,incident_polar_angle_deg,t_matrix_convergence_epsilon, & near_field_plane_position,near_field_plane_vertices(2,2),spacial_step_size, & polarization_angle_deg,plane_wave_epsilon,gaussian_beam_constant, & gaussian_beam_focal_point(3),real_chiral_factor,imag_chiral_factor character*30, optional :: sphere_position_file,output_file,near_field_output_file, & t_matrix_file,run_print_file,scattering_coefficient_file if(present(number_spheres)) numberspheres =number_spheres if(present(sphere_position_file)) positionfile =sphere_position_file if(present(output_file)) outputfile =output_file if(present(length_scale_factor)) lengthscalefactor =length_scale_factor if(present(real_ref_index_scale_factor)) realriscalefactor =real_ref_index_scale_factor if(present(imag_ref_index_scale_factor)) imriscalefactor =imag_ref_index_scale_factor if(present(mie_epsilon)) epsmie =mie_epsilon if(present(translation_epsilon)) epstran =translation_epsilon if(present(solution_epsilon)) epssoln =solution_epsilon if(present(max_number_iterations)) numberiterations =max_number_iterations if(present(track_iterations)) trackiterations =track_iterations if(present(max_memory_per_processor)) maxmemperproc =max_memory_per_processor if(present(fixed_or_random_orientation)) fixedorrandom =fixed_or_random_orientation if(present(scattering_plane_angle_deg)) phideg =scattering_plane_angle_deg if(present(min_scattering_angle_deg)) thetamindeg =min_scattering_angle_deg if(present(max_scattering_angle_deg)) thetamaxdeg =max_scattering_angle_deg if(present(number_scattering_angles)) numbertheta =number_scattering_angles if(present(incident_azimuth_angle_deg)) alphadeg =incident_azimuth_angle_deg if(present(incident_polar_angle_deg)) betadeg =incident_polar_angle_deg if(present(t_matrix_convergence_epsilon)) epstcon =t_matrix_convergence_epsilon if(present(calculate_near_field)) calcnf =calculate_near_field if(present(near_field_plane_coord)) nfplane =near_field_plane_coord if(present(near_field_plane_position)) nfplanepos =near_field_plane_position if(present(near_field_plane_vertices)) nfplanevert =near_field_plane_vertices if(present(spacial_step_size)) deltax =spacial_step_size if(present(polarization_angle_deg)) gammadeg =polarization_angle_deg if(present(near_field_output_file)) nfoutputfile =near_field_output_file if(present(near_field_output_data)) nfoutdata =near_field_output_data if(present(plane_wave_epsilon)) epspw =plane_wave_epsilon if(present(gaussian_beam_constant)) cgaussbeam =gaussian_beam_constant if(present(gaussian_beam_focal_point)) gaussbeamfocus =gaussian_beam_focal_point if(present(t_matrix_file)) tmatrixfile =t_matrix_file if(present(calculate_t_matrix)) calctmatrix =calculate_t_matrix if(present(run_print_file)) printfile =run_print_file if(present(run_print_unit)) runprintunit =run_print_unit if(present(calculate_scattering_coefficients)) calcamn =calculate_scattering_coefficients if(present(scattering_coefficient_file)) amnfile =scattering_coefficient_file if(present(real_chiral_factor)) realchiralfactor =real_chiral_factor if(present(imag_chiral_factor)) imchiralfactor =imag_chiral_factor if(present(store_translation_matrix)) storetranmat =store_translation_matrix if(present(near_field_distance)) nfdistance =near_field_distance if(present(iterations_per_correction)) niterstep =iterations_per_correction end subroutine setrunparameters end module spheredata ! ! module miecoefdata: used to 1) calculate single sphere mie coefficient values, ! 2) store values in an allocated array, 3) provide common access to values, and ! 4) perform multiplication of coefficient values with vectors containing VWH scattering ! coefficients. ! ! ! last revised: 15 January 2011 ! 30 March 2011: added optical activity ! module miecoefdata implicit none integer, private :: numeqns,maxorder integer, allocatable, private :: nodr(:),nodroffset(:),nblk(:),nblkoffset(:) real(8), allocatable, private :: qextmie(:),qabsmie(:) complex(8), allocatable, private :: anmie(:,:,:),cnmie(:,:,:) interface getmiedata module procedure getmiedataall, getmiedataone end interface getmiedata contains ! ! calculation of the max order of sphere expansions and storage of mie coefficients ! ! ! last revised: 15 January 2011 ! 30 March 2011: added optical activity ! subroutine miecoefcalc(nsphere,xsp,ri,qeps) implicit none integer :: n,nodrn,nsphere,nodrtot,ierr,nblktot real(8) :: qext,qabs,qsca,qeps,xsp(nsphere) complex(8) :: ri(2,nsphere) complex(8), allocatable :: anp(:,:,:),cnp(:,:,:) if(allocated(nodr)) deallocate(nodr,nodroffset,nblk, & nblkoffset,qextmie,qabsmie) allocate(nodr(nsphere),nodroffset(nsphere+1), & nblk(nsphere),nblkoffset(nsphere+1), & qextmie(nsphere),qabsmie(nsphere),stat=ierr) if(ierr.ne.0) then write(*,'('' bad allocation in nodr: stat:'',i4)') ierr endif nodrtot=0 nblktot=0 maxorder=0 ! ! calculate the order limits and efficiencies ! do n=1,nsphere call mieoa(xsp(n),ri(1,n),nodrn,qeps,qext,qsca) nodroffset(n)=nodrtot nblkoffset(n)=nblktot nodr(n)=nodrn maxorder=max(maxorder,nodrn) nblk(n)=nodrn*(nodrn+2)*2 nodrtot=nodrtot+nodrn nblktot=nblktot+nblk(n) qextmie(n)=qext qabsmie(n)=qext-qsca enddo nodroffset(nsphere+1)=nodrtot nblkoffset(nsphere+1)=nblktot numeqns=nblktot ! ! calculate the mie coefficients, and store in memory ! if(allocated(anmie)) deallocate(anmie,cnmie) allocate(anmie(2,2,nodrtot),cnmie(2,2,nodrtot),stat=ierr) if(ierr.ne.0) then write(*,'('' bad allocation in anmie: stat:'',i4)') ierr endif do n=1,nsphere if(abs(ri(1,n)-ri(2,n)).eq.0) then allocate(anp(2,1,nodr(n)),cnp(2,1,nodr(n))) call mieregular(xsp(n),ri(1,n),nodrn,qeps,qext,qsca,anp_mie=anp,cnp_mie=cnp) anmie(1,1,nodroffset(n)+1:nodroffset(n+1))=anp(1,1,1:nodr(n)) anmie(2,2,nodroffset(n)+1:nodroffset(n+1))=anp(2,1,1:nodr(n)) anmie(1,2,nodroffset(n)+1:nodroffset(n+1))=0.d0 anmie(2,1,nodroffset(n)+1:nodroffset(n+1))=0.d0 cnmie(1,1,nodroffset(n)+1:nodroffset(n+1))=cnp(1,1,1:nodr(n)) cnmie(2,2,nodroffset(n)+1:nodroffset(n+1))=cnp(2,1,1:nodr(n)) cnmie(1,2,nodroffset(n)+1:nodroffset(n+1))=0.d0 cnmie(2,1,nodroffset(n)+1:nodroffset(n+1))=0.d0 deallocate(anp,cnp) else allocate(anp(2,2,nodr(n)),cnp(2,2,nodr(n))) call mieoa(xsp(n),ri(1,n),nodrn,qeps,qext,qsca,anp_mie=anp,cnp_mie=cnp) anmie(1:2,1:2,nodroffset(n)+1:nodroffset(n+1))=anp(1:2,1:2,1:nodr(n)) cnmie(1:2,1:2,nodroffset(n)+1:nodroffset(n+1))=cnp(1:2,1:2,1:nodr(n)) deallocate(anp,cnp) endif enddo end subroutine miecoefcalc ! ! retrieve the array of mie data ! ! ! last revised: 15 January 2011 ! 30 March 2011: added optical activity ! subroutine getmiedataall(sphere_order, sphere_block, & sphere_order_offset, sphere_block_offset, sphere_qext, & sphere_qabs, sphere_mie_coefficients, sphere_int_mie_coefficients, & number_equations, max_order) use spheredata implicit none integer, optional :: sphere_order(:), sphere_block(:), sphere_order_offset(:), & sphere_block_offset(:),number_equations, max_order integer :: i,nsphere real(8), optional :: sphere_qext(:), sphere_qabs(:) complex(8), optional :: sphere_mie_coefficients(:,:,:,:), & sphere_int_mie_coefficients(:,:,:,:) call getspheredata(number_spheres=nsphere) if(present(sphere_order)) sphere_order=nodr if(present(sphere_block)) sphere_block=nblk if(present(sphere_order_offset)) sphere_order_offset=nodroffset if(present(sphere_block_offset)) sphere_block_offset=nblkoffset if(present(sphere_qext)) sphere_qext=qextmie if(present(sphere_qabs)) sphere_qabs=qabsmie if(present(number_equations)) number_equations=numeqns if(present(max_order)) max_order=maxorder if(present(sphere_mie_coefficients)) then do i=1,nsphere sphere_mie_coefficients(1:2,1:2,1:nodr(i),i) & =anmie(1:2,1:2,nodroffset(i)+1:nodroffset(i+1)) enddo endif if(present(sphere_int_mie_coefficients)) then do i=1,nsphere sphere_int_mie_coefficients(1:2,1:2,1:nodr(i),i) & =cnmie(1:2,1:2,nodroffset(i)+1:nodroffset(i+1)) enddo endif end subroutine getmiedataall ! ! retrieve mie data for a single sphere ! ! ! last revised: 15 January 2011 ! 30 March 2011: added optical activity ! subroutine getmiedataone(which_sphere, sphere_order, sphere_block, & sphere_order_offset, sphere_block_offset, sphere_qext, & sphere_qabs, sphere_mie_coefficients, sphere_int_mie_coefficients, & number_equations, max_order) use spheredata implicit none integer, optional :: sphere_order, sphere_block, sphere_order_offset, & sphere_block_offset, number_equations, max_order integer :: which_sphere integer :: i,nsphere real(8), optional :: sphere_qext, sphere_qabs complex(8), optional :: sphere_mie_coefficients(:,:,:), sphere_int_mie_coefficients(:,:,:) i=which_sphere if(present(sphere_order)) sphere_order=nodr(i) if(present(sphere_block)) sphere_block=nblk(i) if(present(sphere_order_offset)) sphere_order_offset=nodroffset(i) if(present(sphere_block_offset)) sphere_block_offset=nblkoffset(i) if(present(sphere_qext)) sphere_qext=qextmie(i) if(present(sphere_qabs)) sphere_qabs=qabsmie(i) if(present(number_equations)) number_equations=numeqns if(present(max_order)) max_order=maxorder if(present(sphere_mie_coefficients)) & sphere_mie_coefficients(1:2,1:2,1:nodr(i)) & =anmie(1:2,1:2,nodroffset(i)+1:nodroffset(i+1)) if(present(sphere_int_mie_coefficients)) & sphere_int_mie_coefficients(1:2,1:2,1:nodr(i)) & =cnmie(1:2,1:2,nodroffset(i)+1:nodroffset(i+1)) end subroutine getmiedataone ! ! retrieve mie coefficients for sphere n ! 30 March 2011: added optical activity ! function miecoef(n) implicit none integer :: n complex(8), dimension(2,2,nodr(n)) :: miecoef miecoef=anmie(1:2,1:2,nodroffset(n)+1:nodroffset(n+1)) end function miecoef function internalmiecoef(n) implicit none integer :: n complex(8), dimension(2,2,nodr(n)) :: internalmiecoef internalmiecoef=cnmie(1:2,1:2,nodroffset(n)+1:nodroffset(n+1)) end function internalmiecoef ! ! multiples the solution vector cx by mie coefficients and returns in y ! i1: starting sphere, i2: ending sphere ! ! ! last revised: 15 January 2011 ! 30 March 2011: added optical activity ! subroutine miecoeffmult(i1,i2,neqns,cx,cy) implicit none integer :: i1,i2,neqns,i,n,p,nodrvec(3),nodri,nblki,noffi,icon,q complex(8) :: cx(neqns),cy(neqns) complex(8), allocatable :: cxt(:,:,:),an1(:,:,:),cxtt(:,:,:) do i=i1,i2 nodri=nodr(i) nblki=nblk(i) noffi=nblkoffset(i) allocate(cxt(0:nodri+1,nodri,2),cxtt(0:nodri+1,nodri,2),an1(2,2,nodri)) cxt=reshape(cx(noffi+1:noffi+nblki),(/nodri+2,nodri,2/)) cxtt=0.d0 an1=miecoef(i) do n=1,nodri do p=1,2 cxtt(n+1,n:1:-1,p)=an1(p,1,n)*cxt(n+1,n:1:-1,1)+an1(p,2,n)*cxt(n+1,n:1:-1,2) cxtt(0:n,n,p)=an1(p,1,n)*cxt(0:n,n,1)+an1(p,2,n)*cxt(0:n,n,2) enddo enddo cy(noffi+1:noffi+nblki)=reshape(cxtt(0:nodri+1,1:nodri,1:2),(/nblki/)) deallocate(cxt,cxtt,an1) enddo end subroutine miecoeffmult subroutine internalmiecoeffmult(i1,i2,neqns,cx,cy) implicit none integer :: i1,i2,neqns,i,n,p,nodrvec(3),nodri,nblki,noffi,icon,q complex(8) :: cx(neqns),cy(neqns) complex(8), allocatable :: cxt(:,:,:),an1(:,:,:),cxtt(:,:,:) do i=i1,i2 nodri=nodr(i) nblki=nblk(i) noffi=nblkoffset(i) allocate(cxt(0:nodri+1,nodri,2),cxtt(0:nodri+1,nodri,2),an1(2,2,nodri)) cxt=reshape(cx(noffi+1:noffi+nblki),(/nodri+2,nodri,2/)) cxtt=0.d0 an1=internalmiecoef(n) do n=1,nodri do p=1,2 cxtt(n+1,n:1:-1,p)=an1(p,1,n)*cxt(n+1,n:1:-1,1)+an1(p,2,n)*cxt(n+1,n:1:-1,2) cxtt(0:n,n,p)=an1(p,1,n)*cxt(0:n,n,1)+an1(p,2,n)*cxt(0:n,n,2) enddo enddo cy(noffi+1:noffi+nblki)=reshape(cxtt(0:nodri+1,1:nodri,1:2),(/nblki/)) deallocate(cxt,cxtt,an1) enddo end subroutine internalmiecoeffmult ! ! single-sphere lorenz/mie coefficients ! ! ! last revised: 15 January 2011 ! subroutine mieregular(x,ri,nstop,qeps,qext,qsca,anp_mie,cnp_mie) use specialfuncs implicit none integer :: nstop,n,iancalc real(8) :: x,qeps,qext,qsca,prn,prp,qext1,err complex(8), optional :: anp_mie(2,*), cnp_mie(2,*) complex(8) :: ri,y,pcp,xip,da,db,na,nb,an1,an2,cn1,cn2 complex(8), allocatable :: pc(:),xi(:) ! ! modified LM criterion ! if(qeps.gt.0.) nstop=nint(x+4.*x**(1./3.))+15 ! ! user-set order limit ! if(qeps.lt.0) nstop=-qeps ! ! basic calculations follow ! allocate(pc(0:nstop),xi(0:nstop)) y=x*ri call cricbessel(nstop,y,pc) call richankel(nstop,x,xi) qsca=0.0 qext=0.0 do n=1,nstop prn=dble(xi(n)) pcp=pc(n-1)-n*pc(n)/y xip=xi(n-1)-n*xi(n)/x prp=dble(xip) da=ri*xip*pc(n)-xi(n)*pcp db=ri*xi(n)*pcp-xip*pc(n) na=ri*prp*pc(n)-prn*pcp nb=ri*prn*pcp-prp*pc(n) an1=-na/da an2=-nb/db cn1=-dcmplx(0.d0,1.d0)*ri/na cn2=dcmplx(0.d0,1.d0)*ri/nb if(present(anp_mie)) then anp_mie(1,n)=an1 anp_mie(2,n)=an2 endif if(present(cnp_mie)) then cnp_mie(1,n)=cn1 cnp_mie(2,n)=cn2 endif qsca=qsca+(n+n+1)*(cdabs(an1)*cdabs(an1) & +cdabs(an2)*cdabs(an2)) qext1=-(n+n+1)*dble(an1+an2) qext=qext+qext1 err=abs(qext1)/abs(qext) if(err.lt.qeps.or.n.eq.nstop) exit enddo nstop=n qsca=2./x/x*qsca qext=2./x/x*qext deallocate(pc,xi) end subroutine mieregular ! ! optically active lorenz/mie coefficients ! 30 March 2011 ! subroutine mieoa(x,ri,nstop,qeps,qext,qsca,anp_mie,cnp_mie) use specialfuncs implicit none integer :: nstop real(8) :: x,qeps,qext,qsca,fn1,err complex(8) :: ri(2) complex(8), optional :: anp_mie(2,2,*),cnp_mie(2,2,*) integer :: n,i,p,q real(8) :: psi,psip,qext1 complex (8) :: xri(2),xip,psicp,psic,wn(2),vn(2),an(2),bn(2), & den,xi,anct(2,2),cnct(2,2),ri0,ci complex(8), allocatable :: psicn(:,:),xin(:) data ci/(0.d0,1.d0)/ ri0=2.d0/(1.d0/ri(1)+1.d0/ri(2)) if(qeps.ge.0.) then nstop=nint(x+4.*x**(1./3.))+5. else nstop=-qeps endif allocate(psicn(0:nstop+1,2),xin(0:nstop+1)) do i=1,2 xri(i)=x*ri(i) call cricbessel(nstop+1,xri(i),psicn(0,i)) enddo call richankel(nstop+1,x,xin) qsca=0.0 qext=0.0 do n=1,nstop do i=1,2 psic=psicn(n,i) psicp=psicn(n-1,i)-dble(n)*psic/xri(i) xi=xin(n) xip=xin(n-1)-dble(n)*xi/x psi=dble(xi) psip=dble(xip) wn(i)=ri0*psic*xip-xi*psicp vn(i)=psic*xip-ri0*xi*psicp an(i)=ri0*psic*psip-psi*psicp bn(i)=psic*psip-ri0*psi*psicp enddo den=wn(1)*vn(2)+wn(2)*vn(1) anct(1,1)=-(vn(1)*an(2)+vn(2)*an(1))/den anct(2,2)=-(wn(1)*bn(2)+wn(2)*bn(1))/den anct(1,2)=(wn(1)*an(2)-wn(2)*an(1))/den anct(2,1)=anct(1,2) den=an(1)*bn(2)+an(2)*bn(1) cnct(1,1)=-ci*ri(1)*bn(2)/den cnct(1,2)=-ci*ri(1)*an(2)/den cnct(2,1)=ri(2)*ri0*bn(1)/den cnct(2,2)=-ri(2)*ri0*an(1)/den if(present(anp_mie)) then do p=1,2 do q=1,2 anp_mie(p,q,n)=anct(p,q) cnp_mie(p,q,n)=cnct(p,q) enddo enddo endif qext1=0.d0 fn1=n+n+1 do p=1,2 do q=1,2 qsca=qsca+fn1*cdabs(anct(p,q))*cdabs(anct(p,q)) enddo qext1=qext1-fn1*dble(anct(p,p)) enddo qext=qext+qext1 err=abs(qext1)/abs(qext) if(err.lt.qeps.or.n.eq.nstop) exit enddo nstop=min(n,nstop) qsca=2./x/x*qsca qext=2./x/x*qext return end subroutine mieoa end module miecoefdata ! ! module translation contains subroutines for VSWF translation and rotation ! ! ! last revised: 15 January 2011 ! module translation implicit none integer, private :: stored_max_order,store_tran_mat integer, allocatable, private :: nsizerot(:,:),nsizetran(:,:),nsizeephi(:,:), & noffrot(:,:),nofftran(:,:),noffephi(:,:) real(8), private :: near_field_distance real(8), allocatable, private :: sphere_position(:,:) real(8), target, allocatable, private :: rotmatstore(:) complex(8), target, allocatable, private :: tranmatstore(:), ephimatstore(:) complex(8), allocatable, private :: rvec_temp(:,:),tvec_temp(:,:),c_temp(:,:,:), & ct_temp(:,:,:),rvec2_temp(:,:),tvec2_temp(:,:),c2_temp(:,:,:), & ct2_temp(:,:,:) contains ! ! rotation of expansion coefficients amn by euler angles alpha,beta,gamma ! idir=1: forward rotation, idir=-1, reverse rotation. ! ! ! last revised: 15 January 2011 ! subroutine rotvec(alpha,beta,gamma,nmax,mmax,amn,idir) use numconstants use specialfuncs implicit none integer :: nmax,mmax,idir,k,n,m,in,kmax,kn,ka,na,p,im,m1 real(8) :: dc(-nmax-1:nmax+1,-nmax-1:nmax+1),dk0(-nmax-1:nmax+1), & dk01(-nmax-1:nmax+1),sbe,cbe,sbe2,cbe2,sben,dkt, & fmn,dkm0,dkm1,alpha,beta,gamma complex(8) :: ealpha,amn(0:nmax+1,nmax,2),ealpham(-nmax:nmax), & amnt(2,-nmax:nmax),a,b,ci,egamma,egammam(-nmax:nmax) data ci/(0.d0,1.d0)/ call init(nmax) dc=0.d0 dk01=0.d0 dk0=0.d0 ealpha=cdexp(ci*alpha) egamma=cdexp(ci*gamma) cbe=cos(beta) sbe=sqrt((1.d0+cbe)*(1.d0-cbe)) cbe2=.5d0*(1.d0+cbe) sbe2=.5d0*(1.d0-cbe) call ephicoef(ealpha,nmax,ealpham) call ephicoef(egamma,nmax,egammam) in=1 dk0(0)=1.d0 sben=1.d0 dk01(0)=0.d0 do n=1,nmax kmax=min(n,mmax) do k=-kmax,kmax if(k.le.-1) then ka=n+1 na=-k else ka=k na=n endif if(idir.eq.1) then amnt(1,k)=amn(ka,na,1)*ealpham(k) amnt(2,k)=amn(ka,na,2)*ealpham(k) else amnt(1,-k)=amn(ka,na,1)*egammam(k) amnt(2,-k)=amn(ka,na,2)*egammam(k) endif enddo in=-in sben=sben*sbe/2.d0 dk0(n)=in*sben*bcof(n,n) dk0(-n)=in*dk0(n) dk01(n)=0.d0 dk01(-n)=0.d0 dc(0,n)=dk0(n) dc(0,-n)=dk0(-n) do k=-n+1,n-1 dkt=dk01(k) dk01(k)=dk0(k) dk0(k)=(cbe*(n+n-1)*dk01(k)-fnr(n-k-1)*fnr(n+k-1)*dkt) & /(fnr(n+k)*fnr(n-k)) dc(0,k)=dk0(k) enddo im=1 do m=1,kmax im=-im fmn=1./fnr(n-m+1)/fnr(n+m) m1=m-1 dkm0=0. do k=-n,n dkm1=dkm0 dkm0=dc(m1,k) dc(m,k)=(fnr(n+k)*fnr(n-k+1)*cbe2*dkm1 & -fnr(n-k)*fnr(n+k+1)*sbe2*dc(m1,k+1) & -k*sbe*dc(m1,k))*fmn dc(-m,-k)=dc(m,k)*(-1)**(k)*im enddo enddo do m=-n,n if(m.le.-1) then ka=n+1 na=-m else ka=m na=n endif a=0. b=0. do k=-kmax,kmax a=a+dc(-k,-m)*amnt(1,k) b=b+dc(-k,-m)*amnt(2,k) enddo if(idir.eq.1) then amn(ka,na,1)=a*egammam(m) amn(ka,na,2)=b*egammam(m) else amn(ka,na,1)=a*ealpham(m) amn(ka,na,2)=b*ealpham(m) endif enddo enddo end subroutine rotvec ! ! sets up the stored translation matrices for mpi ! ! ! last revised: 15 January 2011 ! november 2011: added near and far field translation ! subroutine mpirottranmtrxsetup(nsphere,nodr,rpos,ri,istore,nfdistance,& runprintunit) use mpidefs use mpidata use intrinsics use numconstants use specialfuncs implicit none integer :: nsphere,nodr(nsphere),i,j,nodrmax,nodrmin,n,ntotrot,ntottran,ntotephi, & ierr,n1,n2,nt,rank,nsrank,runprintunit,isendok,tag,sendrank,numprocs,brank, & nsend,istore real(8) :: rpos(3,nsphere),xij(3),r,ct,memused(1),memusedmax(1),memusedmin(1), & nfdistance,nfdistancei real(8), allocatable :: rotmat(:,:) complex(8) :: ri,ephi complex(8), allocatable :: tranmat(:,:,:),ephimat(:),pivec(:,:,:) data isendok,tag/0,1/ numprocs=proc_per_group rank=group_rank brank=base_rank nsrank=mpi_sphere_number(rank) nodrmax=maxval(nodr) call init(nodrmax) store_tran_mat=istore near_field_distance=nfdistance if(allocated(sphere_position)) deallocate(sphere_position) allocate(sphere_position(3,nsphere)) sphere_position=rpos if(istore.eq.0) then return endif if(allocated(nsizerot)) deallocate(nsizerot,nsizetran,nsizeephi,noffrot,nofftran,noffephi) allocate(nsizerot(nsphere,nsphere),nsizetran(nsphere,nsphere),nsizeephi(nsphere,nsphere), & noffrot(nsphere,nsphere),nofftran(nsphere,nsphere),noffephi(nsphere,nsphere)) if(allocated(rvec_temp)) deallocate(rvec_temp,tvec_temp,c_temp,ct_temp, & rvec2_temp,tvec2_temp,c2_temp,ct2_temp) allocate(rvec_temp(-nodrmax:nodrmax,2),tvec_temp(nodrmax,2), & c_temp(-nodrmax:nodrmax,nodrmax,2),ct_temp(nodrmax,2,2), & rvec2_temp(-nodrmax:nodrmax,2),tvec2_temp(nodrmax,2), & c2_temp(-nodrmax:nodrmax,nodrmax,2),ct2_temp(nodrmax,2,2)) stored_max_order=nodrmax ! ! determine the memory requirements ! ntotrot=0 ntottran=0 ntotephi=0 do i=mpi_sphere_index(rank)+1,mpi_sphere_index(rank)+nsrank do j=1,nsphere xij(:)=rpos(:,i)-rpos(:,j) if(j.ne.i) then if(nfdistance.lt.0.) then nfdistancei=(.5*dble(nodr(i)+nodr(j)))**2. else nfdistancei=nfdistance endif r=sqrt(dot_product(xij,xij)) if(r.le.nfdistancei) then nodrmax=max(nodr(j),nodr(i)) nodrmin=min(nodr(j),nodr(i)) noffrot(i,j)=ntotrot nofftran(i,j)=ntottran noffephi(i,j)=ntotephi nsizerot(i,j)=(2*nodrmin+1)*(1+nodrmax*(nodrmax+2)) nsizetran(i,j)=nodr(i)*nodr(j)*(nodr(j)+3) nsizeephi(i,j)=2*nodrmax+1 ntotrot=ntotrot+nsizerot(i,j) ntottran=ntottran+nsizetran(i,j) ntotephi=ntotephi+nsizeephi(i,j) endif if(r.gt.nfdistancei.and.istore.eq.2) then nodrmax=max(nodr(j),nodr(i)) nofftran(i,j)=ntottran nsizetran(i,j)=2*nodrmax*(nodrmax+2) ntottran=ntottran+nsizetran(i,j) endif endif enddo enddo memused(1)=dble(8*ntotrot+16*(ntottran+ntotephi))*1.d-6 nsend=1 call ms_mpi(mpi_command='reduce',mpi_send_buf_dp=memused,mpi_recv_buf_dp=memusedmax,& mpi_number=1,mpi_rank=0,mpi_operation=ms_mpi_max) call ms_mpi(mpi_command='reduce',mpi_send_buf_dp=memused,mpi_recv_buf_dp=memusedmin,& mpi_number=1,mpi_rank=0,mpi_operation=ms_mpi_min) call ms_mpi(mpi_command='barrier') if(brank.eq.0) then write(runprintunit,'('' maximum translation matrix storage:'',f9.4,'' MB'')') memusedmax write(runprintunit,'('' minimum translation matrix storage:'',f9.4,'' MB'')') memusedmin call flush(runprintunit) endif ! ! calculate the matrices and store in memory ! if(allocated(rotmatstore)) deallocate(rotmatstore,tranmatstore,ephimatstore) allocate(rotmatstore(ntotrot),stat=ierr) allocate(tranmatstore(ntottran),stat=ierr) allocate(ephimatstore(ntotephi),stat=ierr) do i=mpi_sphere_index(rank)+1,mpi_sphere_index(rank)+nsrank do j=1,nsphere if(j.ne.i) then nodrmax=max(nodr(j),nodr(i)) nodrmin=min(nodr(j),nodr(i)) xij=rpos(:,i)-rpos(:,j) call cartosphere(xij,r,ct,ephi) if(nfdistance.lt.0.) then nfdistancei=(.5*dble(nodr(i)+nodr(j)))**2. else nfdistancei=nfdistance endif if(r.le.nfdistancei) then ! ! rotation matrix ! n1=noffrot(i,j)+1 nt=nsizerot(i,j) n2=n1+nt-1 allocate(rotmat(-nodrmin:nodrmin,0:nodrmax*(nodrmax+2))) call rotcoef(ct,nodrmin,nodrmax,rotmat) rotmatstore(n1:n2)=reshape(rotmat,(/nt/)) deallocate(rotmat) ! ! axial translation matrix ! n1=nofftran(i,j)+1 nt=nsizetran(i,j) n2=n1+nt-1 allocate(tranmat(nodr(i),nodr(j)*(nodr(j)+3)/2,2)) call axialtrancoef(3,r,ri,nodr(i),nodr(j),tranmat) tranmatstore(n1:n2)=reshape(tranmat,(/nt/)) deallocate(tranmat) ! ! ephi matrix ! n1=noffephi(i,j)+1 nt=nsizeephi(i,j) n2=n1+nt-1 allocate(ephimat(-nodrmax:nodrmax)) call ephicoef(ephi,nodrmax,ephimat) ephimatstore(n1:n2)=ephimat(-nodrmax:nodrmax) deallocate(ephimat) ! ! ff translation matrix storage ! elseif(istore.eq.2) then n1=nofftran(i,j)+1 nt=nsizetran(i,j) n2=n1+nt-1 nodrmax=max(nodr(j),nodr(i)) allocate(pivec(0:nodrmax+1,nodrmax,2)) call pifunc(ct,ephi,nodrmax,nodrmax,pivec) tranmatstore(n1:n2)=reshape(pivec,(/nt/)) deallocate(pivec) endif endif enddo enddo end subroutine mpirottranmtrxsetup ! ! clear the stored translation matrices ! ! ! last revised: 15 January 2011 ! subroutine rottranmtrxclear() implicit none if(allocated(rotmatstore)) deallocate(rotmatstore,tranmatstore,ephimatstore) if(allocated(sphere_position)) deallocate(sphere_position) end subroutine rottranmtrxclear ! ! translation coefficient vector cx by xij in medium with ri by rotation-translation ! itype: 1 or 3 ! icalc: =1, calculate matrices; = 0, use stored matrix ! idir: =1, translation of xij, =-1, -xij (reverse) ! itran=1, A(i-j) a(j), = -1, a(j) A(i-j) ! ! ! last revised: 15 January 2011 ! subroutine rottran(cx,cy,xij,ri,nodrx,nodry,itype,icalc,idir,itran) use numconstants use specialfuncs implicit none integer :: nodrx,nodry,itype,icalc,idir,itran,nmax,nmin,n,m,p,nblk real(8) :: xij(3),r,ct complex(8) :: ri,ephi,cx(0:nodrx+1,nodrx,2),cy(0:nodry+1,nodry,2) real(8), allocatable, save :: rotmat(:,:) complex(8), allocatable, save :: ephimat(:), tranmat(:,:,:) if(icalc.eq.1) then nmax=max(nodrx,nodry) nmin=min(nodrx,nodry) call cartosphere(xij,r,ct,ephi) if(r.lt.1.d-4) then do p=1,2 do n=1,nmin do m=0,nmin+1 cy(m,n,p)=cy(m,n,p)+cx(m,n,p) enddo enddo enddo return endif if(allocated(ephimat)) deallocate(rotmat,ephimat,tranmat) if(nmax.gt.stored_max_order) then if(allocated(rvec_temp)) deallocate(rvec_temp,tvec_temp, & c_temp,ct_temp,rvec2_temp,tvec2_temp, & c2_temp,ct2_temp) allocate(rvec_temp(-nmax:nmax,2),tvec_temp(nmax,2), & c_temp(-nmax:nmax,nmax,2),ct_temp(nmax,2,2), & rvec2_temp(-nmax:nmax,2),tvec2_temp(nmax,2), & c2_temp(-nmax:nmax,nmax,2),ct2_temp(nmax,2,2)) stored_max_order=nmax endif nblk=(nodrx*(nodrx+3))/2 allocate(rotmat(-nmin:nmin,0:nmax*(nmax+2))) allocate(ephimat(-nmax:nmax)) allocate(tranmat(1:nodry,1:nblk,1:2)) call rotcoef(ct,nmin,nmax,rotmat) ! call axialtrancoef(itype,r,ri,nodry,nodrx,tranmat) call axialtrancoefrecurrence(itype,r,ri,nodry,nodrx,tranmat) call ephicoef(ephi,nmax,ephimat) endif call rottranmtrx(cx,cy,idir,itran,nodrx,nodry,ephimat,rotmat,tranmat) return end subroutine rottran ! ! far field formula for outgoing SVWF translation ! October 2011 ! subroutine farfieldtranslation(cx,cy,xij,ri,nodrx,nodry,icase, & stored_pivec_matrix) use numconstants use specialfuncs implicit none integer :: nodrx,nodry,itype,icalc,icase,nmax,nmin,n,m,p,nblk,im real(8) :: xij(3),r,ct,xijt(3) complex(8) :: ri,ephi,cx(0:nodrx+1,nodrx,2),cy(0:nodry+1,nodry,2), & cxt(0:nodrx+1,nodrx,2),cyt(0:nodry+1,nodry,2), & sumx(2),c1,pivec(0:max(nodrx,nodry)+1,max(nodrx,nodry),2) complex(8), optional :: stored_pivec_matrix(0:max(nodrx,nodry)+1,max(nodrx,nodry),2) call cartosphere(xij,r,ct,ephi) nmax=max(nodrx,nodry) if(present(stored_pivec_matrix)) then pivec=stored_pivec_matrix else call pifunc(ct,ephi,nmax,nmax,pivec) endif if(icase.eq.1) then sumx(1)=sum(pivec(0:nodrx+1,1:nodrx,1:2)*cx(0:nodrx+1,1:nodrx,1:2)) sumx(2)=sum(pivec(0:nodrx+1,1:nodrx,2:1:-1)*cx(0:nodrx+1,1:nodrx,1:2)) sumx=sumx*cdexp((0.d0,1.d0)*ri*r)/((0.d0,1.d0)*ri*r)*8.d0 cyt(0:nodry+1,1:nodry,1) = conjg(pivec(0:nodry+1,1:nodry,1))*sumx(1) & +conjg(pivec(0:nodry+1,1:nodry,2))*sumx(2) cyt(0:nodry+1,1:nodry,2) = conjg(pivec(0:nodry+1,1:nodry,2))*sumx(1) & +conjg(pivec(0:nodry+1,1:nodry,1))*sumx(2) else do n=1,nodrx do p=1,2 im=(-1)**(n+p) cxt(n+1,1:n,p)=im*cx(n+1,1:n,p) cxt(0:n,n,p)=im*cx(0:n,n,p) enddo enddo sumx(1)=sum(conjg(pivec(0:nodrx+1,1:nodrx,1:2))*cxt(0:nodrx+1,1:nodrx,1:2)) sumx(2)=sum(conjg(pivec(0:nodrx+1,1:nodrx,2:1:-1))*cxt(0:nodrx+1,1:nodrx,1:2)) sumx=sumx*cdexp((0.d0,1.d0)*ri*r)/((0.d0,1.d0)*ri*r)*8.d0 cyt(0:nodry+1,1:nodry,1) = pivec(0:nodry+1,1:nodry,1)*sumx(1) & +pivec(0:nodry+1,1:nodry,2)*sumx(2) cyt(0:nodry+1,1:nodry,2) = pivec(0:nodry+1,1:nodry,2)*sumx(1) & +pivec(0:nodry+1,1:nodry,1)*sumx(2) do n=1,nodry do p=1,2 im=(-1)**(n+p) cyt(n+1,1:n,p)=im*cyt(n+1,1:n,p) cyt(0:n,n,p)=im*cyt(0:n,n,p) enddo enddo endif cy=cy+cyt end subroutine farfieldtranslation ! ! far field translation: normal and transpose, for bcgm solution ! october 2011 ! subroutine farfieldtranslationtwovec(cx1,cx2,cy1,cy2,xij,ri,nodrx,nodry, & stored_pivec_matrix) use numconstants use specialfuncs implicit none integer :: nodrx,nodry,itype,icalc,icase,nmax,nmin,n,m,p,nblk,im real(8) :: xij(3),r,ct,xijt(3) complex(8) :: ri,ephi,cx1(0:nodrx+1,nodrx,2),cy1(0:nodry+1,nodry,2), & cx2(0:nodrx+1,nodrx,2),cy2(0:nodry+1,nodry,2), & cxt(0:nodrx+1,nodrx,2),cyt1(0:nodry+1,nodry,2), & cyt2(0:nodry+1,nodry,2), & sumx(2),c1,phasefunc, & pivec(0:max(nodrx,nodry)+1,max(nodrx,nodry),2) complex(8), optional :: stored_pivec_matrix(0:max(nodrx,nodry)+1,max(nodrx,nodry),2) call cartosphere(xij,r,ct,ephi) nmax=max(nodrx,nodry) if(present(stored_pivec_matrix)) then pivec=stored_pivec_matrix else call pifunc(ct,ephi,nmax,nmax,pivec) endif phasefunc=cdexp((0.d0,1.d0)*ri*r)/((0.d0,1.d0)*ri*r)*8.d0 sumx(1)=sum(pivec(0:nodrx+1,1:nodrx,1:2)*cx1(0:nodrx+1,1:nodrx,1:2)) sumx(2)=sum(pivec(0:nodrx+1,1:nodrx,2:1:-1)*cx1(0:nodrx+1,1:nodrx,1:2)) sumx=sumx*phasefunc cyt1(0:nodry+1,1:nodry,1) = conjg(pivec(0:nodry+1,1:nodry,1))*sumx(1) & +conjg(pivec(0:nodry+1,1:nodry,2))*sumx(2) cyt1(0:nodry+1,1:nodry,2) = conjg(pivec(0:nodry+1,1:nodry,2))*sumx(1) & +conjg(pivec(0:nodry+1,1:nodry,1))*sumx(2) do n=1,nodrx do p=1,2 im=(-1)**(n+p) cxt(n+1,1:n,p)=im*cx2(n+1,1:n,p) cxt(0:n,n,p)=im*cx2(0:n,n,p) enddo enddo sumx(1)=sum(conjg(pivec(0:nodrx+1,1:nodrx,1:2))*cxt(0:nodrx+1,1:nodrx,1:2)) sumx(2)=sum(conjg(pivec(0:nodrx+1,1:nodrx,2:1:-1))*cxt(0:nodrx+1,1:nodrx,1:2)) sumx=sumx*phasefunc cyt2(0:nodry+1,1:nodry,1) = pivec(0:nodry+1,1:nodry,1)*sumx(1) & +pivec(0:nodry+1,1:nodry,2)*sumx(2) cyt2(0:nodry+1,1:nodry,2) = pivec(0:nodry+1,1:nodry,2)*sumx(1) & +pivec(0:nodry+1,1:nodry,1)*sumx(2) do n=1,nodry do p=1,2 im=(-1)**(n+p) cyt2(n+1,1:n,p)=im*cyt2(n+1,1:n,p) cyt2(0:n,n,p)=im*cyt2(0:n,n,p) enddo enddo cy1=cy1+cyt1 cy2=cy2+cyt2 end subroutine farfieldtranslationtwovec ! ! correction term for hybrid bcgm solution: difference between exact and ! ff translation field ! november 2011 ! subroutine fftranslationerror(cx,cy,jx,iy,nodrx,nodry) use numconstants use specialfuncs implicit none integer :: nodrx,nodry,idir,itran,iy,jx,istore integer :: nr1,nr2,nt1,nt2,ne1,ne2 real(8) :: xj(3),xi(3),xij(3),rij,nfdist complex(8) :: cx(0:nodrx+1,nodrx,2),cy(0:nodry+1,nodry,2), & cyt(0:nodry+1,nodry,2) xj(:)=sphere_position(:,jx) xi(:)=sphere_position(:,iy) xij=xi-xj rij=sqrt(dot_product(xij,xij)) if(near_field_distance.lt.0.) then nfdist=(.5*(nodrx+nodry))**2. else nfdist=near_field_distance endif if(rij.gt.nfdist) then cyt=0.d0 call farfieldtranslation(cx,cyt,xij,(1.d0,0.d0),nodrx,nodry,1) cyt=-cyt call rottran(cx,cyt,xij,(1.d0,0.d0),nodrx,nodry,3,1,1,1) cy=cy+cyt endif end subroutine fftranslationerror ! ! translation via stored or calculated matrices (replaces rottranstoredmatrix) ! ! 12 October 2011. ! if rij> near_field_distance, the far field formula is ! applied. ! subroutine rottranjtoi(cx,cy,jx,iy,nodrx,nodry,idir,itran) use numconstants use specialfuncs implicit none integer :: nodrx,nodry,idir,itran,iy,jx,istore integer :: nr1,nr2,nt1,nt2,ne1,ne2 real(8) :: xj(3),xi(3),xij(3),rij,nfdist complex(8) :: cx(0:nodrx+1,nodrx,2),cy(0:nodry+1,nodry,2) xj(:)=sphere_position(:,jx) xi(:)=sphere_position(:,iy) xij=xi-xj rij=sqrt(dot_product(xij,xij)) if(near_field_distance.lt.0.) then nfdist=(.5*(nodrx+nodry))**2. else nfdist=near_field_distance endif if(rij.gt.nfdist) then if(store_tran_mat.eq.2) then nt1=nofftran(iy,jx)+1 nt2=nt1+nsizetran(iy,jx)-1 call farfieldtranslation(cx,cy,xij,(1.d0,0.d0),nodrx,nodry,itran, & stored_pivec_matrix=tranmatstore(nt1:nt2)) else call farfieldtranslation(cx,cy,xij,(1.d0,0.d0),nodrx,nodry,itran) endif else if(store_tran_mat.eq.0) then call rottran(cx,cy,xij,(1.d0,0.d0),nodrx,nodry,3,1,idir,itran) else nr1=noffrot(iy,jx)+1 nr2=nr1+nsizerot(iy,jx)-1 nt1=nofftran(iy,jx)+1 nt2=nt1+nsizetran(iy,jx)-1 ne1=noffephi(iy,jx)+1 ne2=ne1+nsizeephi(iy,jx)-1 call rottranmtrx(cx,cy,idir,itran,nodrx,nodry,ephimatstore(ne1:ne2), & rotmatstore(nr1:nr2),tranmatstore(nt1:nt2)) endif endif end subroutine rottranjtoi ! ! normal and transpose translation, for bcgm ! november 2011 ! subroutine rottrantwojtoi(cx1,cx2,cy1,cy2,jx,iy,nodrx,nodry) use numconstants use specialfuncs implicit none integer :: nodrx,nodry,idir,itran,iy,jx,istore integer :: nr1,nr2,nt1,nt2,ne1,ne2 real(8) :: xj(3),xi(3),xij(3),rij,nfdist complex(8) :: cx1(0:nodrx+1,nodrx,2),cy1(0:nodry+1,nodry,2), & cx2(0:nodrx+1,nodrx,2),cy2(0:nodry+1,nodry,2) xj(:)=sphere_position(:,jx) xi(:)=sphere_position(:,iy) xij=xi-xj rij=sqrt(dot_product(xij,xij)) if(near_field_distance.lt.0.) then nfdist=(.5*(nodrx+nodry))**2. else nfdist=near_field_distance endif if(rij.gt.nfdist) then if(store_tran_mat.eq.2) then nt1=nofftran(iy,jx)+1 nt2=nt1+nsizetran(iy,jx)-1 call farfieldtranslationtwovec(cx1,cx2,cy1,cy2,xij,(1.d0,0.d0),nodrx,nodry, & stored_pivec_matrix=tranmatstore(nt1:nt2)) else call farfieldtranslationtwovec(cx1,cx2,cy1,cy2,xij,(1.d0,0.d0),nodrx,nodry) endif else if(store_tran_mat.eq.0) then call rottran(cx1,cy1,xij,(1.d0,0.d0),nodrx,nodry,3,1,1,1) call rottran(cx2,cy2,xij,(1.d0,0.d0),nodrx,nodry,3,0,-1,-1) else nr1=noffrot(iy,jx)+1 nr2=nr1+nsizerot(iy,jx)-1 nt1=nofftran(iy,jx)+1 nt2=nt1+nsizetran(iy,jx)-1 ne1=noffephi(iy,jx)+1 ne2=ne1+nsizeephi(iy,jx)-1 call rottranmtrxtwovec(cx1,cx2,cy1,cy2,nodrx,nodry,ephimatstore(ne1:ne2), & rotmatstore(nr1:nr2),tranmatstore(nt1:nt2)) endif endif end subroutine rottrantwojtoi ! ! the vectorized rotation-translation-rotation operation ! ! ! last revised: 15 January 2011 ! subroutine rottranmtrx(cx,cy,idir,itran,nodrx,nodry,ephimat,rotmat,tranmat) use numconstants use specialfuncs implicit none integer :: nodrx,nodry,itype,icalc,idir,itran,nmax,nmin integer :: m,n,k,l,nn1,nn2,ll1,mn,kl,m1,p,n1,addr(2) real(8), target :: rotmat(-min(nodrx,nodry):min(nodrx,nodry), & 0:max(nodrx,nodry)*(max(nodrx,nodry)+2)) real(8), pointer :: rmat(:,:) complex(8) :: cx(0:nodrx+1,nodrx,2),cy(0:nodry+1,nodry,2), & ephimat(-max(nodrx,nodry):max(nodrx,nodry)) complex(8), target :: tranmat(nodry,nodrx*(nodrx+3)/2,2) complex(8), pointer :: tmat1(:,:),tmat2(:,:) c_temp=(0.d0,0.d0) nmin=min(nodrx,nodry) ! ! rotation to origin of target ! do n=1,nodrx nn1=n*(n+1)-n nn2=nn1+(2*n+1)-1 n1=min(n,nodry) rmat=>rotmat(-n1:n1,nn1:nn2) do p=1,2 rvec_temp(-n:-1,p)=cx(n+1,n:1:-1,p) rvec_temp(0:n,p)=cx(0:n,n,p) if(itran.eq.1) then rvec_temp(-n:n,p)=rvec_temp(-n:n,p)*ephimat(-n:n) else rvec_temp(-n:n,p)=rvec_temp(-n:n,p)*conjg(ephimat(-n:n)) endif enddo c_temp(-n1:n1,n,1:2)=matmul(rmat,rvec_temp(-n:n,1:2)) enddo ! ! axial translation to target ! do m=0,nmin m1=max(1,m) nn1=atcadd(m,m1,nodrx) nn2=atcadd(m,nodrx,nodrx) tmat1=>tranmat(m1:nodry,nn1:nn2,1) tmat2=>tranmat(m1:nodry,nn1:nn2,2) tvec_temp(m1:nodrx,1)=idir*c_temp(m,m1:nodrx,1) tvec_temp(m1:nodrx,2)=c_temp(m,m1:nodrx,2) if(itran*idir.eq.-1) then tvec_temp(m1:nodrx,1)=tvec_temp(m1:nodrx,1)*monen(m1:nodrx) tvec_temp(m1:nodrx,2)=tvec_temp(m1:nodrx,2)*monen(m1:nodrx) endif ct_temp=(0.d0,0.d0) ct_temp(m1:nodry,1,1:2)=matmul(tmat1,tvec_temp(m1:nodrx,1:2)) ct_temp(m1:nodry,2,1:2)=matmul(tmat2,tvec_temp(m1:nodrx,1:2)) c_temp(m,m1:nodry,1)=idir*(ct_temp(m1:nodry,1,1)+ct_temp(m1:nodry,2,2)) c_temp(m,m1:nodry,2)=ct_temp(m1:nodry,2,1)+ct_temp(m1:nodry,1,2) if(itran*idir.eq.-1) then c_temp(m,m1:nodry,1)=c_temp(m,m1:nodry,1)*monen(m1:nodry) c_temp(m,m1:nodry,2)=c_temp(m,m1:nodry,2)*monen(m1:nodry) endif if(m.gt.0) then tvec_temp(m1:nodrx,1)=idir*c_temp(-m,m1:nodrx,1) tvec_temp(m1:nodrx,2)=c_temp(-m,m1:nodrx,2) if(itran*idir.eq.-1) then tvec_temp(m1:nodrx,1)=tvec_temp(m1:nodrx,1)*monen(m1:nodrx) tvec_temp(m1:nodrx,2)=tvec_temp(m1:nodrx,2)*monen(m1:nodrx) endif ct_temp=(0.d0,0.d0) ct_temp(m1:nodry,1,1:2)=matmul(tmat1,tvec_temp(m1:nodrx,1:2)) ct_temp(m1:nodry,2,1:2)=matmul(tmat2,tvec_temp(m1:nodrx,1:2)) c_temp(-m,m1:nodry,1)=idir*(ct_temp(m1:nodry,1,1)-ct_temp(m1:nodry,2,2)) c_temp(-m,m1:nodry,2)=-ct_temp(m1:nodry,2,1)+ct_temp(m1:nodry,1,2) if(itran*idir.eq.-1) then c_temp(-m,m1:nodry,1)=c_temp(-m,m1:nodry,1)*monen(m1:nodry) c_temp(-m,m1:nodry,2)=c_temp(-m,m1:nodry,2)*monen(m1:nodry) endif endif enddo ! ! rotation back to original frame ! do n=1,nodry rvec_temp=(0.d0,0.d0) m1=min(n,nmin) nn1=n*(n+1)-n nn2=n*(n+1)+n rmat=>rotmat(-m1:m1,nn1:nn2) do p=1,2 if(itran.eq.1) then rvec_temp(-n:n,p)=matmul(c_temp(-m1:m1,n,p),rmat)*conjg(ephimat(-n:n)) else rvec_temp(-n:n,p)=matmul(c_temp(-m1:m1,n,p),rmat)*ephimat(-n:n) endif cy(n+1,n:1:-1,p)=cy(n+1,n:1:-1,p)+rvec_temp(-n:-1,p) cy(0:n,n,p)=cy(0:n,n,p)+rvec_temp(0:n,p) enddo enddo end subroutine rottranmtrx ! ! two vector rotation: normal and transpose ! november 2011 ! subroutine rottranmtrxtwovec(cx1,cx2,cy1,cy2,nodrx,nodry, & ephimat,rotmat,tranmat) use numconstants use specialfuncs implicit none integer :: nodrx,nodry,itype,icalc,idir,itran,nmax,nmin integer :: m,n,k,l,nn1,nn2,ll1,mn,kl,m1,p,n1,addr(2) real(8), target :: rotmat(-min(nodrx,nodry):min(nodrx,nodry), & 0:max(nodrx,nodry)*(max(nodrx,nodry)+2)) real(8), pointer :: rmat(:,:) complex(8) :: cx1(0:nodrx+1,nodrx,2),cy1(0:nodry+1,nodry,2), & cx2(0:nodrx+1,nodrx,2),cy2(0:nodry+1,nodry,2), & ephimat(-max(nodrx,nodry):max(nodrx,nodry)) complex(8), target :: tranmat(nodry,nodrx*(nodrx+3)/2,2) complex(8), pointer :: tmat1(:,:),tmat2(:,:) c_temp=(0.d0,0.d0) nmin=min(nodrx,nodry) ! ! rotation to origin of target ! do n=1,nodrx nn1=n*(n+1)-n nn2=nn1+(2*n+1)-1 n1=min(n,nodry) rmat=>rotmat(-n1:n1,nn1:nn2) do p=1,2 rvec_temp(-n:-1,p)=cx1(n+1,n:1:-1,p) rvec_temp(0:n,p)=cx1(0:n,n,p) rvec2_temp(-n:-1,p)=cx2(n+1,n:1:-1,p) rvec2_temp(0:n,p)=cx2(0:n,n,p) rvec_temp(-n:n,p)=rvec_temp(-n:n,p)*ephimat(-n:n) rvec2_temp(-n:n,p)=rvec2_temp(-n:n,p)*conjg(ephimat(-n:n)) enddo c_temp(-n1:n1,n,1:2)=matmul(rmat,rvec_temp(-n:n,1:2)) c2_temp(-n1:n1,n,1:2)=matmul(rmat,rvec2_temp(-n:n,1:2)) enddo ! ! axial translation to target ! do m=0,nmin m1=max(1,m) nn1=atcadd(m,m1,nodrx) nn2=atcadd(m,nodrx,nodrx) tmat1=>tranmat(m1:nodry,nn1:nn2,1) tmat2=>tranmat(m1:nodry,nn1:nn2,2) tvec_temp(m1:nodrx,1)=c_temp(m,m1:nodrx,1) tvec_temp(m1:nodrx,2)=c_temp(m,m1:nodrx,2) tvec2_temp(m1:nodrx,1)=-c2_temp(m,m1:nodrx,1) tvec2_temp(m1:nodrx,2)=c2_temp(m,m1:nodrx,2) ct_temp=(0.d0,0.d0) ct2_temp=(0.d0,0.d0) ct_temp(m1:nodry,1,1:2)=matmul(tmat1,tvec_temp(m1:nodrx,1:2)) ct_temp(m1:nodry,2,1:2)=matmul(tmat2,tvec_temp(m1:nodrx,1:2)) ct2_temp(m1:nodry,1,1:2)=matmul(tmat1,tvec2_temp(m1:nodrx,1:2)) ct2_temp(m1:nodry,2,1:2)=matmul(tmat2,tvec2_temp(m1:nodrx,1:2)) c_temp(m,m1:nodry,1)=(ct_temp(m1:nodry,1,1)+ct_temp(m1:nodry,2,2)) c_temp(m,m1:nodry,2)=ct_temp(m1:nodry,2,1)+ct_temp(m1:nodry,1,2) c2_temp(m,m1:nodry,1)=-(ct2_temp(m1:nodry,1,1)+ct2_temp(m1:nodry,2,2)) c2_temp(m,m1:nodry,2)=ct2_temp(m1:nodry,2,1)+ct2_temp(m1:nodry,1,2) if(m.gt.0) then tvec_temp(m1:nodrx,1)=c_temp(-m,m1:nodrx,1) tvec_temp(m1:nodrx,2)=c_temp(-m,m1:nodrx,2) tvec2_temp(m1:nodrx,1)=-c2_temp(-m,m1:nodrx,1) tvec2_temp(m1:nodrx,2)=c2_temp(-m,m1:nodrx,2) ct_temp=(0.d0,0.d0) ct2_temp=(0.d0,0.d0) ct_temp(m1:nodry,1,1:2)=matmul(tmat1,tvec_temp(m1:nodrx,1:2)) ct_temp(m1:nodry,2,1:2)=matmul(tmat2,tvec_temp(m1:nodrx,1:2)) ct2_temp(m1:nodry,1,1:2)=matmul(tmat1,tvec2_temp(m1:nodrx,1:2)) ct2_temp(m1:nodry,2,1:2)=matmul(tmat2,tvec2_temp(m1:nodrx,1:2)) c_temp(-m,m1:nodry,1)=(ct_temp(m1:nodry,1,1)-ct_temp(m1:nodry,2,2)) c_temp(-m,m1:nodry,2)=-ct_temp(m1:nodry,2,1)+ct_temp(m1:nodry,1,2) c2_temp(-m,m1:nodry,1)=-(ct2_temp(m1:nodry,1,1)-ct2_temp(m1:nodry,2,2)) c2_temp(-m,m1:nodry,2)=-ct2_temp(m1:nodry,2,1)+ct2_temp(m1:nodry,1,2) endif enddo ! ! rotation back to original frame ! do n=1,nodry rvec_temp=(0.d0,0.d0) rvec2_temp=(0.d0,0.d0) m1=min(n,nmin) nn1=n*(n+1)-n nn2=n*(n+1)+n rmat=>rotmat(-m1:m1,nn1:nn2) do p=1,2 rvec_temp(-n:n,p)=matmul(c_temp(-m1:m1,n,p),rmat)*conjg(ephimat(-n:n)) rvec2_temp(-n:n,p)=matmul(c2_temp(-m1:m1,n,p),rmat)*ephimat(-n:n) cy1(n+1,n:1:-1,p)=cy1(n+1,n:1:-1,p)+rvec_temp(-n:-1,p) cy1(0:n,n,p)=cy1(0:n,n,p)+rvec_temp(0:n,p) cy2(n+1,n:1:-1,p)=cy2(n+1,n:1:-1,p)+rvec2_temp(-n:-1,p) cy2(0:n,n,p)=cy2(0:n,n,p)+rvec2_temp(0:n,p) enddo enddo end subroutine rottranmtrxtwovec ! ! GB coefficients for sphere-centered expansions, obtained via translation ! ! last revised: 15 January 2011 ! subroutine spheregaussianbeamcoef(nsphere,neqns,nodr,alpha,beta,cbeam, & rpos,rbeam,epstran,pmnp) use specialfuncs implicit none integer :: m,n,p,nsphere,i,l,nodr(nsphere),nblk,noff,nodrgb,neqns,k real(8) :: alpha,beta,cb,sb,ca,sa,rpos(3,nsphere),rmax,rbeam(3),xib(3),rib, & cbeam,epstran complex(8) :: pmnp(neqns,2) complex(8), allocatable :: pmnp0(:,:,:,:) nodrgb=0 rmax=0.d0 do i=1,nsphere xib(:)=rpos(:,i)-rbeam(:) rib=sqrt(dot_product(xib,xib)) rmax=max(rmax,rib) call tranordertest(rib,(1.d0,0.d0),nodr(i),epstran,n) nodrgb=max(n,nodrgb) enddo allocate(pmnp0(0:nodrgb+1,nodrgb,2,2)) call gaussianbeamcoef(alpha,beta,cbeam,nodrgb,pmnp0) pmnp=0.d0 noff=0 do i=1,nsphere nblk=2*nodr(i)*(nodr(i)+2) xib(:)=rpos(:,i)-rbeam(:) do k=1,2 call rottran(pmnp0(0:nodrgb+1,1:nodrgb,1:2,k),pmnp(noff+1:noff+nblk,k),xib, & (1.d0,0.d0),nodrgb,nodr(i),1,1,1,1) enddo noff=noff+nblk enddo deallocate(pmnp0) end subroutine spheregaussianbeamcoef end module translation ! ! scatprops module: various subroutines for calculation of observables from the solution ! ! ! last revised: 15 January 2011 ! module scatprops implicit none contains ! ! determination of maximum orders for target--based expansions ! ! ! last revised: 15 January 2011 ! subroutine tranorders(nsphere,nodr,rpos,eps,ntran,nodrt) use numconstants use specialfuncs use translation implicit none integer :: nsphere,nodr(nsphere),nodrt,ntran(nsphere),i real(8) :: rpos(3,nsphere),r,eps nodrt=0 do i=1,nsphere r=sqrt(dot_product(rpos(:,i),rpos(:,i))) call tranordertest(r,(1.d0,0.d0),nodr(i),eps,ntran(i)) if(print_intermediate_results.eq.1) & write(*,'('' i, nodr, ntran:'',3i7)') i,nodr(i),ntran(i) nodrt=max(nodrt,ntran(i)) enddo end subroutine tranorders ! ! translation of sphere-based expansions to common target origin ! ! ! last revised: 15 January 2011 ! subroutine amncommonorigin(neqns,nsphere,nodr,ntran,nodrt,rpos,amnp,amnp0) use specialfuncs use translation implicit none integer :: neqns,nsphere,nodr(nsphere),nodrt,i,m,n,p,nblk,ntran(nsphere),noff real(8) :: rpos(3,nsphere),r,eps,xij(3) complex(8) :: amnp(neqns),amnp0(0:nodrt+1,nodrt,2) complex(8), allocatable :: amnpt(:,:,:) amnp0=(0.d0,0.d0) noff=0 do i=1,nsphere allocate(amnpt(0:ntran(i)+1,ntran(i),2)) amnpt=(0.d0,0.d0) nblk=nodr(i)*(nodr(i)+2)*2 xij=-rpos(:,i) call rottran(amnp(noff+1:noff+nblk),amnpt,xij,(1.d0,0.d0), & nodr(i),ntran(i),1,1,1,1) do p=1,2 do n=1,ntran(i) do m=0,ntran(i)+1 amnp0(m,n,p)=amnp0(m,n,p)+amnpt(m,n,p) enddo enddo enddo deallocate(amnpt) noff=noff+nblk enddo end subroutine amncommonorigin ! ! sphereqeff computes the efficiency factors for the sphere, given an1: mie coefficients, ! anp: scattering coefficients, pnp: incident field coefficients. ! ! This subroutine is specific to the OA model for the sphere. ! ! ! original: 15 January 2011 ! revised: 21 February 2011: polarized and cross-polarized efficiency calculation ! 30 March 2011: added optical activity ! subroutine sphereqeff(nsphere,neqns,nodr,nodrmax,xsp,anp1,anp2,& pnp1,pnp2,qext,qabs,qsca) use miecoefdata use spheredata implicit none integer :: nsphere,m,n,p,i,nodr(nsphere),nblk,noff,neqns,nodrmax real(8) :: xsp(nsphere),qext(nsphere),qabs(nsphere),qsca(nsphere), & qe,qa,qs complex(8) :: anp1(neqns),pnp1(neqns),anp2(neqns),pnp2(neqns) complex(8) :: anmie(2,2,nodrmax) qext=0.d0 qabs=0.d0 qsca=0.d0 noff=0 do i=1,nsphere nblk=nodr(i)*(nodr(i)+2)*2 call getmiedata(which_sphere=i,sphere_mie_coefficients=anmie) call qeffcalc(nodr(i),anp1(noff+1:noff+nblk),anp2(noff+1:noff+nblk), & pnp1(noff+1:noff+nblk),pnp2(noff+1:noff+nblk),anmie,qe,qa,qs) noff=noff+nblk qext(i)=2.d0*qe/xsp(i)/xsp(i) qabs(i)=2.d0*qa/xsp(i)/xsp(i) qsca(i)=2.d0*qs/xsp(i)/xsp(i) enddo end subroutine sphereqeff ! ! calculation of sphere efficiency factors for scattered and incident field ! coefficient anp1, pnp1, anp2, pnp2 and mie coefficients anmie ! ! original: 15 January 2011 ! revised: 21 February 2011: polarized and cross-polarized efficiency calculation ! 30 March 2011: added optical activity ! subroutine qeffcalc(nodr,anp1,anp2,pnp1,pnp2,anmie,qe,qa,qs) implicit none integer :: nodr,m,n,p,q real(8) :: qe,qa,qs,babs,aninv(2,2) complex(8) :: anp1(0:nodr+1,nodr,2),pnp1(0:nodr+1,nodr,2), & anp2(0:nodr+1,nodr,2),pnp2(0:nodr+1,nodr,2),anmie(2,2,nodr), & a qe=0.d0 qa=0.d0 qs=0.d0 do n=1,nodr a=anmie(1,1,n)*anmie(2,2,n)-anmie(1,2,n)*anmie(1,2,n) do p=1,2 do q=1,2 aninv(p,q)=(-1)**(p+q)*anmie(3-p,3-q,n)/a enddo aninv(p,p)=aninv(p,p)+1.d0 enddo do p=1,2 ! babs=-(1.d0/anmie(p,n)+1.d0) do m=-n,-1 qe=qe-(anp1(n+1,-m,p)*conjg(pnp2(n+1,-m,p)) & + anp2(n+1,-m,p)*conjg(pnp1(n+1,-m,p)))*.5d0 qs=qs+anp1(n+1,-m,p)*conjg(anp2(n+1,-m,p)) do q=1,2 qa=qa-conjg(anp1(n+1,-m,p))*aninv(p,q)*anp2(n+1,-m,q) enddo enddo do m=0,n qe=qe-(anp1(m,n,p)*conjg(pnp2(m,n,p)) & +anp2(m,n,p)*conjg(pnp1(m,n,p)))*.5d0 qs=qs+anp1(m,n,p)*conjg(anp2(m,n,p)) do q=1,2 qa=qa-conjg(anp1(m,n,p))*aninv(p,q)*anp2(m,n,q) enddo enddo enddo enddo end subroutine qeffcalc ! ! scattering amplitude sa and matrix sm calculation ! ! original: 15 January 2011 ! revised: 21 February 2011: S11 normalization changed ! subroutine scatteringmatrix(amn0,nodrt,xv,ct,phi,sa,sm) use specialfuncs use numconstants implicit none integer :: nodrt,m,n,p,m1,n1,i,j real(8) :: xv,ct,phi,sm(4,4),tau(0:nodrt+1,nodrt,2),cphi,sphi,qsca complex(8) :: amn0(0:nodrt+1,nodrt,2,2),sa(4),ephi,ephim(-nodrt:nodrt), & ci,cin,a,b,sp(4,4) data ci/(0.d0,1.d0)/ call taufunc(ct,nodrt,tau) cphi=cos(phi) sphi=sin(phi) ephi=dcmplx(cphi,sphi) call ephicoef(ephi,nodrt,ephim) sa=(0.d0,0.d0) qsca=0.d0 do n=1,nodrt cin=(-ci)**n do m=-n,n if(m.le.-1) then m1=n+1 n1=-m else m1=m n1=n endif do p=1,2 qsca=qsca+amn0(m1,n1,p,1)*dconjg(amn0(m1,n1,p,1)) & + amn0(m1,n1,p,2)*dconjg(amn0(m1,n1,p,2)) a=amn0(m1,n1,p,1)*cphi+amn0(m1,n1,p,2)*sphi b=amn0(m1,n1,p,1)*sphi-amn0(m1,n1,p,2)*cphi sa(1)=sa(1)+cin*tau(m1,n1,3-p)*b*ephim(m) sa(2)=sa(2)+ci*cin*tau(m1,n1,p)*a*ephim(m) sa(3)=sa(3)+ci*cin*tau(m1,n1,p)*b*ephim(m) sa(4)=sa(4)+cin*tau(m1,n1,3-p)*a*ephim(m) enddo enddo enddo qsca=qsca*2.d0 do i=1,4 do j=1,4 sp(i,j)=sa(i)*dconjg(sa(j))*16.d0/qsca enddo enddo sm(1,1)=sp(1,1)+sp(2,2)+sp(3,3)+sp(4,4) sm(1,2)=-sp(1,1)+sp(2,2)-sp(3,3)+sp(4,4) sm(2,1)=-sp(1,1)+sp(2,2)+sp(3,3)-sp(4,4) sm(2,2)=sp(1,1)+sp(2,2)-sp(3,3)-sp(4,4) sm(3,3)=2.*(sp(1,2)+sp(3,4)) sm(3,4)=-2.*dimag(sp(1,2)+sp(3,4)) sm(4,3)=2.*dimag(sp(1,2)-sp(3,4)) sm(4,4)=2.*(sp(1,2)-sp(3,4)) sm(1,3)=2.*(sp(2,3)+sp(1,4)) sm(3,1)=2.*(sp(2,4)+sp(1,3)) sm(1,4)=2.*dimag(sp(2,3)-sp(1,4)) sm(4,1)=-2.*dimag(sp(2,4)+sp(1,3)) sm(2,3)=2.*(sp(2,3)-sp(1,4)) sm(3,2)=2.*(sp(2,4)-sp(1,3)) sm(2,4)=2.*dimag(sp(2,3)+sp(1,4)) sm(4,2)=-2.*dimag(sp(2,4)-sp(1,3)) ! do i=1,4 ! do j=1,4 ! if(i.ne.1.or.j.ne.1) then ! sm(i,j)=sm(i,j)/sm(1,1) ! endif ! enddo ! enddo end subroutine scatteringmatrix ! c c ! c subroutine scatexp(amn0,nodrt,nodrg,gmn) computes the expansion coefficients c ! c for the spherical harmonic expansion of the scattering phase function from c ! c the scattering coefficients amn0. For a complete expansion, the max. order c ! c of the phase function expansion (nodrg) will be 2*nodrt, where nodrt is c ! c the max. order of the scattered field expansion. In this code nodrg is c ! c typically set to 1, so that the subroutine returns the first moments c ! c of the phase function; gmn(1) and gmn(2). c ! c c ! c The expansion coefficients are normalized so that gmn(0)=1 c ! c c ! c gmn(1)/3 is the asymmetry parameter. c ! c c subroutine s11expansion(amn0,nodrt,mmax,nodrg,gmn) use specialfuncs use numconstants implicit none integer :: nodrt,m,n,p,ma,na,mmax,nodrg,w,w1,w2,u,uw,ww1, & l1,l2,ka,la,k,l,q,ik real(8) :: vc1(0:nodrt*2+1),vc2(0:nodrt*2+1),g0 complex(8) :: amn0(0:nodrt+1,nodrt,2,2),gmn(0:nodrg*(nodrg+3)/2), & a(2,2),c,c2 gmn=(0.d0,0.d0) do n=1,nodrt l1=max(1,n-nodrg) l2=min(nodrt,n+nodrg) do l=l1,l2 c=sqrt(dble((n+n+1)*(l+l+1)))*dcmplx(0.d0,1.d0)**(l-n) w2=min(n+l,nodrg) call vcfunc(-1,l,1,n,vc2) do m=-n,n if(m.le.-1) then ma=n+1 na=-m else ma=m na=n endif do k=-l,min(l,m) if(k.le.-1) then ka=l+1 la=-k else ka=k la=l endif u=m-k if(u.le.mmax) then ik=(-1)**k c2=ik*c do p=1,2 do q=1,2 a(p,q)=c2*(amn0(ma,na,p,1)*conjg(amn0(ka,la,q,1)) & +amn0(ma,na,p,2)*conjg(amn0(ka,la,q,2))) enddo enddo w1=max(abs(n-l),abs(u)) w2=min(n+l,nodrg) call vcfunc(-k,l,m,n,vc1) do w=w1,w2 uw=(w*(w+1))/2+u do p=1,2 if(mod(n+l+w,2).eq.0) then q=p else q=3-p endif gmn(uw)=gmn(uw)-vc1(w)*vc2(w)*a(p,q) enddo enddo endif enddo enddo enddo enddo g0=dble(gmn(0)) gmn(0)=1.d0 do w=1,nodrg ww1=(w*(w+1))/2 gmn(ww1)=dcmplx(dble(gmn(ww1)),0.d0)/g0 do u=1,min(mmax,w) uw=ww1+u gmn(uw)=(-1)**u*2.d0*gmn(uw)/g0 enddo enddo end subroutine s11expansion ! ! calculate azimuth--averaged scattering matrix from expansion, for cos(theta) = ct ! ! ! original: 15 January 2011 ! revised: 21 February 2011: changed normalization on S11 ! subroutine fosmcalc(ntot,s00,s02,sp22,sm22,ct,sm) use numconstants use specialfuncs integer :: ntot,w,i,j,ww1 real(8) :: s00(4,4,0:ntot*2),s02(4,4,0:ntot*2),sp22(4,4,0:ntot*2),sm22(4,4,0:ntot*2), & sm(4,4),dc(-2:2,0:2*ntot*(2*ntot+2)),ct call rotcoef(ct,2,2*ntot,dc) sm=0.d0 do w=0,2*ntot ww1=w*(w+1) sm(:,:)=sm(:,:)+s00(:,:,w)*dc(0,ww1)+s02(:,:,w)*dc(0,ww1+2) & +sp22(:,:,w)*dc(2,ww1+2)+sm22(:,:,w)*dc(-2,ww1+2) enddo sm=sm/s00(1,1,0) ! do i=1,4 ! do j=1,4 ! if(i.ne.1.or.j.ne.1) then ! sm(i,j)=sm(i,j)/sm(1,1) ! endif ! enddo ! enddo end subroutine fosmcalc ! ! determine the generalized spherical function expansion for the azimuth-averaged scattering matrix ! corresponding to the target-based scattering field expansion of amnp. ! ! ! original: 15 January 2011 ! revised: 21 February 2011: fixed flush call. ! subroutine fosmexpansion(ntot,amnp,s00,s02,sp22,sm22) use mpidefs use mpidata use specialfuncs use numconstants use spheredata integer :: ntot,n,p,m,l,wmin,wmax,m1m,q,m1mq,m1mnpl,w,m1w,fe,fo,i,j,wtot integer :: rank,numprocs,nl,nsend,runprintunit integer, allocatable :: nlindex(:),nlnum(:) real(8) :: s00(4,4,0:ntot*2),s02(4,4,0:ntot*2),sp22(4,4,0:ntot*2),sm22(4,4,0:ntot*2), & cm1p1(0:ntot*2),cm1m1(0:ntot*2),cmmpm(0:ntot*2),cmmm2pm(0:ntot*2), & cmmp2pm(0:ntot*2),sum,nlperproc complex(8) :: amnp(0:ntot+1,ntot,2,2),a1(-ntot-2:ntot+2,ntot,2),a2(-ntot-2:ntot+2,ntot,2), & ci,fnl,a1122,a2112,a1p2,a1m2 data ci/(0.d0,1.d0)/ call init(2*ntot) call getrunparameters(run_print_unit=runprintunit) call ms_mpi(mpi_command='rank',mpi_rank=rank) call ms_mpi(mpi_command='size',mpi_size=numprocs) allocate(nlindex(0:numprocs-1),nlnum(0:numprocs-1)) nlperproc=dble(ntot*ntot)/dble(numprocs) sum=0. do i=0,numprocs-1 nlindex(i)=floor(sum) sum=sum+nlperproc enddo do i=0,numprocs-2 nlnum(i)=nlindex(i+1)-nlindex(i) enddo nlnum(numprocs-1)=ntot*ntot-nlindex(numprocs-1) if(rank.eq.0) then write(runprintunit,'('' SM calc, orders per processor:'',f10.4)') nlperproc call flush(runprintunit) endif a1=(0.d0,0.d0) a2=(0.d0,0.d0) s00=0.d0 s02=0.d0 sp22=0.d0 sm22=0.d0 wtot=ntot+ntot do n=1,ntot do p=1,2 do m=-n,-1 a1(m,n,p)=amnp(n+1,-m,p,1) a2(m,n,p)=amnp(n+1,-m,p,2) enddo do m=0,n a1(m,n,p)=amnp(m,n,p,1) a2(m,n,p)=amnp(m,n,p,2) enddo enddo enddo do nl=nlindex(rank)+1,nlindex(rank)+nlnum(rank) n=floor((nl-1)/dble(ntot))+1 l=mod(nl-1,ntot)+1 wmin=abs(n-l) wmax=n+l fnl=sqrt(dble((n+n+1)*(l+l+1)))*ci**(l-n) call vcfunc(-1,n,1,l,cm1p1) call vcfunc(-1,n,-1,l,cm1m1) do m=-min(n,l+2),min(n,l+2) m1m=(-1)**m if(abs(m).le.l) then call vcfunc(-m,n,m,l,cmmpm) else cmmpm=0.d0 endif if(abs(-2+m).le.l) then call vcfunc(-m,n,-2+m,l,cmmm2pm) else cmmm2pm=0.d0 endif if(abs(2+m).le.l) then call vcfunc(-m,n,2+m,l,cmmp2pm) else cmmp2pm=0.d0 endif do p=1,2 do q=1,2 m1mq=(-1)**(m+q) m1mnpl=(-1)**(m+n+p+l) a1122=(a1(m,n,p)*conjg(a1(m,l,q)) + a2(m,n,p)*conjg(a2(m,l,q))) a2112=(a2(m,n,p)*conjg(a1(m,l,q)) - a1(m,n,p)*conjg(a2(m,l,q))) a1p2=(a1(m,n,p)+ci*a2(m,n,p))*conjg(a1(m-2,l,q)-ci*a2(m-2,l,q)) a1m2=(a1(m,n,p)-ci*a2(m,n,p))*conjg(a1(m+2,l,q)+ci*a2(m+2,l,q)) do w=wmin,wmax m1w=(-1)**w if(mod(n+l+w+p+q,2).eq.0) then fe=1 fo=0 else fe=0 fo=1 endif s00(1,1,w) = s00(1,1,w)-(m1m*fe*fnl*a1122*cm1p1(w)*cmmpm(w))/2. s00(3,2,w) = s00(3,2,w)+ (ci/2.*m1m*fnl*fo*a1122*cm1p1(w)*cmmpm(w)) s00(4,2,w) = s00(4,2,w)+ dimag(-ci/2.*m1m*fnl*fo*a1122*cm1p1(w)*cmmpm(w)) s00(1,4,w) = s00(1,4,w)+ dimag(m1m*fe*fnl*(-a2112)*cm1p1(w)*cmmpm(w))/2. s00(2,3,w) = s00(2,3,w)+ (m1m*fe*fnl*(-a2112)*cm1p1(w)*cmmpm(w))/2. s00(4,3,w) = s00(4,3,w)+ dimag(ci/2.*m1m*fnl*fo*a2112*cm1p1(w)*cmmpm(w)) s00(4,4,w) = s00(4,4,w)+ (ci/2.*m1m*fnl*fo*a2112*cm1p1(w)*cmmpm(w)) if(w.lt.2) cycle s02(2,1,w) = s02(2,1,w)-(m1mq*a1122*fe*fnl*cm1m1(w)*cmmpm(w))/2. s02(3,1,w) = s02(3,1,w)+ (-ci/2.*m1mq*a1122*fnl*fo*cm1m1(w)*cmmpm(w)) s02(4,1,w) = s02(4,1,w)+ dimag(ci/2.*m1mq*a1122*fnl*fo*cm1m1(w)*cmmpm(w)) s02(1,3,w) = s02(1,3,w)-(m1mq*a2112*fe*fnl*cm1m1(w)*cmmpm(w))/2. s02(2,4,w) = s02(2,4,w)-dimag(m1mq*a2112*fe*fnl*cm1m1(w)*cmmpm(w))/2. s02(3,3,w) = s02(3,3,w)+ (ci/2.*m1mq*a2112*fnl*fo*cm1m1(w)*cmmpm(w)) s02(3,4,w) = s02(3,4,w)+ dimag(-ci/2.*m1mq*a2112*fnl*fo*cm1m1(w)*cmmpm(w)) s02(1,2,w) = s02(1,2,w)-(m1m*a1p2*fe*fnl*cm1p1(w)*cmmm2pm(w))/4. s02(1,3,w) = s02(1,3,w)+ (-ci/4.*m1m*a1p2*fe*fnl*cm1p1(w)*cmmm2pm(w)) s02(2,4,w) = s02(2,4,w)+ dimag(-ci/4.*m1m*a1p2*fe*fnl*cm1p1(w)*cmmm2pm(w)) s02(3,1,w) = s02(3,1,w)+ (ci/4.*m1m*a1p2*fnl*fo*cm1p1(w)*cmmm2pm(w)) s02(4,1,w) = s02(4,1,w)+ dimag(-ci/4.*m1m*a1p2*fnl*fo*cm1p1(w)*cmmm2pm(w)) s02(4,3,w) = s02(4,3,w)+ dimag(m1m*a1p2*fnl*fo*cm1p1(w)*cmmm2pm(w))/4. s02(4,4,w) = s02(4,4,w)+ (m1m*a1p2*fnl*fo*cm1p1(w)*cmmm2pm(w))/4. sm22(1,4,w) = sm22(1,4,w)+ dimag(-ci/8.*m1mnpl*m1w*a1p2*fnl*cm1m1(w)*cmmm2pm(w)) sm22(2,2,w) = sm22(2,2,w)-(m1mnpl*m1w*a1p2*fnl*cm1m1(w)*cmmm2pm(w))/8. sm22(2,3,w) = sm22(2,3,w)+ (-ci/8.*m1mnpl*m1w*a1p2*fnl*cm1m1(w)*cmmm2pm(w)) sm22(3,2,w) = sm22(3,2,w)+ (ci/8.*m1mnpl*m1w*a1p2*fnl*cm1m1(w)*cmmm2pm(w)) sm22(3,3,w) = sm22(3,3,w)-(m1mnpl*m1w*a1p2*fnl*cm1m1(w)*cmmm2pm(w))/8. sm22(3,4,w) = sm22(3,4,w)+ dimag(m1mnpl*m1w*a1p2*fnl*cm1m1(w)*cmmm2pm(w))/8. sm22(4,2,w) = sm22(4,2,w)+ dimag(-ci/8.*m1mnpl*m1w*a1p2*fnl*cm1m1(w)*cmmm2pm(w)) sp22(1,4,w) = sp22(1,4,w)+ dimag(-ci/8.*m1mq*a1p2*fnl*cm1m1(w)*cmmm2pm(w)) sp22(2,2,w) = sp22(2,2,w)-(m1mq*a1p2*fnl*cm1m1(w)*cmmm2pm(w))/8. sp22(2,3,w) = sp22(2,3,w)+ (-ci/8.*m1mq*a1p2*fnl*cm1m1(w)*cmmm2pm(w)) sp22(3,2,w) = sp22(3,2,w)+ (-ci/8.*m1mq*a1p2*fnl*cm1m1(w)*cmmm2pm(w)) sp22(3,3,w) = sp22(3,3,w)+ (m1mq*a1p2*fnl*cm1m1(w)*cmmm2pm(w))/8. sp22(3,4,w) = sp22(3,4,w)-dimag(m1mq*a1p2*fnl*cm1m1(w)*cmmm2pm(w))/8. sp22(4,2,w) = sp22(4,2,w)+ dimag(ci/8.*m1mq*a1p2*fnl*cm1m1(w)*cmmm2pm(w)) s02(1,2,w) = s02(1,2,w)-(m1m*a1m2*fe*fnl*cm1p1(w)*cmmp2pm(w))/4. s02(1,3,w) = s02(1,3,w)+ (ci/4.*m1m*a1m2*fe*fnl*cm1p1(w)*cmmp2pm(w)) s02(2,4,w) = s02(2,4,w)+ dimag(ci/4.*m1m*a1m2*fe*fnl*cm1p1(w)*cmmp2pm(w)) s02(3,1,w) = s02(3,1,w)+ (ci/4.*m1m*a1m2*fnl*fo*cm1p1(w)*cmmp2pm(w)) s02(4,1,w) = s02(4,1,w)+ dimag(-ci/4.*m1m*a1m2*fnl*fo*cm1p1(w)*cmmp2pm(w)) s02(4,3,w) = s02(4,3,w)-dimag(m1m*a1m2*fnl*fo*cm1p1(w)*cmmp2pm(w))/4. s02(4,4,w) = s02(4,4,w)-(m1m*a1m2*fnl*fo*cm1p1(w)*cmmp2pm(w))/4. sm22(1,4,w) = sm22(1,4,w)+ dimag(ci/8.*m1mq*a1m2*fnl*cm1m1(w)*cmmp2pm(w)) sm22(2,2,w) = sm22(2,2,w)-(m1mq*a1m2*fnl*cm1m1(w)*cmmp2pm(w))/8. sm22(2,3,w) = sm22(2,3,w)+ (ci/8.*m1mq*a1m2*fnl*cm1m1(w)*cmmp2pm(w)) sm22(3,2,w) = sm22(3,2,w)+ (-ci/8.*m1mq*a1m2*fnl*cm1m1(w)*cmmp2pm(w)) sm22(3,3,w) = sm22(3,3,w)-(m1mq*a1m2*fnl*cm1m1(w)*cmmp2pm(w))/8. sm22(3,4,w) = sm22(3,4,w)+ dimag(m1mq*a1m2*fnl*cm1m1(w)*cmmp2pm(w))/8. sm22(4,2,w) = sm22(4,2,w)+ dimag(ci/8.*m1mq*a1m2*fnl*cm1m1(w)*cmmp2pm(w)) sp22(1,4,w) = sp22(1,4,w)+ dimag(ci/8.*m1mnpl*m1w*a1m2*fnl*cm1m1(w)*cmmp2pm(w)) sp22(2,2,w) = sp22(2,2,w)-(m1mnpl*m1w*a1m2*fnl*cm1m1(w)*cmmp2pm(w))/8. sp22(2,3,w) = sp22(2,3,w)+ (ci/8.*m1mnpl*m1w*a1m2*fnl*cm1m1(w)*cmmp2pm(w)) sp22(3,2,w) = sp22(3,2,w)+ (ci/8.*m1mnpl*m1w*a1m2*fnl*cm1m1(w)*cmmp2pm(w)) sp22(3,3,w) = sp22(3,3,w)+ (m1mnpl*m1w*a1m2*fnl*cm1m1(w)*cmmp2pm(w))/8. sp22(3,4,w) = sp22(3,4,w)-dimag(m1mnpl*m1w*a1m2*fnl*cm1m1(w)*cmmp2pm(w))/8. sp22(4,2,w) = sp22(4,2,w)+ dimag(-ci/8.*m1mnpl*m1w*a1m2*fnl*cm1m1(w)*cmmp2pm(w)) enddo enddo enddo enddo enddo call ms_mpi(mpi_command='barrier') nsend=4*4*(2*ntot+1) call ms_mpi(mpi_command='allreduce',mpi_recv_buf_dp=s00,& mpi_number=nsend,mpi_operation=ms_mpi_sum) call ms_mpi(mpi_command='allreduce',mpi_recv_buf_dp=s02,& mpi_number=nsend,mpi_operation=ms_mpi_sum) call ms_mpi(mpi_command='allreduce',mpi_recv_buf_dp=sp22,& mpi_number=nsend,mpi_operation=ms_mpi_sum) call ms_mpi(mpi_command='allreduce',mpi_recv_buf_dp=sm22,& mpi_number=nsend,mpi_operation=ms_mpi_sum) ! ! a patch ! do i=3,4 do j=1,i s00(j,i,0:wtot)=-s00(j,i,0:wtot) s02(j,i,0:wtot)=-s02(j,i,0:wtot) sm22(j,i,0:wtot)=-sm22(j,i,0:wtot) sp22(j,i,0:wtot)=-sp22(j,i,0:wtot) enddo enddo deallocate(nlindex,nlnum) end subroutine fosmexpansion ! ! compute the coefficients for the GSF expansion of the random orientation ! scattering matrix. ! ! ! original: 15 January 2011 ! revised: 21 February 2011: changed normalization on S11 ! subroutine ranorientscatmatrix(xv,nsphere,nodr,nodrw,cbeam,tmatrixfile,& sm,qext,qabs,qsca) use mpidefs use mpidata use intrinsics use specialfuncs use spheredata use numconstants implicit none integer :: nodr,nodrw,nodr2,m,n,p,k,l,q,s,t,v,u,w,nblk,kl,mn,nn1,tn, & lmax,ll1,tvl,ku,k1,ns,ik,ik1,m1,nu,n1s,n1e,nu1,p1,n1max, & in,n1,i,lt,kt,qt,nt,mt,ikm,klm,mnm,nodrt,nsphere, & rank,iunit,numprocs real(8) :: sm(4,4,0:nodrw),fl,vc(0:4*nodr+2),xv,fl2,fc1,fc2,fc3,fc4, & cbeam,gbn,qext(nsphere),qabs(nsphere),qsca(nsphere),qel, & qal,qsl,fc(4),time1,time2,qsca0 complex(8) :: ci,cin,a complex(8) :: aw(0:2,-1:1,0:nodrw),bw(0:2,-1:1,0:nodrw),cw(0:nodrw), & dw(0:nodrw),pp(nodr,2,2), & bm(2,nodr*(nodr+2),2),am(2,nodr+1,2),fm(3,nodr,2,nodr,2) complex(8), allocatable :: dm(:,:,:,:,:,:) complex(4), allocatable :: tc(:,:,:,:) integer :: nblkw,wv,sizedm,ierr,sizetm,nsend integer, allocatable :: windex(:),vindex(:),wvindex(:),wvnum(:) real(8) :: wvperproc,sum character*30 :: tmatrixfile data ci/(0.d0,1.d0)/ call ms_mpi(mpi_command='rank',mpi_rank=rank) call ms_mpi(mpi_command='size',mpi_size=numprocs) call getrunparameters(run_print_unit=iunit) if(rank.eq.0) time1=mytime() ! ! read the T matrix from the file ! if(rank.eq.0) then open(3,file=tmatrixfile) read(3,*) nodrt endif nodrt=nodr nblk=nodr*(nodr+2) sizetm=4*nblk*nblk allocate(tc(2,nblk,2,nblk)) tc=(0.,0.) if(rank.eq.0) then qext=0.d0 qabs=0.d0 qsca=0.d0 do l=1,nodr gbn=dexp(-((dble(l)+.5d0)*cbeam)**2.) do k=-l,l kl=l*(l+1)+k klm=l*(l+1)-k do q=1,2 read(3,*) lt,kt,qt do n=1,l do m=-n,n mn=n*(n+1)+m mnm=n*(n+1)-m read(3,*) nt,mt,fc tc(1,mn,q,kl)=cmplx(fc(1),fc(2)) tc(2,mn,q,kl)=cmplx(fc(3),fc(4)) if(n.lt.l) then ikm=(-1)**(m+k) do p=1,2 tc(q,klm,p,mnm)=tc(p,mn,q,kl)*ikm enddo endif enddo enddo enddo enddo do i=1,nsphere read(3,*) n,qel,qal,qsl qext(i)=qext(i)+qel*gbn*gbn qabs(i)=qabs(i)+qal*gbn*gbn qsca(i)=qsca(i)+qsl*gbn*gbn enddo enddo close(3) endif ! ! send to the other processors ! if(numprocs.gt.1) then call ms_mpi(mpi_command='barrier') call ms_mpi(mpi_command='bcast',mpi_send_buf_dp=qext,mpi_number=nsphere,mpi_rank=0) call ms_mpi(mpi_command='bcast',mpi_send_buf_dp=qabs,mpi_number=nsphere,mpi_rank=0) call ms_mpi(mpi_command='bcast',mpi_send_buf_dp=qsca,mpi_number=nsphere,mpi_rank=0) call ms_mpi(mpi_command='bcast',mpi_send_buf_c=tc,mpi_number=sizetm,mpi_rank=0) endif allocate(dm(-nodr-1:nodr+1,3,nodr,2,nodr,2)) if(rank.eq.0) then time2=mytime()-time1 call timewrite(iunit,' t matrix read time:',time2) time1=mytime() endif nodr2=nodr+nodr nblk=nodr*(nodr+2) dm=(0.d0,0.d0) sizedm=size(dm) call init(nodr2) ! ! compute the GB modified T matrix ! do n=1,nodr gbn=dexp(-((dble(n)+.5d0)*cbeam)**2.) cin=ci**(n+1) pp(n,1,1) =-.5d0*cin*fnr(n+n+1)*gbn pp(n,2,1) =-pp(n,1,1) pp(n,1,2)=-pp(n,1,1) pp(n,2,2)=pp(n,2,1) enddo do n=1,nodr nn1=n*(n+1) do m=-n,n mn=nn1+m do p=1,2 do l=1,nodr do k=-l,l kl=l*(l+1)+k a=tc(p,mn,1,kl) tc(p,mn,1,kl)=tc(p,mn,1,kl)*pp(l,1,1)& +tc(p,mn,2,kl)*pp(l,1,2) tc(p,mn,2,kl)=a*pp(l,2,1)+tc(p,mn,2,kl)*pp(l,2,2) enddo enddo enddo enddo enddo ! ! determine the distribution of work load among the processors ! nblkw=nodr2*(nodr2+2)+1 allocate(windex(nblkw),vindex(nblkw),wvindex(0:numprocs-1),wvnum(0:numprocs-1)) w=0 do n=0,nodr2 do m=-n,n w=w+1 windex(w)=n vindex(w)=m enddo enddo wvperproc=dble(nblkw)/dble(numprocs) sum=0. do i=0,numprocs-1 wvindex(i)=floor(sum) sum=sum+wvperproc enddo do i=0,numprocs-2 wvnum(i)=wvindex(i+1)-wvindex(i) enddo wvnum(numprocs-1)=nblkw-wvindex(numprocs-1) if(rank.eq.0) then write(iunit,'('' d matrix calculation, order+degree per proc.:'',f9.2)') & wvperproc call flush(iunit) endif ! ! the big loop ! do wv=wvindex(rank)+1,wvindex(rank)+wvnum(rank) w=windex(wv) v=vindex(wv) bm=(0.d0,0.d0) do n=1,nodr nn1=n*(n+1) do l=max(1,abs(w-n)),min(nodr,w+n) am(1,l,1)=0.d0 am(1,l,2)=0.d0 am(2,l,1)=0.d0 am(2,l,2)=0.d0 enddo do t=-n,n tn=nn1+t lmax=min(nodr,w+n) call vcfunc(v,w,-t,n,vc) do l=max(1,abs(v-t),abs(n-w)),lmax ll1=l*(l+1) tvl=ll1+t-v do k=1,2 do p=1,2 am(k,l,p)=am(k,l,p)+vc(l)*tc(p,tn,k,tvl) enddo enddo enddo enddo do m=-n,n mn=nn1+m do k=1,2 u=m-(-3+2*k) if(abs(u).le.w) then lmax=min(nodr,w+n) call vcfunc(-u,w,m,n,vc) do l=max(1,abs(w-n)),lmax fl=-(-1)**l*vc(l)/dble(l+l+1) do p=1,2 bm(k,mn,p)=bm(k,mn,p)+am(k,l,p)*fl enddo enddo endif enddo enddo enddo do u=-min(w,nodr+1),min(w,nodr+1) do ku=1,3 if(ku.eq.1) then k=-1 k1=-1 elseif(ku.eq.2) then k=1 k1=1 else k=1 k1=-1 endif m=u+k ns=max(1,abs(m)) ik=(k+1)/2+1 ik1=(k1+1)/2+1 m1=u+k1 do n=ns,nodr nu=n*(n+1)+m n1s=max(1,abs(m1),n-nodrw) n1e=min(nodr,n+nodrw) do n1=n1s,n1e cin=ci**(n-n1) nu1=n1*(n1+1)+m1 fl=-fnr(n+n+1)*fnr(n1+n1+1)*dble(w+w+1) do p=1,2 do p1=1,2 a=bm(ik,nu,p)*cin*fl*conjg(bm(ik1,nu1,p1)) dm(u,ku,n,p,n1,p1)=dm(u,ku,n,p,n1,p1)+a enddo enddo enddo enddo enddo enddo enddo deallocate(tc) call ms_mpi(mpi_command='barrier') call ms_mpi(mpi_command='allreduce',mpi_recv_buf_dc=dm,& mpi_number=sizedm,mpi_operation=ms_mpi_sum) if(rank.eq.0) then time2=mytime()-time1 call timewrite(iunit,' d matrix time:',time2) time1=mytime() endif ! ! compute the expansion coefficients ! aw=0.d0 bw=0.d0 cw=0.d0 dw=0.d0 do w=0,nodrw do n=1,nodr n1s=max(1,abs(n-w)) n1e=min(nodr,n+w) do n1=n1s,n1e do k=1,3 do p=1,2 do p1=1,2 fm(k,n,p,n1,p1)=0. enddo enddo enddo enddo enddo do u=-nodr-1,nodr+1 do k=-1,1,2 m=u+k ik=(k+1)/2+1 ns=max(1,abs(m)) do n=ns,nodr n1max=min(w+n,nodr) call vcfunc(m,n,0,w,vc) do n1=ns,nodr if((n+n1.lt.w).or.(abs(n-n1).gt.w)) cycle fl=-(-1)**n*vc(n1)*fnr(w+w+1)/fnr(n1+n1+1) do p=1,2 do p1=1,2 fm(ik,n,p,n1,p1)=fm(ik,n,p,n1,p1)+dm(u,ik,n,p,n1,p1)*fl enddo enddo enddo enddo enddo if(w.lt.2) cycle m=u+1 m1=u-1 ns=max(1,abs(m)) n1s=max(1,abs(m1)) do n=ns,nodr n1max=min(w+n,nodr) call vcfunc(m,n,-2,w,vc) do n1=n1s,nodr if((n+n1.lt.w).or.(abs(n-n1).gt.w)) cycle fl=-(-1)**n*vc(n1)*fnr(w+w+1)/fnr(n1+n1+1) do p=1,2 do p1=1,2 fm(3,n,p,n1,p1)=fm(3,n,p,n1,p1)+dm(u,3,n,p,n1,p1)*fl enddo enddo enddo enddo enddo do n=1,nodr n1s=max(1,abs(n-w)) n1e=min(nodr,n+w) in=(-1)**n n1max=min(w+n,nodr) call vcfunc(1,n,0,w,vc) do n1=n1s,n1e fl=2.d0*in*vc(n1)*fnr(w+w+1)/fnr(n1+n1+1) i=mod(n+n1-w,2)+1 do p=1,2 p1=(2-i)*p+(i-1)*(3-p) do k=-1,1,2 ik=(k+1)/2+1 aw(0,k,w)=aw(0,k,w)+fm(ik,n,p,n1,p1)*fl bw(0,k,w)=bw(0,k,w)+fm(ik,n,p,n1,3-p1)*fl enddo bw(2,0,w)=bw(2,0,w)+fm(3,n,p,n1,3-p1)*fl aw(2,0,w)=aw(2,0,w)+fm(3,n,p,n1,p1)*fl enddo enddo if(w.lt.2) cycle call vcfunc(1,n,-2,w,vc) do n1=n1s,n1e fl=2.d0*in*vc(n1)*fnr(w+w+1)/fnr(n1+n1+1) i=mod(n+n1-w,2)+1 do p=1,2 p1=(2-i)*p+(i-1)*(3-p) do k=-1,1,2 ik=(k+1)/2+1 aw(2,k,w)=aw(2,k,w)+fm(ik,n,p,n1,p1)*fl*(-1)**p1 bw(2,k,w)=bw(2,k,w)+fm(ik,n,p,n1,3-p1)*fl*(-1)**(3-p1) enddo enddo fl2=2.*(-1)**(n1+w)*vc(n1)*fnr(w+w+1)/fnr(n1+n1+1) do p=1,2 do p1=1,2 cw(w)=cw(w)+fm(3,n,p,n1,p1)*fl*(-1)**p1 dw(w)=dw(w)+fm(3,n,p,n1,p1)*fl2*(-1)**p enddo enddo enddo enddo enddo do w=0,nodrw do k=-1,1 do i=0,2 aw(i,k,w)=aw(i,k,w)*2./xv/xv bw(i,k,w)=bw(i,k,w)*2./xv/xv enddo enddo cw(w)=cw(w)*2./xv/xv dw(w)=dw(w)*2./xv/xv enddo do n=0,nodrw sm(1,1,n)=aw(0,-1,n)+aw(0,1,n) sm(1,2,n)=aw(2,-1,n)+aw(2,1,n) sm(1,3,n)=2.d0*dimag(aw(2,0,n)) sm(1,4,n)=aw(0,1,n)-aw(0,-1,n) sm(2,2,n)=dw(n) sm(2,3,n)=dimag(dw(n)) sm(2,4,n)=aw(2,1,n)-aw(2,-1,n) sm(3,2,n)=dimag(cw(n)) sm(3,3,n)=cw(n) sm(3,4,n)=dimag(bw(2,-1,n)-bw(2,1,n)) sm(4,4,n)=bw(0,1,n)-bw(0,-1,n) enddo ! ! normalization ! qsca0=sm(1,1,0) do n=0,nodrw sm(1,1,n)=sm(1,1,n)/qsca0 sm(1,2,n)=sm(1,2,n)/qsca0 sm(1,3,n)=sm(1,3,n)/qsca0 sm(1,4,n)=sm(1,4,n)/qsca0 sm(2,2,n)=sm(2,2,n)/qsca0 sm(2,3,n)=sm(2,3,n)/qsca0 sm(2,4,n)=sm(2,4,n)/qsca0 sm(3,2,n)=sm(3,2,n)/qsca0 sm(3,3,n)=sm(3,3,n)/qsca0 sm(3,4,n)=sm(3,4,n)/qsca0 sm(4,4,n)=sm(4,4,n)/qsca0 enddo call ms_mpi(mpi_command='barrier') if(rank.eq.0) then time2=mytime()-time1 call timewrite(iunit,' scat matrix coef time:',time2) endif deallocate(windex,vindex,wvindex,wvnum,dm) end subroutine ranorientscatmatrix ! ! calculation of the RO scattering matrix from the GSF expansion ! ! ! original: 15 January 2011 ! revised: 21 February 2011: changed normalization on S11 ! subroutine ranorienscatmatrixcalc(nt,tmin,tmax,iscale,smc,nodrexp,sm) use specialfuncs use numconstants implicit none integer :: nt,iscale,nodrexp,i,j,k,n,nn0,nnp2,nnm2 real(8) :: tmin,tmax,smc(4,4,0:nodrexp),sm(4,4,nt),dc(-2:2,0:nodrexp*(nodrexp+2)), & ct,th,qsca do i=1,nt if(nt.eq.1) then th=tmin else th=tmin+(tmax-tmin)*dble(i-1)/dble(nt-1) endif ct=cos(th*pi/180.d0) ! ! dc is the normalized generalized spherical function ! dc(k,n*(n+1)+m) = ((n-k)!(n+m)!/(n+k)!/(n-m)!)^(1/2) D^k_{mn}, ! where D^k_{mn} is defined in M&M JOSA 96 ! call rotcoef(ct,2,nodrexp,dc) do j=1,4 do k=j,4 sm(j,k,i)=0.d0 enddo enddo do n=0,nodrexp nn0=n*(n+1) nnp2=nn0+2 nnm2=nn0-2 sm(1,1,i)=sm(1,1,i)+dc(0,nn0)*smc(1,1,n) sm(1,4,i)=sm(1,4,i)+dc(0,nn0)*smc(1,4,n) sm(4,4,i)=sm(4,4,i)+dc(0,nn0)*smc(4,4,n) if(n.ge.2) then sm(1,2,i)=sm(1,2,i)+dc(2,nn0)*smc(1,2,n) sm(2,4,i)=sm(2,4,i)+dc(2,nn0)*smc(2,4,n) sm(3,4,i)=sm(3,4,i)+dc(2,nn0)*smc(3,4,n) sm(1,3,i)=sm(1,3,i)+dc(2,nn0)*smc(1,3,n) sm(2,2,i)=sm(2,2,i)+dc(2,nnm2)*smc(2,2,n)+dc(2,nnp2)*smc(3,3,n) sm(2,3,i)=sm(2,3,i)+dc(2,nnp2)*smc(2,3,n)+dc(2,nnp2)*smc(3,2,n) sm(3,3,i)=sm(3,3,i)-dc(2,nnm2)*smc(2,2,n)+dc(2,nnp2)*smc(3,3,n) endif enddo ! ! discontiued scaling option: now done in main program ! ! if(iscale.eq.1) then ! do j=1,4 ! do k=j,4 ! if(j.ne.1.or.k.ne.1) then ! sm(j,k,i)=sm(j,k,i)/sm(1,1,i) ! endif ! enddo ! enddo ! endif ! ! here are the VV and HH differential cross sections ! ! gvv=.25*(sm(1,1)+sm(2,2)-2.*sm(1,2)) ! ghh=.25*(sm(1,1)+sm(2,2)+2.*sm(1,2)) ! enddo return end subroutine ranorienscatmatrixcalc end module scatprops ! ! module nearfield contains local data and subroutines for near field calculation ! ! ! last revised: 15 January 2011 ! 30 March 2011: added optical activity ! module nearfield implicit none integer, private :: axialinc,ndimpw,nodrpwmax integer, private, allocatable :: nblk_nf(:),noff_nf(:) real(8), private :: rplotmax complex(8), allocatable, private :: amnp_nf(:),cmnp_nf(:),pmnp_nf(:), & amn3mp_nf(:),cmn3mp_nf(:),pmn3mp_nf(:) contains ! ! nearfieldspherepart calculates the field at point xg generated by the spheres ! ! ! last revised: 15 January 2011 ! 30 March 2011: added optical activity ! subroutine nearfieldspherepart(xg,nsphere,xsp,rpos,ri,& nodr,neqns,insphere,efield,hfield) use specialfuncs use numconstants implicit none integer :: nsphere,nodr(nsphere),neqns,i,insphere,nblki,n real(8) :: xg(3),xsp(nsphere),rpos(3,nsphere),x(3),r complex(8) :: ri(2,nsphere),vwh(3,neqns),efield(3),hfield(3),ri0,cn1,cn2 complex(8), allocatable :: vwhleft(:,:,:),vwhright(:,:,:) ! ! find if the point is inside a sphere ! insphere=0 do i=1,nsphere x=xg(:)-rpos(:,i) r=sqrt(dot_product(x,x)) if(r.le.xsp(i)) then insphere=i exit endif enddo ! ! do the calculations ! if(insphere.eq.0) then ! ! outside a sphere: field = scattered ! do i=1,nsphere x=xg(:)-rpos(:,i) ri0=(1.d0,0.d0) call vwhcalc(x,ri0,nodr(i),3,vwh(1:3,noff_nf(i)+1:noff_nf(i)+nblk_nf(i))) enddo efield(:)=matmul(vwh(:,1:neqns),amnp_nf(1:neqns)) hfield(:)=matmul(vwh(:,1:neqns),amn3mp_nf(1:neqns))/dcmplx(0.d0,1.d0) else ! ! inside a sphere: field = internal ! i=insphere if(abs(ri(1,i)-ri(2,i)).eq.0) then x=xg(:)-rpos(:,i) call vwhcalc(x,ri(1,i),nodr(i),1,vwh) efield(:)=matmul(vwh(:,1:nblk_nf(i)), & cmnp_nf(noff_nf(i)+1:noff_nf(i)+nblk_nf(i))) hfield(:)=matmul(vwh(:,1:nblk_nf(i)), & cmn3mp_nf(noff_nf(i)+1:noff_nf(i)+nblk_nf(i)))*ri(1,i)/dcmplx(0.d0,1.d0) else nblki=nodr(i)*(nodr(i)+2) allocate(vwhleft(3,2,nblki),vwhright(3,2,nblki)) x=xg(:)-rpos(:,i) call vwhcalc(x,ri(1,i),nodr(i),1,vwhleft) call vwhcalc(x,ri(2,i),nodr(i),1,vwhright) efield=0.d0 hfield=0.d0 do n=1,nblki cn1=cmnp_nf(noff_nf(i)+2*n-1) cn2=cmnp_nf(noff_nf(i)+2*n) efield(:)=efield(:)+(vwhleft(:,1,n)+vwhleft(:,2,n))*cn1 & +(vwhright(:,1,n)-vwhright(:,2,n))*cn2 hfield(:)=hfield(:)+((vwhleft(:,2,n)+vwhleft(:,1,n))*cn1*ri(1,i) & +(vwhright(:,2,n)-vwhright(:,1,n))*cn2*ri(2,i))/dcmplx(0.d0,1.d0) enddo deallocate(vwhleft,vwhright) endif endif end subroutine nearfieldspherepart ! ! nearfieldincidentpart calculates the incident field at point xg using a regular ! vswh expansion ! ! ! last revised: 15 January 2011 ! subroutine nearfieldincidentpart(xg,nodrpw,efield,hfield) use specialfuncs use numconstants implicit none integer :: nblkpw,nodrpw real(8) :: xg(3),r,epspw complex(8) :: vwhpw(3,nodrpw*(nodrpw+2)*2),vwhpwaxial(3,4*nodrpw), & efield(3),hfield(3) ! ! oblique incidence: use the full expansion ! if(axialinc.eq.0) then call vwhcalc(xg,(1.d0,0.d0),nodrpw,1,vwhpw) nblkpw=nodrpw*(nodrpw+2)*2 efield(:)=matmul(vwhpw(:,1:nblkpw),pmnp_nf(1:nblkpw)) hfield(:)=matmul(vwhpw(:,1:nblkpw),pmn3mp_nf(1:nblkpw))/dcmplx(0.d0,1.d0) else ! ! axial incidence: use the shortcut ! call vwhaxialcalc(xg,(1.d0,0.d0),nodrpw,1,vwhpwaxial) nblkpw=4*nodrpw efield(:)=matmul(vwhpwaxial(:,1:nblkpw),pmnp_nf(1:nblkpw)) hfield(:)=matmul(vwhpwaxial(:,1:nblkpw),pmn3mp_nf(1:nblkpw))/dcmplx(0.d0,1.d0) endif end subroutine nearfieldincidentpart ! ! nearfieldincidentcoef generates the reshaped array of incident field coefficients ! ! ! last revised: 15 January 2011 ! subroutine nearfieldincidentcoef(nodrpw,alpha,beta,gamma,cbeam,epspw) use specialfuncs use spheredata use miecoefdata use numconstants implicit none integer :: m,n,p,nn1,mn,mnp,nodrpw real (8) :: alpha,beta,cbeam,gamma,epspw,cgamma,sgamma complex(8), allocatable :: pmnp0(:,:,:,:) allocate(pmnp0(0:nodrpw+1,nodrpw,2,2)) if(allocated(pmnp_nf)) deallocate(pmnp_nf,pmn3mp_nf) if(beta.ne.0.d0) then axialinc=0 ndimpw=2*nodrpw*(nodrpw+2) else axialinc=1 ndimpw=4*nodrpw endif allocate(pmnp_nf(ndimpw),pmn3mp_nf(ndimpw)) if(cbeam.eq.0.d0) then call planewavecoef(alpha,beta,nodrpw,pmnp0) else call gaussianbeamcoef(alpha,beta,cbeam,nodrpw,pmnp0) endif cgamma=cos(gamma) sgamma=sin(gamma) if(axialinc.eq.0) then do n=1,nodrpw nn1=n*(n+1) do p=1,2 do m=-n,-1 mn=nn1+m mnp=2*(mn-1)+p pmnp_nf(mnp)=pmnp0(n+1,-m,p,1)*cgamma+pmnp0(n+1,-m,p,2)*sgamma pmn3mp_nf(mnp)=pmnp0(n+1,-m,3-p,1)*cgamma+pmnp0(n+1,-m,3-p,2)*sgamma enddo do m=0,n mn=nn1+m mnp=2*(mn-1)+p pmnp_nf(mnp)=pmnp0(m,n,p,1)*cgamma+pmnp0(m,n,p,2)*sgamma pmn3mp_nf(mnp)=pmnp0(m,n,3-p,1)*cgamma+pmnp0(m,n,3-p,2)*sgamma enddo enddo enddo else do n=1,nodrpw do p=1,2 mnp=4*(n-1)+p pmnp_nf(mnp)=pmnp0(n+1,1,p,1)*cgamma+pmnp0(n+1,1,p,2)*sgamma pmn3mp_nf(mnp)=pmnp0(n+1,1,3-p,1)*cgamma+pmnp0(n+1,1,3-p,2)*sgamma pmnp_nf(mnp+2)=pmnp0(1,n,p,1)*cgamma+pmnp0(1,n,p,2)*sgamma pmn3mp_nf(mnp+2)=pmnp0(1,n,3-p,1)*cgamma+pmnp0(1,n,3-p,2)*sgamma enddo enddo endif deallocate(pmnp0) end subroutine nearfieldincidentcoef ! ! nearfieldpointcalc: if newcalc = 1, generates the reshaped incident, scattered, and ! internal field coefficients, and returns with newcalc=0 ! if newcalc = 0, generates the field at point xg ! ! ! last revised: 15 January 2011 ! 30 March 2011: added optical activity ! subroutine nearfieldpointcalc(neqns,nsphere,nodr,alpha,beta,cbeam,xsp,rpos,ri,amnp, & gamma,epspw,xg,newcalc,efield,hfield) use specialfuncs use spheredata use miecoefdata use numconstants implicit none integer :: nsphere,neqns,nodr(nsphere),i,j,k,m,n,p,nn1,mn,nodrpw,newcalc, & insphere real (8) :: alpha,beta,cbeam,xsp(nsphere),rpos(3,nsphere),xg(3),xgp(3), & gamma,epspw,rplot,cgamma,sgamma complex(8) :: amnp(neqns,2),ri(2,nsphere),efield(3),hfield(3),einc(3),hinc(3),ri0, & ct1,ct2 complex(8), allocatable :: pmnp0(:,:,:,:),cnmie(:,:,:),amnpt(:,:),cmnpt(:,:) ! ! initialization operations: newcalc=1 ! if(newcalc.eq.1) then if(allocated(amnp_nf)) deallocate(amnp_nf,cmnp_nf,amn3mp_nf,cmn3mp_nf, & noff_nf,nblk_nf) allocate(amnp_nf(neqns),cmnp_nf(neqns),amn3mp_nf(neqns),cmn3mp_nf(neqns), & noff_nf(nsphere),nblk_nf(nsphere)) noff_nf(1)=0 do i=1,nsphere nblk_nf(i)=2*nodr(i)*(nodr(i)+2) if(i.lt.nsphere) noff_nf(i+1)=noff_nf(i)+nblk_nf(i) enddo cgamma=cos(gamma) sgamma=sin(gamma) do i=1,nsphere ri0=2.d0/(1.d0/ri(1,i)+1.d0/ri(2,i)) allocate(pmnp0(0:nodr(i)+1,nodr(i),2,2),cnmie(2,2,nodr(i)),& amnpt(2,nodr(i)*(nodr(i)+2)),cmnpt(2,nodr(i)*(nodr(i)+2))) call getmiedata(which_sphere=i,sphere_int_mie_coefficients=cnmie) do p=1,2 pmnp0(0:nodr(i)+1,1:nodr(i),1:2,p) & =reshape(amnp(noff_nf(i)+1:noff_nf(i)+nblk_nf(i),p),(/nodr(i)+2,nodr(i),2/)) enddo if(abs(ri(1,i)-ri(2,i)).gt.1.d-10) then do n=1,nodr(i) nn1=n*(n+1) do p=1,2 do m=-n,-1 mn=nn1+m amnpt(p,mn)=pmnp0(n+1,-m,p,1)*cgamma+pmnp0(n+1,-m,p,2)*sgamma enddo do m=0,n mn=nn1+m amnpt(p,mn)=pmnp0(m,n,p,1)*cgamma+pmnp0(m,n,p,2)*sgamma enddo enddo do m=-n,n mn=nn1+m ct1=amnpt(1,mn)*cnmie(1,1,n)+amnpt(2,mn)*cnmie(1,2,n) ct2=amnpt(1,mn)*cnmie(2,1,n)+amnpt(2,mn)*cnmie(2,2,n) cmnpt(1,mn)=ct1 cmnpt(2,mn)=-(0.d0,1.d0)/ri0*ct2 enddo enddo else do n=1,nodr(i) nn1=n*(n+1) do p=1,2 do m=-n,-1 mn=nn1+m amnpt(p,mn)=pmnp0(n+1,-m,p,1)*cgamma+pmnp0(n+1,-m,p,2)*sgamma cmnpt(p,mn)=amnpt(p,mn)*cnmie(p,p,n) enddo do m=0,n mn=nn1+m amnpt(p,mn)=pmnp0(m,n,p,1)*cgamma+pmnp0(m,n,p,2)*sgamma cmnpt(p,mn)=amnpt(p,mn)*cnmie(p,p,n) enddo enddo enddo endif amnp_nf(noff_nf(i)+1:noff_nf(i)+nblk_nf(i))= & reshape(amnpt(1:2,1:nodr(i)*(nodr(i)+2)),(/nblk_nf(i)/)) cmnp_nf(noff_nf(i)+1:noff_nf(i)+nblk_nf(i))= & reshape(cmnpt(1:2,1:nodr(i)*(nodr(i)+2)),(/nblk_nf(i)/)) amn3mp_nf(noff_nf(i)+1:noff_nf(i)+nblk_nf(i))= & reshape(amnpt(2:1:-1,1:nodr(i)*(nodr(i)+2)),(/nblk_nf(i)/)) cmn3mp_nf(noff_nf(i)+1:noff_nf(i)+nblk_nf(i))= & reshape(cmnpt(2:1:-1,1:nodr(i)*(nodr(i)+2)),(/nblk_nf(i)/)) deallocate(pmnp0,cnmie,amnpt,cmnpt) enddo rplot=sqrt(dot_product(xg,xg)) rplotmax=rplot call planewavetruncationorder(rplot,epspw,nodrpw) nodrpwmax=nodrpw call nearfieldincidentcoef(nodrpw,alpha,beta,gamma,cbeam,epspw) newcalc=0 return endif ! ! point calculation operations: newcalc=0 ! first determine the required order of the incident field expansion, and recalculate ! field coefficients, if necessary ! rplot=sqrt(dot_product(xg,xg)) rplotmax=max(rplot,rplotmax) call planewavetruncationorder(rplot,epspw,nodrpw) if(nodrpw.gt.nodrpwmax) then nodrpwmax=nodrpw call nearfieldincidentcoef(nodrpw,alpha,beta,gamma,cbeam,epspw) endif efield=0.d0 hfield=0.d0 ! ! calculate the sphere contribution to the field ! call nearfieldspherepart(xg,nsphere,xsp,rpos,ri,& nodr,neqns,insphere,efield,hfield) ! ! if the point is external to the spheres, calculate the incident field ! if(insphere.eq.0) then call nearfieldincidentpart(xg,nodrpw,einc,hinc) efield=efield+einc hfield=hfield+hinc endif end subroutine nearfieldpointcalc ! ! nearfieldgridcalc is an MPI--enabled subroutine for calculating field points on a ! rectangular grid. Writes the data to nfoutunit. ! ! ! last revised: 15 January 2011 ! 30 March 2011: added optical activity ! changed so that input positions are defined relative to sphere position file origin, and ! not the gb focal point. ! subroutine nearfieldgridcalc(neqns,nsphere,nodr,alpha,beta,cbeam,xsp,rpos,ri,amnp, & nfplane,nfplanepos0,nfplanevert0,gbfocus,deltax,gamma,nfoutunit,epspw, & nfoutdata,runprintunit) use mpidefs use mpidata use intrinsics use specialfuncs use spheredata use miecoefdata use numconstants implicit none integer :: nsphere,neqns,nodr(nsphere),nfplane,runprintunit,npoints1,npoints2, & npoints,i,j,k,np23,gcoord(3),rank,numprocs,nrowperproc,nrowrem, & npoints1by5,plotincfield,nsp,nfoutunit,nfoutdata,newcalc real (8) :: alpha,beta,cbeam,xsp(nsphere),rpos(3,nsphere),nfplanepos,& nfplanevert(2,2),frowperproc,rowsum,xg(3),xgp(3),deltax,gamma,epspw, & time1,time2,xplot(3,nsphere),xi0,ri0,esquare,xgpmax(3),rplot, & gbfocus(3),nfplanepos0,nfplanevert0(2,2) complex(8) :: amnp(neqns,2),ri(2,nsphere),efield(3),hfield(3) integer, allocatable :: efindex(:),efnum(:) complex(8), allocatable :: efieldrow(:,:),efieldrowt(:,:),hfieldrow(:,:),hfieldrowt(:,:) call ms_mpi(mpi_command='size',mpi_size=numprocs) call ms_mpi(mpi_command='rank',mpi_rank=rank) ! ! determine the plane ! if(nfplane.eq.1) then gcoord=(/2,3,1/) elseif(nfplane.eq.2) then gcoord=(/3,1,2/) else gcoord=(/1,2,3/) endif ! ! shift the coordinates to gb focal origin ! nfplanevert(1,1)=nfplanevert0(1,1)-gbfocus(gcoord(1)) nfplanevert(1,2)=nfplanevert0(1,2)-gbfocus(gcoord(1)) nfplanevert(2,1)=nfplanevert0(2,1)-gbfocus(gcoord(2)) nfplanevert(2,2)=nfplanevert0(2,2)-gbfocus(gcoord(2)) nfplanepos=nfplanepos0-gbfocus(gcoord(3)) xg(gcoord(3))=nfplanepos ! ! determine the number of points ! npoints1=nint((nfplanevert(1,2)-nfplanevert(1,1))/deltax)+1 npoints2=nint((nfplanevert(2,2)-nfplanevert(2,1))/deltax)+1 npoints=npoints1*npoints2 ! ! find the maximum point-to-target origin distance and initialize the field calculation ! xgp(3)=nfplanepos rplotmax=0.d0 xgpmax=0.d0 do i=1,npoints1 xgp(1)=nfplanevert(1,1)+deltax*dble(i-1) do j=1,npoints2 xgp(2)=nfplanevert(2,1)+deltax*dble(j-1) rplot=sqrt(dot_product(xgp,xgp)) if(rplot.gt.rplotmax) then rplotmax=rplot xgpmax=xgp endif enddo enddo newcalc=1 call nearfieldpointcalc(neqns,nsphere,nodr,alpha,beta,cbeam,xsp,rpos,ri,amnp, & gamma,epspw,xgpmax,newcalc,efield,hfield) ! ! determine the intersecting spheres ! nsp=0 do i=1,nsphere xi0=abs(rpos(gcoord(3),i)-xg(gcoord(3))) if(xi0.le.xsp(i)) then nsp=nsp+1 xplot(1,nsp)=rpos(gcoord(1),i)+gbfocus(gcoord(1)) xplot(2,nsp)=rpos(gcoord(2),i)+gbfocus(gcoord(2)) ri0=xsp(i)*xsp(i)-xi0*xi0 if(ri0.ne.0.) ri0=sqrt(ri0) xplot(3,nsp)=ri0 endif enddo ! ! report to runprintunit ! if(rank.eq.0) then write(runprintunit,'('' near field calculations'')') write(runprintunit,'('' plane, position:'',i5,f9.3)') nfplane,nfplanepos0 write(runprintunit,'('' rectangular plot vertices:'')') write(runprintunit,'('' min:'',3f9.3)') nfplanevert0(1:2,1) write(runprintunit,'('' max:'',3f9.3)') nfplanevert0(1:2,2) write(runprintunit,'('' number of plotting points, step size:'',i8,f8.3)') npoints, deltax write(runprintunit,'('' max plane wave order:'',i5)') nodrpwmax endif ! ! determine the distribution of work among the processors ! allocate(efindex(0:numprocs-1),efnum(0:numprocs-1), & efieldrow(3,npoints2),efieldrowt(3,npoints2), & hfieldrow(3,npoints2),hfieldrowt(3,npoints2)) np23=3*npoints2 frowperproc=dble(npoints2)/dble(numprocs) rowsum=0. do i=0,numprocs-1 efindex(i)=floor(rowsum) rowsum=rowsum+frowperproc enddo do i=0,numprocs-2 efnum(i)=efindex(i+1)-efindex(i) enddo efnum(numprocs-1)=npoints2-efindex(numprocs-1) npoints1by5=int(npoints1/5.+.5) ! ! do the calculations and write the results to the file ! if(rank.eq.0) then write(nfoutunit,*) npoints1,npoints2 write(nfoutunit,*) nsp do i=1,nsp write(nfoutunit,'(3e13.5)') xplot(1,i),xplot(2,i),xplot(3,i) enddo time1=mytime() endif xg(gcoord(3))=nfplanepos newcalc=0 do i=1,npoints1 xg(gcoord(1))=nfplanevert(1,1)+deltax*dble(i-1) xgp(gcoord(1))=nfplanevert0(1,1)+deltax*dble(i-1) efieldrowt=0.d0 efieldrow=0.d0 hfieldrowt=0.d0 hfieldrow=0.d0 do j=efindex(rank)+1,efindex(rank)+efnum(rank) xg(gcoord(2))=nfplanevert(2,1)+deltax*dble(j-1) call nearfieldpointcalc(neqns,nsphere,nodr,alpha,beta,cbeam,xsp,rpos,ri,amnp, & gamma,epspw,xg,newcalc,efieldrowt(:,j),hfieldrowt(:,j)) enddo call ms_mpi(mpi_command='barrier') call ms_mpi(mpi_command='reduce',mpi_send_buf_dc=efieldrowt,mpi_recv_buf_dc=efieldrow,& mpi_number=np23,mpi_rank=0,mpi_operation=ms_mpi_sum) if(nfoutdata.ge.2) then call ms_mpi(mpi_command='reduce',mpi_send_buf_dc=hfieldrowt,mpi_recv_buf_dc=hfieldrow,& mpi_number=np23,mpi_rank=0,mpi_operation=ms_mpi_sum) endif if(rank.eq.0) then do j=1,npoints2 xgp(gcoord(2))=nfplanevert0(2,1)+deltax*dble(j-1) if(nfoutdata.eq.0) then esquare=dot_product(efieldrow(:,j),efieldrow(:,j)) write(nfoutunit,'(2f9.4,e13.5)') xgp(gcoord(1)),xgp(gcoord(2)),esquare elseif(nfoutdata.eq.1) then write(nfoutunit,'(2f9.4,6e13.5)') xgp(gcoord(1)),xgp(gcoord(2)),efieldrow(:,j) else write(nfoutunit,'(2f9.4,12e13.5)') xgp(gcoord(1)),xgp(gcoord(2)),efieldrow(:,j), & hfieldrow(:,j) endif enddo if(mod(i,npoints1by5).eq.0) then k=i/npoints1by5 time2=(mytime()-time1)*dble(5-k)/dble(k) call timewrite(runprintunit,' estimated time remaining:',time2) endif endif enddo deallocate(efindex,efnum,efieldrow,efieldrowt,hfieldrow,hfieldrowt) end subroutine nearfieldgridcalc ! ! nearfieldaverage is an MPI--enabled subroutine for calculating average field ! values along a line. ! subroutine nearfieldaverage(neqns,nsphere,nodr,alpha,beta,cbeam,xsp,rpos,ri,amnp, & nfplane,nfplaneposstart0,nfplaneposend0,numberplanes,nfplanevert0,gbfocus, & deltax,gamma,epspw,runprintunit,efieldavez,hfieldavez,svecavez) use mpidefs use mpidata use intrinsics use specialfuncs use spheredata use miecoefdata use numconstants implicit none integer :: nsphere,neqns,nodr(nsphere),nfplane,runprintunit,npoints1,npoints2, & npoints,i,j,k,np23,gcoord(3),rank,numprocs,nrowperproc,nrowrem, & npoints1by5,plotincfield,nsp,nfoutunit,nfoutdata,newcalc,nsend integer :: numberplanes real (8) :: alpha,beta,cbeam,xsp(nsphere),rpos(3,nsphere),nfplanepos,& nfplanevert(2,2),frowperproc,rowsum,xg(3),xgp(3),deltax,gamma,epspw, & time1,time2,xplot(3,nsphere),xi0,ri0,esquare,xgpmax(3),rplot, & gbfocus(3),nfplanepos0,nfplanevert0(2,2),nfplaneposstart,nfplaneposend, & deltanfplane,nfplaneposstart0,nfplaneposend0,svec(3),svecave(3) real (8) :: svecavez(3,numberplanes) real(8), allocatable :: xypoint(:,:) complex(8) :: amnp(neqns,2),ri(2,nsphere),efield(3),hfield(3),hfieldave(3),efieldave(3) complex(8) :: efieldavez(3,numberplanes),hfieldavez(3,numberplanes) integer, allocatable :: efindex(:),efnum(:) call ms_mpi(mpi_command='size',mpi_size=numprocs) call ms_mpi(mpi_command='rank',mpi_rank=rank) ! ! determine the plane ! if(nfplane.eq.1) then gcoord=(/2,3,1/) elseif(nfplane.eq.2) then gcoord=(/3,1,2/) else gcoord=(/1,2,3/) endif ! ! shift the coordinates to gb focal origin ! nfplanevert(1,1)=nfplanevert0(1,1)-gbfocus(gcoord(1)) nfplanevert(1,2)=nfplanevert0(1,2)-gbfocus(gcoord(1)) nfplanevert(2,1)=nfplanevert0(2,1)-gbfocus(gcoord(2)) nfplanevert(2,2)=nfplanevert0(2,2)-gbfocus(gcoord(2)) nfplaneposstart=nfplaneposstart0-gbfocus(gcoord(3)) nfplaneposend=nfplaneposend0-gbfocus(gcoord(3)) xg(gcoord(3))=nfplanepos ! ! determine the number of points ! npoints1=nint((nfplanevert(1,2)-nfplanevert(1,1))/deltax)+1 npoints2=nint((nfplanevert(2,2)-nfplanevert(2,1))/deltax)+1 npoints=npoints1*npoints2 allocate(efindex(0:numprocs-1),efnum(0:numprocs-1),xypoint(2,npoints)) frowperproc=dble(npoints)/dble(numprocs) rowsum=0. do i=0,numprocs-1 efindex(i)=floor(rowsum) rowsum=rowsum+frowperproc enddo do i=0,numprocs-2 efnum(i)=efindex(i+1)-efindex(i) enddo efnum(numprocs-1)=npoints-efindex(numprocs-1) xgp(gcoord(3))=max(abs(nfplaneposstart),abs(nfplaneposend)) rplotmax=0.d0 xgpmax=0.d0 k=0 do i=1,npoints1 xgp(gcoord(1))=nfplanevert(1,1)+deltax*dble(i-1) do j=1,npoints2 xgp(gcoord(2))=nfplanevert(2,1)+deltax*dble(j-1) k=k+1 xypoint(1,k)=xgp(gcoord(1)) xypoint(2,k)=xgp(gcoord(2)) rplot=sqrt(dot_product(xgp,xgp)) if(rplot.gt.rplotmax) then rplotmax=rplot xgpmax=xgp endif enddo enddo newcalc=1 call nearfieldpointcalc(neqns,nsphere,nodr,alpha,beta,cbeam,xsp,rpos,ri,amnp, & gamma,epspw,xgpmax,newcalc,efield,hfield) newcalc=0 do k=1,numberplanes if(numberplanes.eq.1) then nfplanepos=nfplaneposstart else nfplanepos=nfplaneposstart+(nfplaneposend-nfplaneposstart)*(k-1)/dble(numberplanes-1) endif ! ! find the maximum point-to-target origin distance and initialize the field calculation ! xg(3)=nfplanepos efieldave=0. hfieldave=0. svecave=0. do i=efindex(rank)+1,efindex(rank)+efnum(rank) xg(gcoord(1))=xypoint(1,i) xg(gcoord(2))=xypoint(2,i) efield=0.d0 hfield=0.d0 call nearfieldpointcalc(neqns,nsphere,nodr,alpha,beta,cbeam,xsp,rpos,ri,amnp, & gamma,epspw,xg,newcalc,efield,hfield) efieldave=efieldave+efield hfieldave=hfieldave+hfield hfield=conjg(hfield)/2.d0 svec(1)=-efield(3)*hfield(2)+efield(2)*hfield(3) svec(2)=efield(3)*hfield(1)-efield(1)*hfield(3) svec(3)=-efield(2)*hfield(1)+efield(1)*hfield(2) svecave=svecave+svec enddo call ms_mpi(mpi_command='barrier') efield=0.d0 hfield=0.d0 svec=0.d0 call ms_mpi(mpi_command='reduce',mpi_send_buf_dc=efieldave,mpi_recv_buf_dc=efield,& mpi_number=3,mpi_rank=0,mpi_operation=ms_mpi_sum) call ms_mpi(mpi_command='reduce',mpi_send_buf_dc=hfieldave,mpi_recv_buf_dc=hfield,& mpi_number=3,mpi_rank=0,mpi_operation=ms_mpi_sum) call ms_mpi(mpi_command='reduce',mpi_send_buf_dp=svecave,mpi_recv_buf_dp=svec,& mpi_number=3,mpi_rank=0,mpi_operation=ms_mpi_sum) call ms_mpi(mpi_command='barrier') if(rank.eq.0) then efield=efield/dble(npoints) hfield=hfield/dble(npoints) svec=svec/dble(npoints) efieldavez(:,k)=efield hfieldavez(:,k)=hfield svecavez(:,k)=svec i=gcoord(3) write(runprintunit,'('' plane:'',i5,f9.3,2e12.4)') k,nfplanepos, & svec(i) call flush(runprintunit) endif enddo nsend=3*numberplanes call ms_mpi(mpi_command='bcast',mpi_send_buf_dc=efieldavez, & mpi_number=nsend,mpi_rank=0) call ms_mpi(mpi_command='bcast',mpi_send_buf_dc=hfieldavez, & mpi_number=nsend,mpi_rank=0) call ms_mpi(mpi_command='bcast',mpi_send_buf_dp=svecavez, & mpi_number=nsend,mpi_rank=0) call ms_mpi(mpi_command='barrier') deallocate(efindex,efnum,xypoint) end subroutine nearfieldaverage end module nearfield ! ! module solver: subroutines for solving interaction equations for fixed orientation ! and T matrix problems ! ! ! last revised: 15 January 2011 ! module solver implicit none contains ! ! tmatrixsoln: calculation of T matrix via solution of interaction equations for ! a generalized plane wave expansion ! ! ! original: 15 January 2011 ! revised: 21 February 2011: call for sphereqeff changed ! subroutine tmatrixsoln(neqns,nsphere,nodr,nodrt,xsp,rpos,epssoln,epscon,niter,& calctmatrix,tmatrixfile,fftranpresent,niterstep,qext,qabs,qsca,istat) use mpidefs use mpidata use intrinsics use numconstants use specialfuncs use miecoefdata use spheredata use translation use scatprops implicit none integer :: iter,niter,neqns,nsphere,nodr(nsphere),ntran(nsphere),nodrt, & nodrmax,i,ierr,istat,m,n,p,k,l,q,noff,nblk,ma,na,ka,la,mn,istart,iunit, & rank,iexit(1),calctmatrix,lt,kt,qt,nt,mt,it,nodrtt,lstart(1),numsolns,isoln, & isolnstart,igroup,ngroup,rgrank,lold,grank,nsend,nodrta(1), & fftranpresent,niterstep integer, allocatable :: lindex(:),kindex(:) real(8) :: eps,err,qext(nsphere),qabs(nsphere),qsca(nsphere),xsp(nsphere),xv, & rpos(3,nsphere),xij(3),qabsklq(nsphere),qscaklq(nsphere),qextklq(nsphere), & f2,qexttot,qabstot,qscatot,qextold(1),qscaold(1),errqe,errqs, & timetran,timesolve,time1,time2,epssoln,epscon,dtemp(4),timeorder,& at1,at2,at3,at4 real(8) :: qextl(nsphere),qabsl(nsphere),qscal(nsphere) real(8), allocatable :: qextgroup(:,:),qabsgroup(:,:),qscagroup(:,:) complex(8) :: amnp(neqns),pmnp(neqns),pmnpan(neqns) complex(8), allocatable :: pmnp0(:,:,:),ac(:,:,:,:),pmnpt(:,:,:),amnp0(:,:,:), & amnp0group(:,:,:) character*30 :: tmatrixfile character*4 :: timeunit data istart,iexit/1,0/ rank=base_rank rgrank=root_group_rank grank=group_rank ngroup=number_groups call getrunparameters(run_print_unit=iunit) xv=(sum(xsp**3.d0))**(1.d0/3.d0) nodrmax=maxval(nodr) qext=0.d0 qabs=0.d0 qsca=0.d0 qextold=0.d0 qscaold=0.d0 ! ! perform T matrix file operations as needed ! if(rank.eq.0) then if(calctmatrix.eq.1) then open(3,file=tmatrixfile) write(3,'(i4)') nodrt lstart(1)=1 else open(3,file=tmatrixfile) write(iunit,'('' finding end of record to file '',a)') tmatrixfile read(3,*) nodrtt do l=1,nodrt do k=-l,l do q=1,2 read(3,'(3i5)',end=20,err=20) lt,kt,qt do n=1,l do m=-n,n read(3,'(2i5,4e17.9)',end=20,err=20) nt,mt,at1,at2,at3,at4 enddo enddo enddo enddo do i=1,nsphere read(3,'(i5,3e17.9)',end=20,err=20) it,qextl(i),qabsl(i),qscal(i) enddo qext=qext+qextl qabs=qabs+qabsl qsca=qsca+qscal enddo 20 close(3) open(3,file=tmatrixfile) qextold(1)=0.d0 qabstot=0.d0 do i=1,nsphere qextold=qextold+qext(i)*xsp(i)*xsp(i)/xv/xv qabstot=qabstot+qabs(i)*xsp(i)*xsp(i)/xv/xv enddo qscaold(1)=qextold(1)-qabstot lstart(1)=lt read(3,*) nodrtt write(iunit,'('' calculations begin with order '',i5)') lstart(1) do l=1,lstart(1)-1 do k=-l,l do q=1,2 read(3,'(3i5)') lt,kt,qt do n=1,l do m=-n,n read(3,'(2i5,4e17.9)') nt,mt,at1,at2,at3,at4 enddo enddo enddo enddo do i=1,nsphere read(3,'(i5,3e17.9)') it,at1,at2,at3 enddo enddo endif endif call ms_mpi(mpi_command='bcast',mpi_send_buf_i=lstart,mpi_number=1,mpi_rank=0) call ms_mpi(mpi_command='bcast',mpi_send_buf_dp=qextold,mpi_number=1,mpi_rank=0) call ms_mpi(mpi_command='bcast',mpi_send_buf_dp=qscaold,mpi_number=1,mpi_rank=0) call ms_mpi(mpi_command='bcast',mpi_send_buf_dp=qext,mpi_number=nsphere,mpi_rank=0) call ms_mpi(mpi_command='bcast',mpi_send_buf_dp=qabs,mpi_number=nsphere,mpi_rank=0) call ms_mpi(mpi_command='bcast',mpi_send_buf_dp=qsca,mpi_number=nsphere,mpi_rank=0) call ms_mpi(mpi_command='barrier') allocate(amnp0group(2,nodrt*(nodrt+2),0:ngroup-1),qextgroup(nsphere,0:ngroup-1), & qabsgroup(nsphere,0:ngroup-1),qscagroup(nsphere,0:ngroup-1)) numsolns=2*nodrt*(nodrt+2) allocate(lindex(numsolns),kindex(numsolns)) ! ! find the starting point ! i=0 do l=1,nodrt do k=-l,l do q=1,2 i=i+1 lindex(i)=l kindex(i)=k enddo enddo enddo do i=1,numsolns if(lindex(i).eq.lstart(1).and.kindex(i).eq.-lstart(1)) exit enddo isolnstart=i qextl=0.d0 qabsl=0.d0 qscal=0.d0 lold=0 ! ! begin the loop over RHS of the interaction equations. The solutions are distributed ! among ngroup groups of processors ! do isoln=isolnstart,numsolns,ngroup if(rank.eq.0) timeorder=mytime() do igroup=0,ngroup-1 l=lindex(isoln+igroup) k=kindex(isoln+igroup) q=mod(isoln+igroup-1,2)+1 ! ! calculate the RHS ! if(l.eq.-k.and.q.eq.1) then if(allocated(ac)) deallocate(ac,amnp0) allocate(ac(2,nodrmax*(nodrmax+2),-l:l,nsphere),amnp0(0:l+1,l,2)) do i=1,nsphere xij=rpos(:,i) call gentrancoef(1,xij,(1.d0,0.d0),1,nodrmax,l,l,0,0,ac(1,1,-l,i)) enddo endif if(igroup.eq.rgrank) then if(k.le.-1) then ka=l+1 la=-k else ka=k la=l endif noff=0 do i=1,nsphere nblk=nodr(i)*(nodr(i)+2)*2 allocate(pmnp0(0:nodr(i)+1,nodr(i),2)) do p=1,2 do n=1,nodr(i) do m=-n,n mn=n*(n+1)+m if(m.le.-1) then ma=n+1 na=-m else ma=m na=n endif pmnp0(ma,na,p)=ac(abs(p-q)+1,mn,k,i) enddo enddo enddo pmnp(noff+1:noff+nblk)=reshape(pmnp0,(/nblk/)) deallocate(pmnp0) noff=noff+nblk enddo ! ! multiply RHS by mie coefficients ! call miecoeffmult(1,nsphere,neqns,pmnp,pmnpan) amnp=pmnpan ! ! call the solver ! if(fftranpresent.eq.1) then call cbicgff(neqns,nsphere,niter,epssoln,pmnpan,amnp,0, & niterstep,iter,err) else call cbicg(neqns,nsphere,niter,epssoln,pmnpan,amnp,0,iter,err) endif if(iter.gt.niter.or.err.gt.epssoln) istat=1 call sphereqeff(nsphere,neqns,nodr,nodrmax,xsp,amnp,amnp, & pmnp,pmnp,qextklq,qabsklq,qscaklq) qextgroup(1:nsphere,igroup)=qextklq(1:nsphere) qabsgroup(1:nsphere,igroup)=qabsklq(1:nsphere) qscagroup(1:nsphere,igroup)=qscaklq(1:nsphere) ! ! compute the target-based expansion ! ntran=l amnp0=0.d0 call amncommonorigin(neqns,nsphere,nodr,ntran,l,rpos, & amnp,amnp0) do n=1,l do m=-n,n if(m.le.-1) then ma=n+1 na=-m else ma=m na=n endif mn=n*(n+1)+m do p=1,2 amnp0group(p,mn,igroup)=amnp0(ma,na,p) enddo enddo enddo endif enddo ! ! send the solutions to the rank 0 processor ! call ms_mpi(mpi_command='barrier') if(grank.eq.0) then if(rank.ne.0) then l=lindex(isoln+rgrank) nblk=l*(l+2) nsend=2*nblk call ms_mpi(mpi_command='send',mpi_send_buf_dc=amnp0group(1,1,rgrank),& mpi_number=nsend,mpi_rank=0,mpi_comm=root_group_comm) call ms_mpi(mpi_command='send',mpi_send_buf_dp=qextgroup(1,rgrank),& mpi_number=nsphere,mpi_rank=0,mpi_comm=root_group_comm) call ms_mpi(mpi_command='send',mpi_send_buf_dp=qabsgroup(1,rgrank),& mpi_number=nsphere,mpi_rank=0,mpi_comm=root_group_comm) call ms_mpi(mpi_command='send',mpi_send_buf_dp=qscagroup(1,rgrank),& mpi_number=nsphere,mpi_rank=0,mpi_comm=root_group_comm) else do igroup=1,ngroup-1 l=lindex(isoln+igroup) nblk=l*(l+2) nsend=2*nblk call ms_mpi(mpi_command='recv',mpi_recv_buf_dc=amnp0group(1,1,igroup),& mpi_number=nsend,mpi_rank=igroup,mpi_comm=root_group_comm) call ms_mpi(mpi_command='recv',mpi_recv_buf_dp=qextgroup(1,igroup),& mpi_number=nsphere,mpi_rank=igroup,mpi_comm=root_group_comm) call ms_mpi(mpi_command='recv',mpi_recv_buf_dp=qabsgroup(1,igroup),& mpi_number=nsphere,mpi_rank=igroup,mpi_comm=root_group_comm) call ms_mpi(mpi_command='recv',mpi_recv_buf_dp=qscagroup(1,igroup),& mpi_number=nsphere,mpi_rank=igroup,mpi_comm=root_group_comm) enddo endif endif call ms_mpi(mpi_command='barrier') ! ! write results, check for convergence ! if(rank.eq.0) then do igroup=0,ngroup-1 l=lindex(isoln+igroup) k=kindex(isoln+igroup) q=mod(isoln+igroup-1,2)+1 qextl=qextl+qextgroup(1:nsphere,igroup) qabsl=qabsl+qabsgroup(1:nsphere,igroup) qscal=qscal+qscagroup(1:nsphere,igroup) qext=qext+qextgroup(1:nsphere,igroup) qabs=qabs+qabsgroup(1:nsphere,igroup) qsca=qsca+qscagroup(1:nsphere,igroup) write(3,'(3i5)') l,k,q do n=1,l do m=-n,n mn=n*(n+1)+m write(3,'(2i5,4e17.9)') n,m,amnp0group(1,mn,igroup), & amnp0group(2,mn,igroup) enddo enddo if(istart.eq.1.and.igroup.eq.0) then time1=mytime()-timeorder call timewrite(iunit,' time per group solution:',time1) time2=time1*dble(numsolns-isolnstart)/dble(ngroup) call timewrite(iunit,' estimated t matrix calcuation time:',time2) write(iunit,'('' n # its qext qabs'',& &'' qsca error est. time rem.'')') call flush(iunit) istart=0 endif if(igroup.eq.0) then timeorder=mytime()-timeorder time2=timeorder*dble(numsolns-isoln)/dble(ngroup) if(time2.gt.3600.d0) then time2=time2/3600.d0 timeunit=' hrs' elseif(time2.gt.60.d0) then time2=time2/60.d0 timeunit=' min' else timeunit=' sec' endif endif iexit(1)=0 if(k.eq.l.and.q.eq.2) then qexttot=0.d0 qabstot=0.d0 do i=1,nsphere qexttot=qexttot+qext(i)*xsp(i)*xsp(i)/xv/xv qabstot=qabstot+qabs(i)*xsp(i)*xsp(i)/xv/xv write(3,'(i5,3e17.9)') i,qextl(i),qabsl(i),qscal(i) enddo qextl=0.d0 qabsl=0.d0 qscal=0.d0 qscatot=qexttot-qabstot errqe=qexttot-qextold(1) errqs=qscatot-qscaold(1) err=max(errqe,errqs) write(iunit,'(i4,i5,4e13.5,f8.2,a4)') l,iter,qexttot,qabstot, & qscatot,err,time2,timeunit call flush(iunit) qextold(1)=qexttot qscaold(1)=qscatot if(err.le.epscon) iexit(1)=1 endif if(iexit(1).eq.1) then nodrt=l exit endif enddo endif call ms_mpi(mpi_command='bcast',mpi_send_buf_i=iexit,mpi_number=1,mpi_rank=0) call ms_mpi(mpi_command='barrier') if(iexit(1).eq.1) then ! ! solution has converged ! deallocate(amnp0group,qextgroup,qabsgroup,qscagroup,lindex,kindex,ac,amnp0) nodrta(1)=nodrt call ms_mpi(mpi_command='bcast',mpi_send_buf_i=nodrta,mpi_number=1,mpi_rank=0) nodrt=nodrta(1) if(rank.eq.0) then write(iunit,'('' T matrix converged, order:'',i5)') nodrt close(3) open(3,file=tmatrixfile,form='formatted',access='direct',recl=4) write(3,'(i4)',rec=1) nodrt close(3) endif return endif enddo deallocate(amnp0group,qextgroup,qabsgroup,qscagroup,lindex,kindex,ac,amnp0) if(rank.eq.0) then write(*,'('' T matrix did not converge to set epsilon'')') close(3) endif end subroutine tmatrixsoln ! ! solution of interaction equations for a fixed orientation ! ! ! original: 15 January 2011 ! revised: 21 February 2011: modification of efficiency calculation, to calculate ! polarized components ! 30 March 2011: took out gbfocus argument: this is not needed since positions are defined ! relative to the gb focus. ! 20 April 2011: used 2-group MPI formulation ! subroutine fixedorsoln(neqns,nsphere,nodr,alpha,beta,cbeam,xsp,rpos,& eps,epstran,niter,amnp,qext,qabs,qsca,maxerr,maxiter,iterwrite, & fftranpresent,niterstep,istat) use mpidefs use mpidata use intrinsics use numconstants use specialfuncs use miecoefdata use translation use scatprops implicit none integer :: iter,niter,neqns,nodrmax,k,nsphere,i,ierr,istat,rank,maxiter,iterwrite integer :: nodr(nsphere),m1,n1,p,rgrank,grank,ngroup,sendrank,numprocs, & fftranpresent,niterstep real(8) :: alpha,beta,eps,err,qext(nsphere,3),maxerr,& qabs(nsphere,3),qsca(nsphere,3),cbeam,gbfocus(3),epstran real(8) :: xsp(nsphere), rpos(3,nsphere),maxerra(1) complex(8) :: amnp(neqns,2) complex(8), allocatable :: pmnp(:,:),pmnpan(:) rank=base_rank rgrank=root_group_rank grank=group_rank ngroup=number_groups numprocs=number_proc sendrank=numprocs/2 nodrmax=maxval(nodr) allocate(pmnp(neqns,2)) gbfocus=0.d0 if(cbeam.eq.0.d0) then call sphereplanewavecoef(nsphere,neqns,nodr,nodrmax,alpha,beta,rpos,pmnp) else call spheregaussianbeamcoef(nsphere,neqns,nodr,alpha,beta,cbeam, & rpos,gbfocus,epstran,pmnp) endif istat=0 maxiter=0 maxerr=0. ! ! calculate the two solutions ! allocate(pmnpan(neqns)) if(ngroup.eq.1) then do k=1,2 call miecoeffmult(1,nsphere,neqns,pmnp(1,k),pmnpan) amnp(1:neqns,k)=pmnpan(1:neqns) if(fftranpresent.eq.1) then call cbicgff(neqns,nsphere,niter,eps,pmnpan,amnp(1,k),iterwrite, & niterstep,iter,err) else call cbicg(neqns,nsphere,niter,eps,pmnpan,amnp(1,k),iterwrite, & iter,err) endif maxiter=max(iter,maxiter) maxerr=max(err,maxerr) if(iter.gt.niter.or.err.gt.eps) istat=1 call ms_mpi(mpi_command='barrier') enddo else k=rgrank+1 call miecoeffmult(1,nsphere,neqns,pmnp(1,k),pmnpan) amnp(1:neqns,k)=pmnpan(1:neqns) if(fftranpresent.eq.1) then call cbicgff(neqns,nsphere,niter,eps,pmnpan,amnp(1,k),iterwrite, & niterstep,iter,err) else call cbicg(neqns,nsphere,niter,eps,pmnpan,amnp(1,k),iterwrite, & iter,err) endif maxiter=max(iter,maxiter) maxerr=max(err,maxerr) call ms_mpi(mpi_command='barrier') call ms_mpi(mpi_command='bcast',mpi_send_buf_dc=amnp(1,2),& mpi_number=neqns,mpi_rank=sendrank) endif deallocate(pmnpan) ! ! efficiency factor calculations ! call sphereqeff(nsphere,neqns,nodr,nodrmax,xsp,amnp(1,1),amnp(1,1), & pmnp(1,1),pmnp(1,1),qext(1,1),qabs(1,1),qsca(1,1)) call sphereqeff(nsphere,neqns,nodr,nodrmax,xsp,amnp(1,2),amnp(1,2), & pmnp(1,2),pmnp(1,2),qext(1,2),qabs(1,2),qsca(1,2)) call sphereqeff(nsphere,neqns,nodr,nodrmax,xsp,amnp(1,1),amnp(1,2), & pmnp(1,1),pmnp(1,2),qext(1,3),qabs(1,3),qsca(1,3)) call ms_mpi(mpi_command='barrier') deallocate(pmnp) end subroutine fixedorsoln ! ! hybrid bcgm, using far field translation ! november 2011 ! subroutine cbicgff(neqns,nsphere,niter,eps,pnp,anp,iterwrite,niterstep,iter,err) use mpidefs use mpidata use intrinsics use spheredata use miecoefdata use numconstants use specialfuncs use translation implicit none integer :: neqns,niter,iter,nsphere,writetime,nodr(nsphere),noff(nsphere+1),& nblk(nsphere),rank,ierr,iexit,i,j,m,n,p,iunit,iterwrite,ip1,ip2, & np1,np2,nsend,numprocs,grank,istore,itermax,istep,niterstep real(8) :: eps,err,erra(1),enorm,time1,time2,epsstep,errstep complex(8) :: pnp(neqns),anp(neqns),gnp(neqns),gnpold(neqns),pgnp(neqns), & cr(neqns) data writetime/0/ rank=base_rank grank=group_rank numprocs=proc_per_group call getmiedata(sphere_order=nodr,sphere_block=nblk,sphere_block_offset=noff) if(rank.eq.0) then call getrunparameters(run_print_unit=iunit) endif err=0.d0 iter=0 enorm=dot_product(pnp,pnp) gnpold=0.d0 if(enorm.eq.0.d0) return gnp=0.d0 ip1=mpi_sphere_index(grank)+1 ip2=mpi_sphere_index(grank)+mpi_sphere_number(grank) do i=ip1,ip2 do j=1,nsphere if(i.ne.j) then call fftranslationerror(anp(noff(j)+1:noff(j)+nblk(j)), & gnp(noff(i)+1:noff(i)+nblk(i)),j,i,nodr(j),nodr(i)) endif enddo enddo do i=0,numprocs-1 ip1=mpi_sphere_index(i)+1 ip2=mpi_sphere_index(i)+mpi_sphere_number(i) np1=noff(ip1)+1 np2=noff(ip2)+nblk(ip2) nsend=np2-np1+1 call ms_mpi(mpi_command='bcast',mpi_send_buf_dc=gnp(np1:np2),mpi_number=nsend, & mpi_rank=i,mpi_comm=group_comm) enddo call miecoeffmult(1,nsphere,neqns,gnp,gnp) iter=0 epsstep=eps istep=0 do istep=istep+1 gnpold=gnp pgnp=pnp+gnp call cbicg(neqns,nsphere,niterstep,epsstep,pgnp,anp,0,itermax,errstep) iter=iter+min(itermax,niterstep) gnp=0.d0 ip1=mpi_sphere_index(grank)+1 ip2=mpi_sphere_index(grank)+mpi_sphere_number(grank) do i=ip1,ip2 do j=1,nsphere if(i.ne.j) then call fftranslationerror(anp(noff(j)+1:noff(j)+nblk(j)), & gnp(noff(i)+1:noff(i)+nblk(i)),j,i,nodr(j),nodr(i)) endif enddo enddo do i=0,numprocs-1 ip1=mpi_sphere_index(i)+1 ip2=mpi_sphere_index(i)+mpi_sphere_number(i) np1=noff(ip1)+1 np2=noff(ip2)+nblk(ip2) nsend=np2-np1+1 call ms_mpi(mpi_command='bcast',mpi_send_buf_dc=gnp(np1:np2),mpi_number=nsend, & mpi_rank=i,mpi_comm=group_comm) enddo call miecoeffmult(1,nsphere,neqns,gnp,gnp) err=dot_product(gnp-gnpold,gnp-gnpold)/enorm if(rank.eq.0.and.iterwrite.eq.1) then write(iunit,'('' step,iteration,bcgm err,correc err:'',2i5,2e13.5)') & istep,iter,errstep,err call flush(iunit) endif epsstep=eps err=max(err,errstep) if((err.lt.eps).or.iter.gt.niter) exit enddo end subroutine cbicgff ! ! iteration solver ! generalized complex biconjugate gradient method ! original code: Piotr Flatau, although not much remains. ! specialized to the multiple sphere problem ! ! ! last revised: 15 January 2011 ! october 2011: translation calls modified ! subroutine cbicg(neqns,nsphere,niter,eps,pnp,anp,iterwrite,iter,err) use mpidefs use mpidata use intrinsics use spheredata use miecoefdata use numconstants use specialfuncs use translation implicit none integer :: neqns,niter,iter,nsphere,writetime,nodr(nsphere),noff(nsphere+1),& nblk(nsphere),rank,ierr,iexit,i,j,m,n,p,iunit,iterwrite,ip1,ip2, & np1,np2,nsend,numprocs,grank,istore real(8) :: eps,err,erra(1),enorm,time1,time2 complex(8) :: pnp(neqns),anp(neqns) complex(8) :: cak(1),csk,cbk,csk2(1) complex(8) :: cr(neqns),cp(neqns),cw(neqns),cq(neqns),cap(neqns),caw(neqns), & crt(neqns),capt(neqns),cawt(neqns),ccw(neqns) data writetime/0/ rank=base_rank grank=group_rank numprocs=proc_per_group call getmiedata(sphere_order=nodr,sphere_block=nblk,sphere_block_offset=noff) ip1=mpi_sphere_index(grank)+1 ip2=mpi_sphere_index(grank)+mpi_sphere_number(grank) np1=noff(ip1)+1 np2=noff(ip2)+nblk(ip2) nsend=np2-np1+1 crt=0.d0 iexit=0 if(rank.eq.0) then call getrunparameters(run_print_unit=iunit) endif err=0.d0 iter=0 enorm=dot_product(pnp,pnp) cr=0.d0 if(enorm.eq.0.d0) return do i=ip1,ip2 do j=1,nsphere if(i.ne.j) then call rottranjtoi(anp(noff(j)+1:noff(j)+nblk(j)), & cr(noff(i)+1:noff(i)+nblk(i)),j,i,nodr(j),nodr(i),1,1) endif enddo enddo do i=0,numprocs-1 ip1=mpi_sphere_index(i)+1 ip2=mpi_sphere_index(i)+mpi_sphere_number(i) np1=noff(ip1)+1 np2=noff(ip2)+nblk(ip2) nsend=np2-np1+1 call ms_mpi(mpi_command='bcast',mpi_send_buf_dc=cr(np1:np2),mpi_number=nsend, & mpi_rank=i,mpi_comm=group_comm) enddo call miecoeffmult(1,nsphere,neqns,cr,cr) cr=pnp-anp+cr cq=conjg(cr) cw=cq cp=cr csk=dot_product(conjg(cr),cr) if(cdabs(csk).eq.0.d0) return ! ! here starts the main iteration loop ! do iter=1,niter ip1=mpi_sphere_index(grank)+1 ip2=mpi_sphere_index(grank)+mpi_sphere_number(grank) np1=noff(ip1)+1 np2=noff(ip2)+nblk(ip2) nsend=np2-np1+1 cak(1)=(0.d0,0.d0) cawt=(0.d0,0.d0) capt=(0.d0,0.d0) if(rank.eq.0) then if(writetime.eq.0) time1=mytime() endif ccw=conjg(cw) cap=0.d0 caw=0.d0 call miecoeffmult(1,nsphere,neqns,ccw,ccw) do i=ip1,ip2 do j=1,nsphere if(i.ne.j) then call rottrantwojtoi(cp(noff(j)+1:noff(j)+nblk(j)), & ccw(noff(j)+1:noff(j)+nblk(j)), & cap(noff(i)+1:noff(i)+nblk(i)), & caw(noff(i)+1:noff(i)+nblk(i)), & j,i,nodr(j),nodr(i)) ! call rottranjtoi(cp(noff(j)+1:noff(j)+nblk(j)), & ! cap(noff(i)+1:noff(i)+nblk(i)),j,i,nodr(j),nodr(i),1,1) ! call rottranjtoi(ccw(noff(j)+1:noff(j)+nblk(j)), & ! caw(noff(i)+1:noff(i)+nblk(i)),j,i,nodr(j),nodr(i),-1,-1) endif enddo enddo call miecoeffmult(ip1,ip2,neqns,cap,cap) cap(np1:np2)=cp(np1:np2)-cap(np1:np2) caw(np1:np2)=cw(np1:np2)-conjg(caw(np1:np2)) cak(1)=dot_product(cw(np1:np2),cap(np1:np2)) call ms_mpi(mpi_command='allreduce',mpi_recv_buf_dc=cak,mpi_number=1, & mpi_operation=ms_mpi_sum,mpi_comm=group_comm) cak(1)=csk/cak(1) anp(np1:np2)=anp(np1:np2)+cak(1)*cp(np1:np2) cr(np1:np2)=cr(np1:np2)-cak(1)*cap(np1:np2) cq(np1:np2)=cq(np1:np2)-conjg(cak(1))*caw(np1:np2) csk2(1)=dot_product(cq(np1:np2),cr(np1:np2)) err=dot_product(cr(np1:np2),cr(np1:np2)) erra(1)=err call ms_mpi(mpi_command='allreduce',mpi_recv_buf_dc=csk2,mpi_number=1, & mpi_operation=ms_mpi_sum,mpi_comm=group_comm) call ms_mpi(mpi_command='allreduce',mpi_recv_buf_dp=erra,mpi_number=1, & mpi_operation=ms_mpi_sum,mpi_comm=group_comm) err=erra(1) err=err/enorm if(err.lt. eps) exit cbk=csk2(1)/csk cp(np1:np2)=cr(np1:np2)+cbk*cp(np1:np2) cw(np1:np2)=cq(np1:np2)+conjg(cbk)*cw(np1:np2) csk=csk2(1) do i=0,numprocs-1 ip1=mpi_sphere_index(i)+1 ip2=mpi_sphere_index(i)+mpi_sphere_number(i) np1=noff(ip1)+1 np2=noff(ip2)+nblk(ip2) nsend=np2-np1+1 call ms_mpi(mpi_command='bcast',mpi_send_buf_dc=cp(np1:np2),mpi_number=nsend, & mpi_rank=i,mpi_comm=group_comm) call ms_mpi(mpi_command='bcast',mpi_send_buf_dc=cw(np1:np2),mpi_number=nsend, & mpi_rank=i,mpi_comm=group_comm) enddo if(rank.eq.0.and.iter.eq.1.and.writetime.eq.0) then time2=mytime()-time1 call timewrite(iunit,' time per iteration:',time2) writetime=1 endif if(rank.eq.0.and.iterwrite.eq.1) then write(iunit,'('' iter, err:'',i5,e13.5)') iter,err call flush(iunit) endif enddo ! ! arrive here with a converged solution ! do i=0,numprocs-1 ip1=mpi_sphere_index(i)+1 ip2=mpi_sphere_index(i)+mpi_sphere_number(i) np1=noff(ip1)+1 np2=noff(ip2)+nblk(ip2) nsend=np2-np1+1 call ms_mpi(mpi_command='bcast',mpi_send_buf_dc=anp(np1:np2),mpi_number=nsend, & mpi_rank=i,mpi_comm=group_comm) enddo end subroutine cbicg end module solver