www.pudn.com > Miescatter.rar > BNMIE2.F, change:2003-06-02,size:31192b
c ========================================================================
c This subroutine gives scattering matrix or its inverse
c for an anisotropic magnetic sphere
c in form of Ha and Hb, scattering matrices for each m
c Parameter:
c nc: maximum value of n input
c khab: input
c khab = 1, give scattering matirx for each m
c khab =-1, give inverse of sct matrix for each m
c xg: xg=kv*r0*dsqrt(em*um), with input
c r0: radius of sphere
c kv: wave vector in vacuum
c ms: relative refractive index of sphere ms=sqrt(us*es) input
c em: permittivity of medium -----
c um: permeability of medium -----
c es: permittivity of sphere / em -----
c us: \mu_s / um [cf. Eq.(4) in note] input
c urb: \mu'_r-1, [see Eq.(4) in kkrbn.tex for \mu'_r] input
c ukp: \mu'_k [cf. Eq.(4) in kkrbn.tex] input
c for permeability tensor see Eq.(3) and (4) in note
c ........................................................................
subroutine BnMie2(nc,khab,xg,ms,us,urb,ukp,Ha,Hb)
implicit none
integer*4 nc0, na0, nd0, nc02, nr0, nk0, nn0
include 'akkran2.g'
parameter (nd0=2*nc0*(nc0+2),nc02=nc0+nc0)
parameter (nn0=nc0*(nc02*nc0+6*nc0+1)/3)
integer*4 nc, khab, m, n1, nda, nm, i, j, ij, ns
integer*4 kTm(nc02)
real*8 D1x(nc0), xg
real*8 Gmn(nd0/2)
complex*16 z0, z1
complex*16 urb, ukp, us, ms
complex*16 D1z(nc0), D3x(nc0)
complex*16 Am(nc02), Ap(nc02)
complex*16 Hma(nc0,nc0), Hmb(nc0,nc0)
complex*16 Ham(nc0,nc0), Hbm(nc0,nc0)
complex*16 Snx(nc0), Snbx(nc0), Tnx(nc0), Tnbx(nc0)
complex*16 Ha(nn0), Hb(nn0)
complex*16 va(nc0,nc0), vb(nc0,nc0), wa(nc0), wb(nc0)
common/cGmn/Gmn
c ...........
c complex*16 ua(nc0,nc0), ub(nc0,nc0)
c equivalence (ua,Hma), (ub,Hmb)
c ...........
z0=(0.0d0,0.0d0)
z1=(1.0d0,0.0d0)
nda=2*nc*(nc+2)
call getD1x(nc,xg,D1x)
call getD3x(nc,xg,D3x)
call getSnx(nc,xg,D1x,D3x,Snx,Snbx)
c ---------------------------------------------------------
m=0
nm=nc
n1=nm+1
call getkTm(nm,kTm)
call getAmAp(nc,m,nm,Snx,Snbx,kTm,Am,Ap)
call getHab(nc,m,nm,urb,ukp,kTm,Hma,Hmb)
call getUmVm(nc,m,nm,xg,ms,us,D1x,kTm,Hma,Hmb,wa,wb,va,vb)
call getSctm(nm,khab,Am, Ap, Hma,va,Ham)
call getSctm(nm,khab,Am(n1),Ap(n1),Hmb,vb,Hbm)
ij=0
do i=1,nm
do j=1,nm
ij=ij+1
Ha(ij)=Ham(i,j)
Hb(ij)=Hbm(i,j)
enddo
enddo
ns=nm*nm
do m=1,nc
nm=nc+1-m
n1=nm+1
call getkTm(nm,kTm)
call getAmAp(nc,m,nm,Snx,Snbx,kTm,Am,Ap)
c .............
call getHab(nc,-m,nm,urb,ukp,kTm,Hma,Hmb)
call getUmVm(nc,-m,nm,xg,ms,us,D1x,kTm,Hma,Hmb,wa,wb,va,vb)
call getSctm(nm,khab,Am, Ap, Hma,va,Ham)
call getSctm(nm,khab,Am(n1),Ap(n1),Hmb,vb,Hbm)
ij=ns
do i=1,nm
do j=1,nm
ij=ij+1
Ha(ij)=Ham(i,j)
Hb(ij)=Hbm(i,j)
enddo
enddo
ns=ns+nm*nm
c .............
call getHab(nc,m,nm,urb,ukp,kTm,Hma,Hmb)
call getUmVm(nc,m,nm,xg,ms,us,D1x,kTm,Hma,Hmb,wa,wb,va,vb)
call getSctm(nm,khab,Am, Ap, Hma,va,Ham)
call getSctm(nm,khab,Am(n1),Ap(n1),Hmb,vb,Hbm)
ij=ns
do i=1,nm
do j=1,nm
ij=ij+1
Ha(ij)=Ham(i,j)
Hb(ij)=Hbm(i,j)
enddo
enddo
ns=ns+nm*nm
enddo
c ---------------------------------------------------------
return
end
c ========================================================================
c This sub gives scattering matrix for a given m
c ........................................................................
subroutine getSctm(nm,khab,Am,Ap,ua,va,Ha)
implicit none
integer*4 nc0, na0, nr0, nk0
include 'akkran2.g'
integer*4 ipiv(nc0), lwork, info
parameter (lwork=64*nc0)
integer*4 i, j, nm, khab
complex*16 Am(nm), Ap(nm), z0, z1, zv
complex*16 va(nc0,nc0), ua(nc0,nc0), Ha(nc0,nc0)
complex*16 Zm(nc0,nc0), Rm(nc0,nc0)
complex*16 work(lwork)
common/cwork/work,Zm,Rm
z0=(0.0d0,0.0d0)
z1=(1.0d0,0.0d0)
c ....................................................
c Zm = Z = (A-A')^{-1} (V-U)
c
do i=1,nm
zv=z1/(Am(i)-Ap(i))
do j=1,nm
Zm(i,j)=zv*(va(i,j)-ua(i,j))
enddo
enddo
c ....................................................
c ....................................................
c Rm = R^{-1} = U + A Z
c
do i=1,nm
zv=Am(i)
do j=1,nm
Rm(i,j)=ua(i,j) + zv*Zm(i,j)
enddo
enddo
c ....................................................
if(khab.gt.0) then
c ....................................................
c get inverse of matrix Rm such that Rm = R
c
call zgetrf(nm,nm,Rm,nc0,ipiv,info)
call zgetri(nm,Rm,nc0,ipiv,work,lwork,info)
c ....................................................
c scattering matrix Ha = S = Z R
c
call zgemm('N','N',nm,nm,nm,z1,Zm,nc0,Rm,nc0,z0,Ha,nc0)
c ....................................................
elseif(khab.lt.0) then
c ....................................................
c get inverse of matrix Zm such that Zm = Z^{-1}
c
call zgetrf(nm,nm,Zm,nc0,ipiv,info)
call zgetri(nm,Zm,nc0,ipiv,work,lwork,info)
c ....................................................
c Inverse of scattering matrix: Ha = S^{-1} = R^{-1} Z^{-1}
c
call zgemm('N','N',nm,nm,nm,z1,Rm,nc0,Zm,nc0,z0,Ha,nc0)
c ....................................................
endif
return
end
c ========================================================================
c This sub gives diagonal matrices Am, Ap
c ........................................................................
subroutine getAmAp(nc,m,nm,Snx,Snbx,kTm,Am,Ap)
implicit none
integer*4 nc0, na0, nr0, nk0, nd0
include 'akkran2.g'
parameter (nd0=nc0+nc0)
integer*4 j, n, nm, nc, m, ma, nda
integer*4 kTm(nd0)
complex*16 Am(nd0), Ap(nd0), At(nd0)
complex*16 Snx(nc0), Snbx(nc0)
ma=max(1,iabs(m))
nda=nm+nm
do n=ma,nc
j=n+1-ma
At(j )=Snbx(n)
At(j+nm)=Snx(n)
enddo
do j=1,nda
Am(j)=At(kTm(j))
enddo
do n=ma,nc
j=n+1-ma
At(j )=Snx(n)
At(j+nm)=Snbx(n)
enddo
do j=1,nda
Ap(j)=At(kTm(j))
enddo
return
end
c ========================================================================
c This sub gives U, U', V, and V' matrices (see Eqs.(39-42) in note)
c from eigensystem of Ha and Hb for a given m
c Parameter
c wa: eigenvalues of matirx Ha input
c va: eigenvectors of matirx Ha input
c wb: eigenvalues of matirx Hb input
c vb: eigenvectors of matirx Hb input
c
c ms: relative refractive index of sphere: ms=sqrt(us*es) input
c es: permittivity of sphere / em -----
c us: \mu_s / um [cf. Eq.(4) in note] input
c em: permittivity of medium -----
c um: permeability of medium -----
c xg: xg=kv*rs*dsqrt(em*um), with input
c rs: radius of sphere
c kv: wave vector in vacuum
c ........................................................................
subroutine getUmVm(nc,m,nm,xg,ms,us,D1x,kTm,ua,ub,wa,wb,va,vb)
implicit none
integer*4 nc0, na0, nr0, nk0, nd0
include 'akkran2.g'
parameter (nd0=nc0+nc0)
integer*4 ifail, lwork
parameter (lwork=64*nc0)
integer*4 i, j, nc, m, ma, n, nm, nda, i1
integer*4 kTm(nd0)
real*8 xg
real*8 D1x(nc0), rwork(nd0)
complex*16 us, ms, mskb, za, zb, zv, zg, z0, z1
complex*16 va(nc0,nc0), vb(nc0,nc0), wa(nc0), wb(nc0)
complex*16 ua(nc0,nc0), ub(nc0,nc0), D1z(nc0)
complex*16 Tnx(nc0), Tnbx(nc0)
complex*16 Rw(nd0,nd0)
complex*16 Um(nd0,nd0), Vm(nd0,nd0), Wm(nd0,nd0)
complex*16 Gtld(nc0,nc0), Etld(nc0,nc0), Ftld(nc0,nc0)
complex*16 Gbar(nc0,nc0), Ebar(nc0,nc0), Fbar(nc0,nc0)
complex*16 work(lwork)
common/defmtrx/Gtld,Gbar,Etld,Ebar,Ftld,Fbar
common/cwork/work,Rw,Um,Vm,Wm
z0=(0.0d0,0.0d0)
z1=(1.0d0,0.0d0)
ma=max(1,iabs(m))
do j=1,nm+nm
do i=1,nm+nm
Vm(i,j)=z0
Um(i,j)=z0
enddo
enddo
c ...................................................
c get eigensystems using NAG-18 subroutine f02gbf
c
ifail=0
call f02gbf('V',nm,ua,nc0,wa,va,nc0,rwork,work,lwork,ifail)
ifail=0
call f02gbf('V',nm,ub,nc0,wb,vb,nc0,rwork,work,lwork,ifail)
c ...................................................
c ...................................................
c get w_{mn,l}, [cf. equation after Eq.(32) in note]
c
do j=1,nm
do i=1,nm
Vm(i, j )=va(i,j)
Vm(i+nm,j+nm)=vb(i,j)
Um(i, j )=Ftld(i,j)
Um(i, j+nm)=Fbar(i,j)
enddo
enddo
nda=nm+nm
do i=1,nda
i1=kTm(i)
do j=1,nda
Rw(i1,j)=Vm(i,j)
enddo
enddo
call zgemm('N','N',nm, nda,nda,z1,Um,nd0,Rw,nd0,z0,Wm,nd0)
c ...................................................
c ...................................................
c get U, V and W matrices: Um = U, Vm = V, Wm = W
c
do i=1,nm
c ..............................
c for first nm eigenvalues
c
mskb=ms/cdsqrt(wa(i))
if(dimag(mskb).lt.0.0d0) mskb=-mskb
zg=xg*mskb
call getD1z(nc,zg,D1z)
call getTnx(nc,xg,zg,D1x,D1z,Tnx,Tnbx)
za=1.0d0/mskb
zb=wa(i)/us
do n=ma,nc
j=n+1-ma
zv=Tnbx(n)/(D1z(n)*zg*us)
Wm(j,i)=zv *Wm(j,i)
zv=Rw(j,i)
Um(j, i)=za*Tnbx(n)*zv
Vm(j, i)=zb*Tnx(n) *zv
zv=Rw(j+nm,i)
Um(j+nm,i)=za*Tnx(n) *zv
Vm(j+nm,i)=zb*Tnbx(n)*zv
enddo
c ..............................
c ..............................
c for last nm eigenvalues
c
mskb=ms/cdsqrt(wb(i))
if(dimag(mskb).lt.0.0d0) mskb=-mskb
zg=xg*mskb
call getD1z(nc,zg,D1z)
call getTnx(nc,xg,zg,D1x,D1z,Tnx,Tnbx)
za=1.0d0/mskb
zb=wb(i)/us
do n=ma,nc
j=n+1-ma
zv=Tnbx(n)/(D1z(n)*zg*us)
Wm(j,i+nm)=zv *Wm(j,i+nm)
zv=Rw(j,i+nm)
Um(j, i+nm)=za*Tnbx(n)*zv
Vm(j, i+nm)=zb*Tnx(n) *zv
zv=Rw(j+nm,i+nm)
Um(j+nm,i+nm)=za*Tnx(n) *zv
Vm(j+nm,i+nm)=zb*Tnbx(n)*zv
enddo
c ..............................
enddo
c ...................................................
c transform and decompose into two submatrices
c
do i=1,nda
i1=kTm(i)
do j=1,nda
Rw(i,j)=Vm(i1,j)
enddo
enddo
do j=1,nm
do i=1,nm
va(i,j)=Rw(i,j)
vb(i,j)=Rw(i+nm,j+nm)
enddo
enddo
do i=1,nda
i1=kTm(i)
do j=1,nda
Rw(i,j)=Um(i1,j)
enddo
enddo
do j=1,nm
do i=1,nm
ua(i,j)=Rw(i,j)
ub(i,j)=Rw(i+nm,j+nm)
enddo
enddo
do i=1,nm
do j=1,nm+nm
Wm(i+nm,j)=Wm(i,j)
Wm(i, j)=z0
enddo
enddo
do i=1,nda
i1=kTm(i)
do j=1,nda
Rw(i,j)=Wm(i1,j)
enddo
enddo
do j=1,nm
do i=1,nm
va(i,j)=va(i,j)+Rw(i,j)
vb(i,j)=vb(i,j)+Rw(i+nm,j+nm)
enddo
enddo
c ...................................................
return
end
c ========================================================================
c this sub decomposes a (2 nm) x (2 nm) eignenmatrix for a given m
c into two nm x nm eigenmatrices
c by similar transformation: kTm
c ........................................................................
subroutine getHab(nc,m,nm,urb,ukp,kTm,Ha,Hb)
implicit none
integer*4 nc0, na0, nr0, nk0, nd0
include 'akkran2.g'
parameter (nd0=nc0+nc0)
integer*4 nc, m, nm, i, j, i1, j1, nda
integer*4 kTm(nd0)
complex*16 urb, ukp
complex*16 Ha(nc0,nc0), Hb(nc0,nc0)
complex*16 Rw(nd0,nd0), Hm(nd0,nd0)
complex*16 Gtld(nc0,nc0), Etld(nc0,nc0), Ftld(nc0,nc0)
complex*16 Gbar(nc0,nc0), Ebar(nc0,nc0), Fbar(nc0,nc0)
common/defmtrx/Gtld,Gbar,Etld,Ebar,Ftld,Fbar
common/cwork/Rw,Hm
call gefset(nc,m,urb,ukp)
c .......................................................
c construct a (2 nm) x (2 nm) eignenmatrix
c from four matrices: Gtld, Gbar, Etld, Ebar
c
do j=1,nm
do i=1,nm
Hm(i, j )=Gtld(i,j)
Hm(i, j+nm)=Gbar(i,j)
Hm(i+nm,j )=Etld(i,j)
Hm(i+nm,j+nm)=Ebar(i,j)
enddo
enddo
c .......................................................
c .......................................................
c decompose into two nm x nm eigenmatrices
c
nda=nm+nm
do i=1,nda
i1=kTm(i)
do j=1,nda
Rw(i,j)=Hm(i1,j)
enddo
enddo
do j=1,nda
j1=kTm(j)
do i=1,nda
Hm(i,j)=Rw(i,j1)
enddo
enddo
do j=1,nm
do i=1,nm
Ha(i,j)=Hm(i,j)
Hb(i,j)=Hm(i+nm,j+nm)
enddo
enddo
c .......................................................
return
end
c ========================================================================
c this subroutine gives a transformation matrix for a given m such that
c a (2 nm) x (2 nm) eignensystem decomposed into two nm x nm ones
c ........................................................................
subroutine getkTm(nm,kTm)
implicit none
integer*4 nc0, na0, nr0, nk0
include 'akkran2.g'
integer*4 i,j,nm
integer*4 kTm(2*nc0)
do i=1,nm-1,2
kTm(i )=i
kTm(i+1 )=i+nm+1
kTm(i+nm )=i+nm
kTm(i+nm+1)=i+1
enddo
if(mod(nm,2).ne.0) then
kTm(nm )=nm
kTm(nm+nm)=nm+nm
endif
return
end
c ========================================================================
c This subroutine gives 6 matrices for a given m
c based on equations (18 - 19) and (28) in note kkrbn.tex
c Parameter:
c m: given value of m input
c urb, ukp: u'_r-1 and ukppa' in note kkrbn.tex input
c Gmn(j): Gmn(j)=sqrt[n(n+1)(2n+1)(n-m)!/(n+m)!] input
c with j=n*(n+1)+m
c Gtld etc: 6 matrices given in (28) for a given m output
c ........................................................................
subroutine gefset(nc,m,urb,ukp)
implicit none
integer*4 nc0, na0, nr0, nk0, nd0
include 'akkran2.g'
parameter (nd0=nc0*(nc0+2))
integer*4 m, n, i, j, n1, ma, nc, imn, iuv
real*8 Gmn(nd0), sg
complex*16 urb, ukp, zi, z0, tmp
complex*16 Gtld(nc0,nc0), Etld(nc0,nc0), Ftld(nc0,nc0)
complex*16 Gbar(nc0,nc0), Ebar(nc0,nc0), Fbar(nc0,nc0)
common/cGmn/Gmn
common/defmtrx/Gtld,Gbar,Etld,Ebar,Ftld,Fbar
zi=(0.0d0,1.0d0)
z0=(0.0d0,0.0d0)
ma=max(1,iabs(m))
do j=1,nc0
do i=1,nc0
Gtld(i,j)=z0
Etld(i,j)=z0
Ftld(i,j)=z0
Gbar(i,j)=z0
Ebar(i,j)=z0
Fbar(i,j)=z0
enddo
enddo
do n=ma,nc
n1=n+1
j=n1-ma
i=j
Gtld(i,j)=1.0d0+((n1*n-m*m)*urb+m*ukp)/(n*n1)
enddo
do n=ma+1,nc
n1=n+1
j=n1-ma
i=j-1
imn=n*(n+1)+m
iuv=(n-1)*n+m
sg=(iuv-m)*Gmn(imn)/((imn-m)*Gmn(iuv))
Etld(i,j)=-sg*(n+m)*(m*urb-n1*ukp)/(n*(n+n1))
enddo
do n=ma,nc-1
n1=n+1
j=n1-ma
i=j+1
imn=n*(n+1)+m
iuv=(n+1)*(n+2)+m
sg=(iuv-m)*Gmn(imn)/((imn-m)*Gmn(iuv))
Etld(i,j)=sg*(n1-m)*(m*urb+n*ukp)/(n1*(n+n1))
enddo
do n=ma+1,nc
n1=n+1
j=n1-ma
i=j-1
imn=n*(n+1)+m
iuv=(n-1)*n+m
sg=(iuv-m)*Gmn(imn)/((imn-m)*Gmn(iuv))
Ftld(i,j)=sg*(n+m)*(m*urb-n1*ukp)/(n+n1)
enddo
do n=ma,nc-1
n1=n+1
j=n1-ma
i=j+1
imn=n*(n+1)+m
iuv=(n+1)*(n+2)+m
sg=(iuv-m)*Gmn(imn)/((imn-m)*Gmn(iuv))
Ftld(i,j)=sg*(n1-m)*(m*urb+n*ukp)/(n+n1)
enddo
do n=ma+1,nc
n1=n+1
j=n1-ma
i=j-1
imn=n*(n+1)+m
iuv=(n-1)*n+m
sg=(iuv-m)*Gmn(imn)/((imn-m)*Gmn(iuv))
Gbar(i,j)=sg*(n+m)*n1*(m*urb+(n-1)*ukp) /(n*(n-1)*(n+n1))
enddo
do n=ma,nc-1
n1=n+1
j=n1-ma
i=j+1
imn=n*(n+1)+m
iuv=(n+1)*(n+2)+m
sg=(iuv-m)*Gmn(imn)/((imn-m)*Gmn(iuv))
Gbar(i,j)=-sg*(n1-m)*n*(m*urb-(n+2)*ukp) /(n1*(n+2)*(n+n1))
enddo
do n=ma,nc
n1=n+1
j=n1-ma
i=j
tmp=(m*m*(2*n*n1+3)+n*n1*(2*n*n1-3))*urb+(m*(4*n*n1-3))*ukp
Ebar(i,j)=1.0d0+tmp/(n*n1*(n+n-1)*(n+n+3))
enddo
do n=ma+2,nc
n1=n+1
j=n1-ma
i=j-2
imn=n*(n+1)+m
iuv=(n-2)*(n-1)+m
sg=(iuv-m)*Gmn(imn)/((imn-m)*Gmn(iuv))
Ebar(i,j)= sg*(n+m-1)*(n+m)*n1*urb/((n-1)*(n+n-1)*(n+n1))
enddo
do n=ma,nc-2
n1=n+1
j=n1-ma
i=j+2
imn=n*(n+1)+m
iuv=(n+2)*(n+3)+m
sg=(iuv-m)*Gmn(imn)/((imn-m)*Gmn(iuv))
Ebar(i,j)= sg*(n1-m)*(n+2-m)*n*urb/((n+2)*(n+n+3)*(n+n1))
enddo
do n=ma,nc
n1=n+1
j=n1-ma
i=j
tmp=(n*n1-3*m*m)*urb-m*(n+n-1)*(n+n+3)*ukp
Fbar(i,j)=-tmp/((n+n-1)*(n+n+3))
enddo
do n=ma+2,nc
n1=n+1
j=n1-ma
i=j-2
imn=n*(n+1)+m
iuv=(n-2)*(n-1)+m
sg=(iuv-m)*Gmn(imn)/((imn-m)*Gmn(iuv))
Fbar(i,j)=-sg*(n+m-1)*(n+m)*n1*urb/((n+n-1)*(n+n1))
enddo
do n=ma,nc-2
n1=n+1
j=n1-ma
i=j+2
imn=n*(n+1)+m
iuv=(n+2)*(n+3)+m
sg=(iuv-m)*Gmn(imn)/((imn-m)*Gmn(iuv))
Fbar(i,j)= sg*(n1-m)*(n+2-m)*n*urb/((n+n+3)*(n+n1))
enddo
return
end
c ========================================================================
c The sub gives logrithmic derivative of Ricatti-Bessel function x j_n(x)
c with real argument
c Parameter:
c nv: maximum order of n input
c x: real argument input
c D1x: logrithmic derivatives of Ricattic-Bessel output
c D1x(n) = [x j_{n}(x)]'/[x j_{n}(x)]
c ........................................................................
subroutine getD1x(nv,x,D1x)
implicit none
integer*4 nv,n,nst
real*8 x, ax, D1x(nv), tmp
c ...................................................
c Downward recurrence for D1x
c
ax = dabs(x)
nst = int( dsqrt( 101.0d0 + ax ) )
nst = nst + ax + 4.05*ax**(1.0d0/3.0d0) + 2
nst = max(nst,nv) + 15
tmp=0.0d0
do n=nst,nv+1,-1
tmp=tmp+dble(n)/x
tmp=dble(n)/x-1.0d0/tmp
enddo
D1x(nv)=tmp
do n=nv,2,-1
tmp=tmp+dble(n)/x
tmp=dble(n)/x-1.0d0/tmp
D1x(n-1)=tmp
enddo
c ...................................................
return
end
c ========================================================================
c The sub gives logrithmic derivative of Ricatti-Bessel function x h_n(x)
c with real argument, where h_n = j_n + i y_n
c Parameter:
c nv: maximum order of n input
c x: real argument input
c D3x: logrithmic derivatives of Ricattic-Bessel output
c D3x(n) = [x h_{n}(x)]'/[x h_{n}(x)]
c ........................................................................
subroutine getD3x(nv,x,D3x)
implicit none
integer*4 nv,n
real*8 x
complex*16 D3x(nv), ztmp
c ...................................................
c Upward recurrence for D3x
c
ztmp = (0.0d0,1.0d0)
do n=1,nv
ztmp=dble(n)/x - ztmp
ztmp=1.0d0/ztmp-dble(n)/x
D3x(n)=ztmp
enddo
c ...................................................
return
end
c ========================================================================
c The sub gives logrithmic derivative of Ricatti-Bessel function z j_n(z)
c with complex argumentc
c Parameter:
c nv: maximum order of n input
c zx: complex argument input
c D1z: logrithmic derivatives of Ricattic-Bessel output
c D1z(n) = [z j_{n}(z)]'/[z j_{n}(z)]
c ........................................................................
subroutine getD1z(nv,zx,D1z)
implicit none
integer*4 nv,n,nst
real*8 ax
complex*16 zx, D1z(nv), ztmp
c ...................................................
c Downward recurrence for D1z
c
ax = cdabs(zx)
nst = int( dsqrt( 101.0d0 + ax ) )
nst = nst + ax + 4.05*ax**(1.0d0/3.0d0) + 2
nst = max(nst,nv) + 15
ztmp=(0.0d0,0.0d0)
do n=nst,nv+1,-1
ztmp=ztmp+dble(n)/zx
ztmp=dble(n)/zx-1.0d0/ztmp
enddo
D1z(nv)=ztmp
do n=nv,2,-1
ztmp=ztmp+dble(n)/zx
ztmp=dble(n)/zx-1.0d0/ztmp
D1z(n-1)=ztmp
enddo
c ...................................................
return
end
c ========================================================================
c The sub gives logrithmic derivative of Ricatti-Bessel function z h_n(z)
c with complex argument, where h_n = j_n + i y_n
c Parameter:
c nv: maximum order of n input
c zx: complex argument input
c D3z: logrithmic derivatives of Ricattic-Bessel output
c D3z(n) = [z h_{n}(z)]'/[z h_{n}(z)]
c ........................................................................
subroutine getD3z(nv,zx,D3z)
implicit none
integer*4 nv,n
complex*16 zx, D3z(nv), ztmp
c ...................................................
c Upward recurrence for D3z
c
ztmp = (0.0d0,1.0d0)
do n=1,nv
ztmp=dble(n)/zx - ztmp
ztmp=1.0d0/ztmp-dble(n)/zx
D3z(n)=ztmp
enddo
c ...................................................
return
end
c ========================================================================
c The sub gives Sn(x) and Snb(x) in equation (48) in note
c Parameter:
c nv: maximum order of n input
c D1x(n): D1x(n) = [x j_{n}(x)]'/[x j_{n}(x)] input
c D3x(n): D3x(n) = [x h_{n}(x)]'/[x h_{n}(x)] input
c Snx(n): Snx(n) = [x h_{n}(x)] /[x j_{n}(x)] output
c Snbx(n): Snbx(n) = [x h_{n}(x)]'/[x j_{n}(x)]' output
c ........................................................................
subroutine getSnx(nv,x,D1x,D3x,Snx,Snbx)
implicit none
integer*4 nv, n
real*8 x, D1x(nv)
complex*16 zs, zi
complex*16 Snx(nv), Snbx(nv), D3x(nv)
zi =(0.0d0,1.0d0)
zs=cdexp(zi*(x+x))
zs=(zs+zs)/(zs-1.0d0)
do n=1,nv
zs=zs*(D1x(n)+dble(n)/x)/(D3x(n)+dble(n)/x)
Snx(n)=zs
Snbx(n)=zs*D3x(n)/D1x(n)
enddo
return
end
c ========================================================================
c The sub gives Tn(x,z) and Tnb(x,z) in equation (48) in note
c Parameter:
c nv: maximum order of n input
c D1x(n): D1x(n) = [x j_{n}(x)]'/[x j_{n}(x)] input
c D1z(n): D1z(n) = [z j_{n}(z)]'/[z j_{n}(z)] input
c Tnx(n): Tnx(n) = [z j_{n}(z)] /[x j_{n}(x)] output
c Tnbx(n): Tnbx(n) = [z j_{n}(z)]'/[x j_{n}(x)]' output
c ........................................................................
subroutine getTnx(nv,x,z,D1x,D1z,Tnx,Tnbx)
implicit none
integer*4 nv, n
real*8 x, D1x(nv)
complex*16 zt, zs, zi, z
complex*16 Tnx(nv), Tnbx(nv), D1z(nv)
zi =(0.0d0,1.0d0)
zs=cdexp(zi*x)
zt=cdexp(zi*z)
c ........................................................
c From the definition of Tnx and Tnbx
c one should start with the following initial value
c
c zt=(zt*zt-1.0d0)/(zs*zs-1.0d0) * zs/zt
c ........................................................
c However, for light scattering by gyrotropic sphere
c one can start with the following initial value
c without changing the results of field OUTSIDE sphere
c and avoiding overfloat for metal sphere
c Lin: May 20, 2003
c
zt=(zt*zt-1.0d0)/(zs*zs-1.0d0) * zs
c ..............................................
do n=1,nv
zt=zt*(D1x(n)+dble(n)/x)/(D1z(n)+dble(n)/z)
Tnx(n)=zt
Tnbx(n)=zt*D1z(n)/D1x(n)
enddo
return
end
c ========================================================================
c this subroutine gives a transformation matrix that transforms inverse
c of scattering matrix of (m,n) form into conventional (n,m) form
c
c (m,n) form: two indices (m,n) are maped to a single index by imn(m,n,nc)
c m=0,-1,1,-2,2,...,-nc,nc, while
c n=max(1,|m|), ..., nc for a given value of m
c scattering matrix (and its inverse) is block-diagonal in (m,n) form
c
c (n,m) form: two indices (n,m) are maped to a single index by jmn(m,n)
c n=1,2,...,nc, while
c m=-1,-n+1, ..., n-1,n for a given value of n
c
c see also integer functions imn(m,n,nc) and jmn(m,n)
c ........................................................................
subroutine getkTt(nc,ndh,kTt)
implicit none
integer*4 nc0, na0, nd0, nr0, nk0
include 'akkran2.g'
parameter (nd0=2*nc0*(nc0+2))
integer*4 i,j,m,n,nm,nc,ndh,ns
integer*4 kTt(nd0)
integer*4 imn, jmn
external imn, jmn
m=0
nm=nc
do n=1,nc,2
i=imn(m,n,nc)
j=jmn(m,n)
kTt(i )=j
kTt(i+nm)=j+ndh
enddo
do n=2,nc,2
i=imn(m,n,nc)
j=jmn(m,n)
kTt(i )=j+ndh
kTt(i+nm)=j
enddo
ns=nm
do m=1,nc
nm=nc+1-m
do n=m,nc,2
i=imn(-m,n,nc)
j=jmn(-m,n)
kTt(i+ns )=j
kTt(i+ns+nm)=j+ndh
enddo
do n=m+1,nc,2
i=imn(-m,n,nc)
j=jmn(-m,n)
kTt(i+ns )=j+ndh
kTt(i+ns+nm)=j
enddo
ns=ns+nm
do n=m,nc,2
i=imn(m,n,nc)
j=jmn(m,n)
kTt(i+ns )=j
kTt(i+ns+nm)=j+ndh
enddo
do n=m+1,nc,2
i=imn(m,n,nc)
j=jmn(m,n)
kTt(i+ns )=j+ndh
kTt(i+ns+nm)=j
enddo
ns=ns+nm
enddo
return
end
c ========================================================================
integer*4 function isgn(m)
integer*4 m
if(m.gt.0) then
isgn=1
else
isgn=0
endif
return
end
c ========================================================================
c two indices (m,n) are maped to a single index by imn(m,n,nc) such that
c m=0,-1,1,-2,2,...,-nc,nc, while
c n=max(1,|m|), ..., nc for a given value of m
c scattering matrix (and its inverse) is block-diagonal in this form
c ........................................................................
integer*4 function imn(m,n,nc)
integer*4 m,n,ma,nc
ma=iabs(m)
imn=ma*(nc+nc+3-ma)+isgn(m)*(nc-ma+1)+n-ma-isgn(ma)*(nc+1)
return
end
c ========================================================================
c two indices (n,m) are maped to a single index by jmn(m,n) such that
c n=1,2,...,nc, while
c m=-1,-n+1, ..., n-1,n for a given value of n
c ........................................................................
integer*4 function jmn(m,n)
integer*4 m,n
jmn=n*(n+1)+m
return
end