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 2005 C C This is CHEWFS, a computer code to calculate the normalized spontaneous C decay rates and frequency shifts for an atom interacting C with a complex spherical nanoshell. CHEWFS calculates the total decay C rates W^t and frequency shifts for two principal (tangential and radial) C and averaged dipole orientations. The quantities are normalized with C respect to the radiative decay rates of the same dipole in a free-space C filled in with a homogeneous medium which is identical to that at C the dipole position. C C CHEWFS code is based on the article C C [1] A. Moroz, C A recursive transfer-matrix solution for a dipole radiating inside and C outside a stratified sphere, C Annals of Physics (NY) 315, 352-418 (2005) C (http://dx.doi.org/10.1016/j.aop.2004.07.002) C C and its follow up C C [2] A. Moroz, C Spectroscopic properties of a two-level atom interacting with C a complex spherical nanoshell, C Chem. Phys. 317(1), 1-15 (2005) (published online on 9 August 2005) C (http://dx.doi.org/10.1016/j.chemphys.2005.05.003) C [preprint quant-ph/0412094 is available from C http://xxx.lanl.gov/abs/quant-ph/0412094] C C The program default options correspond to the cases considered in the latter C article. Default option is to have four concentric shell. C The second and fourth shells have the dielectric constant of bulk gold C (material data supplied by Audat.dat). The executable will ask you to supply: C C 1) the respective four spherical shell radii in nm (separated by empty space); C C 2) the real dielectric constants of the nanoshell core and 3rd shell, and C of ambient (separated by empty space); C C 3) the vacuum (detection) wavelength lambda. C C In order to modify the default option, read through the comments C in the program and the description of parameters and variables. C C Some additional plots obtained with the code have been published in a follow up: C C [3] A. Moroz, C Non-radiative decay of a dipole emitter close to a metallic C nanoparticle: Importance of higher-order multipole contributions, C Opt. Commun. 283 (10), 2277-2287 (2010). C http://dx.doi.org/10.1016/j.optcom.2010.01.061 C C If you want to run the code for a homogeneous sphere, follow C the instructions available at C http://www.wave-scattering.com/chew-man.pdf C C ================== C C IMPORTANT!!! C C Before comparing the simulations with your experiment, read Sec. 12.1 of my [1] C ======================================================================== C OUTPUT: C C CHEWFS returns files C C dpfshiftm.dat: Frequency shift in the presence of a single C & coated sphere as a function of the distance normalized to the C & decay rates in a free space filled in with the local medium' C C dipqmratesm.dat: Total decay rate W^t in the presence of a single C & coated sphere as a function of the distance normalized to the C & decay rates in a free space filled in with the local medium' C C Each output file has data organized in four columns: C normalized radial position (in the units of the bead radius) C and the corresponding quantities for the C respective tangential, radial, and averaged dipole orientations C (from left to right). Note in passing that the fourth column is C not a simple average of the second and third columns but is calculated C according to Eq. (123) of Ref. [1]. C ======================================================================== C C TESTS: C C 1) For nanoshell A and wavelength of 595 nm: C C the first lines of dpfshiftm.dat should read as follows: C C .010000 -.19402060E+01 -.19408732E+01 -.19404284E+01 C .019950 -.19417693E+01 -.19444298E+01 -.19426561E+01 C .029900 -.19443989E+01 -.19503934E+01 -.19463971E+01 C .039851 -.19481362E+01 -.19588303E+01 -.19517009E+01 C .049801 -.19530402E+01 -.19698350E+01 -.19586385E+01 C .059751 -.19591885E+01 -.19835316E+01 -.19673028E+01 C .069701 -.19666781E+01 -.20000763E+01 -.19778109E+01 C .079652 -.19756280E+01 -.20196602E+01 -.19903054E+01 C .089602 -.19861804E+01 -.20425129E+01 -.20049579E+01 C .099552 -.19985038E+01 -.20689074E+01 -.20219717E+01 C .109502 -.20127966E+01 -.20991655E+01 -.20415862E+01 C .119453 -.20292905E+01 -.21336649E+01 -.20640820E+01 C C and the first lines of dipqmratesm.dat should read as follows: C C .010000 .87510855E+00 .87533539E+00 .87518416E+00 C .019950 .87581103E+00 .87671599E+00 .87611268E+00 C .029900 .87699284E+00 .87903345E+00 .87767305E+00 C .039851 .87867296E+00 .88231726E+00 .87988773E+00 C .049801 .88087821E+00 .88660921E+00 .88278854E+00 C .059751 .88364361E+00 .89196393E+00 .88641705E+00 C .069701 .88701281E+00 .89844965E+00 .89082509E+00 C .079652 .89103858E+00 .90614907E+00 .89607541E+00 C .089602 .89578348E+00 .91516059E+00 .90224252E+00 C .099552 .90132071E+00 .92559968E+00 .90941370E+00 C .109502 .90773506E+00 .93760060E+00 .91769024E+00 C .119453 .91512410E+00 .95131846E+00 .92718889E+00 C C 2) IN THE CASE OF A HOMOGENEOUS NONABSORBING SPHERE, C C - THE FREQUENCY SHIFT RESULTS SHOULD COINCIDE WITH C THOSE OBTAINED BY C V. V. Klimov, M. Ducloy, and V. S. Letokhov, C ``Spontaneous emission rate and level shift of an atom inside C a dielectric microsphere", C J. Mod. Opt. 43, 549-563 (1996) (see, e.g., Fig. 3 therein) C C - THE DECAY RATE RESULTS SHOULD COINCIDE WITH THOSE OBTAINED BY C C H. Chew, Transition rates of atoms near spherical surfaces, C J. Chem. Phys. 87, 1355-1360 (1987) C (see, e.g., Figs. 2,3 for dipole outside a C dielectric particle in the ambient normalization C for the size parameter x=ka=5 and the refractive C index n=1.3 and 1.5, respectively) C H. Chew, Radiation and lifetimes of atoms inside dielectric particles, C Phys. Rev. A 38, 3410-3416 (1988) C (see Fig.1 for dipole inside a C dielectric particle in the ambient normalization) C C 3) You can also download related executable nanoshell.exe from C C http://www.wave-scattering.com/nanoshell.zip C C and test your results against the executable output. C======================================================================== C CONVERGENCE: C C You receive no warning unless convergence is worse than up to C 3 exact digits. This should only happen in the case of frequency C shifts in the immediate proximity of metal shell boundaries. The C respective values of tol1 and tol2 parameters will tell you with C which precision the result was obtained. C======================================================================== C C REQUIREMENTS & NOT DISTRIBUTED DEPENDENT ROUTINES: C C 1) You should download data file "Audat.dat" from C http://www.wave-scattering.com/Audat.dat.txt C C 2) You should supply your own routine C C GNZBESS(RX(I),LMAX,jl,drjl,nl,drnl) C C which generates an array of the Bessel functions C JL(0:LMAX),NL(0:LMAX) of the first and second order C and their derivatives DJL(0:LMAX),DNL(0:LMAX) for a C complex argument RX(I) C C and C C CALL GNRICBESSH(CQEPS(I),OMEGA*rmf(j),LMX,hl,drhl) C C which returns an array of Riccati-Hankel functions C zeta=x*hl(x), and C dzeta=d[x*hl(x)]/dx of the complex argument C RX(i)=CQEPS(i)*OMEGA*rmff(j) C C My program was tested together with the AMOS package of C complex Bessel functions routines which C you can, together with their dependencies, C freely download at Netlib (http://www.netlib.org/amos). C Before using the AMOS package, machine constants C in the real function D1MACH(I) and in the integer function I1MACH(I) C have to be adjusted to your computer. When using Microsoft Windows, C you can activate the same machine constants as listed therein under C the option of Silicon Graphics IRIS. C C My Bessel function package bessel.zip is available at C http://www.wave-scattering.com/bessel.zip 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. 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: 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 chewfs C--------/---------/---------/---------/---------/---------/---------/-- C C GIVEN INTERIOR AND EXTERIOR PERMITTIVITIES EPS1 AND EPS0, C RESPECTIVELY, RADIUS RMUF AND FREQUENCY OMEGA, C THIS SUBROUTINE CALCULATES THE FREQUENCY SHIFTS AND TOTAL C DECAY RATES NORMALIZED TO THE DEACY RATES IN THE VACUUM FILLED C WITH THE LOCAL DIPOLE MEDIUM. C * * * * * * C ============================== C ZEPS0 ... background permittivity C CSEPS ... in the case of coated sphere, the outer shell permittivity C C INTERNAL FUNCTIONS: C generic SQRT,INT,ABS,SIGN C >>> specific CMPLX(any)=complex ----> DCMPLX C REAL(integer)=real C C >>> ROUTINES CALLED: C C GNZBESS to generate arrays of spheical Bessel functions c QROM0 ... NUMERICAL RECIPES ROUTINE WHICH PERFORMS A 1D INTEGRAL C If you do not have this routine, substitute it with C an alternative one C-------------------------------------------------------------------- C >>> INPUT: <<< C C LMAX ... the maximal value of the orbital angular momentum C for which vector structure constants are calculated C EPS0 ... background permeability C EPS1 ... permeability of the exterior layer of the coated sphere C C RMUF ... sphere radius in the units of a length scale (RMUF=1 by default) C OMEGA ... frequency of photons in dimensionless units 2*PI*RSNM/(LAMBDA*RMUF), C where RSNM is the sphere radius and LAMBDA the wavelength C of light in the vacuum C--------/---------/---------/---------/---------/---------/---------/-- implicit none integer NOUT,NX,LMAX,NMAT,LCS,NFIN,NRS,NLD,ilcs,iabs real*8 TOL,OMEGA,RINT,pi,rmuf complex*16 ZEPS0,CCEPS,CSEPS,cp character*1 ync logical ynperfcon,ynbrug C--------/---------/---------/---------/---------/---------/---------/-- C INPUT REQUIRED: C-------------------------------------------------------------------- C PARAMETER INPUT : C-------------------------------------------------------------------- C ::: number of the output unit PARAMETER (NOUT=60) C >>> ANGULAR MOMENTUM CUTOFF ON ARRAY DIMENSION C (The actual angular-momentum cutoff on summation is specified C by the value of variable LMX below) PARAMETER (LMAX=60) * Specify the upper bound on the number of additional wavelengths * used in simulations. If NLD=0 then simulation is performed * for a single wavelength. If you chose NLD different from zero, * go below to "DO 400 ILD=0,NLD" and specify elementary increment * in the wavelengths. PARAMETER (NLD=0) * Specify the upper bound on the number of additional sphere radii * used in simulations. If NRS=0 then simulation is performed * for a single sphere radius. If you chose NRS different from zero, * go below to "300 IRS=0,NRS" and specify elementary increment * in the sphere radius. PARAMETER (NRS=0) C ::: relative error allowed for the LDOS. If the convergence * within TOL is not reached, program issues warning PARAMETER (TOL=1.d-3) c: frequency c PARAMETER (omega=0.2D0) C ::: number of steps on direct space (NX.GE.1): PARAMETER (NX=201) C ::: interval on which the LDOS is calculated * RMUF=1 is the sphere radius PARAMETER (RINT=2.d0) * * If sphere is coated, ync=y, otherwise ync=n parameter (ync='y') * * ynperfcon=.true. if core is a perfect conductor, otherwise * ynperfcon=.false. * PARAMETER (ynperfcon=.false.) !The .true. option has not yet been !incorporated in the code * * coated sphere parameters: * * ratio of the interior/exterior radius of the coated sphere * c PARAMETER (rff=0.2d0) * * number of layers of the coated sphere. If lcs=1 - homogeneous sphere * PARAMETER (lcs=4) * c The coating layer to which material data are read in parameter (ilcs=2) c The absorbing layer. If there is no absorbing layer, or you do not want C to calculate the nonradiative decay rates, then set IABS to a negative number parameter (iabs=2) C >>> BACKGROUND PERMITTIVITY <<< c PARAMETER (ZEPS0=1.7689D0) !2.281D0; 1.7689D0=1.33**2=H2O c sphere core dielectric constant PARAMETER (CCEPS=(2.1025d0,0.d0)) !glycerol=2.1609d0 c PARAMETER (CCEPS=(6.d0,0.d0)) !KDL C >>> SPHERE (OUTER SHELL SCATTERER) PERMITTIVITY <<< * n(silica)=1.45 <---> EPS(1)=2.1025D0 * n(ZnS)=2. <---> EPS(1)=4.D0 PARAMETER (CSEPS=2.1025D0) !2.1025=1.45**2=SiO2 c Temporarily option for reading of the real data for the dielectric constant c The number of the entries in the data file AGC.DAT to be read below cx PARAMETER (NFIN=73) * material code number c NMAT=0 dispersionless dielectric c NMAT=1 Drude metal c NMAT=2 Ag c NMAT=3 Au c NMAT=4 ZnS c NMAT=5 Cu c NMAT=6 Al c NMAT=7 Pt c NMAT=8 Si * PARAMETER(NMAT=3) * 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. * Data files for Au,Cu,Al,Pt should be ordered with the decreased wavelength * (omega increases in the loop and is oriented along the data file) * c AGC.DAT NFIN=73 ! from Palik c Audat.dat NFIN=66 ! from Palik c Au_2dat.dat NFIN=76 ! from JAW c Au*new.dat NFIN=142 c Cudat.dat NFIN=47 ! from Palik c Aldat.dat NFIN=80 ! from Palik c Ptdat.dat NFIN=53 ! from Palik c Nidat.dat NFIN=68 ! from Palik c Sidat.dat NFIN=291 c sieps.dat NFIN=117 * PARAMETER (NFIN=66) c If ynbrug=.true., performs Bruggeman approximation for ZEPS1. Otherwise c ynbrug=false. parameter (ynbrug=.false.) C*********************************************************************** C >>> DECLARATIONS: <<< C*********************************************************************** integer lmx,ieps,ikl,irs,ild integer i,j,l,il,ml real*8 xd1,xd2,tol1,tol2,gim REAL*8 rsnm,rsnm0,lambda,lambda0 REAL*8 rmff,delo,gperp,gpar,gqmpar,gqmperp real*8 omxf,reepsz,plasma,omxp,ff,filfrac !ff for Bruggemann;filfrac for ZnS complex*16 ci,cone,czero,zdet1,zdet2,zeps1,Z1,Z2,zz,zgperp,zgpar * Array declarations: REAL*8 RMF(lcs),rff(lcs) REAL*8 dos(nx+1),dosperp(nx+1),dospar(nx+1) COMPLEX*16 RX(2),SG(2),CQEPS(2),zeps(lcs+1) real*8 omf(NFIN) complex*16 ceps1(NFIN) c COMPLEX*16 KMT(2,lmx),alpha(lmx),beta(lmx) * * coated sphere declarations: * moving lmax ===> * COMPLEX*16 cxm(1:lmax),dxm(1:lmax) COMPLEX*16 cxe(1:lmax),dxe(1:lmax) COMPLEX*16 cfel(1:lmax),cfml(1:lmax),cdrfel(1:lmax), & cfel1(1:lmax),cfml1(1:lmax),cdrfel1(1:lmax), & cfel2(1:lmax),cfml2(1:lmax),cdrfel2(1:lmax) COMPLEX*16 cmx11(1:lmax),cmx12(1:lmax),cmx21(1:lmax),cmx22(1:lmax) COMPLEX*16 tt1(2,2,1:lmax,2,lcs),tt2(2,2,1:lmax,2,lcs) * C--------/---------/---------/---------/---------/---------/---------/-- C Bessel functions declaration COMPLEX*16 JL(0:lmax), NL(0:lmax) COMPLEX*16 DRJL(0:lmax), DRNL(0:lmax) COMPLEX*16 UL(2,0:lmax), DRUL(2,0:lmax) COMPLEX*16 WL(2,0:lmax), DRWL(2,0:lmax) COMPLEX*16 HL(0:lmax),DRHL(0:lmax) C******************************************************************** C READING THE DATA : C DATA PI/3.141592653589793d0/ DATA ci/(0.d0,1.d0)/,cone/(1.d0,0.d0)/,czero/(0.d0,0.d0)/ DATA RMUF/1.0D0/ DATA LMX/60/ * * security trap - remainder * write(6,*)'YOU ARE CALCULATING NORMALIZED DECAY RATES AND' write(6,*)'FREQUENCY SHIFTS FOR AN ATOM DIPOLE INTERACTING WITH A' write(6,*)'COMPLEX NANOSHELL' write(6,*) * if ((ync.eq.'n').and.(ilcs.gt.1)) then write(6,*)'YNC=', ync,' and ILCS=',ilcs,' params incompatible' stop end if * if (nmat.gt.1) then write(6,*)'(Real material data are to be provided.)' * if (ynbrug) write(6,*)'Bruggeman approx. used!' end if * if (ync.eq.'y' .and. lcs.eq.1) then write(6,*)'Check compatibility of YNC and LCS' stop end if * * Reading in the input data: * radii of the coated sphere: * 3 write(6,*)'Read the respective four spherical shell radii in nm' write(6,*)'(separated by empty space)' read(5,*) rff(1),rff(2),rff(3),rsnm * check entry: if ((rff(1).eq.0).or.((rff(1).gt.rff(2)).or.(rff(2).gt.rff(3)) & .or.(rff(3).gt.rsnm))) & then write(6,*)'Shell radii supplied incorrectly!' write(6,*)'Please, supply them again.' pause go to 3 end if * rff(1)=rff(1)/rsnm rmf(1)=rff(1)*rmuf * rmf(lcs)=rmuf * * n=2.35 <---> eps = 5.5225 * * n(silica)=1.45 <---> EPS(1)=2.1025D0 * n(ZnS)=2. <---> EPS(1)=4.D0 * 2.1025 * * dielectric constants of the coated sphere: * * 4 write(6,*)'Read the real dielectric constants of the nanoshell' write(6,*)'core, 3rd shell, and ambient' write(6,*)'(separated by empty space):' read(5,*) xd1,xd2,gim ZEPS(1)=xd1 ZEPS(3)=xd2 ZEPS(lcs+1)=gim zeps0=gim * check entry: if ((xd1.le.0).or.(xd2.le.0).or.(gim.le.0)) & then write(6,*)'One of the dielectric constants is negative.' write(6,*)'Please, supply them again!' pause go to 4 end if * * radii of the coated sphere: * if (lcs.ge.3) then !multicoated particle * do ikl=2,lcs-1 rff(ikl)=rff(ikl)/rsnm rmf(ikl)=rff(ikl)*rmuf * enddo * end if write(6,*)'Give the vacuum (detection) wavelength =' read (5,*) lambda ct lambda=1800.d0 C--------/---------/---------/---------/---------/---------/---------/-- * The header of diprates*.dat output files: OPEN(UNIT=NOUT,FILE='dpfshiftm.dat') rewind(NOUT) WRITE(NOUT,*)'#Frequency shift in the presence of a single & coated sphere as a function of the distance normalized to the & decay rates in a free space filled in with the local medium' * write(NOUT,*)'#Angular momentum cut-off LMAX=',LMX write(nout,*)'#Material number =',NMAT if (ynbrug) write(nout,*)'#Bruggeman approximation performed' if (ync .eq.'n') then write(nout,*)'# Homogeneous sphere' WRITE(NOUT,*) '#Set ync .eq. y if you want coated sphere!' else if (ync.eq.'y') then write(nout,*)'# Coated sphere' WRITE(NOUT,*) '#Set ync .eq. n if you want homogeneous sphere!' end if * c WRITE(NOUT,*) write(nout,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT,*) * WRITE(NOUT,*)'#the vacuum (detection) wavelength =', lambda write(nout,*) '# rsnm/wavelength=', rsnm/lambda * write(nout,*)'#Numbers of layers lcs=', lcs write(nout,*)'#host dielectric constant=', dble(zeps(lcs+1)) write(nout,*)'#Diel. constants of diel. shells 1 and 3:', & dble(zeps(1)),dble(zeps(3)) write(nout,*)'#Nanoshell radii:' write(nout,*)rff(1)*rsnm,'/',rff(2)*rsnm,'/',rff(3)*rsnm,'/',rsnm C--------/---------/---------/---------/---------/---------/---------/-- write(nout,*)'#Length of the searched interval rint=', rint write(nout,*)'#Number of the points on the uniform mesh nx=',nx c WRITE(NOUT,*)'#The LDOS in the free case (in background zeps0)', c 1 freeld WRITE(NOUT,*)'#In columns:' WRITE(NOUT,*)'#position, the dipole transition rates' write(nout,*)'#for tangentially and radially oriented dipole,' write(nout,*)'#and their local average, respectively' WRITE(NOUT,*) cxxxxxxxxxxxx C--------/---------/---------/---------/---------/---------/---------/-- * The header of dipqmrates*.dat output files: OPEN(UNIT=nout+7,FILE='dipqmratesm.dat') rewind(nout+7) WRITE(nout+7,*)'#Total decay rates in the presence of a single & coated sphere as a function of the distance normalized to the & decay rates in a free space filled in with the local medium' * write(nout+7,*)'#Angular momentum cut-off LMAX=',LMX write(nout+7,*)'#Material number =',NMAT if (ynbrug) write(nout+7,*)'#Bruggeman approximation performed' if (ync .eq.'n') then write(nout+7,*)'# Homogeneous sphere' WRITE(nout+7,*) '#Set ync .eq. y if you want coated sphere!' else if (ync.eq.'y') then write(nout+7,*)'# Coated sphere' WRITE(nout+7,*) '#Set ync .eq. n if you want homogeneous sphere!' end if * c WRITE(nout+7,*) write(nout+7,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(nout+7,*) * WRITE(nout+7,*)'#the vacuum (detection) wavelength =', lambda write(nout+7,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+7,*)'#Numbers of layers lcs=', lcs write(nout+7,*)'#host dielectric constant=',dble(zeps(lcs+1)) write(nout+7,*)'#Diel. constants of diel. shells 1 and 3:', & dble(zeps(1)),dble(zeps(3)) write(nout+7,*)'#Nanoshell radii:' write(nout+7,*) & rff(1)*rsnm,'/',rff(2)*rsnm,'/',rff(3)*rsnm,'/',rsnm C--------/---------/---------/---------/---------/---------/---------/-- write(nout+7,*)'#Length of the searched interval rint=', rint write(nout+7,*)'#Number of the points on the uniform mesh nx=',nx c WRITE(nout+7,*)'#Decay rates in the free case (in background zeps0)', c 1 freeld WRITE(nout+7,*)'#In columns:' WRITE(nout+7,*)'#position, the dipole decay rates' write(nout+7,*)'#for tangentially and radially oriented dipole,' write(nout+7,*)'#and their local average, respectively' WRITE(NOUT+7,*) cxxxxxxxxxxxx C--------/---------/---------/---------/---------/---------/---------/-- rsnm0=rsnm !initial radius lambda0=lambda !initial wavelength DO 400 ILD=0,NLD lambda=lambda0+dble(ild)*5 DO 300 IRS=0,NRS rsnm=rsnm0+dble(irs) * if (irs.ne.0) then write(nout,*)'#Sphere radius in nm=', rsnm write(nout+7,*)'#Sphere radius in nm=', rsnm WRITE(NOUT,*) WRITE(NOUT+7,*) end if * omega=2.d0*pi*rsnm/lambda !=kr_s or size parameter * C--------/---------/---------/---------/---------/---------/---------/-- * By activating line below you can feed in a desired sphere size parameter * to check program results against earlier calculations (for instance those * of Chew), or you can use an output of other programs to feed in here * the size parameter corresponding to a Mie resonance cc omega=2.51022320078408d0 !temporarily check of Chew's results cc omega=5.d0 !temporarily check of Chew's results outside cc omega=3.3d0 !temporarily check of KDL's results C--------/---------/---------/---------/---------/---------/---------/-- * If nmat.ge.1, one of the ZEPS entries is overwritten below: * Sphere optical constants in the case of a dispersion: if (nmat.le.1) goto 1 ! goto frequency loop C else > READING IN MATERIAL DATA: * Reading real material data, e.g., 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 * material constant reading: * if (nmat.eq.2) then ! silver data OPEN(UNIT=30,FILE='agc.dat') rewind(30) do ieps=1,nfin read(30,*) omf(ieps),ceps1(ieps) enddo close(30) else if (nmat.eq.3) then ! Gold data c OPEN(UNIT=30,FILE='Au293Knew.dat') !Gold data for different T OPEN(UNIT=30,FILE='Audat.dat') !Gold data in nm write(6,*)'Gold particles' 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) cc else if (nmat.eq.4) then else if (nmat.eq.5) then ! Copper data OPEN(UNIT=30,FILE='Cudat.dat') !Copper data in nm write(6,*)'Copper particles' rewind(30) do ieps=1, nfin read(30,*) omf(ieps),ceps1(ieps) omf(ieps)=2.d0*pi*rsnm/(omf(ieps)*rmuf) enddo close(30) else if (nmat.eq.6) then ! Aluminium data OPEN(UNIT=30,FILE='Aldat.dat') !Aluminium data in nm write(6,*)'Aluminum particles' rewind(30) do ieps=1, nfin read(30,*) omf(ieps),ceps1(ieps) omf(ieps)=2.d0*pi*rsnm/(omf(ieps)*rmuf) enddo close(30) else if (nmat.eq.7) then ! Platinum data OPEN(UNIT=30,FILE='Ptdat.dat') !Platinum data in nm write(6,*)'Platinum particles' rewind(30) do ieps=1, nfin read(30,*) omf(ieps),ceps1(ieps) omf(ieps)=2.d0*pi*rsnm/(omf(ieps)*rmuf) enddo close(30) else if (nmat.eq.8) then ! Silicon data c OPEN(UNIT=30,FILE='sieps.dat') !Silicon data in nm OPEN(UNIT=30,FILE='Sidat.dat') !Silicon data in nm for larger interval write(6,*)'Silicon particles' rewind(30) do ieps=1, nfin read(30,*) omf(ieps),ceps1(ieps) omf(ieps)=2.d0*pi*rsnm/(omf(ieps)*rmuf) enddo close(30) end if ! material constant reading ********************* * The sphere size parameter=(2*pi*rmuf/lambda)*eps * the vacuum sphere size parameter= xs= 2*pi*rmuf/lambda) * Here you can opt for either c write(6,*)'Give the vacuum sphere size parameter=' c read (5,*) omega c irx=1 * or *xxx---------------- c 1 write(6,*)'Give the sphere radius to vacuum wavelength ratio =' c read (5,*) omega 1 continue *----------------------------------------------------------------- if ((nmat.eq.0).or.(ynperfcon)) go to 8 !dispersionless dielectric !or ideal metal * In case of a dispersion, ZEPS1 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. reepsz=2.d0*pi*rsnm/(324.269d0*rmuf) IF (NMAT.EQ.1) THEN !Material decision IF - Drude metal plasma=reepsz omxp=plasma/omega zeps1=1.d0-omxp**2/(1.d0+ci*plasma/(144.d0*omega)) go to 5 * ELSE IF (nmat.eq.4) then !Material decision IF - ZnS * filfrac=0.62d0 ! filfrac of ZnS in ZnS core call znsrefind(LAMBDA,FILFRAC,zeps1) go to 5 * ELSE IF (NMAT.EQ.2) THEN !Material decision IF - Ag c >>> real material data: !silver * lambda_z=323.83d0 * 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 5 else if (omxf.eq.omf(1)) then zeps1=ceps1(1) go to 5 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 5 end if enddo end if * ELSE IF ((NMAT.GE.3).or.((nmat.ge.5).and.(nmat.le.7))) then !Material decision IF !Au,Cu,Al,Pt c >>> * 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 cc write(6,*)'Material data not available for this wavelength' cc stop * call sordalc(NMAT,lambda,ZEPS1) go to 5 * end if * if (omega.eq.omf(nfin)) then zeps1=ceps1(nfin) go to 5 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 5 end if enddo end if ELSE IF (NMAT.EQ.8) then !Material decision IF - Silicon c >>> * 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 5 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 5 end if enddo end if END IF ! END of Material decision IF * The end of reading real data according to Palik's book *_____________________________________ * activate Bruggeman: 5 if (ynbrug) then ff=0.8d0 z1 = (3.d0*ff-1.d0)*zeps1+(2.d0 - 3.d0*ff)*zeps0 z2 = sqrt(z1*z1 + 8.d0*zeps1*zeps0) * if (IMAG(z2).GE.0.0) then zeps1= (z1 + z2)/4.d0 else zeps1= (z1 - z2)/4.d0 end if end if zeps(ilcs)=zeps1 zeps(ilcs+2)=zeps1 !tempor for Naomi *______________________________________ 8 continue * Check iabs vs Im zeps.ne.0 consistency: if (iabs.ge.1) then if (imag(zeps(iabs)).eq.0.d0) then write(6,*)'The iabs=',iabs,'shell is not absorbing!' write(6,*)'Check the set-up!' pause end if end if * 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 >>> * CALL GNZBESS(RX(I),LMAX,jl,drjl,nl,drnl) * CALL GNRICBESSH(CQEPS(I),OMEGA*rmf(j),LMX,hl,drhl) C--------/---------/---------/---------/---------/---------/---------/-- C Returns Riccati-Bessel functions zeta,dzeta of the argument C RX(i)=CQEPS(i)*OMEGA*rmff(j) C--------/---------/---------/---------/---------/---------/---------/-- * DO 15 L=1,LMX WL(I,L)=HL(L)/SG(I) DRWL(I,L)=DRHL(L) C >>> (AS 10.1.22): UL(I,L)=RMF(j)*JL(L) DRJL(L)=SG(I)*DRJL(L) DRUL(I,L)=JL(L)+RMF(j)*DRJL(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 Calculation of the regular and asymptotic solutions. C TT1 is the transfer matrix which translates the asymptotic solution c from outside to inside the sphere C TT2 is the transfer matrix which translates the regular solution c from inside to outside the sphere * do 27 l=1,lmx * * magnetic part * C--------/---------/---------/---------/---------/---------/---------/-- tt1(1,1,l,1,j)= -ci*sg(1)*UL(2,L)*WL(1,L)*(DRWL(1,L)/WL(1,L) & - DRUL(2,L)/UL(2,L)) tt1(1,2,l,1,j)= -ci*sg(1)*WL(2,L)*WL(1,L)*(DRWL(1,L)/WL(1,L) & - DRWL(2,L)/WL(2,L)) tt1(2,1,l,1,j)= -ci*sg(1)*UL(2,L)*UL(1,L)*(-DRUL(1,L)/UL(1,L) & +DRUL(2,L)/UL(2,L)) tt1(2,2,l,1,j)= -ci*sg(1)*WL(2,L)*UL(1,L)*(-DRUL(1,L)/UL(1,L) & + DRWL(2,L)/WL(2,L)) * tt2(1,1,l,1,j)=-ci*sg(2)*UL(1,L)*WL(2,L)*(DRWL(2,L)/WL(2,L) & - DRUL(1,L)/UL(1,L)) tt2(1,2,l,1,j)=-ci*sg(2)*WL(1,L)*WL(2,L)*(DRWL(2,L)/WL(2,L) & - DRWL(1,L)/WL(1,L)) tt2(2,1,l,1,j)=-ci*sg(2)*UL(1,L)*UL(2,L)*(-DRUL(2,L)/UL(2,L) & +DRUL(1,L)/UL(1,L)) tt2(2,2,l,1,j)=-ci*sg(2)*WL(1,L)*UL(2,L)*(-DRUL(2,L)/UL(2,L) & +DRWL(1,L)/WL(1,L)) * * electric part * c if ( ((eps(j) .ge. 0) .and. (eps(j+1) .gt. 0)) .or. c 1 ((eps(j) .le. 0) .and. (eps(j+1) .lt. 0)) ) then c cp=dsqrt(eps(j)/eps(j+1)) c else if ( (eps(j) .lt. 0) .and. (eps(j+1) .gt. 0) ) then c cp=ci*dsqrt(abs(eps(j))/eps(j+1)) c else if ( (eps(j) .gt. 0) .and. (eps(j+1) .lt. 0) ) then c cp=- ci*dsqrt(eps(j)/abs(eps(j+1))) c end if cp=sqrt(zeps(j)/zeps(j+1)) C--------/---------/---------/---------/---------/---------/---------/-- tt1(1,1,l,2,j)=-ci*sg(1)*UL(2,L)*WL(1,L)* & (DRWL(1,L)/(cp*WL(1,L)) - cp*DRUL(2,L)/UL(2,L)) tt1(1,2,l,2,j)=-ci*sg(1)*WL(2,L)*WL(1,L)* & (DRWL(1,L)/(cp*WL(1,L)) - cp*DRWL(2,L)/WL(2,L)) tt1(2,1,l,2,j)=-ci*sg(1)*UL(2,L)*UL(1,L)* & (-DRUL(1,L)/(cp*UL(1,L)) + cp*DRUL(2,L)/UL(2,L)) tt1(2,2,l,2,j)=-ci*sg(1)*WL(2,L)*UL(1,L)* & (-DRUL(1,L)/(cp*UL(1,L)) + cp*DRWL(2,L)/WL(2,L)) * tt2(1,1,l,2,j)= -ci*sg(2)*UL(1,L)*WL(2,L)* & (cp*DRWL(2,L)/WL(2,L) - DRUL(1,L)/(cp*UL(1,L))) tt2(1,2,l,2,j)= -ci*sg(2)*WL(1,L)*WL(2,L)* & (cp*DRWL(2,L)/WL(2,L) - DRWL(1,L)/(cp*WL(1,L))) tt2(2,1,l,2,j)= -ci*sg(2)*UL(1,L)*UL(2,L)* & (-cp*DRUL(2,L)/UL(2,L) +DRUL(1,L)/(cp*UL(1,L))) tt2(2,2,l,2,j)= -ci*sg(2)*WL(1,L)*UL(2,L)* & (-cp*DRUL(2,L)/UL(2,L) +DRWL(1,L)/(cp*WL(1,L))) C--------/---------/---------/---------/---------/---------/---------/-- * 27 CONTINUE ! over l * 28 CONTINUE ! over shells * * Forward and backward transfer matrices * if (lcs.eq.1) go to 40 * **** DO 35 j=2,lcs DO 30 l=1,lmx * ============ TT2 (FORWARD) PART <==> TAU-MATRICES ============== * * m-part * CMX11(l)=tt2(1,1,l,1,j)*tt2(1,1,l,1,j-1) + & tt2(1,2,l,1,j)*tt2(2,1,l,1,j-1) CMX12(l)=tt2(1,1,l,1,j)*tt2(1,2,l,1,j-1) + & tt2(1,2,l,1,j)*tt2(2,2,l,1,j-1) CMX21(l)=tt2(2,1,l,1,j)*tt2(1,1,l,1,j-1) + & tt2(2,2,l,1,j)*tt2(2,1,l,1,j-1) CMX22(l)=tt2(2,1,l,1,j)*tt2(1,2,l,1,j-1) + & tt2(2,2,l,1,j)*tt2(2,2,l,1,j-1) * tt2(1,1,l,1,j)= CMX11(l) tt2(1,2,l,1,j)= CMX12(l) tt2(2,1,l,1,j)= CMX21(l) tt2(2,2,l,1,j)= CMX22(l) * * e-part * CMX11(l)=tt2(1,1,l,2,j)*tt2(1,1,l,2,j-1) + & tt2(1,2,l,2,j)*tt2(2,1,l,2,j-1) CMX12(l)=tt2(1,1,l,2,j)*tt2(1,2,l,2,j-1) + & tt2(1,2,l,2,j)*tt2(2,2,l,2,j-1) CMX21(l)=tt2(2,1,l,2,j)*tt2(1,1,l,2,j-1) + & tt2(2,2,l,2,j)*tt2(2,1,l,2,j-1) CMX22(l)=tt2(2,1,l,2,j)*tt2(1,2,l,2,j-1) + & tt2(2,2,l,2,j)*tt2(2,2,l,2,j-1) * tt2(1,1,l,2,j)= CMX11(l) tt2(1,2,l,2,j)= CMX12(l) tt2(2,1,l,2,j)= CMX21(l) tt2(2,2,l,2,j)= CMX22(l) * * ============ TT1 (BACKWARD) PART <==> CAL M-MATRICES ============== * * m-part * CMX11(l)=tt1(1,1,l,1,lcs-j+1)*tt1(1,1,l,1,lcs-j+2) + & tt1(1,2,l,1,lcs-j+1)*tt1(2,1,l,1,lcs-j+2) CMX12(l)=tt1(1,1,l,1,lcs-j+1)*tt1(1,2,l,1,lcs-j+2) + & tt1(1,2,l,1,lcs-j+1)*tt1(2,2,l,1,lcs-j+2) CMX21(l)=tt1(2,1,l,1,lcs-j+1)*tt1(1,1,l,1,lcs-j+2) + & tt1(2,2,l,1,lcs-j+1)*tt1(2,1,l,1,lcs-j+2) CMX22(l)=tt1(2,1,l,1,lcs-j+1)*tt1(1,2,l,1,lcs-j+2) + & tt1(2,2,l,1,lcs-j+1)*tt1(2,2,l,1,lcs-j+2) * tt1(1,1,l,1,lcs-j+1)= CMX11(l) tt1(1,2,l,1,lcs-j+1)= CMX12(l) tt1(2,1,l,1,lcs-j+1)= CMX21(l) tt1(2,2,l,1,lcs-j+1)= CMX22(l) * * e-part * CMX11(l)=tt1(1,1,l,2,lcs-j+1)*tt1(1,1,l,2,lcs-j+2) + & tt1(1,2,l,2,lcs-j+1)*tt1(2,1,l,2,lcs-j+2) CMX12(l)=tt1(1,1,l,2,lcs-j+1)*tt1(1,2,l,2,lcs-j+2) + & tt1(1,2,l,2,lcs-j+1)*tt1(2,2,l,2,lcs-j+2) CMX21(l)=tt1(2,1,l,2,lcs-j+1)*tt1(1,1,l,2,lcs-j+2) + & tt1(2,2,l,2,lcs-j+1)*tt1(2,1,l,2,lcs-j+2) CMX22(l)=tt1(2,1,l,2,lcs-j+1)*tt1(1,2,l,2,lcs-j+2) + & tt1(2,2,l,2,lcs-j+1)*tt1(2,2,l,2,lcs-j+2) * tt1(1,1,l,2,lcs-j+1)= CMX11(l) tt1(1,2,l,2,lcs-j+1)= CMX12(l) tt1(2,1,l,2,lcs-j+1)= CMX21(l) tt1(2,2,l,2,lcs-j+1)= CMX22(l) * 30 CONTINUE ! over l * 35 CONTINUE ! over shells * 40 CONTINUE * ************************* RADIAL SCAN ************************* *** In this part, DRHL, DRJL are derivatives with respect to *** the Bessel function arguments *>>> scanning the interval RINT: do 200 il=0,nx !nx !over the interval RINT cx if (nx.eq.0) then cx delo=1.d0 cx else if (nx.ne.0) then delo=rint/dble(nx) cx end if rmff= 1.d-2 + dble(il)*delo cd rmff=(rsnm+1.d0)/rsnm cx rmff= 0.8771d0 + dble(il)*delo * identifying the dipole layer for a given rmff do 45 ml=1,lcs if (rmff.lt.rmf(ml)) then j=ml go to 50 end if 45 continue 50 if (rmff.gt.rmf(lcs)) j=lcs+1 * * Test if an absorption layer is present: * ct if (j.eq.iabs) go to 200 if ((j.eq.iabs).or.(j.eq.4)) go to 200 !tempor Naomi !Otherwise iabs shell is absorptive, !or Dipole layers is different from absorptive layer' ****** * Initialize frequency shifts: gperp=0.d0 gpar= 0.d0 * Initialize total decay rates: gqmpar= 0.d0 gqmperp= 0.d0 * Assign Bessel functions: CQEPS(1)=SQRT(ZEPS(j)) SG(1)=OMEGA*CQEPS(1) !in the dipole medium RX(1)=SG(1)*rmff * CALL GNZBESS(RX(1),LMAX,jl,drjl,nl,drnl) * DO 150 L=1,LMX !summation over l HL(L)=JL(L)+ci*NL(L) DRHL(L)=HL(L)+RX(1)*(DRJL(L)+ci*DRNL(L)) DRJL(L)=JL(L)+RX(1)*DRJL(L) * * Setting up the P and Q coefficients of the F multipoles for * the D expansion coefficient: * * [Follows Sec. 4.1.4 of my Ann. Phys. (NY) 315, 352-418 (2005)] * * TT1 <==> CAL M-matrices, TT2 <==> TAU-matrices if (j.eq.1) then ! dipole in the core ! (see Eqs. (63),(71) of [1]) cxm(l)=cone/tt1(2,2,l,1,j) ! F-multipole P-coefficient dxm(l)=czero ! F-multipole Q-coefficient cxe(l)=cone/tt1(2,2,l,2,j) dxe(l)=czero else if (j.eq.lcs+1) then ! dipole outside the sphere for r>r_d ! (see Eqs. (66),(72) of [1]) cxm(l)=cone ! F-multipole P-coefficient dxm(l)=tt2(2,1,l,1,j-1)/tt2(1,1,l,1,j-1) ! F-multipole Q-coefficient cxe(l)=cone dxe(l)=tt2(2,1,l,2,j-1)/tt2(1,1,l,2,j-1) else ! dipole in the generic shell ! (see Eqs. (60),(70) of [1]) zdet1=tt1(2,2,l,1,j)*tt2(1,1,l,1,j-1) - & tt1(1,2,l,1,j)*tt2(2,1,l,1,j-1) cxm(l)=tt2(1,1,l,1,j-1)/zdet1 ! F-multipole P-coefficient dxm(l)=tt2(2,1,l,1,j-1)/zdet1 ! F-multipole Q-coefficient zdet2=tt1(2,2,l,2,j)*tt2(1,1,l,2,j-1) - & tt1(1,2,l,2,j)*tt2(2,1,l,2,j-1) C--------/---------/---------/---------/---------/---------/---------/-- cxe(l)=tt2(1,1,l,2,j-1)/zdet2 ! F-multipole P-coefficient dxe(l)=tt2(2,1,l,2,j-1)/zdet2 ! F-multipole Q-coefficient end if C A-coefficient: if (j.eq.lcs+1) then !Eq. (65) cfml1(l)=czero cfel1(l)=czero cdrfel1(l)=czero else !Eq. (86) cfml1(l)=tt1(1,2,l,1,j)*(cxm(l)*jl(l) + dxm(l)*hl(l)) cfel1(l)=tt1(1,2,l,2,j)*(cxe(l)*jl(l) + dxe(l)*hl(l)) cdrfel1(l)=tt1(1,2,l,2,j)*(cxe(l)*drjl(l) + dxe(l)*drhl(l)) cc else if (j.eq.1) then !2nd variant of the A-coefficient cc cc cfml1(l)=-tt2(1,2,l,1,lcs)*jl(l)/tt2(1,1,l,1,lcs) cc cfel1(l)=-tt2(1,2,l,2,lcs)*jl(l)/tt2(1,1,l,2,lcs) cc cdrfel1(l)=-tt2(1,2,l,2,lcs)*drjl(l)/tt2(1,1,l,2,lcs) end if C B-coefficient: if (j.eq.1) then !Eq. (144) cfml2(l) = czero cfel2(l) = czero cdrfel2(l)= czero else if (j.eq.lcs+1) then !Eq. (65) cfml2(l)=dxm(l)*hl(l) cfel2(l)=dxe(l)*hl(l) cdrfel2(l)=dxe(l)*drhl(l) else !Eqs. (87,145) cfml2(l)=tt2(2,1,l,1,j-1)*(tt1(1,2,l,1,j)*jl(l) & + tt1(2,2,l,1,j)*hl(l))/zdet1 cfel2(l)=tt2(2,1,l,2,j-1)*(tt1(1,2,l,2,j)*jl(l) & + tt1(2,2,l,2,j)*hl(l))/zdet2 cdrfel2(l)=tt2(2,1,l,2,j-1)*(tt1(1,2,l,2,j)*drjl(l) & + tt1(2,2,l,2,j)*drhl(l))/zdet2 end if * * Eqs. (146): * cfml(l)=cfml1(l)*jl(l) + cfml2(l)*hl(l) cfel(l)=cfel1(l)*jl(l) + cfel2(l)*hl(l) cdrfel(l)=cdrfel1(l)*drjl(l) + cdrfel2(l)*drhl(l) * *================= Trace over angular momenta =================== * Final multiplication by dipole Bessel functions * FREQUENCY SHIFTS & TOTAL DECAY RATES * Perpendicular (radially oscillating) dipole if (l.gt.1) then xd1=gperp xd2=gpar end if zgperp=dble(l*(l+1)*(2*l+1))*cfel(l)/rx(1)**2 gperp=gperp+3.d0*dimag(zgperp)/4.d0 gqmperp=gqmperp+3.d0*dble(zgperp)/2.d0 * Parallel (tangentially oscillating) dipole zgpar=dble(2*l+1)*(cfml(l) + cdrfel(l)/rx(1)**2) gpar=gpar+3.d0*dimag(zgpar)/8.d0 gqmpar=gqmpar+3.d0*dble(zgpar)/4.d0 * Convergence discussion: if (l.gt.1) then if (gperp.ne.0.) tol1=abs((gperp-xd1)/gperp) if (gpar.ne.0.) tol2=abs((gpar-xd2)/gpar) if ((tol1.le.1.d-8).and.(tol2.le.1.d-8)) goto 160 end if 150 CONTINUE ! over l's * END of trace over the angular-momentum indices ********************************************************************* 160 CONTINUE !desired convegence has been reached * Final convergence discussion: if (lmx.gt.1) then if (gperp.ne.0.) tol1=abs((gperp-xd1)/gperp) if (gpar.ne.0.) tol2=abs((gpar-xd2)/gpar) if (tol1.gt.tol) WRITE(NOUT,*)'#Warning: tol1.gt.tol,tol1=',tol1 if (tol2.gt.tol) WRITE(NOUT,*)'#Warning: tol2.gt.tol,tol2=',tol2 if (tol1.gt.tol) & WRITE(NOUT,*)'#Warning: tol1.gt.tol,tol1=',tol1 if (tol2.gt.tol) & WRITE(NOUT,*)'#Warning: tol2.gt.tol,tol2=',tol2 end if C--------/---------/---------/---------/---------/---------/---------/-- gim=gperp+2.d0*gpar !See Eq. (7) of PRA34_3410 dospar(il+1) = gpar !tangential oscillations dosperp(il+1)= gperp !radial oscillations dos(il+1) = gim/3.d0 !polarization average: there is !one radial and two indep. tangential !orientations possible WRITE(NOUT,999) rmff, & dospar(il+1), & dosperp(il+1), & dos(il+1) c >>> QM decay rates directly from Green's function gim=3.d0+gqmperp+2.d0*gqmpar !See Eq. (7) of PRA34_3410 dospar(il+1) = 1.d0+gqmpar !tangential oscillations dosperp(il+1)= 1.d0+gqmperp !radial oscillations dos(il+1) = gim/3.d0 !polarization average: there is !one radial and two indep. tangential !orientations possible WRITE(NOUT+7,999) rmff, & dospar(il+1), & dosperp(il+1), & dos(il+1) 200 CONTINUE !over r 300 CONTINUE !over rsnm 400 CONTINUE !over lambda close(nout) close(nout+7) write(6,*)'The frequency shift in dpfshiftm.dat' write(6,*)'The total decay rates in dipqmratesm.dat' write(6,*) write(6,*)'Warnings are issued when' write(6,*)'the convergence is worse than up to 3 exact digits.' write(6,*)'This should only happen in dpfshiftm.dat in the' write(6,*)'immediate proximity of metal shell boundaries.' write(6,*)'The respective values of tol1, tol2 parameters will' write(6,*)'tell you which precision the results' write(6,*)'were obtained with.' write(6,*) write(6,*)'Improvements, comments, additions: all welcome at' write(6,*)'wavescattering@yahoo.com' 999 format(f9.6,3X,3(e15.8,3x)) c 999 format(f8.6,3X,e23.16,3x,2(e15.8,3x)) END subroutine sordalc(NMAT,LAMBDA,ZEPS) C---------------------------------------------------------------------- C SUBROUTINE TO CALCULATE THE DIELECTRIC CONSTANT OF METALS C ACCORDING TO AN ARTICLE BY M. A. Ordal et al, C Optical properties of fourteen metals in the infrared and far infrared: C Al, Co, Cu, Au, Fe, Pb, Mo, Ni, Pd, Pt, Ag, Ti, V, and W, C Appl. Opt. {\bf 22}, 1099 (1983); ibid. {\bf 24}, 4493 (1985) C C f77 -g -check_bounds ordalc.f -o rnordalc C C omega=1/\lambda [cm^{-1}] in Ordal [spectroscopic convention] C C Re eps = - \fr{\om_p^2}{\om^2+\om_\tau^2} C Im eps = \fr{\om_p^2 \om_\tau}{\om^3+\om \om_\tau^2} C 1eV=1.24 \mu m C---------------------------------------------------------------------- IMPLICIT NONE INTEGER NMAT REAL*8 lambda,plasma,tau COMPLEX*16 zeps C ------------------------------- REAL*8 pi,x,y,omega DATA PI/3.141592653589793d0/ C --------- C ::: speed of light in vacuum in nm/s C PARAMETER (c0=2.9927925d17) C C According to Table I of Ordal et al, Appl. Opt. {\bf 24}, 4493 (1985): C (plasma=omega_plasma and tau=omega_tau below) C C plasma[THz]/plasma[cm-1] tau[THz]/tau[cm-1] C Al 3570/119000 19.4/660 C Cu 1914/59600 8.34/73.2 C Au 2175/72800 6.5/215 C Ag 2175/72700 4.35/145 C Pt /41500 /558 C C ::: conversion factor between normal angular frequency and eV: c PARAMETER(XCN=4.13566727d-15) if (nmat.eq.3) then !Au C ::: Convert the plasma frequency from Table I of Ordal et al C ::: from [cm-1] to [nm-1] ===> conversion factor 10^{-7} PLASMA=72800.d-7 !d12 in Hz/d-7 in [nm-1] C ::: Convert the tau frequency from Table I of Ordal et al C ::: from [cm-1] to [nm-1] ===> conversion factor 10^{-7} TAU=215.d-7 else if (nmat.eq.5) then !Cu PLASMA=59600d-7 TAU=73.2d-7 else if (nmat.eq.6) then !Al PLASMA=119000d-7 TAU=660.d-7 else if (nmat.eq.7) then !Pt PLASMA=41150d-7 TAU=558.d-7 end if C ------------------------------- c write(6,*)'Read in wavelength in nm' c read(5,*) lambda omega=1.d0/lambda * X=-plasma**2/(omega**2+tau**2) Y=tau*plasma**2/(omega**3+omega*tau**2) * zeps=dcmplx(X,Y) * END * SUBROUTINE znsrefind(LAMBDA,FILFRAC,zeps) C--------/---------/---------/---------/---------/---------/---------/-- C FILFRAC ... filfrac of ZnS in ZnS core C--------/---------/---------/---------/---------/---------/---------/-- IMPLICIT NONE COMPLEX*16 zeps REAL*8 F,FILFRAC,LAMBDA,EPSHOST REAL*8 EPSZnS REAL*8 eMG1,EPSPAR * c if the host is different from the medium (silica n=1.45) EPSHOST = 1.45d0*1.45d0 c ZnS filling fraction f f = FILFRAC c Particle material: ZnS bulk dielectric constant 350 nm - 900 nm EPSZnS = 5.164d0 + 1.208d+7/(lambda*lambda*100 - 0.732d+7) * particle material EPSPAR = EPSZnS c Bruggeman (effective medium host) c wor = (3.d0*F-1.d0)*EPSPAR+(2.d0 - 3.d0*F)*EPSHOST c wor1 = sqrt(wor*wor + 8.d0*EPSPAR*EPSHOST) c if (AIMAG(wor1).GE.0.0) then c eBr = (wor + wor1 )/4.d0 c else c e_Br = (wor - wor1 )/4.d0 c end if c Maxwell-Garnett MG1 (medium material host) eMG1=EPSHOST*(2.d0*EPSHOST+EPSPAR + 2.d0*F*(EPSPAR - EPSHOST))/ & (2.d0*EPSHOST + EPSPAR - F*(EPSPAR - EPSHOST)) C--------/---------/---------/---------/---------/---------/---------/-- c Maxwell-Garnett MG2 (particle material host) C eMG2 = EPSPAR*(2.*EPSPAR + EPSHOST + 2.*(1.- F)*(EPSHOST - EPSPAR))/ C & (2.*EPSPAR + EPSHOST -(1.- F)*(EPSHOST - EPSPAR)) c WRITE(*,*) LAMBDA, eMG1 * zeps = eMG1 * RETURN END * C (C) Copr. 8/2004 Alexander Moroz