C (C) Copr. 06/2009 Alexander Moroz C Wave-scattering.com C http://www.wave-scattering.com/ C wavescattering@yahoo.com C C Last update: June 2009 C C This is CRMNT, a computer code to calculate normalized spontaneous C emission decay rates and the Ohmic loss contribution to the c nonradiative decay rates in the dipole-dipole interaction approximation C for various dipole polarizabilities C C The selection of polarizability employed in the calculation value C is determines by the value of integer variable imod: C C if imod.eq.0 ... Rayleigh polarizability C if imod.eq.1 ... Carminati et al polarizability C if imod.eq.2 ... T-asymptotic polarizability (e.g. Eq. (5) of [2] C in the dipolar case; the r.h.s. of Eq. (42) C of [3] in the general case). C if imod.eq.3 ... polarizability from the exact expression from the C dipolar T-matrix C if imod.eq.4 ... Meier&Wokaun polarizability C C The PROGRAM HAS BEEN WRITTEN ACCORDING TO THEORY DESCRIBED IN C C [1] R. Carminati, J.-J. Greffet, C. Henkel, J.M. Vigoureux, C Radiative and non-radiative decay of a single molecule close to C a metallic nanoparticle, Opt. Commun. 261 (2006) 368-375. C http://dx.doi.org/10.1016/j.optcom.2005.12.009 C C and in my article C C [2] 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 For the nonradiative rates, see also a recent follow up C C [3] A. Moroz, A superconvergent representation of the C Gersten-Nitzan and Ford-Weber nonradiative rates C J. Phys. Chem. C 115(40), 19546–19556 (2011). c http://dx.doi.org/10.1021/jp2057833 C C The program default options correspond to the cases considered in the above C article. C C GIVEN INTERIOR AND EXTERIOR PERMITTIVITIES CCEPS AND ZEPS0, C RESPECTIVELY, RADIUS RMUF AND FREQUENCY OMEGA, C THIS PROGRAM CALCULATES EITHER THE C SPONTANEOUS EMISSION DECAY RATES NORMALIZED EITHER TO THE RATES IN A VACCUM C OR IN THE VACUUM FILLED WITH THE LOCAL DIPOLE MEDIUM. C BOTH RADIATIVE AND NONRADIATIVE (DUE TO OHMIC LOSS) DECAY RATES ARE C CALCULATED. C C OUTPUT: C Returns files *.dat, *.dat containing the radiative decay C rates and files *.dat and *.dat'.dat, which contain the C relevant nonradiative decay rates and *.dat, which contains the C relevant total decay rate. C======================================================================== C C REQUIREMENTS & DEPENDENT ROUTINES: C C 1) In the case you want to run the code for real C material data you should download data file "Audat.dat" from C http://www.wave-scattering.com/Audat.dat.txt in the case of gold C or to supply your own material data. C C 2) In the case you want to run the "imod.eq.3" option, C 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 medium.f, sordalc.f, znsrefind.f are available as C http://www.wave-scattering.com/medium.f C http://www.wave-scattering.com/sordalc.f C http://www.wave-scattering.com/znsrefind.f C======================================================================== C TESTS: C C IN THE CASE OF A HOMOGENEOUS SPHERE THE RESULTS SHOULD COINCIDE WITH C THOSE OBTAINED BY C C R. Carminati, J.-J. Greffet, C. Henkel, and J.M. Vigoureux, C Radiative and non-radiative decay of a single molecule close to C a metallic nanoparticle, C Optics Commun. 261(2) (2006) 368-375. C http://dx.doi.org/10.1016/j.optcom.2005.12.009 C C In order to test your compilation, I have also provided an accompanying C limited MS Windows executable at C C http://www.wave-scattering.com/crmnt.exe C C which reads in the particle radius, dielectric constants, and C the wavelength and outputs the relevant rates. C C C (C) Copr. 6/2009 Alexander Moroz C * * * * * * C ============================== C ZEPS0 ... background permittivity C C-------------------------------------------------------------------- C >>> INPUT: <<< C 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======================================================================== 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 AUTHORS DISCLAIM ALL LIABILITY FOR C ANY DAMAGES THAT MAY RESULT 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 at yahoo.com C program crmnt C--------/---------/---------/---------/---------/---------/---------/-- implicit none integer NOUT,NX,NMAT,LCS,NFIN,NRS,NLD,nabs,iabs real*8 TOL,OMEGA,RINT,pi,rmuf complex*16 ZEPS0,CCEPS,cp logical ynbrug C--------/---------/---------/---------/---------/---------/---------/-- C INPUT REQUIRED: C-------------------------------------------------------------------- C PARAMETER INPUT : C-------------------------------------------------------------------- C ::: number of the output unit PARAMETER (NOUT=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=200) C ::: radial interval on which the decay rates are calculated in the units C of the sphere radius (RMUF=1) PARAMETER (RINT=20.d0) 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=1) C >>> BACKGROUND PERMITTIVITY <<< PARAMETER (ZEPS0=1.0D0) !2.281D0; 1.7689D0=1.33**2=H2O c sphere core dielectric constant PARAMETER (CCEPS=(- 2.03d0, 0.6d0)) !ld=354 for Carminati c PARAMETER (CCEPS=(- 15.04d0, 1.02d0)) !lambda=612 for Carminati c PARAMETER (CCEPS=(6.d0,0.d0)) !glycerol=2.1609d0 c Ag (612 nm) = (- 15.04d0, 1.02d0) c Ag (354 nm) = (- 2.03d0, 0.6d0) C >>> SPHERE (OUTER SHELL SCATTERER) PERMITTIVITY <<< * n(silica)=1.45 <---> EPS(1)=2.1025D0 * n(ZnS)=2. <---> EPS(1)=4.D0 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=0) * 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=73) c 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,ij1,il,ml,irx,lmx,ieps,l1,l2,ikl,irs,idt,ild, & imod REAL*8 xpar,xperp,xx,xxinv, & rsnm,rsnm0,lambda,lambda0, & gqmpar,gqmperp,gaspar,gasperp,gperp,gpar,gnrperp,gnrpar,gnrim, & gim,rmff,delo, & omxf,reepsz,plasma,omxp,ff,filfrac !ff for Bruggemann;filfrac for ZnS complex*16 ci,cone,czero,zdet1,zdet2,zeps1,Z1,Z2,zz,zz1,zk, & zpar,zperp,ztmas,zalph0,zalph COMPLEX*16 ZSIGMA,ZSIGME1,ZSIGME2 * Array declarations: REAL*8 RMF(1),rff(1),xim(1,nabs),xie(1,nabs) REAL*8 dos(nx+1),dosperp(nx+1),dospar(nx+1) COMPLEX*16 RX(2),SG(2),CQEPS(2),zeps(2) real*8 omf(NFIN) complex*16 ceps1(NFIN) * * coated sphere declarations: * moving lmax ===> * COMPLEX*16 cxm(1),dxm(1) COMPLEX*16 cxe(1),dxe(1) COMPLEX*16 exm,fxm,exe1,fxe1,exe2,fxe2 COMPLEX*16 cfel(1),cfml(1),cdrfel(1) COMPLEX*16 cfnrel(1),cfnrml(1),cdrfnrel(1) COMPLEX*16 cmx11(1),cmx12(1),cmx21(1),cmx22(1) COMPLEX*16 tt1(2,2,1,2,1), tt2(2,2,1,2,1) * C--------/---------/---------/---------/---------/---------/---------/-- COMPLEX*16 JL(0:1), NL(0:1) COMPLEX*16 DRJL(0:1), DRNL(0:1) COMPLEX*16 UL(2,0:1),DRUL(2,0:1) COMPLEX*16 WL(2,0:1), DRWL(2,0:1) COMPLEX*16 HL(0:1),DRHL(0:1) C ===================== 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/1/ * * security trap - remainder * if (nmat.gt.1) then write(6,*)'Real material data are to be provided' 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=5.d0 * rmf(1)=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 zeps(2)=zeps0 write(6,*)'Give the vacuum (detection) wavelength =' read (5,*) lambda cz lambda=354.d0 C--------/---------/---------/---------/---------/---------/---------/-- * The header of diprates*.dat output files: OPEN(UNIT=NOUT,FILE='qmrates.dat') rewind(NOUT) WRITE(NOUT,*)'#QM-Decay rates as a function of the distance & normalized to the decay rates in a host medium' * write(nout,*)'#Material number =',NMAT write(nout,*)'# Homogeneous sphere' * c WRITE(NOUT,*) write(nout,*)'#Sphere radius rmuf in relative units=', rmuf write(nout,*)'#Sphere radius in nm=', rsnm c WRITE(NOUT,*) * WRITE(NOUT,*)'#the vacuum (detection) wavelength =', lambda write(nout,*) '# rsnm/wavelength=', rsnm/lambda * write(nout,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout,*)'#sphere diel. const.=', zeps(1) 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,*)'#Angular-momentum cutoff LMX=', LMX WRITE(NOUT,*)'#In columns:' WRITE(NOUT,*)'#radial position in the units of sphere radius,' WRITE(NOUT,*)'#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='rrates.dat') rewind(NOUT+5) WRITE(NOUT+5,*)'#Radiative rates as a function of the distance & normalized to the decay rates in a free space filled in with & the host medium' WRITE(NOUT+5,*)'#the vacuum (detection) wavelength =', lambda write(nout+5,*) '# rsnm/wavelength=', rsnm/lambda write(nout+5,*)'#Sphere radius rmuf in relative units=', rmuf write(nout+5,*)'#Sphere radius in nm=', rsnm write(nout+5,*)'#Material number =',NMAT write(nout+5,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+5,*)'#sphere diel. const.=', zeps(1) write(nout+5,*)'#Angular-momentum cutoff LMX=', LMX WRITE(NOUT+5,*)'#In columns:' WRITE(NOUT+5,*)'#radial position in the units of sphere radius,' WRITE(NOUT+5,*)'#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='nrrates.dat') rewind(NOUT+6) write(nout+6,*)'#Nonradiative rates as a function of the distance & normalized to radiative decay rates in a host medium' write(nout+6,*) * 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 write(nout+6,*)'#Sphere radius rmuf in relative units=', rmuf write(nout+6,*)'#Sphere radius in nm=', rsnm c WRITE(NOUT+6,*) write(nout+6,*)'#Material number =',NMAT write(nout+6,*)'# Homogeneous sphere' * write(nout+6,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+6,*)'#sphere diel. const.=', zeps(1) 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,*)'#Angular-momentum cutoff LMX=', LMX WRITE(NOUT+6,*)'#In columns:' WRITE(NOUT+6,*)'#radial position in the units of sphere radius,' WRITE(NOUT+6,*)'#the 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='asrates.dat') rewind(NOUT+7) WRITE(NOUT,*)'#Asymptotic QM-Decay rates as a function of & the distance normalized to the decay rates in a 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 write(nout+7,*)'#Sphere radius rmuf in relative units=', rmuf write(nout+7,*)'#Sphere radius in nm=', rsnm c WRITE(NOUT+7,*) write(nout+7,*)'#Material number =',NMAT write(nout+7,*)'# Homogeneous sphere' write(nout+7,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+7,*)'#sphere diel. const.=', zeps(1) 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,*)'#Angular-momentum cutoff LMX=', LMX WRITE(NOUT+7,*)'#In columns:' WRITE(NOUT+7,*)'#radial position in the units of sphere radius,' WRITE(NOUT+7,*)'#the 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 zk=2.d0*pi*sqrt(zeps0)/lambda !wave vector in the host DO 300 IRS=0,NRS rsnm=rsnm0+dble(irs)*5 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 omega=2.d0*pi*rsnm/(lambda*rmuf) !size parameter for default !value of rmuf=1 * 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 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) go to 8 !dispersionless dielectric !or ideal metal * In case of a dispersion, EPSSPH is modified. * For ideal Drude metal * plasma=2.d0*pi*sphere radius in nm/(lambda_z in nm*rmuf) * where lambda_z is the wavelength for which Re eps_s=0. call medium(ynbrug,nmat,nfin,omega,lambda,rmuf, & rsnm,omf,ceps1,zeps0,zeps1) 7 zeps(1)=zeps1 if (lambda.eq.612) then zeps1=(- 15.04d0, 1.02d0) !lambda=612 for Carminati else if ((lambda.ne.354).and.(lambda.ne.612)) then write(6,*)'Read in complex dielectric constant' read(5,*) zeps1 end if *______________________________________ 8 continue 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,1,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)!! DO 28 j=1,1 CQEPS(1)=SQRT(ZEPS(j)) SG(1)=OMEGA*CQEPS(1) !in the sphere CQEPS(2)=SQRT(ZEPS(j+1)) SG(2)=OMEGA*CQEPS(2) !in the host * C Cc Chew test: x=ka=5; n=1.3, 1.5 cc zk=cone cc zeps(1)=1.3d0 * DO 25 I=1,2 * RX(I)=SG(I)*RMF(j) !corresponding size param. * CALL GNZBESS(RX(I),1,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 (see Eqs. (17,25)) C TT2 is the transfer matrix which translates the regular solution c from inside to outside the sphere (see Eqs. (22,26)) * do 27 l=1,lmx * * magnetic part * C--------/---------/---------/---------/---------/---------/---------/-- * Eq. (17): * 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)) * Eq. (22): * 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--------/---------/---------/---------/---------/---------/---------/-- * Eq. (25): * 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)) * Eq. (26): * 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 determined 40 CONTINUE * * Initialize total, radiative, and nonradiative decay rates: gqmperp= 0.d0 gqmpar= 0.d0 gasperp=0.d0 gaspar=0.d0 gperp= 0.d0 gpar= 0.d0 gnrperp=0.d0 gnrpar=0.d0 * * Assign polarizability: *==== 33 write(6,*)'imod=0 - Rayleigh polarizability' write(6,*)'imod=1 - Carminati et al polarizability' write(6,*)'imod=2 - T-matrix asymptotic polarizability' write(6,*)'imod=3 - T-matrix' write(6,*)'imod=4 - Meier&Wokaun polarizability' write(6,*)'Read in the value of imod' read(5,*) imod if ((imod.lt.0).or.(imod.gt.4)) then write(6,*)'The value of imod has to be between 0 and 4' write(6,*)'Read in imod again' go to 33 end if if (imod.eq.0) then write(6,*)'Results for Rayleigh polarizability' write(nout,*)'#Rayleigh polarizability' write(nout+5,*)'#Rayleigh polarizability' write(nout+6,*)'#Rayleigh polarizability' write(nout+7,*)'#Rayleigh polarizability' * Rayleigh option zalph=4.d0*pi*rsnm**3*(zeps(1)-cone)/(zeps(1)+2.d0*cone) *==== else if (imod.eq.1) then write(6,*)'Results for Carminati et al polarizability' write(nout,*)'#Carminati et al polarizability' write(nout+5,*)'#Carminati et al polarizability' write(nout+6,*)'#Carminati et al polarizability' write(nout+7,*)'#Carminati et al polarizability' * Carminati et al option without a depolarization term zalph0=4.d0*pi*rsnm**3*(zeps(1)-cone)/(zeps(1)+2.d0*cone) zalph=zalph0/(cone-ci*zk**3*zalph0/(6.d0*pi)) *==== else if (imod.eq.2) then write(6,*)'Results for T-matrix asymptotic polarizability' write(nout,*)'#Exact asymptotic polarizability' write(nout+5,*)'#Exact asymptotic polarizability' write(nout+6,*)'#Exact asymptotic polarizability' write(nout+7,*)'#Exact asymptotic polarizability' * Exact asymptotic with a depolarization term ztmas=ci*2.d0*rx(2)**3*(zeps(1)-cone)* & (cone-(zeps(1)+cone)*rx(2)**2/10.d0)/3.d0 ztmas=ztmas/(zeps(1)+2.d0*cone & -(zeps(1)-cone)*(zeps(1)+10.d0*cone)*rx(2)**2/10.d0 & - ci*2.d0*(zeps(1)-cone)*rx(2)**3/3.d0) zalph=-ci*6.d0*pi*ztmas/zk**3 *==== else if (imod.eq.3) then write(6,*)'Results for polarizability directly from the T-matrix' write(nout,*)'#Polarizability directly from the T-matrix' write(nout+5,*)'#Polarizability directly from the T-matrix' write(nout+6,*)'#Polarizability directly from the T-matrix' write(nout+7,*)'#Polarizability directly from the T-matrix' * Exact expression from the T-matrix ztmas=tt2(2,1,1,2,1)/tt2(1,1,1,2,1) zalph=-ci*6.d0*pi*ztmas/zk**3 else if (imod.eq.4) then write(6,*)'Results for the Meier&Wokaun polarizability' write(nout,*)'#Polarizability directly from the Meier&Wokaun' write(nout+5,*)'#Polarizability directly from the Meier&Wokaun' write(nout+6,*)'#Polarizability directly from the Meier&Wokaun' write(nout+7,*)'#Polarizability directly from the Meier&Wokaun' ztmas=(zeps(1)-cone) ztmas=ztmas/(zeps(1)+2.d0*cone -ztmas*rx(2)**2 & - ci*2.d0*ztmas*rx(2)**3/3.d0) zalph=4.d0*pi*rsnm**3*ztmas end if write(nout,*) write(nout+5,*) write(nout+6,*) write(nout+7,*) *>>> scanning the interval RINT: do 200 il =1,nx !10 !nx !over the interval RINT delo= rint/dble(nx) rmff= 1.d0 + dble(il)*delo cu delo= (rint-1.d0)/dble(nx) ct delo=0.5d0 ct rmff=(rsnm+ dble(il)*delo)/rsnm * identifying the layer for a given rmff do 45 ml=1,1 if (rmff.lt.rmf(ml)) then j=ml go to 50 end if 45 continue 50 if (rmff.gt.rmf(1)) j=2 CQEPS(1)=SQRT(ZEPS(j)) SG(1)=OMEGA*CQEPS(1) !in the dipole medium RX(1)=SG(1)*rmff !size param in the host *%%%%% "Exact" calculations: xperp= 3.d0*zk**3/(2.d0*pi) xpar= 3.d0*zk**3/(8.d0*pi) xx= sg(2)*rmff !"size-like param" in the host xxinv=1.d0/xx zperp=exp(2.d0*ci*xx)*(-xxinv**4 - 2.d0*ci*xxinv**5 +xxinv**6) zpar=exp(2.d0*ci*xx)*(xxinv**2 + 2.d0*ci*xxinv**3 & -3.d0*xxinv**4 - 2.d0*ci*xxinv**5 + xxinv**6) gqmperp= 1.d0 + xperp*imag(zalph*zperp) gqmpar = 1.d0 + xpar*imag(zalph*zpar) *%%%%% Short-distance limit: zperp=2.d0*ci*xxinv**3/(3.d0) + xxinv**4 +xxinv**6 zpar= -4.d0*ci*xxinv**3/(3.d0) -xxinv**4 + xxinv**6 gasperp= 1.d0 + xperp*imag(zalph*zperp) gaspar= 1.d0 + xpar*imag(zalph*zpar) *%%%%% Nonradiative decay rates: zperp= xxinv**4 + xxinv**6 zpar= xxinv**2 - xxinv**4 + xxinv**6 gnrperp= & xperp*zperp*(imag(zalph)-zk**3*abs(zalph)**2/(6.d0*pi)) gnrpar= & xpar*zpar*(imag(zalph)-zk**3*abs(zalph)**2/(6.d0*pi)) *%%%%% Radiative decay rates: xperp= zk**6/(4.d0*pi**2) xpar= zk**6/(16.d0*pi**2) zperp= xxinv**4 + xxinv**6 zpar= - xxinv**4 + xxinv**6 gperp= 1.d0+ xperp*zperp*abs(zalph)**2 & + dble(zalph)*xxinv**3*zk**3/pi gpar= 1.d0+ xpar*zpar*abs(zalph)**2 & - dble(zalph)*xxinv**3*zk**3/(2.d0*pi) C--------/---------/---------/---------/---------/---------/---------/-- * OUTPUT WRITING: * * QM rates: gim=gqmperp+2.d0*gqmpar !See Eq. (7) of PRA34_3410 dospar(il+1) = gqmpar !tangential oscillations dosperp(il+1)= gqmperp !radial oscillations dos(il+1) = gim/3.d0 !polarizability average: there is !one radial and two indep. tangential !orientations possible WRITE(NOUT,999) rmff, & dospar(il+1), & dosperp(il+1), & dos(il+1) * Radiative rates: gim=gperp+2.d0*gpar dospar(il+1) = gpar dosperp(il+1)= gperp dos(il+1) = 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 rates: gnrim=gnrperp+2.d0*gnrpar dospar(il+1) = gnrpar !tangential oscillations dosperp(il+1)= gnrperp !radial oscillations dos(il+1) = gnrim/3.d0 !polarizability 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) * Asymptotic QM rates: gnrim=gasperp+2.d0*gaspar dospar(il+1) = gaspar dosperp(il+1)= gasperp dos(il+1) = 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,*)'QM rates in qmrates.dat' write(6,*)'Nonradiative rates in nrrates.dat' write(6,*)'Radiative rates in rrates.dat' write(6,*)'Asymptotic QM rates in asrates.dat' 999 format(f9.6,3X,3(e15.8,3x)) c 999 format(f8.6,3X,e23.16,3x,2(e15.8,3x)) END * C (C) Copr. 08/2008 Alexander Moroz