C (C) Copr. 06/2004 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 CHEW, a computer code to calculate normalized spontaneous C emission decay rates and the Ohmic loss contribution to the c nonradiative decay rates. CHEW calculates these quantities for C two principal (tangential and radial) and averaged dipole orientations. C The quantities are normalized with respect to the radiative decay C rates of the same dipole either in a free-space (vacuum) or in a free-space C filled in with a homogeneous medium which is identical to that at C the dipole position. C C CHEW code 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 The program default options correspond to the cases considered in the C article. Default option is to have three concentric shell. C The first and third shells have the dielectric constant of silica, C whereas the second has 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 three spherical shell radii in nm; C C 2) the vacuum (detection) wavelength lambda. C C The real dielectric constants of the nanoshell core and 3rd shell, and C of ambient are specified as parameters. 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 have been published in follow up's: 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 [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 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 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 ======================================================================== C OUTPUT: C C CHEW returns files C C dipratesm.dat, dipratesv.dat containing the radiative decay rates C C and C C dipnrratesm.dat, dipnrratesv.dat.dat, which contain the C relevant Ohmic loss contribution to the nonradiative decay rates. C C Each output file has data organized in four columns: C normalized radial position (in the units of the bead radius) and C 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) IN THE CASE OF A HOMOGENEOUS SPHERE THE RESULTS SHOULD COINCIDE WITH C 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 2) For sphere radius rsnm=300 nm, inner core radius 200 nm, the C first-second shell interface at 248 nm, and wavelength of 595 nm C the first result lines of dipratesm.dat with default options C should read as follows: C C #Sphere radius in nm= 300.000000000000 C 0.010000 0.65263317E-01 0.65680119E-01 0.65402251E-01 C 0.019950 0.68834444E-01 0.70495002E-01 0.69387963E-01 C 0.029900 0.74755963E-01 0.78492106E-01 0.76001344E-01 C 0.039851 0.82994428E-01 0.89646012E-01 0.85211623E-01 C 0.049801 0.93503336E-01 0.10392131E+00 0.96975993E-01 C 0.059751 0.10622341E+00 0.12127273E+00 0.11123985E+00 C 0.069701 0.12108297E+00 0.14164537E+00 0.12793710E+00 C 0.079652 0.13799836E+00 0.16497488E+00 0.14699053E+00 C 0.089602 0.15687452E+00 0.19118776E+00 0.16831227E+00 C 0.099552 0.17760549E+00 0.22020170E+00 0.19180423E+00 C C whereas those of dipnrratesm.dat should read as follows: C C #Sphere radius in nm= 300.000000000000 C 0.010000 0.83692940E-01 0.84234047E-01 0.83873309E-01 C 0.019950 0.88357112E-01 0.90512902E-01 0.89075709E-01 C 0.029900 0.96090989E-01 0.10094135E+00 0.97707775E-01 C 0.039851 0.10685075E+00 0.11548595E+00 0.10972915E+00 C 0.049801 0.12057546E+00 0.13410013E+00 0.12508368E+00 C 0.059751 0.13718743E+00 0.15672431E+00 0.14369972E+00 C 0.069701 0.15659272E+00 0.18328623E+00 0.16549056E+00 C 0.079652 0.17868169E+00 0.21370116E+00 0.19035484E+00 C 0.089602 0.20332966E+00 0.24787228E+00 0.21817720E+00 C 0.099552 0.23039771E+00 0.28569106E+00 0.24882883E+00 C C 3) You can also download executable from C http://www.wave-scattering.com/chew.exe C and test your results against the executable output. 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 3) For running of the program you will also need the following C copyrighted routines from 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 qromo(func,a,b,ss,choose), C which returns as ss the integral of the function func from a to b, C together with its dependencies polint.f and midpnt.f C C which can all be found on C C http://www.nr.com C C The routines have to be transformed into double precision. 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 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 chew C--------/---------/---------/---------/---------/---------/---------/-- C 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,nabs real*8 TOL,OMEGA,RINT,pi,rmuf complex*16 ZEPS0,CCEPS,CSEPS,cp character*1 ync logical ynperfcon,ynbrug,ynint 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. 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: PARAMETER (NX=201) C ::: radial interval on which the decay rates are calculated in the units C of the sphere radius (RMUF=1) 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=3) * c The coating layer to which material data are read in parameter (ilcs=2) c Total number of absorbing layers. Otherwise, if there is no absorbing layer, c set NABS=1 parameter (nabs=1) 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 <<< 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)) !glycerol=2.1609d0 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 i,j,l,im,il,ml,irx,lmx,ieps,l1,l2,ikl,irs,idt,ild REAL*8 xd1,xd2,xnrd1,xnrd2,tol1,tol2,tolnr1,tolnr2,gim, & rsnm,rsnm0,lambda,lambda0,xinf,xsup,xmin, & fintfm,fintfe1,fintfe2,gperp,gpar,gnrperp,gnrpar,gnrim, & rmff,delo,ximl,ximle1,ximle2,xoutf,xinnrf,xoutnrf, & omxf,reepsz,plasma,omxp,ff,filfrac !ff for Bruggemann;filfrac for ZnS complex*16 ci,cone,czero,zdet1,zdet2,zeps1,Z1,Z2,zz,zz1 COMPLEX*16 ZSIGMA,ZSIGME1,ZSIGME2 * Array declarations: REAL*8 RMF(lcs),rff(lcs),xim(lmax,nabs),xie(lmax,nabs) 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 exm,fxm,exe1,fxe1,exe2,fxe2 COMPLEX*16 cfel(1:lmax),cfml(1:lmax),cdrfel(1:lmax) COMPLEX*16 cfnrel(1:lmax),cfnrml(1:lmax),cdrfnrel(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-------------------------------------------------------------------- * COMMON/TOFINTFMI/ L COMMON/TOFINTFEI1/ L1 COMMON/TOFINTFEI2/ l2 COMMON/TOFINTFMC/ ZSIGMA,exm,fxm COMMON/TOFINTFEC1/ ZSIGME1,exe1,fxe1 COMMON/TOFINTFEC2/ ZSIGME2,exe2,fxe2 * external fintfm,fintfe1,fintfe2,midpnt * 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 * * 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: * write(6,*)'Read sphere radius rsnm in nm' read(5,*) rsnm cy rsnm=150.d0 if(lcs.ge.2) then write(6,*)'Read inner core radius in nm' read(5,*) rff(1) rff(1)=rff(1)/rsnm rmf(1)=rff(1)*rmuf end if * 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 ZEPS(1)=cceps if ((lcs.gt.1).and. & ((nmat.eq.0).or.((nmat.ne.0).and.(ilcs.ne.lcs)))) & zeps(lcs)=cseps !tempor Naomi zeps(lcs+1)=zeps0 * * radii of the coated sphere: * if (lcs.ge.3) then !multicoated particle * write(6,*)'Coated sphere core radii r(l) labelled from 1 for the & inner core till LCS for the outer shell & ===> r(lcs) is the sphere radius' C--------/---------/---------/---------/---------/---------/---------/-- do ikl=2,lcs-1,1 * write(6,*)'Read in r(l) for l=',ikl read (5,*) rff(ikl) rff(ikl)=rff(ikl)/rsnm rmf(ikl)=rff(ikl)*rmuf * if ((nmat.eq.0).or.((nmat.gt.0).and.(ikl.ne.ilcs))) then write(6,*)'Read in the lth-sphere layer diel. const. for l=',ikl cy write(6,*)'(In case of dispersive component, give 1. )' write(6,*) & '(In case of an absorbing component, complex number in brackets)' read (5,*) zeps(ikl) end if * enddo * end if write(6,*)'Give the vacuum (detection) wavelength =' read (5,*) lambda cy lambda=595.d0 C--------/---------/---------/---------/---------/---------/---------/-- * The header of diprates*.dat output files: OPEN(UNIT=NOUT,FILE='dipratesm.dat') rewind(NOUT) WRITE(NOUT,*)'#Dipole radiated power 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=', zeps(lcs+1) if (nmat.ge.1) write(nout,*)'#Dispersive layer number=', ilcs if ((lcs.eq.1).and.(nmat.eq.0)) & write(nout,*)'#sphere diel. const.=', zeps(1) if (lcs.ge.2) then write(nout,*)'#core radius/sphere radius =',rff(1) if ((ilcs.ne.1).or.(nmat.eq.0)) & write(nout,*)'#sphere core diel. const.=', zeps(1) if ((ilcs.ne.lcs).or.(nmat.eq.0)) & write(nout,*)'#coating diel. const.=',zeps(lcs) end if 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,*)'#Decay rates in the free case (in background zeps0)', c 1 freeld WRITE(NOUT,*)'#In columns:' WRITE(NOUT,*)'#position, the dipole radiated power' write(nout,*)'#for tangentially and radially oriented dipole,' write(nout,*)'#and their orientational average, respectively' cxxxxxxxxxxxx OPEN(UNIT=NOUT+5,FILE='dipratesv.dat') rewind(NOUT+5) WRITE(NOUT+5,*)'#Dipole radiated power in the presence of a & coated sphere as a function of the distance normalized to the & decay rates in a free space filled in with the surrounding & sphere host medium' write(NOUT+5,*)'#Angular momentum cut-off LMAX=',LMX WRITE(NOUT+5,*)'#the vacuum (detection) wavelength =', lambda write(nout+5,*) '# rsnm/wavelength=', rsnm/lambda write(nout+5,*)'#Material number =',NMAT if (ynbrug) write(nout+5,*)'#Bruggeman approximation performed' WRITE(NOUT+5,*)'#In columns:' WRITE(NOUT+5,*)'#position, the dipole radiated power' write(nout+5,*)'#for tangentially and radially oriented dipole,' write(nout+5,*)'#and their orientational average, respectively' cxxxxxxxxxxxx OPEN(UNIT=NOUT+6,FILE='dipnrratesm.dat') rewind(NOUT+6) write(nout+6,*)'#The power loss due to Ohmic losses in the & presence of a single coated sphere as a function of the distance & normalized to radiative decay rates in a free space filled in & with the local medium' write(nout+6,*)'#coated sphere as a function of the distance' * WRITE(NOUT+6,*)'#the vacuum (detection) wavelength =', lambda write(nout+6,*)'#Sphere radius rmuf in relative units=', rmuf write(nout+6,*) '# rsnm/wavelength=', rsnm/lambda c WRITE(NOUT+6,*) write(NOUT+6,*)'#Angular momentum cut-off LMAX=',LMX write(nout+6,*)'#Material number =',NMAT if (ynbrug) write(nout+6,*)'#Bruggeman approximation performed' if (ync .eq.'n') then write(nout+6,*)'# Homogeneous sphere' WRITE(NOUT+6,*) '#Set ync .eq. y if you want coated sphere!' else if (ync.eq.'y') then write(nout+6,*)'# Coated sphere' WRITE(NOUT+6,*) '#Set ync .eq. n if you want homogeneous sphere!' end if * write(nout+6,*)'#Numbers of layers lcs=', lcs write(nout+6,*)'#host dielectric constant=', zeps(lcs+1) if (nmat.ge.1) write(nout+6,*)'#Dispersive layer number=', ilcs if ((lcs.eq.1).and.(nmat.eq.0)) & write(nout+6,*)'#sphere diel. const.=', zeps(1) if (lcs.ge.2) then write(nout+6,*)'#core radius/sphere radius =',rff(1) if ((ilcs.ne.1).or.(nmat.eq.0)) & write(nout+6,*)'#sphere core diel. const.=', zeps(1) if ((ilcs.ne.lcs).or.(nmat.eq.0)) & write(nout+6,*)'#coating diel. const.=',zeps(lcs) end if C--------/---------/---------/---------/---------/---------/---------/-- write(nout+6,*)'#Length of the searched interval rint=', rint write(nout+6,*)'#Number of the points on the uniform mesh nx=',nx c WRITE(NOUT+6,*)'#Decay rates in the free case (in background zeps0)', c 1 freeld WRITE(NOUT+6,*)'#In columns:' WRITE(NOUT+6,*)'#position, dipole power loss due Ohmic losses' write(nout+6,*)'#for tangentially and radially oriented dipole,' write(nout+6,*)'#and their orientational average, respectively' cxxxxxxxxxxxx OPEN(UNIT=NOUT+7,FILE='dipnrratesv.dat') rewind(NOUT+7) write(nout+7,*)'#The power loss due to Ohmic losses in the & presence of a single coated sphere as a function of the distance & normalized to radiative decay rates in a free space filled in & with the surrounding sphere host medium' WRITE(NOUT+7,*)'#the vacuum (detection) wavelength =', lambda write(nout+7,*)'#Sphere radius rmuf in relative units=', rmuf write(nout+7,*) '# rsnm/wavelength=', rsnm/lambda c WRITE(NOUT+7,*) 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 * write(nout+7,*)'#Numbers of layers lcs=', lcs write(nout+7,*)'#host dielectric constant=', zeps(lcs+1) if (nmat.ge.1) write(nout+7,*)'#Dispersive layer number=', ilcs if ((lcs.eq.1).and.(nmat.eq.0)) & write(nout+7,*)'#sphere diel. const.=', zeps(1) if (lcs.ge.2) then write(nout+7,*)'#core radius/sphere radius =',rff(1) if ((ilcs.ne.1).or.(nmat.eq.0)) & write(nout+7,*)'#sphere core diel. const.=', zeps(1) if ((ilcs.ne.lcs).or.(nmat.eq.0)) & write(nout+7,*)'#coating diel. const.=',zeps(lcs) end if C--------/---------/---------/---------/---------/---------/---------/-- write(nout+7,*)'#Length of the searched interval rint=', rint write(nout+7,*)'#Number of the points on the uniform mesh nx=',nx WRITE(NOUT+7,*)'#In columns:' WRITE(NOUT+7,*)'#position, dipole power loss due Ohmic losses' write(nout+7,*)'#for tangentially and radially oriented dipole,' write(nout+7,*)'#and their orientational average, respectively' 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) write(nout,*)'#Sphere radius in nm=', rsnm write(nout+5,*)'#Sphere radius in nm=', rsnm write(nout+6,*)'#Sphere radius in nm=', rsnm write(nout+7,*)'#Sphere radius in nm=', rsnm WRITE(NOUT,*) WRITE(NOUT+5,*) WRITE(NOUT+6,*) WRITE(NOUT+7,*) omega=2.d0*pi*rsnm/(lambda*rmuf) !size parameter for default ! rmuf=1 * 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 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 7 zeps(ilcs)=zeps1 cn zeps(ilcs+2)=zeps1 !tempor Naomi *______________________________________ 8 continue * Check iabs vs Im zeps.ne.0 consistency: cn iabs=2 !tempor for Naomi 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) DRUL(I,L)=JL(L)+RX(I)*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,1 DO 30 l=1,lmx * ========= TT2 (FORWARD) PART - {cal T}-matrices (Eq. (30)) =========== * {cal T}-matrices below are defined as {cal T}(n) = \prod_{j=1}^n T^+(j), * a slight modification of Eq. (30)) * * 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 (Eq. (30)) =========== * * 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 * *** In this part, DRHL, DRJL are derivatives with respect to *** the Bessel function arguments if (iabs.gt.0) then ZSIGMA=OMEGA*SQRT(ZEPS(iabs)) ZSIGME1=ZSIGMA ZSIGME2=ZSIGMA end if ynint=.true. *>>> scanning the interval RINT: do 200 il =0,nx !nx !over the interval RINT delo=rint/dble(nx) 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 CQEPS(1)=SQRT(ZEPS(j)) SG(1)=OMEGA*CQEPS(1) !in the dipole medium RX(1)=SG(1)*rmff ** Common prefactors for the radiative and nonradiative decay rates xinf=abs(CQEPS(1))/sqrt(ZEPS0) !refractive index contrast - !to dipole medium normalization xoutf=abs(ZEPS(j))/ZEPS0 !diel. constant contrast - !to ambient normalization cn iabs=2 !tempor Naomi if ((iabs.gt.0).and.(iabs.ne.j)) then * xinnrf comes here as the imaginary part of ZEPS(IABS) xinnrf=imag(zeps(iabs))*abs(SG(1))**3/abs(ZEPS(j)) !to dipole medium normalization xoutnrf=xinnrf*abs(CQEPS(1))/abs(sqrt(ZEPS0)) !to ambient normalization if (ynint) idt=j if (idt.ne.j) then ynint=.true. idt=j end if end if ** C In a nonmagnetic medium, xoutf=xinf**2=xinf* differ by the factor * CALL GNZBESS(RX(1),LMAX,jl,drjl,nl,drnl) * * Initialize radiative and nonradiative decay rates: gperp= 0.d0 gpar= 0.d0 gnrperp=0.d0 gnrpar=0.d0 * Assign Bessel functions DO 150 L=1,LMX !summation over l * Initialize common variables L1,L2 l1=l l2=l HL(L)=JL(L)+ci*NL(L) DRHL(L)=DRJL(L)+ci*DRNL(L) * * For the calculation of the dipole transitient rates, only * the D-coefficient is relevant. * Setting up the P and Q coefficients of the F multipoles for * the D expansion coefficient in the radiation zone. * cxm/e ... prefactors of a_{\gamma L}^d * dxm/e ... prefactors of \alpha_{\gamma L}^d * * [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) 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--------/---------/---------/---------/---------/---------/---------/-- C D-coefficients in terms of the F-multipoles (for j=lcs+1 for r>r_d): cfml(l)=cxm(l)*jl(l) + dxm(l)*hl(l) cfel(l)=cxe(l)*jl(l) + dxe(l)*hl(l) cdrfel(l)=cxe(l)*drjl(l) + dxe(l)*drhl(l) * if (iabs.le.0) go to 146 ct if ((iabs.le.0).or.(iabs.eq.j)) go to 146 * do 145 ikl=1,nabs * cn iabs=2*ikl !tempor Naomi if (.not.ynint) go to 140 * Test if an absorption layer is present: * if (j.eq.iabs) go to 140 cn if ((j.eq.2).or.(j.eq.4)) go to 140 !tempor Naomi !Otherwise iabs shell is absorptive, !or Dipole layers is different from absorptive layer' !or shell integration has been performed * * Setting up the coefficients of the G multipole in the absorptive shell * from the known P and Q coefficients of the F multipole in the radiation zone: * * By default, the absorptive shell option is only considered in the absence * of the dipole source * if (j.lt.iabs) then !only D-coef is nonzero outside the sphere ! (see Eq. (74)) * exm/e, fxm/e below are prefactors of D-coef * (backward transfer matrix applied:) exm= tt1(1,2,l,1,iabs) fxm= tt1(2,2,l,1,iabs) exe1= tt1(1,2,l,2,iabs) fxe1= tt1(2,2,l,2,iabs) * else if (j.eq.lcs+1) then !both C and D-coef nonzero outside the sphere * ( C ) ( 1 ) * ( ) = ( ) x dipole alpha coefficient * ( D ) ( dxm/e ) * * exm/e, fxm/e below are prefactors of alpha coefficient * (backward transfer matrix applied (see Eq. (73))) exm= tt1(1,1,l,1,iabs)+tt1(1,2,l,1,iabs)*dxm(l) fxm= tt1(2,1,l,1,iabs)+tt1(2,2,l,1,iabs)*dxm(l) exe1= tt1(1,1,l,2,iabs)+tt1(1,2,l,2,iabs)*dxe(l) fxe1= tt1(2,1,l,2,iabs)+tt1(2,2,l,2,iabs)*dxe(l) * else if ((j.gt.iabs).and.(j.lt.lcs+1)) then !transfer from the core !to iabs * exm/e, fxm/e below are prefactors of A(1) in Eqs. (80) if (iabs.eq.1) then exm= cone fxm= czero exe1= cone fxe1= czero else if (iabs.gt.1) then * forward transfer matrix applied: exm= tt2(1,1,l,1,iabs-1) fxm= tt2(2,1,l,1,iabs-1) exe1= tt2(1,1,l,2,iabs-1) fxe1= tt2(2,1,l,2,iabs-1) end if !iabs * end if !j * if (iabs.eq.1) then !there are no singular modes in the core fxm=czero fxe1=czero end if * exe2= exe1 fxe2= fxe1 C--------/---------/---------/---------/---------/---------/---------/-- * Integrate positive function over the absorptive layer: * if (iabs.gt.1) then xmin=rmf(iabs-1) else xmin=0.d0 end if xsup=rmf(iabs) call qromo(fintfm,xmin,xsup,ximl,midpnt) call qromo(fintfe1,xmin,xsup,ximle1,midpnt) call qromo(fintfe2,xmin,xsup,ximle2,midpnt) C--------/---------/---------/---------/---------/---------/---------/-- C Returns the integral of the function fintf** from xinf to xsup. C Integration is performed by Romberg's method of order 2K, where, C e.g.,K=2 is Simpson's rule. C Internal QROMO Parameters: C EPS ... the fractional accuracy desired, as determined C by the extrapolation error estimate; C JMAX ... limits the total number of steps C K ... the number of points used in the extrapolation. C C Includes trapzd.f and polint.f C--------/---------/---------/---------/---------/---------/---------/-- * * supplying prefactors: xim(l,ikl)= ximl ximle1= ximle1*dble(l*(l+1))/abs(zsigma)**2 xie(l,ikl)= ximle1 + ximle2/abs(zsigma)**2 !xie is now the sum !of integrals I^1_E+I^2_E * if ((l.eq.lmx).and.(ikl.eq.nabs)) ynint=.false. * *************** Trace over angular momenta ********************* 140 IF ((iabs.gt.0).and.(iabs.ne.j)) THEN !iabs cn 140 IF ((j.eq.1).or.(j.eq.3).or.(j.eq.(lcs+1))) THEN !tempor Naomi * NON-RADIATIVE DECAY RATES: * * Final multiplication by dipole Bessel functions * * Forming linear combinations d_{\gamma L} of Bessel functions * according to Eq. (120) if (l.gt.1) then xnrd1=gnrperp xnrd2=gnrpar end if if (j.lt.iabs) then gnrperp=gnrperp+dble(3*l*(l+1)*(2*l+1))*xie(l,ikl)* & abs(cfel(l)/rx(1))**2/2.d0 gnrpar=gnrpar+3.d0*dble(2*l+1)*(xim(l,ikl)*(abs(cfml(l)))**2 + & xie(l,ikl)*abs((cfel(l) + rx(1)*cdrfel(l))/rx(1))**2)/4.d0 else if (j.eq.lcs+1) then gnrperp=gnrperp+dble(3*l*(l+1)*(2*l+1))*xie(l,ikl)* & abs(hl(l)/rx(1))**2/2.d0 gnrpar=gnrpar+3.d0*dble(2*l+1)*(xim(l,ikl)*(abs(hl(l)))**2 + & xie(l,ikl)*abs((hl(l) + rx(1)*drhl(l))/rx(1))**2)/4.d0 else if ((j.gt.iabs).and.(j.lt.lcs+1)) then ! v_{\gamma L} of Eq. (85) C A1-coefficient in terms of the F-multipoles [Eq. (85)]: cxm(l)= tt1(1,2,l,1,j)/zdet1 dxm(l)= tt1(2,2,l,1,j)/zdet1 cxe(l)= tt1(1,2,l,2,j)/zdet2 dxe(l)= tt1(2,2,l,2,j)/zdet2 cfnrml(l)=cxm(l)*jl(l) + dxm(l)*hl(l) cfnrel(l)=cxe(l)*jl(l) + dxe(l)*hl(l) cdrfnrel(l)=cxe(l)*drjl(l) + dxe(l)*drhl(l) * gnrperp=gnrperp+dble(3*l*(l+1)*(2*l+1))*xie(l,ikl)* & abs(cfnrel(l)/rx(1))**2/2.d0 gnrpar=gnrpar+3.d0*dble(2*l+1)*(xim(l,ikl)*(abs(cfnrml(l)))**2 + & xie(l,ikl)*abs((cfnrel(l) + rx(1)*cdrfnrel(l))/rx(1))**2)/4.d0 end if !j END IF !j w iabs * 145 continue !over iabs - tempor Naomi 146 continue * ************************ IF (j.ne.iabs) THEN !iabs cn IF ((j.eq.1).or.(j.eq.3).or.(j.eq.(lcs+1))) THEN !tempor Naomi * RADIATIVE DECAY RATES: * Perpendicular (radially oscillating) dipole if (l.gt.1) then xd1=gperp xd2=gpar end if gperp=gperp+dble(3*l*(l+1)*(2*l+1))*abs(cfel(l)/rx(1))**2/2.d0 * Parallel (tangentially oscillating) dipole gpar=gpar+3.d0*dble(2*l+1)*(abs(cfml(l))**2 + & abs((cfel(l) + rx(1)*cdrfel(l))/rx(1))**2)/4.d0 END IF !j w iabs if (l.gt.1) then if (gperp.ne.0.) tol1=abs((gperp-xd1)/gperp) if (gpar.ne.0.) tol2=abs((gpar-xd2)/gpar) if (iabs.gt.0) then if (gnrperp.ne.0.) tolnr1=abs((gnrperp-xnrd1)/gnrperp) if (gnrpar.ne.0.) tolnr2=abs((gnrpar-xnrd2)/gnrpar) if ((tol1.le.1.d-8).and.(tol2.le.1.d-8).and.(tolnr1.le.1.d-8) & .and.(tolnr2.le.1.d-8)) goto 160 else if ((tol1.le.1.d-8).and.(tol2.le.1.d-8)) goto 160 end if end if 150 CONTINUE ! over l's * END of trace over the angular-momentum indices ********************************************************************* 160 ynint=.false. !desired convegence has been reached * Convergence discussion: if (lmx.gt.1) then ct if (gperp.ne.0.) tol1=abs((gperp-xd1)/gperp) ct if (gpar.ne.0.) tol2=abs((gpar-xd2)/gpar) ct if (gnrperp.ne.0.) tolnr1=abs((gnrperp-xnrd1)/gnrperp) ct if (gnrpar.ne.0.) tolnr2=abs((gnrpar-xnrd2)/gnrpar) 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+5,*)'#Warning: tol1.gt.tol,tol1=',tol1 if (tol2.gt.tol) & WRITE(NOUT+5,*)'#Warning: tol2.gt.tol,tol2=',tol2 * if (iabs.gt.0) then if (tolnr1.gt.tol) & WRITE(NOUT+6,*)'#Warning: tolnr1.gt.tol,tolnr1=',tolnr1 if (tolnr2.gt.tol) & WRITE(NOUT+6,*)'#Warning: tolnr2.gt.tol,tolnr2=',tolnr2 if (tolnr1.gt.tol) & WRITE(NOUT+7,*)'#Warning: tolnr1.gt.tol,tolnr1=',tolnr1 if (tolnr2.gt.tol) & WRITE(NOUT+7,*)'#Warning: tolnr2.gt.tol,tolnr2=',tolnr2 end if * end if C--------/---------/---------/---------/---------/---------/---------/-- * OUTPUT WRITING: * * Radiative part: gim=gperp+2.d0*gpar !See Eq. (7) of PRA34_3410 dospar(il+1) = xinf*gpar !tangential oscillations dosperp(il+1)= xinf*gperp !radial oscillations dos(il+1) = xinf*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) dospar(il+1) = xoutf*gpar dosperp(il+1)= xoutf*gperp dos(il+1) = xoutf*gim/3.d0 WRITE(NOUT+5,999) rmff, & dospar(il+1), & dosperp(il+1), & dos(il+1) ***************************************** IF ((iabs.gt.0).and.(iabs.ne.j)) THEN !iabs * Nonradiative part: gnrim=gnrperp+2.d0*gnrpar dospar(il+1) = xinnrf*gnrpar !tangential oscillations dosperp(il+1)= xinnrf*gnrperp !radial oscillations dos(il+1) = xinnrf*gnrim/3.d0 !polarization average: there is !one radial and two indep. tangential !orientations possible WRITE(NOUT+6,999) rmff, & dospar(il+1), & dosperp(il+1), & dos(il+1) dospar(il+1) = xoutnrf*gnrpar dosperp(il+1)= xoutnrf*gnrperp dos(il+1) = xoutnrf*gnrim/3.d0 WRITE(NOUT+7,999) rmff, & dospar(il+1), & dosperp(il+1), & dos(il+1) END IF 200 CONTINUE !over r 300 CONTINUE !over rsnm 400 CONTINUE !over lambda close(nout) close(nout+5) close(nout+6) close(nout+7) write(6,*)'Dipole radiated power in dipratesm.dat and & dipratesv.dat' write(6,*)'Dipole power loss due to Ohmic losses in & dipnrratesm.dat and dipnrratesv.dat' 999 format(f9.6,3X,3(e15.8,3x)) c 999 format(f8.6,3X,e23.16,3x,2(e15.8,3x)) END * DOUBLEPRECISION FUNCTION fintfm(x) C--------/---------/---------/---------/---------/---------/---------/-- implicit none integer LMAX,LM PARAMETER (LMAX=60) real*8 x COMPLEX*16 ZX,ZSIGMA,exm,fxm,ci,zhl COMPLEX*16 JL(0:lmax), NL(0:lmax) COMPLEX*16 DRJL(0:lmax), DRNL(0:lmax) * COMMON/TOFINTFMI/ lm COMMON/TOFINTFMC/ ZSIGMA,exm,fxm * DATA ci/(0.d0,1.d0)/ * ZX=ZSIGMA*X CALL GNZBESS(ZX,LMAX,jl,drjl,nl,drnl) * CALL GNRICBESSH(ZSIGMA,X,LMAX,nl,drnl) * zhl=nl(lm)/zx fintfm=(exm*jl(lm)+fxm*zhl)*dconjg(exm*jl(lm)+fxm*zhl) fintfm=fintfm*x**2 return END * DOUBLEPRECISION FUNCTION fintfe1(x) C--------/---------/---------/---------/---------/---------/---------/-- implicit none integer LMAX,LM PARAMETER (LMAX=60) real*8 x COMPLEX*16 ZX,ZSIGMA,exe,fxe,ci,zhl COMPLEX*16 JL(0:lmax), NL(0:lmax) COMPLEX*16 DRJL(0:lmax), DRNL(0:lmax) COMMON/TOFINTFEI1/ lm COMMON/TOFINTFEC1/ ZSIGMA,exe,fxe * DATA ci/(0.d0,1.d0)/ * ZX=ZSIGMA*X CALL GNZBESS(ZX,LMAX,jl,drjl,nl,drnl) * CALL GNRICBESSH(ZSIGMA,X,LMAX,nl,drnl) * zhl=nl(lm)/zx fintfe1=(exe*jl(lm)+fxe*zhl)*dconjg(exe*jl(lm)+fxe*zhl) * return END * DOUBLEPRECISION FUNCTION fintfe2(x) C--------/---------/---------/---------/---------/---------/---------/-- implicit none integer LMAX,LM PARAMETER (LMAX=60) real*8 x COMPLEX*16 ZX,ZSIGMA,exe,fxe,ci,drhl,drj COMPLEX*16 JL(0:lmax), NL(0:lmax) COMPLEX*16 DRJL(0:lmax), DRNL(0:lmax) * COMMON/TOFINTFEI2/ lm COMMON/TOFINTFEC2/ ZSIGMA,exe,fxe * DATA ci/(0.d0,1.d0)/ ZX=ZSIGMA*X CALL GNZBESS(ZX,LMAX,jl,drjl,nl,drnl) * CALL GNRICBESSH(ZSIGMA,X,LMAX,nl,drnl) * drj=zx*drjl(lm)+jl(lm) drhl=drnl(lm) c3 drhl=drjl(lm)+ci*drnl(lm) c3 drhl=zx*drhl+jl(lm)+ci*nl(lm) fintfe2=(exe*drj+fxe*drhl)* & dconjg(exe*drj+fxe*drhl) * return 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. 06/2004 Alexander Moroz