C (C) Copr. 06/2005 Alexander Moroz C Wave-scattering.com C http://www.wave-scattering.com/ C wavescattering@yahoo.com C C Last update: January 2002 C C This is OPTTRAP, a F77 code which calculates the optical axial trapping C efficiency Q for a multilayered/coated micro/nano/sphere. C The program was made according to the theory of C P. A. Maia Neto and H. M. Nussenzveig, Theory of optical tweezers, C Europhys. Lett. 50, 702-708 (2000) by extending it to the case of C (multi)coated spheres. C C The axial trapping efficiency Q, or the dimensionless normalized C axial force, when multiplied by the focused laser beam power P, determines C the resulting axial force F on a sphere immersed in a homogeneous medium C with refractive index n according to the formula C F=(n/c) PQ C For the results obtained by this program, see C C A. Moroz, Metallo-dielectric diamond and zinc-blende photonic structures, C Phys. Rev. B 66, 115109 (2002) C (http://www.arXiv.org/abs/cond-mat/0209188) C (see discussion subsection "Fabrication" with Fig. 23 therein) C C and C C A. van der Horst, P. D. J. van Oostrum, A. Moroz, A. van Blaaderen, and C M. C. Dogterom, C High-refractive index particles trapped in dynamic arrays of dual-beam C optical tweezers, Appl. Opt. 47(17), 3196-3202 (2008) C (http://ao.osa.org/abstract.cfm?uri=ao-47-17-3196) C C ======================================================================== C DEPENDENCES: c C Requires in total the following routines: C C opttrap.f zkmatl.f plgndr.f qromb.f polint.f trapzd.f C gnzbess.f znsrefind.f zbessf.f bessjy.f bessik.f beschb.f chebev.f C C See, e.g., a makefile "maketrapdv" for UNIX/LINUX compilers, which C for some computer security reasons, which do not allow to distribute C *.dat files and files without an extenmsion, can be downloaded as C http://www.wave-scattering.com/maketrapdv.txt C (After a download it has to be renamed to maketrapdv.) C Check your local compiler options and, if necessary, amend the makefile C appropriately. C Save all the required files within the same directory and issue C the command C C make -f maketrapdv C C to compile the source code. C C The present file comprises opttrap.f and zkmatl.f C C The supplementary routines C C gnzbess.f znsrefind.f zbessf.f bessjy.f bessik.f beschb.f chebev.f C C are a part of the Bessel function package bessel.zip available at C http://www.wave-scattering.com/bessel.zip C C The gold data file "Audat.dat" (the option NMAT=3) can be downloaded at C http://www.wave-scattering.com/Audat.dat.txt C For security reasons, *.dat files cannot be posted; therefore C the additional *.txt extension. You have to rename the data file C back to Audat.dat after the download. C C The supplementary routines C C plgndr.f qromb.f polint.f trapzd.f C C are not distributed here. They are copyrighted routines from C Numerical Recipes Software (I have not C received any reply yet from Numerical Recipes Software on my request C for a permission to distribute the routines): C C qromb(func,a,b,ss), C which returns as ss the integral of the function func from a to b, C together with its dependencies polint.f and trapzd.f and C C plgndr(l,m,x) C which returns associated Legendre polynomials C (See http://www.numerical-recipes.com/recipe-list.html C for the listing of all the Numerical Recipes routines) C C All the routines can be found on C http://www.nr.com C C IMPORTANT!!! C The Numerical Recipes routines have to be transformed into C double precision, i.e. you have to change C change declaration real to real*8 and complex to complex*16. C ======================================================================== C OUTPUT: C C OPTTRAP returns file C C optforce.dat containing the axial force on a single (coated) sphere C======================================================================== C C This code is made available on the following conditions: C C - You can use and modify it,as long as acknowledgment is given C to the author. I'd appreciate a reprint or e-copy C of any paper where my code has been used. C C - You can not sell it or use it in any way to get money out of it C (save for the research grants you might get thanks to it) C C - I'm making this code available to the best of my knowledge. C It has been tested as throughly as possible. C It work for Windows and Linux/Unix platforms. However, it's C your responsibility to check for its accuracy in the environment C and size/shape ranges you use. C C DISCLAIMER - The usual stuff: C C WHILE THE COMPUTER PROGRAMS HAVE BEEN TESTED FOR A VARIETY OF CASES, C IT IS NOT INCONCEIVABLE THAT THEY MAY CONTAIN UNDETECTED ERRORS. ALSO, C INPUT PARAMETERS CAN BE USED WHICH ARE OUTSIDE THE ENVELOPE OF C VALUES FOR WHICH RESULTS ARE COMPUTED ACCURATELY. FOR THIS REASON, C THE AUTHOR DISCLAIM ALL LIABILITY FOR ANY DAMAGES THAT MAY RESULT C FROM THE USE OF THE PROGRAMS. C C I would highly appreciate informing me of any problems encountered C with this code. Improvements, comments, additions: all welcome at C C wavescattering@yahoo.com C C ===================================================================== program opttrap C--------/---------/---------/---------/---------/---------/---------/-- C This routines calculates the resulting axial trapping force excerted C on a (including coated) sphere in a focused laser beam C C make -f maketrap C C gamma the ratio of the objective focal length to the beam waist C gamma2= gamma**2 C qpar with origin at the sphere center, - qpar*z is the focal C point position C theta0 opening angle of a focused beam beyond the objective C wvno 2*pi/lambda laser wavenumber C cnh the host refractive index C cns the sphere refractive index C force the resulting axial trapping force excerted on a sphere C--------/---------/---------/---------/---------/---------/---------/-- implicit none integer LMAXV,LCS,ILCS,NMAT integer NOUT,NFIN real*8 rmuf,tol,sol,pi COMPLEX*16 ZEPS0,CCEPS,CSEPS character*1 ync c Parameters: C ::: number of the output unit PARAMETER (NOUT=35) PARAMETER (PI=3.141592653589793d0) c the speed of light in vacuum in m/s PARAMETER (sol=2.9979245812d8) c number of spherical harmonics used PARAMETER (lmaxv=45) c material code number PARAMETER(NMAT=2) c NMAT=0 dispersionless dielectric c NMAT=1 Drude metal c NMAT=2 Ag c NMAT=3 Au c NMAT=4 ZnS c Temporarily option for reading of the real data for the dielectric constant c The number of the entries in a material data file to be read below c AGC.DAT NFIN=73 ! from Palik c Audat.dat NFIN=65 ! from Palik c Au_2dat.dat NFIN=76 ! from JAW c Au*new.dat NFIN=142 c PARAMETER (NFIN=73) * c sphere radius PARAMETER (rmuf=1.d0) c If sphere is coated, ync=y, otherwise ync=n parameter (ync='y') c number of coatings parameter (lcs=2) c The coating layer to which material data are read in parameter (ilcs=1) c background (host) dielectric constant PARAMETER (ZEPS0=(1.33D0,0.d0)**2) c sphere/core dielectric constant PARAMETER (CCEPS=(2.00D0,0.d0)**2) c shell dielectric constant PARAMETER (CSEPS=(1.45D0,0.d0)**2) * n(silica)=1.45 <---> EPS(1)=2.1025D0 * n(ZnS)=2. <---> EPS(1)=4.D0 C ::: relative error allowed for the TCS. If the convergence * within TOL is not reached, program issues warning PARAMETER (TOL=1.d-6) * integer nstep,l,ipar,npar,jq,ieps real*8 rsnm,lambda,rff,rmf(lcs) real*8 enw,xstep real*8 qpar,xs real*8 omf(NFIN),omxf,reepsz,plasma,omxp,filfrac real*8 delo,omega0,omega,xres1,xres2,xdel real*8 force,gamma2,theta0,wvno,pow,ama,qe,qs complex*16 ci,CQE1,CQE2,CQS,zeps(lcs) COMPLEX*16 KMT(2,lmaxv),cgl(lmaxv),cglp(lmaxv) COMPLEX*16 alpha(lmaxv),beta(lmaxv) c Quantities for the use of the material data COMPLEX*16 CEPS1(NFIN),zeps1 external ZKMATL,funki1,funki2,funki3,funki4 COMMON/BDD7/gamma2,xdel COMMON/IDD7/l C--------/---------/---------/---------/---------/---------/---------/-- * Checking set up: if (nmat.gt.1) write(6,*)'Real material data are to be provided' if (ync.eq.'y' .and. lcs.eq.1) then write(6,*)'Check compatibility of YNC and LCS' stop end if C-------------------------------------------------------------------- * Reading in the input data * DATA ci/(0.d0,1.d0)/ * zeps(1)=cceps if (lcs.gt.1) then zeps(lcs)=cseps c write(6,*)'Read core to sphere radius ratio' c read(5,*) rmf(1) end if * rmf(lcs)=rmuf * interactively: write(6,*)' BEAM PARAMETERS' c write(6,*)'The focused laser beam power=' c read(5,*) pow pow=1.d0 c write(6,*)'Beam opening angle theta0 (in degrees)=' c read(5,*) theta0 theta0=78.d0 theta0=theta0*pi/180.d0 c write(6,*) 'sin(theta0)**2=',sin(theta0)**2 c write(6,*) 'sin(theta0)**2=',(sin(theta0))**2 C--------/---------/---------/---------/---------/---------/---------/-- c write(6,*)'The ratio of the objective focal length to the beam c & waist=' c read(5,*) gamma2 c in the following only gamma2=gamma**2 is used c gamma2=gamma2**2 gamma2=0.99d0 * write(6,*)' SPHERE PARAMETERS' write(6,*)'Read sphere radius rsnm in nm' read(5,*) rsnm c rsnm=700 if(lcs.eq.2) then c write(6,*)'Read core radius in nm' c read(5,*) rff c rff=rff/rsnm rff=0.3d0 rmf(1)=rff*rmuf end if * c write(6,*)' SCANNING' c write(6,*)'Read initial (maximal) lambda (in the vacuum) in nm' c read(5,*) lambda lambda=1000.0d0 * * size parameter is customarily defined as the ratio of * circumference of sphere to the wavelength in the host medium * in which the sphere is embedded * x=kr=\sg a=2*pi*r/lambda * xs=2.d0*pi*rsnm*dble(sqrt(zeps0))/lambda * convert lambda to the lambda in vacuum: c lambda=lambda*dble(sqrt(zeps0)) c omega=2.d0*pi*rsnm/(lambda*rmuf)=xs/(rmuf*dble(sqrt(zeps0))), c where rsnm is the sphere radius (in nm) and lambda is the wavelengths (in nm) c in the vacuum: omega=2.d0*pi*rsnm/(lambda*rmuf) omega0=omega wvno=2.d0*pi*dble(sqrt(zeps0))/lambda WRITE(6,*)'omega=', omega * * Option for omega input: c write(6,*)'Read omega =' c read(5,*) omega * xs=RMUF*omega*dble(sqrt(zeps0)) * write(6,*)'Size parameter x=2*pi*rsnm*n_0/lambda=', xs if(xs.gt.30.d0) then write(6,*)'Convergence ensured only for size parameters <=30!!!' write(6,*)'Raise LMAXV so that LMAXV > xs+4.*xs^{1/3}+2' stop 777 end if * c write(6,*)'Scan down to wavelength (in nm)' c read(5,*) enw enw=1000.0d0 c write(6,*)'Scanning step in nm' c read(5,*) xstep xstep=1 C ::: number of steps on frequency interval: if(lambda.eq.enw) then NSTEP=0 delo=0.d0 else NSTEP=(lambda-enw)/xstep C ::: width of the searched frequency interval ENW=2.d0*pi*rsnm/(enw*rmuf) enw=enw-omega0 delo=enw/dble(nstep) end if c write(6,*)'Scan axial force profile till nm' c read(5,*) npar npar=5000 c write(6,*)'with step in nm' c read(5,*) ipar ipar=100 * -------------------------------- * output initial statements OPEN(UNIT=NOUT,FILE='optforce.dat') rewind(NOUT) WRITE(NOUT,*)'#Axial force on a single (coated) sphere' write(nout,*)'#Sphere radius in nm=', rsnm write(nout,*) if (ync .eq.'n') then write(nout,*)'#Homogeneous sphere' else if (ync.eq.'y') then write(nout,*)'#coated sphere' end if write(nout,*) write(nout,*)'#host dielectric constant=', zeps0 if (ync.eq.'n') write(nout,*)'#sphere diel. constant=', zeps(1) if (ync.eq.'y') write(nout,*)'#core diel. constant=', zeps(1) if (ync.eq.'y') write(nout,*)'#coating diel. constant=',zeps(lcs) if (ync.eq.'y') write(nout,*)'#core radius/sphere radius =', rff write(nout,*) write(nout,*)'#In columns, the center offset from the focus in' write(nout,*)'#units q/rsnm and the axial force per power of 1W' *--------/---------/---------/---------/---------/---------/---------/-- * Sphere optical constants in the case of a dispersion: if (nmat.eq.0) goto 20 ! dispersionless dielectric C else > READING IN MATERIAL DATA: * Reading real silver data according to Palik's book * requires reading data files OMF and CEPS1 of dimension NFIN * OMF is reepsz/omega and CEPS1 contains the sphere EPS * if (nmat.eq.2) then write(6,*)'Silver particles' OPEN(UNIT=30,FILE='agc.dat') !silver data write(nout,*)'#Palik data used' write(6,*)'Palik data used' rewind(30) do ieps=1,nfin read(30,*) omf(ieps),ceps1(ieps) enddo close(30) * else if (nmat.eq.3) then * write(6,*)'Gold particles' c OPEN(UNIT=30,FILE='Au293Knew.dat') !Gold data for different T c OPEN(UNIT=30,FILE='Audat.dat') !Gold data (Palik) in nm c write(nout,*)'#Palik data used' c write(6,*)'Palik data used' OPEN(UNIT=30,FILE='Au_2dat.dat') !Gold data (JAW) in nm write(nout,*)'#JAW data used' write(6,*)'JAW data used' rewind(30) do ieps=1, nfin read(30,*) omf(ieps),ceps1(ieps) c omf(ieps)=2.d0*pi*rsnm*omf(ieps)/(1240.d0*rmuf) omf(ieps)=2.d0*pi*rsnm/(omf(ieps)*rmuf) enddo close(30) end if * * In case of a dispersion, EPSSPH is modified. * For ideal Drude metal * plasma=2.d0*pi*sphere radius in nm/(lambda_z in nm*rmuf) * where lambda_z is the wavelength for which Re eps_s=0. if ((nmat.eq.1).or.(nmat.eq.2)) & reepsz=2.d0*pi*RSNM/(324.269d0*rmuf) *--------/---------/---------/---------/---------/---------/---------/-- if (nmat.eq.1) then reepsz=2.d0*pi*RSNM/(324.269d0*rmuf) plasma=reepsz omxp=plasma/omega zeps1=1.d0-omxp**2/(1.d0+ci*plasma/(144.d0*omega)) go to 18 * else if (nmat.eq.2) then !silver * c >>> real material data: !silver * lambda_z=324.2693391740981d0 * lambda_p=164.d0 * When real material data are used, * reepsz differs from plasma!!! The plasma wavelength is * calculated below: plasma=reepsz*7.2d0/3.8291d0 * security trap - remainder (not optimized!) omxf=omega/reepsz if (omxf.gt.omf(1)) then write(6,*)'Calculation of has to stop with' write(6,*)' OMF(1)' write(6,*)' OMXF=', omxf stop end if if (omxf.lt.omf(nfin)) then omxp=plasma/omega zeps1=1.d0-omxp**2/(1.d0+ci*plasma/(144.d0*omega)) * damping coefficient for silver is plasma/144 where plasma is different from * the Re eps zero crossing at 3.8291 eV according to Palik!!! go to 18 else if (omxf.eq.omf(1)) then zeps1=ceps1(1) go to 18 else do ieps=2,nfin * data file ordered with the increased wavelength * omxf increases in the loop and is oriented opposite to the data file if (omxf.gt.omf(ieps)) then ! linear interpolation zeps1=ceps1(ieps)+(omxf-omf(ieps))*(ceps1(ieps-1)-ceps1(ieps)) 1 /(omf(ieps-1)-omf(ieps)) go to 18 end if enddo end if * else if (nmat.eq.3) then ! Au * c >>> real gold data: * data file ordered with the decreased wavelength * omega increases in the loop and is oriented along the data file * if ( (omega.lt.omf(1)).or.(omega.gt.omf(nfin)) ) then write(6,*)'Material data not available for this wavelength' stop end if * if (omega.eq.omf(nfin)) then zeps1=ceps1(nfin) go to 18 else do ieps=1,nfin-1 if (omega.lt.omf(ieps+1)) then ! linear interpolation zeps1=ceps1(ieps)+(omega-omf(ieps))*(ceps1(ieps+1)-ceps1(ieps)) 1 /(omf(ieps+1)-omf(ieps)) go to 18 end if enddo end if * else if (nmat.eq.4) then write(nout,*)'#Palik data used' write(6,*)'Palik data used' * filfrac=0.62d0 ! filfrac of ZnS in ZnS core call znsrefind(LAMBDA,FILFRAC,zeps1) go to 18 * end if * *The end of reading real data according to Palik's book * 18 zeps(ilcs)=zeps1 !ovewriting zeps(ilcs) with material data * *--------/---------/---------/---------/---------/---------/---------/-- * Execution: * * Calculating phase-shifts: * ZKMATL for coated and * KPCSMAT for perfectly conducting sphere !!! * 20 CALL ZKMATL(lcs,lmaxv,omega,zeps,zeps0,rmf,kmt) * * CALL KPCSMAT(lmaxv,rmuf(IB),omega,eps0,rmf,kmt) C-------------------------------------------------------------------- C Field KMT contains - tan(phase shift) C-------------------------------------------------------------------- ama=1.d0-exp(-2.d0*gamma2*sin(theta0)**2) write(6,*)'ama=',ama ama=4.d0*gamma2/ama * do 100 jq=-npar,npar,ipar qpar=dble(jq) xdel=qpar*wvno * * summation over partial wave contributions: * do 40 l=1,lmaxv * alpha(l)=-kmt(1,l)/(ci-kmt(1,l)) beta(l)=-kmt(2,l)/(ci-kmt(2,l)) * Calculating gl fields: call qromb(funki1,0.d0,theta0,xres1) call qromb(funki2,0.d0,theta0,xres2) cgl(l)=dcmplx(xres1,xres2) call qromb(funki3,0.d0,theta0,xres1) call qromb(funki4,0.d0,theta0,xres2) cglp(l)=dcmplx(xres1,xres2) * 40 continue * * cqs=dcmplx(0.d0,0.d0) cqe1=dcmplx(0.d0,0.d0) cqe2=dcmplx(0.d0,0.d0) * do 50 l=1,lmaxv cqs=cqs+dble(2*l+1)*(alpha(l)+beta(l))*cgl(l)*dconjg(cglp(l)) cqe2=cqe2+dble(2*l+1)*alpha(l)*dconjg(beta(l))*cgl(l)* & dconjg(cgl(l))/dble(l*(l+1)) if(l.eq.lmaxv) go to 50 cqe1=cqe1+dble(l*(l+2))*(alpha(l)*dconjg(alpha(l+1))+ & beta(l)*dconjg(beta(l+1)))*cgl(l)*dconjg(cgl(l+1))/dble(l+1) 50 continue qs=ama*dble(cqs) qe=-2.d0*ama*dble(cqe1+cqe2) force=sqrt(zeps0)*pow*(qs+qe)/sol write(nout,*) qpar/rsnm, qs+qe 100 continue * write(6,*)'Collect results from the file optforce.dat' * *--------/---------/---------/---------/---------/---------/---------/-- end * function funki1(x) implicit none integer l real*8 funki1,x,xdel,gamma2,ds1 external ds1 COMMON/BDD7/gamma2,xdel COMMON/IDD7/l funki1=sin(x)*sqrt(cos(x))*exp(-gamma2*sin(x)**2) & *cos(xdel*cos(x))*ds1(l,x) return end * function funki2(x) implicit none integer l real*8 funki2,x,xdel,gamma2,ds1 external ds1 COMMON/BDD7/gamma2,xdel COMMON/IDD7/l funki2=sin(x)*sqrt(cos(x))*exp(-gamma2*sin(x)**2) & *sin(xdel*cos(x))*ds1(l,x) return end * function funki3(x) implicit none integer l real*8 funki3,x,xdel,gamma2,ds1 external ds1 COMMON/BDD7/gamma2,xdel COMMON/IDD7/l funki3=sin(2.d0*x)*sqrt(cos(x))*exp(-gamma2*sin(x)**2) & *cos(xdel*cos(x))*ds1(l,x)/2.d0 return end * function funki4(x) implicit none integer l real*8 funki4,x,xdel,gamma2,ds1 external ds1 COMMON/BDD7/gamma2,xdel COMMON/IDD7/l funki4=sin(2.d0*x)*sqrt(cos(x))*exp(-gamma2*sin(x)**2) & *sin(xdel*cos(x))*ds1(l,x)/2.d0 return end * DOUBLE PRECISION FUNCTION DS1(J,BETA) c---------------------------------------------------------------------------------- C C SMALL D. RETURNS REDUCED ROTATION MATRIX ELEMENTS D^J_{11} C C J ... the total angular momentum J c BETA ... the Euler angle of rotation (in radians) c---------------------------------------------------------------------------------- implicit none integer j real*8 plgndr,xt,beta external plgndr xt=cos(beta) ds1=(plgndr(j-1,0,xt)-plgndr(j+1,0,xt))/(1.d0+xt) ds1=ds1/dble(2*j+1)+plgndr(j,0,xt) return end * C (C) Copr. 1/1999 Alexander Moroz C============================================================================= subroutine ZKMATL(lcs,lmax,omega,zeps1,zeps0,rmf,kmt) C-------------------------------------------------------------------- C >>> lmax,omega,zeps1,zeps0,rmf C <<< kmt C Field KMT contains -tan(phase shift) C =================================== C GIVEN INTERIOR AND EXTERIOR PERMITTIVITIES EPS1 AND EPS0, RESPECTIVELY, C RADIUS RMUF AND FREQUENCY OMEGA, THIS SUBROUTINE CALCULATES THE ELEMENTS C OF THE REACTION MATRIX K FOR A DIELECTRIC SPHERE C C A VARIANT OF KMATL FOR METALLIC SPHERES C (KMT IS COMPLEX FIELD HERE) C C CAN BE ALSO USED FOR A COATED SPHERE C C * * * * * * C ============================== C EPS0 ... background permittivity C EPS1 ... in the case of a coated sphere, the core permittivity C KMT ... K MATRIX FOR A SPHERICAL SCATTERER C DKL ... (SIGN OF THE) DETERMINANT OF THE K MATRIX FOR A C SPHERICAL SCATTERER C C INTERNAL FUNCTIONS: C generic SQRT,INT,ABS,SIGN C >>> specific CMPLX(any)=complex ----> DCMPLX C REAL(integer)=real C C >>> CALLS BESSEL to generate the Bessel functions C-------------------------------------------------------------------- C >>> INPUT: <<< C C LMAX ... (LMAXV in main) maximal value of the orbital angular momentum C for which vector structure constants are calculated. C Partial wave expansion is used, which is badly convergent C for large size parameters $x> 100$. In numerical applications, the C series is to be cut off after C LMAX (LMX parameter here) \approx x+4x^{1/3}+2$. C In the case of LMX=50 this means that x <= 35. C If one wants to observe ripples, the cutoff for a given x has to C be even larger C C NLB ... (LMAXV+1)**2 in main - max. value of aligned (lm) index. C Specifies the dimension of vector structure constants VEC, kmt, C tin, cin, AMIN, B, BP C C EPS0 ... background permeability C EPS1 ... permeability of the exterior layer of the coated sphere C C RMUF ... muffin tin radius (in the natural units) C OMEGA ... frequency of photons. Its relation to EPS is according to C >>> OMEGA= PI*SQRT(EPS/EPS0) <<< C Conventional plot is for dimensionless SQRT(EPS(K)), C C SQRT(EPS)= (OMEGA*A)/(2*PI*C) C C NB ... labels differences of basis vectors of direct lattice C C >>> OUTPUT : <<< C C kmt ... - tan (phase shift) = AC*KMAT, C stored in a complex field C C TT1, TT2 ... transfer matrices for a coated sphere C C CM, DM, CE, DE ... nominator and denominator of the K-matrix C for electric and magnetic modes C C RFF ... radii of the coatings in the units of the radius of C the whole sphere C ============================== C RAISING LMAXV IN MAIN : C C Modify dims of JL, NL, PSI, PHI, UL, VL, etc C-------------------------------------------------------------------- implicit none integer LMM,LMAX C INPUT REQUIRED: C-------------------------------------------------------------------- C PARAMETER INPUT : C-------------------------------------------------------------------- * LMM is the local LMAX * PARAMETER (LMM=45) * C PARAMETER (NLB=(LMAX+1)**2) C*********************************************************************** C >>> DECLARATIONS: <<< C*********************************************************************** integer LCS,i,j,l cl integer lt1,lt2 REAL*8 RMF(lcs),omega * * RFF the ratio of the interior/exterior radius of the coated sphere * c PARAMETER (rff=0.1d0) * * LCS the number of layers of the coated sphere. * lcs=1 - homogeneous sphere * * declarations: COMPLEX*16 cqeps(2),zeps0,zeps1(lcs),zeps(lcs+1) COMPLEX*16 RX(2),SG(2) COMPLEX*16 KMT(2,lmax) cl COMPLEX*16 KMT(2,nlb) * * coated sphere declarations: * moving lmax ===> * cm, dm, ce, de(lcs,lmax) etc. * am, bm, ae, be(lmax) * tt1, tt2(2,2,lmax,2) * COMPLEX*16 cm(lcs,lmm),dm(lcs,lmm),ce(lcs,lmm),de(lcs,lmm) COMPLEX*16 tt1(2,2,lmm,2),tt2(2,2,lmm,2) COMPLEX*16 AM(lmm),AE(lmm),BM(lmm),BE(lmm) * C------------------------- C KMT,NML,DML,NEL,DEL,EPS,RX,SG, and Bessel functions below must be C declared complex in the case of absorptiom and/or amplification C-------------------------------------------------------------------- C ==================== C SUBR MODIF: C ===================== C REAL*4 JL(2,0:LMAX+1), NL(2,0:LMAX) C REAL*4 PSI(LMAX+1), PHI(LMAX+1) C REAL*4 DRJL(0:LMAX), DRNL(0:LMAX) C REAL*4 UL(2,0:LMAX), VL(2,0:LMAX) C REAL*4 DRUL(2,0:LMAX), DRVL(2,0:LMAX) c REAL*4 JL(2,0:*), NL(2,0:*) c REAL*4 PSI(LMAX+1), PHI(LMAX+1) c REAL*4 DRJL(0:*), DRNL(0:*) c REAL*4 UL(2,0:*), VL(2,0:*) c REAL*4 DRUL(2,0:*), DRVL(2,0:*) C COMPLEX*16 JL(0:lmm),NL(0:lmm) COMPLEX*16 DRJL(0:lmm),DRNL(0:lmm) COMPLEX*16 UL(2,0:lmm),VL(2,0:lmm) COMPLEX*16 DRUL(2,0:lmm),DRVL(2,0:lmm) C ===================== C-------------------------------------------------------------------- C******************************************************************** C READING THE DATA : C c DATA PI/3.141592653589793d0/ * * security trap - remainder * * if (lcs.eq.1) write(6,*)'Homogeneous sphere in ZKMATL' if (lcs.gt.1) write(6,*)'Coated sphere in ZKMATL' * if (lmax.gt.lmm) then write(6,*)' EXECUTION STOPPING IN ZKMATL' write(6,*)'Modify dims of JL, NL, UL, VL, etc to LMM=',LMAX stop end if * zeps(1)=zeps1(1) if (lcs.gt.1) zeps(lcs)=zeps1(lcs) zeps(lcs+1)=zeps0 * do l=1,lmax AM(l)=dcmplx(1.d0,0.d0) AE(l)=dcmplx(1.d0,0.d0) BM(l)=dcmplx(0.d0,0.d0) BE(l)=dcmplx(0.d0,0.d0) enddo * C******************************************************************** C >>> ASSIGNING THE VALUES FOR BESSEL FUNCTIONS and their C derivatives inside and outside the atom boundary. C THE DATA ARE PRODUCED BY THE SUBROUTINE BESSEL(Y,LMAX,PHI,PSI) C RECURSION RELATIONS ARE USED ACCORDING (AS 10.1.21-22) C NL(=YL) ARE CONSTRUCTED FROM JL BY USING (AS 10.1.15) C THE PREFIX DR MEANS THE DERIVATIVE WITH RESPECT TO RMUF AND C NOT RX(I)!! C * DO 28 j=1,lcs CQEPS(1)=SQRT(ZEPS(j)) SG(1)=omega*CQEPS(1) * CQEPS(2)=SQRT(ZEPS(j+1)) SG(2)=omega*CQEPS(2) * DO 25 I=1,2 * RX(I)=SG(I)*RMF(j) c WRITE(6,*)'i, rx(i)=', i, rx(i) C >>> * call gnzbess(RX(I),lmax,jl,drjl,nl,drnl) * c write(6,*)'jl=', jl * DO 15 L=1,LMAX C >>> (AS 10.1.22): UL(I,L)=RMF(j)*JL(L) VL(I,L)=RMF(j)*NL(L) DRJL(L)=SG(I)*DRJL(L) DRNL(L)=SG(I)*DRNL(L) DRUL(I,L)=JL(L)+RMF(j)*DRJL(L) DRVL(I,L)=NL(L)+RMF(j)*DRNL(L) 15 continue 25 CONTINUE * C >>> END OF THE LOOP TO ASSIGN VALUES OF BESSEL FUNCTIONS C JL and NL start to oscillate after RX.GT. approx 2.5 C******************************************************************** C Transfer matrix for a layered (coated) sphere C******************************************************************** * do l=1,lmax * * magnetic part * tt1(1,1,l,1)= UL(1,L) tt1(1,2,l,1)= VL(1,L) tt1(2,1,l,1)= DRUL(1,L) tt1(2,2,l,1)= DRVL(1,L) * tt2(1,1,l,1)= sg(2)*DRVL(2,L) tt2(1,2,l,1)= - sg(2)*VL(2,L) tt2(2,1,l,1)= - sg(2)*DRUL(2,L) tt2(2,2,l,1)= sg(2)*UL(2,L) * * electric part * tt1(1,1,l,2)=cqeps(1)*UL(1,L) tt1(1,2,l,2)=cqeps(1)*VL(1,L) tt1(2,1,l,2)=DRUL(1,L)/cqeps(1) tt1(2,2,l,2)= DRVL(1,L)/cqeps(1) * tt2(1,1,l,2)= sg(2)*DRVL(2,L)/cqeps(2) tt2(1,2,l,2)= -sg(2)*cqeps(2)*VL(2,L) tt2(2,1,l,2)= -sg(2)*DRUL(2,L)/cqeps(2) tt2(2,2,l,2)= sg(2)*cqeps(2)*UL(2,L) * * m-part * cm(j,l)=AM(l)*(tt2(1,1,l,1)*tt1(1,1,l,1) 1 +tt2(1,2,l,1)*tt1(2,1,l,1))+BM(l)*( 2 tt2(1,1,l,1)*tt1(1,2,l,1)+tt2(1,2,l,1)*tt1(2,2,l,1)) * dm(j,l)=AM(l)*(tt2(2,1,l,1)*tt1(1,1,l,1) 1 +tt2(2,2,l,1)*tt1(2,1,l,1))+BM(l)*( 2 tt2(2,1,l,1)*tt1(1,2,l,1)+tt2(2,2,l,1)*tt1(2,2,l,1)) * * e-part * ce(j,l)=AE(l)*(tt2(1,1,l,2)*tt1(1,1,l,2) 1 +tt2(1,2,l,2)*tt1(2,1,l,2))+BE(l)*( 2 tt2(1,1,l,2)*tt1(1,2,l,2)+tt2(1,2,l,2)*tt1(2,2,l,2)) * de(j,l)=AE(l)*(tt2(2,1,l,2)*tt1(1,1,l,2) 1 +tt2(2,2,l,2)*tt1(2,1,l,2))+BE(l)*( 2 tt2(2,1,l,2)*tt1(1,2,l,2)+tt2(2,2,l,2)*tt1(2,2,l,2)) * AM(l)=cm(j,l) BM(l)=dm(j,l) AE(l)=ce(j,l) BE(l)=de(j,l) * enddo * 28 CONTINUE C--------/---------/---------/---------/---------/---------/---------/-- C ASSIGNING VALUES TO ELEMENTS OF T AND dT^{-1}/dE MATRIX C TOGETHER WITH ALIGNEMENT OF INDICES ACCORDING TO C (L,M) ---> LIN=L**2 + L+1+M C-------------------------------------------------------------------- C >>> * DO 40 L=1,LMAX * c LT1=L**2+1 c LT2=(L+1)**2 c DO 40 LIN=LT1,LT2 * Field KMT contain -tan(phase shift)=sg*KMAT KMT(1,L)=bm(l)/am(l) KMT(2,L)=be(l)/ae(l) * 40 CONTINUE RETURN END C (C) Copr. 1/1999 Alexander Moroz C============================================================================= subroutine KPCSMAT(lmax,nlb,rmuf,omega,eps0,rff,kmt) C-------------------------------------------------------------------- C >>> lmax,nlb,rmuf,omega,eps1,eps0 C <<< kmt,dkl C =================================== C THIS SUBROUTINE CALCULATES THE ELEMENTS OF THE REACTION MATRIX K AT C FREQUENCY OMEGA FOR A PERFECTLY CONDUCTING SPHERE OF RADIUS RMUF C EMBEDDED IN A HOST WITH PERMITTIVITY EPS0. C C FOR PERFECTLY CONDUCTING SPHERE ELECTRIC FIELD INSIDE VANISHES C (cf. Jackson, p. 571) and the corresponding boundary conditions C for fileds outside the sphere are E_t=0 and B_n=0 C C A VARIANT OF KMATL FOR METALLIC SPHERES C (KMT IS COMPLEX FIELD HERE) C C CAN BE ALSO USED FOR A COATED SPHERE C C * * * * * * C ============================== C EPS0 ... background permittivity C EPS1 ... in the case of a coated sphere, the shell permittivity C KMT ... K MATRIX FOR A SPHERICAL SCATTERER C KMT=-tan(phase shift)/sg C DKL ... (SIGN OF THE) DETERMINANT OF THE K MATRIX FOR A C SPHERICAL SCATTERER C C INTERNAL FUNCTIONS: C generic SQRT,INT,ABS,SIGN C >>> specific CMPLX(any)=complex ----> DCMPLX C REAL(integer)=real C C >>> CALLS BESSEL to generate the Bessel functions C-------------------------------------------------------------------- C >>> INPUT: <<< C C LMAX ... (LMAXV in main) maximal value of the orbital angular momentum C for which vector structure constants are calculated C NLB ... (LMAXV+1)**2 in main - max. value of aligned (lm) index. C Specifies the dimension of vector structure constants VEC, kmt, C tin, cin, AMIN, B, BP C C EPS0 ... background permeability C EPS1 ... permeability of the exterior layer of the coated sphere C C RMUF ... muffin tin radius (in the natural units) C OMEGA ... frequency of photons. Its relation to EPS is according to C >>> OMEGA= PI*SQRT(EPS/EPS0) <<< C Conventional plot is for dimensionless SQRT(EPS(K)), C C SQRT(EPS)= (OMEGA*A)/(2*PI*C) C C NB ... labels differences of basis vectors of direct lattice C C >>> OUTPUT : <<< C C kmt ... t- matrix corresponding to $j_l-i*\sg t_l h_l^+$, C stored in a complex field C C TT1, TT2 ... transfer matrices for a coated sphere C C CM, DM, CE, DE ... nominator and denominator of the K-matrix C for electric and magnetic modes C C ============================== C RAISING LMAXV IN MAIN : C C Modify dims of JL, NL, PSI, PHI, UL, VL, etc C-------------------------------------------------------------------- C INPUT REQUIRED: C-------------------------------------------------------------------- C PARAMETER INPUT : C-------------------------------------------------------------------- * LMM is the local LMAX * PARAMETER (LMM=11) * C PARAMETER (NLB=(LMAX+1)**2) C*********************************************************************** C >>> DECLARATIONS: <<< C*********************************************************************** IMPLICIT REAL*8 (A-H,O-Z) * * coated sphere parameters: * radius of the interior sphere of the coated sphere * c PARAMETER (rmuf1=0.2d0) * * ratio of the interior/exterior radius of the coated sphere * c PARAMETER (rff=0.8d0) * * number of layers of the coated sphere. If lcs=1 - homogeneous sphere * PARAMETER (lcs=1) * * declarations: * c real*8 wt,xt,w,x,z c integer L,M,AL,I,J,K,N,LIN COMPLEX*16 ey,cqeps(2) REAL*8 EPS(lcs+1) REAL*8 RMF(lcs),RMUF,NEL,NML REAL*8 DEL,DML COMPLEX*16 RX(2),SG(2) COMPLEX*16 KMT(2,NLB),KMT1,KMT2 * * coated sphere declarations: * moving lmax ===> * cm, dm, ce, de(lcs,lmax) etc. * am, bm, ae, be(lmax) * tt1, tt2(2,2,lmax,2) * COMPLEX*16 cm(lcs,lmm),dm(lcs,lmm),ce(lcs,lmm),de(lcs,lmm) COMPLEX*16 tt1(2,2,lmm,2), tt2(2,2,lmm,2) COMPLEX*16 AM(lmm), AE(lmm), BM(lmm), BE(lmm) * C------------------------- C KMT,NML,DML,NEL,DEL,EPS,RX,SG, and Bessel functions below must be C declared complex in the case of absorptiom and/or amplification C-------------------------------------------------------------------- C ==================== C SUBR MODIF: C ===================== C REAL*4 JL(2,0:LMAX+1), NL(2,0:LMAX) C REAL*4 PSI(LMAX+1), PHI(LMAX+1) C REAL*4 DRJL(0:LMAX), DRNL(0:LMAX) C REAL*4 UL(2,0:LMAX), VL(2,0:LMAX) C REAL*4 DRUL(2,0:LMAX), DRVL(2,0:LMAX) c REAL*4 JL(2,0:*), NL(2,0:*) c REAL*4 PSI(LMAX+1), PHI(LMAX+1) c REAL*4 DRJL(0:*), DRNL(0:*) c REAL*4 UL(2,0:*), VL(2,0:*) c REAL*4 DRUL(2,0:*), DRVL(2,0:*) C COMPLEX*16 JL(0:lmm), NL(0:lmm) COMPLEX*16 DRJL(0:lmm), DRNL(0:lmm) COMPLEX*16 UL(2,0:lmm), VL(2,0:lmm) COMPLEX*16 DRUL(2,0:lmm), DRVL(2,0:lmm) C ===================== C-------------------------------------------------------------------- C******************************************************************** C READING THE DATA : C c DATA PI/3.141592653589793d0/ DATA ey/(0.d0,1.d0)/ * * security trap - remainder * * if (lcs.eq.1) write(6,*)'Homogeneous perfectly conducting 1 sphere in KPCSMAT' if (lcs.gt.1) write(6,*)'Coated sphere in KPCSMAT' if (lcs.ne.1 .and. rff.eq.2) then write(6,*)' EXECUTION STOPPING IN KPCSMAT' write(6,*)'Verify if LCS should be .gt.1' stop end if if (lmax.gt.lmm) then write(6,*)' EXECUTION STOPPING IN KPCSMAT' write(6,*)'Modify dims of JL, NL, UL, VL, etc to LMM=',LMAX stop end if * * radii of the coated sphere: * rmf(1)=rff*rmuf rmf(lcs)=rmuf om=omega * * n=2.35 <---> eps = 5.5225 * * n(silica)=1.45 <---> EPS(1)=2.1025D0 * n(ZnS)=2. <---> EPS(1)=4.D0 * EPS(1) ... core diel. constant * EPS(2) ... coating shell diel. constant * EPS(0) ... background diel. constant c EPS(2)=2.1025D0 EPS(lcs+1)=EPS0 * c do l=1,lmax c AM(l)=dcmplx(1.d0,0.d0) c AE(l)=dcmplx(1.d0,0.d0) c BM(l)=dcmplx(0.d0,0.d0) c BE(l)=dcmplx(0.d0,0.d0) c enddo * C******************************************************************** C >>> ASSIGNING THE VALUES FOR BESSEL FUNCTIONS and their C derivatives inside and outside the atom boundary. C THE DATA ARE PRODUCED BY THE SUBROUTINE BESSEL(Y,LMAX,PHI,PSI) C RECURSION RELATIONS ARE USED ACCORDING (AS 10.1.21-22) C NL(=YL) ARE CONSTRUCTED FROM JL BY USING (AS 10.1.15) C THE PREFIX DR MEANS THE DERIVATIVE WITH RESPECT TO RMUF AND C NOT RX(I)!! C * if (eps(2).lt.0) CQEPS(2)=ey*SQRT(DABS(EPS(2))) if (eps(2).ge.0) CQEPS(2)=SQRT(EPS(2)) SG(2)=OM*CQEPS(2) RX(1)=SG(2)*RMF(1) * CALL GNCBESS(RX(1),LMAX,jl,drjl,nl,drnl) * DO 10 L=1,LMAX C >>> (AS 10.1.22): UL(1,L)=RMF(1)*JL(L) VL(1,L)=RMF(1)*NL(L) DRJL(L)=SG(2)*DRJL(L) DRNL(L)=SG(2)*DRNL(L) DRUL(1,L)=JL(L)+RMF(1)*DRJL(L) DRVL(1,L)=NL(L)+RMF(1)*DRNL(L) AM(l)= NL(L) ! cm(1,l) BM(l)=-JL(L) ! dm(1,l) AE(l)= DRVL(1,L) ! ce(1,l) BE(l)=-DRUL(1,L) ! de(1,l) * cf. Jackson 1962, p. 571, Eqs. (16.147); * B/A should yield -tan(phase shift) 10 continue if (lcs.eq.1) go to 30 DO 28 j=2,lcs SG(1)=SG(2) if (eps(j+1).lt.0) CQEPS(2)=ey*SQRT(DABS(EPS(j+1))) if (eps(j+1).ge.0) CQEPS(2)=SQRT(EPS(j+1)) SG(2)=OM*CQEPS(2) * DO 25 I=1,2 * RX(I)=SG(I)*RMF(j) C >>> * CALL GNCBESS(RX(I),LMAX,jl,drjl,nl,drnl) * DO 15 L=1,LMAX C >>> (AS 10.1.22): UL(I,L)=RMF(j)*JL(L) VL(I,L)=RMF(j)*NL(L) DRJL(L)=SG(I)*DRJL(L) DRNL(L)=SG(I)*DRNL(L) DRUL(I,L)=JL(L)+RMF(j)*DRJL(L) DRVL(I,L)=NL(L)+RMF(j)*DRNL(L) 15 continue 25 CONTINUE * C >>> END OF THE LOOP TO ASSIGN VALUES OF BESSEL FUNCTIONS C JL and NL start to oscillate after RX.GT. approx 2.5 C******************************************************************** C Transfer matrix for a layered (coated) sphere C******************************************************************** * do l=1,lmax * * magnetic part * tt1(1,1,l,1)= UL(1,L) tt1(1,2,l,1)= VL(1,L) tt1(2,1,l,1)= DRUL(1,L) tt1(2,2,l,1)= DRVL(1,L) * tt2(1,1,l,1)= sg(2)*DRVL(2,L) tt2(1,2,l,1)= - sg(2)*VL(2,L) tt2(2,1,l,1)= - sg(2)*DRUL(2,L) tt2(2,2,l,1)= sg(2)*UL(2,L) * * electric part * tt1(1,1,l,2)=cqeps(1)*UL(1,L) tt1(1,2,l,2)=cqeps(1)*VL(1,L) tt1(2,1,l,2)=DRUL(1,L)/cqeps(1) tt1(2,2,l,2)= DRVL(1,L)/cqeps(1) * tt2(1,1,l,2)= sg(2)*DRVL(2,L)/cqeps(2) tt2(1,2,l,2)= -sg(2)*cqeps(2)*VL(2,L) tt2(2,1,l,2)= -sg(2)*DRUL(2,L)/cqeps(2) tt2(2,2,l,2)= sg(2)*cqeps(2)*UL(2,L) * * m-part * cm(j,l)=AM(l)*(tt2(1,1,l,1)*tt1(1,1,l,1) 1 +tt2(1,2,l,1)*tt1(2,1,l,1))+BM(l)*( 2 tt2(1,1,l,1)*tt1(1,2,l,1)+tt2(1,2,l,1)*tt1(2,2,l,1)) * dm(j,l)=AM(l)*(tt2(2,1,l,1)*tt1(1,1,l,1) 1 +tt2(2,2,l,1)*tt1(2,1,l,1))+BM(l)*( 2 tt2(2,1,l,1)*tt1(1,2,l,1)+tt2(2,2,l,1)*tt1(2,2,l,1)) * * e-part * ce(j,l)=AE(l)*(tt2(1,1,l,2)*tt1(1,1,l,2) 1 +tt2(1,2,l,2)*tt1(2,1,l,2))+BE(l)*( 2 tt2(1,1,l,2)*tt1(1,2,l,2)+tt2(1,2,l,2)*tt1(2,2,l,2)) * de(j,l)=AE(l)*(tt2(2,1,l,2)*tt1(1,1,l,2) 1 +tt2(2,2,l,2)*tt1(2,1,l,2))+BE(l)*( 2 tt2(2,1,l,2)*tt1(1,2,l,2)+tt2(2,2,l,2)*tt1(2,2,l,2)) * AM(l)=cm(j,l) BM(l)=dm(j,l) AE(l)=ce(j,l) BE(l)=de(j,l) * enddo * 28 CONTINUE C--------/---------/---------/---------/---------/---------/---------/-- C ASSIGNING VALUES TO ELEMENTS OF K MATRIX C TOGETHER WITH ALIGNEMENT OF INDICES ACCORDING TO C (L,M) ---> LIN=L**2 + L+1+M C-------------------------------------------------------------------- C >>> 30 CONTINUE * DO 40 L=1,LMAX * LT1=L**2+1 LT2=(L+1)**2 DO 40 LIN=LT1,LT2 * KMT(1,LIN)=bm(l)/(am(l)*sg(2)) KMT(2,LIN)=be(l)/(ae(l)*sg(2)) * * B/A yields -tan(phase shift) ==> KMT=-tan(phase shift)/sg * 40 CONTINUE RETURN END