C (C) Copr. 06/2011 Alexander Moroz C Wave-scattering.com C http://www.wave-scattering.com/ C wavescattering@yahoo.com C C Last update: June 2011 C C This is GN, a computer code to calculate normalized C nonradiative decay rates based on the theory of C C [1] J. Gersten and A. Nitzan, C Spectroscopic properties of molecules interacting with C small dielectric particles, C J. Chem. Phys. 75, 1139-1152 (1981). C http://link.aip.org/link/?JCPSA6/75/1139/1 C C The GN code code has been used to generate the results in my articles C C [2] A. Moroz, Non-radiative decay of a dipole emitter close to C a metallic nanoparticle: Importance of higher-order multipole C contributions, C Opt. Commun. 283(10), 2277-2287 (2010). C http://dx.doi.org/10.1016/j.optcom.2010.01.061 C C and 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 GN theory has been supplemented with two generalizations of the C corresponding l-polar polarizabilities. The respective options could C be selected by chosing an appropriate value of the variable imod. C C if (imod.eq.1) then the standard electrostatic polarizabilities in the C Rayleigh limit are used (e.g. Eq. (3) of [2] in the dipolar C case; the l.h.s. of Eq. (42) of [3] in the general case). C C if (imod.eq.2) then the T-matrix asymptotic up to the order x^3 in the C size parameter has been used for the l-polar polarizabilities C (e.g. Eq. (5) of [2] in the dipolar case; the r.h.s. of Eq. (42) C of [3] in the general case). C if (imod.eq.3) then the so-called T-polarization is used, C wherein the polarizabilities are directly related to the C corresponding T-matrix (Mie coefficients) elements C (e.g. Eq. (4) of [2] in the dipolar case). C C GN calculates these quantities for two principal (tangential and radial) C and averaged dipole orientations. C C In a homogeneous medium the radiative rate is proportional to the medium C refractive index n: C C gamma_r = n*gamma_0r C C where gamma_0r is the radiative rate in vacuum. The code then outputs the ratios C of a given (either radiative or nonradiative) rate to gamma_0r. C Hence for large sphere-dipole seprations the radiative rate ratios will C approach the refractive index n of the host medium. C C The following output files are generated: C C Rayleighgnnr.dat - contains the GN non-radiative rates as a function of the C distance normalized to the decay rates in a host medium for imod=1. C C Tasgnnr.dat - the GN non-radiative rates as a function of the C distance normalized to the decay rates in a host medium for imod=2. C C Tmatgnnr.dat - contains the GN non-radiative rates as a function of the C distance normalized to the decay rates in a host medium for imod=3. C C Ratiosgn.dat - contains the Ratios of non-radiative rates obtained C for T-matrix and static Rayleigh polarization as a function C of the distance. C C Rlr.dat - contains the GN radiative rates as a function of the C distance normalized to the decay rates in a host medium for imod=1. C C Tasr.dat - contains the GN radiative rates as a function of the C distance normalized to the decay rates in a host medium for imod=2. C C Tmatr.dat - contains the GN radiative rates as a function of the C distance normalized to the decay rates in a host medium for imod=3. C C Cnr.dat - contains a dual expression for the C GN non-radiative rates as a function of the distance C C Casmpnr.dat - contains the 5-term analytic expression for the C GN non-radiative rates as a function of the distance C C Cnr4o.dat - contains the 4-term analytic expression for the C GN non-radiative rates as a function of the distance C C Cnr3o.dat - contains the 3-term analytic expression for the C GN non-radiative rates as a function of the distance C C Asmpnrd3.dat - contains the dominant d**(-3) short-distance C GN non-radiative rates as a function of the distance C C Asmp2onr.dat - contains the 2-term short-distance asymptotic of C GN non-radiative rates as a function of the distance C C Asmp3onr.dat - contains the 3-term short-distance asymptotic of C GN non-radiative rates as a function of the distance C C Asmpnr.dat - contains the short-distance asymptotic of GN non-radiative C rates as a function of the distance normalized to C the decay rates in a host medium (up to dilogarithm) C C The data contain in columns the dipole position, and the C corresponding rates for tangentially and radially oriented dipole, and their C orientational average, respectively, normalized to the radiative rates in vacuum. C C Coefasmp.dat - contains a comparison of the left and right hand sides C of the asymptotic expansions according to Eqs. 53-55 of [3] c C rmff ... the dipole position in the units of sphere radius C xd=d/a ... the dipole distance from the sphere surface in the C units of sphere radius C The above data files could have either rmff or xd in C the first column C C In order to read an external file and compare it with just calclated values, C find the code line introduced with the text: "ct for convergence test:". If C the program line containing "dos(il+1)=dos(il+1)/xr(il)" is uncommented, C the averaged nonradiative rates are overwritten by the ratios to the C external data file values. C Only if the line containing "dos(il+1)=dos(il+1)/xr(il)" has been C commented out, the absolute average rates are outputted. C Default is the the absolute average rates on the output. C ================== C C IMPORTANT!!! C C Before comparing the simulations with your experiment, read Sec. 12.1 of C C [4] 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--------/---------/---------/---------/---------/---------/---------/- C C RUNNING THE CODE C C The code uses similar set of parameters and variables as my other codes C available from http://www.wave-scattering.com/codes.html C C C In the case one performs calculations for a given wavelength, one should C chose the NMAT=0 option and specify the relevant dielectric functions C manually as parameter values (default option). C C The dependent routines readmat.f and medium.f here are only to supply C material dielectric function. The routines are bypassed if the NMAT=0 C option is selected (default option). C C The number of multipoles taken into account in the calculations is C controlled by the of variable LMX. This may not be larger than the C parameter value LMAXD, which defines corresponding array dimensions. C If LMAXD.gt.60 one has to comment out the routines gnzbess.f and gnricbessh.f C calling the spherical Besel functions. The latter are only needed if C and only if the option imod.eg.3 is selected/required. The AMOS package of C complex Bessel functions routines you can, together with their dependencies, C freely download at http://www.netlib.org/amos/. C 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 In the case the imod.eq.3 option is not needed (default option), as in C the case of Ref. [3], there is no upper limit on the values of C LMAXD and LMX. In order to bypass any call to Bessel function routines, C I have inserted two program jumps "go to 25" and "go to 40". Obviously, C the 2 program jumps have to be commented out if the imod.eq.3 option is C to be used again. C C--------/---------/---------/---------/---------/---------/---------/- C TESTS C C A number of selected output files is available from C http://www.wave-scattering.com/gn.html C C--------/---------/---------/---------/---------/---------/---------/- C EXTENSIONS c C By selecting nonzero values of NLD and NRS, the respective loops C over the wavelengths (DO 400 ILD=0,NLD) and over the sphere radii C over the sphere radii (DO 300 IRS=0,NRS) could be activated. 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 gn C--------/---------/---------/---------/---------/---------/---------/-- C Gersten&Nitzan nonradiative rates C--------/---------/---------/---------/---------/---------/---------/-- implicit none integer LMAXD,NOUT,NX,NMAT,NFIN,NRS,NLD,NMD real*8 TOL,OMEGA,RINT,pi,rmuf complex*16 ZEPS0,CCEPS,cp logical ynbrug C-------------------------------------------------------------------- C PARAMETER INPUT : C-------------------------------------------------------------------- c Maximal number of spherical harmonics used. The floating number is c specified below by the value of LMAX parameter PARAMETER (LMAXD=500) * 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 ::: The number of different model situations for polarization PARAMETER (NMD=3) C ::: relative error allowed. PARAMETER (TOL=1.d-3) c: frequency c PARAMETER (omega=0.2D0) C ::: number of steps on direct space: PARAMETER (NX=400) C ::: radial interval on which the decay rates are calculated in the units C of the sphere radius (RMUF=1) PARAMETER (RINT=8.d0) C >>> BACKGROUND PERMITTIVITY <<< PARAMETER (ZEPS0=1.0D0) !2.281D0; 1.7689D0=1.33**2=H2O c sphere core dielectric constant PARAMETER (CCEPS=(- 1.0d0, 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 xc,xpar,xperp,xx,xxinv,xxf,xrperp,xrpar,xq,xkoef,xrint, & rsnm,rsnm0,lambda,lambda0,gamma0,xd, & gim,rmff,delo,xu,xf1,xf2,xf3,xfs,x1,x2,x3, & omxf,reepsz,plasma,omxp,ff,filfrac !ff for Bruggemann;filfrac for ZnS real*8 dtperp(nmd),dtpar(nmd),gnrperp(nmd),gnrpar(nmd), & gperp(nmd),gpar(nmd),ximal(lmaxd,nmd),xfac(lmaxd+1), & gnrperps,gnrpars,rff1,rff2,xli2 complex*16 ci,cone,czero,zdet1,zdet2,zeps1,Z1,Z2,zz,zz1,zk, & zp,zx,zdpar,zdperp,zpar,zperp,ztmas,zalph0,zalph,zcrmn, & zff,zv complex*16 zal(lmaxd,nmd),za(4),zb(4) COMPLEX*16 ZSIGMA,ZSIGME1,ZSIGME2 * Array declarations: REAL*8 RMF(1),rff(1) REAL*8 dos(nx+1),dosperp(nx+1),dospar(nx+1),xr(nx+1) COMPLEX*16 RX(2),SG(2),CQEPS(2),zeps(2) real*8 omf(NFIN) complex*16 ceps1(NFIN) * moving lmax ===> * COMPLEX*16 cxm(lmaxd),dxm(lmaxd) COMPLEX*16 cxe(lmaxd),dxe(lmaxd) COMPLEX*16 exm,fxm,exe1,fxe1,exe2,fxe2 COMPLEX*16 cfel(lmaxd),cfml(lmaxd),cdrfel(lmaxd) COMPLEX*16 cfnrel(1),cfnrml(1),cdrfnrel(1) COMPLEX*16 cmx11(1),cmx12(1),cmx21(1),cmx22(1) COMPLEX*16 tt1(2,2,lmaxd,2,1), tt2(2,2,lmaxd,2,1) * C--------/---------/---------/---------/---------/---------/---------/-- COMPLEX*16 JL(0:lmaxd), NL(0:lmaxd) COMPLEX*16 DRJL(0:lmaxd), DRNL(0:lmaxd) COMPLEX*16 UL(2,0:lmaxd),DRUL(2,0:lmaxd) COMPLEX*16 WL(2,0:lmaxd), DRWL(2,0:lmaxd) COMPLEX*16 HL(0:lmaxd),DRHL(0:lmaxd) C-------------------------------------------------------------------- C READING THE DATA : C DATA xc/2.99792458d17/ !speed of light in nm/s 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/500/ * * 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 ct write(6,*)'Give the intrinsic radiative rate of dye in 1/s =' ct read (5,*) xrint C According to Soller et al, Nano Lett. {\bf 7}(7), 1941-1946, 2007: C Ru(pby)${}_3$-Bi has the intrinsic R_{rad}=7.1d4 (1/s) C PdCS3P-Bi has the intrinsic R_{rad}=62 (1/s) xrint=7.1d4 *>>> xfac(1)=1.d0 !(2l-1)!! if (lmx.ge.2) then do il=2,lmx xfac(il)=(2*il-1)*xfac(il-1) enddo end if xfac(lmx+1)=(2*lmx+1)*xfac(lmx) C--------/---------/---------/---------/---------/---------/---------/-- * The header of gn*.dat output files: OPEN(UNIT=NOUT,FILE='Rayleighgnnr.dat') rewind(NOUT) WRITE(NOUT,*)'#GN non-radiative rates as a function of the & distance normalized to the decay rates in a host medium' * write(nout,*)'#Static Rayleigh polarization' * write(nout,*)'#Material number =',NMAT write(nout,*)'# Homogeneous sphere' * 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,*)'#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 write(nout,*)'#Angular-momentum cutoff LMX=', LMX WRITE(NOUT,*)'#In columns:' WRITE(NOUT,*)'#position, the dipole nonradiative rates' write(nout,*)'#for tangentially and radially oriented dipole,' write(nout,*)'#and their orientational average, respectively' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+1,FILE='Tasgnnr.dat') rewind(NOUT+1) WRITE(NOUT+1,*)'#GN non-radiative rates as a function of the & distance normalized to the decay rates in a host medium' * write(nout+1,*)'#T-asymptotic polarization' * write(nout+1,*)'#Material number =',NMAT write(nout+1,*)'# Homogeneous sphere' * c WRITE(NOUT+1,*) write(nout+1,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+1,*) * WRITE(NOUT+1,*)'#the vacuum (detection) wavelength =', lambda write(nout+1,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+1,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+1,*)'#sphere diel. const.=', zeps(1) C--------/---------/---------/---------/---------/---------/---------/-- write(nout+1,*)'#Length of the searched interval rint=', rint write(nout+1,*)'#Number of the points on the uniform mesh nx=',nx write(nout+1,*)'#Angular-momentum cutoff LMX=', LMX WRITE(NOUT+1,*)'#In columns:' WRITE(NOUT+1,*)'#position, the dipole nonradiative rates' write(nout+1,*)'#for tangentially and radially oriented dipole,' write(nout+1,*)'#and their orientational average, respectively' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+2,FILE='Tmatgnnr.dat') rewind(NOUT+2) WRITE(NOUT+2,*)'#GN non-radiative rates as a function of the & distance normalized to the decay rates in a host medium' * write(nout+2,*)'#T-matrix polarization' * write(nout+2,*)'#Material number =',NMAT write(nout+2,*)'# Homogeneous sphere' * c WRITE(NOUT+2,*) write(nout+2,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+2,*) * WRITE(NOUT+2,*)'#the vacuum (detection) wavelength =', lambda write(nout+2,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+2,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+2,*)'#sphere diel. const.=', zeps(1) C--------/---------/---------/---------/---------/---------/---------/-- write(nout+2,*)'#Length of the searched interval rint=', rint write(nout+2,*)'#Number of the points on the uniform mesh nx=',nx write(nout+2,*)'#Angular-momentum cutoff LMX=', LMX WRITE(NOUT+2,*)'#In columns:' WRITE(NOUT+2,*)'#position, the dipole nonradiative rates' write(nout+2,*)'#for tangentially and radially oriented dipole,' write(nout+2,*)'#and their orientational average, respectively' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+3,FILE='Ratiosgn.dat') rewind(NOUT+3) WRITE(NOUT+3,*)'#Ratios of non-radiative rates obtained & for T-matrix and static Rayleigh polarization as a function & of the distance' * write(nout+3,*)'#Material number =',NMAT write(nout+3,*)'# Homogeneous sphere' * c WRITE(NOUT+3,*) write(nout+3,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+3,*) * WRITE(NOUT+3,*)'#the vacuum (detection) wavelength =', lambda write(nout+3,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+3,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+3,*)'#sphere diel. const.=', zeps(1) write(nout+3,*)'#Length of the searched interval rint=', rint write(nout+3,*)'#Number of the points on the uniform mesh nx=',nx write(nout+3,*)'#Angular-momentum cutoff LMX=', LMX WRITE(NOUT+3,*)'#In columns:' WRITE(NOUT+3,*)'#position, the dipole nonradiative rates ratios' write(nout+3,*)'#for tangentially and radially oriented dipole:' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+4,FILE='Rlr.dat') rewind(NOUT+4) WRITE(NOUT+4,*)'#GN radiative rates as a function of the & distance normalized to the decay rates in a host medium' * write(nout+4,*)'#Static Rayleigh polarization' * write(nout+4,*)'#Material number =',NMAT write(nout+4,*)'# Homogeneous sphere' * c WRITE(NOUT+4,*) write(nout+4,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+4,*) * WRITE(NOUT+4,*)'#the vacuum (detection) wavelength =', lambda write(nout+4,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+4,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+4,*)'#sphere diel. const.=', zeps(1) write(nout+4,*)'#Length of the searched interval rint=', rint write(nout+4,*)'#Number of the points on the uniform mesh nx=',nx write(nout+4,*)'#Angular-momentum cutoff LMX=', LMX WRITE(NOUT+4,*)'#In columns:' WRITE(NOUT+4,*)'#position, the dipole radiated power' write(nout+4,*)'#for tangentially and radially oriented dipole,' write(nout+4,*)'#and their orientational average, respectively' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+5,FILE='Tasr.dat') rewind(NOUT+5) WRITE(NOUT+5,*)'#GN radiative rates as a function of the & distance normalized to the decay rates in a host medium' * write(nout+5,*)'#T-asymptotic polarization' * write(nout+5,*)'#Material number =',NMAT write(nout+5,*)'# Homogeneous sphere' * c WRITE(NOUT+5,*) write(nout+5,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+5,*) * WRITE(NOUT+5,*)'#the vacuum (detection) wavelength =', lambda write(nout+5,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+5,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+5,*)'#sphere diel. const.=', zeps(1) write(nout+5,*)'#Length of the searched interval rint=', rint write(nout+5,*)'#Number of the points on the uniform mesh nx=',nx write(nout+5,*)'#Angular-momentum cutoff LMX=', LMX 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' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+6,FILE='Tmatr.dat') rewind(NOUT+6) WRITE(NOUT+6,*)'#GN radiative rates as a function of the & distance normalized to the decay rates in a host medium' * write(nout+6,*)'#T-matrix polarization' * write(nout+6,*)'#Material number =',NMAT write(nout+6,*)'# Homogeneous sphere' * c WRITE(NOUT+6,*) write(nout+6,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+6,*) * WRITE(NOUT+6,*)'#the vacuum (detection) wavelength =', lambda write(nout+6,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+6,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+6,*)'#sphere diel. const.=', zeps(1) write(nout+6,*)'#Length of the searched interval rint=', rint write(nout+6,*)'#Number of the points on the uniform mesh nx=',nx write(nout+6,*)'#Angular-momentum cutoff LMX=', LMX WRITE(NOUT+6,*)'#In columns:' WRITE(NOUT+6,*)'#position, the dipole radiated power' write(nout+6,*)'#for tangentially and radially oriented dipole,' write(nout+6,*)'#and their orientational average, respectively' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+7,FILE='Cnr.dat') rewind(NOUT+7) WRITE(NOUT+7,*)'#Dual expression for the & GN non-radiative rates as a function of the distance as & a function of the distance normalized to the decay rates in a & host medium' write(nout+7,*)'#Material number =',NMAT write(nout+7,*)'# Homogeneous sphere' * 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,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+7,*)'#sphere diel. const.=', zeps(1) 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, the dipole nonradiative rates' write(nout+7,*)'#for tangentially and radially oriented dipole,' write(nout+7,*)'#and their orientational average, respectively' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+8,FILE='Casmpnr.dat') rewind(NOUT+8) WRITE(NOUT+8,*)'#The 5-term analytic expression for the & GN non-radiative rates as a function of the distance as & a function of the distance normalized to the decay rates in a & host medium' write(nout+8,*)'#Material number =',NMAT write(nout+8,*)'# Homogeneous sphere' * c WRITE(NOUT+8,*) write(nout+8,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+8,*) * WRITE(NOUT+8,*)'#the vacuum (detection) wavelength =', lambda write(nout+8,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+8,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+8,*)'#sphere diel. const.=', zeps(1) write(nout+8,*)'#Length of the searched interval rint=', rint write(nout+8,*)'#Number of the points on the uniform mesh nx=',nx WRITE(NOUT+8,*)'#In columns:' WRITE(NOUT+8,*)'#position, the dipole nonradiative rates' write(nout+8,*)'#for tangentially and radially oriented dipole,' write(nout+8,*)'#and their orientational average, respectively' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+9,FILE='Asmpnrd3.dat') rewind(NOUT+9) WRITE(NOUT+9,*)'#The dominant d**(-3) short-distance asymptotic & of GN non-radiative rates as a function of the distance & normalized to the decay rates in a host medium' write(nout+9,*)'#Material number =',NMAT write(nout+9,*)'# Homogeneous sphere' * c WRITE(NOUT+9,*) write(nout+9,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+9,*) * WRITE(NOUT+9,*)'#the vacuum (detection) wavelength =', lambda write(nout+9,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+9,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+9,*)'#sphere diel. const.=', zeps(1) write(nout+9,*)'#Length of the searched interval rint=', rint write(nout+9,*)'#Number of the points on the uniform mesh nx=',nx WRITE(NOUT+9,*)'#In columns:' WRITE(NOUT+9,*)'#position, the dipole nonradiative rates' write(nout+9,*)'#for tangentially and radially oriented dipole,' write(nout+9,*)'#and their orientational average, respectively' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+10,FILE='Asmp2onr.dat') rewind(NOUT+10) WRITE(NOUT+10,*)'#The two-order short-distance asymptotic & of GN non-radiative rates as a function of the distance & normalized to the decay rates in a host medium' write(nout+10,*)'#Material number =',NMAT write(nout+10,*)'# Homogeneous sphere' * c WRITE(NOUT+10,*) write(nout+10,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+10,*) * WRITE(NOUT+10,*)'#the vacuum (detection) wavelength =', lambda write(nout+10,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+10,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+10,*)'#sphere diel. const.=', zeps(1) write(nout+10,*)'#Length of the searched interval rint=', rint write(nout+10,*)'#Number of the points on the uniform mesh nx=',nx WRITE(NOUT+10,*)'#In columns:' WRITE(NOUT+10,*)'#position, the dipole nonradiative rates' write(nout+10,*)'#for tangentially and radially oriented dipole,' write(nout+10,*)'#and their orientational average, respectively' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+11,FILE='Asmp3onr.dat') rewind(NOUT+11) WRITE(NOUT+11,*)'#The three-order short-distance asymptotic & of GN non-radiative rates as a function of the distance & normalized to the decay rates in a host medium' write(nout+11,*)'#Material number =',NMAT write(nout+11,*)'# Homogeneous sphere' * c WRITE(NOUT+11,*) write(nout+11,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+11,*) * WRITE(NOUT+11,*)'#the vacuum (detection) wavelength =', lambda write(nout+11,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+11,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+11,*)'#sphere diel. const.=', zeps(1) write(nout+11,*)'#Length of the searched interval rint=', rint write(nout+11,*)'#Number of the points on the uniform mesh nx=',nx WRITE(NOUT+11,*)'#In columns:' WRITE(NOUT+11,*)'#position, the dipole nonradiative rates' write(nout+11,*)'#for tangentially and radially oriented dipole,' write(nout+11,*)'#and their orientational average, respectively' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+12,FILE='Asmpnr.dat') rewind(NOUT+12) WRITE(NOUT+12,*)'#Short-distance asymptotic of GN non-radiative & rates as a function of the distance normalized to & the decay rates in a host medium (up to F(u,v))' write(nout+12,*)'#Material number =',NMAT write(nout+12,*)'# Homogeneous sphere' * c WRITE(NOUT+12,*) write(nout+12,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+12,*) * WRITE(NOUT+12,*)'#the vacuum (detection) wavelength =', lambda write(nout+12,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+12,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+12,*)'#sphere diel. const.=', zeps(1) write(nout+12,*)'#Length of the searched interval rint=', rint write(nout+12,*)'#Number of the points on the uniform mesh nx=',nx WRITE(NOUT+12,*)'#In columns:' WRITE(NOUT+12,*)'#position, the dipole nonradiative rates' write(nout+12,*)'#for tangentially and radially oriented dipole,' write(nout+12,*)'#and their orientational average, respectively' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+13,FILE='Cnr4o.dat') rewind(NOUT+13) WRITE(NOUT+13,*)'#The 4-term analytic expression for the & GN non-radiative rates as a function of the distance as & a function of the distance normalized to the decay rates in a & host medium' write(nout+13,*)'#Material number =',NMAT write(nout+13,*)'# Homogeneous sphere' * c WRITE(NOUT+13,*) write(nout+13,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+13,*) * WRITE(NOUT+13,*)'#the vacuum (detection) wavelength =', lambda write(nout+13,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+13,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+13,*)'#sphere diel. const.=', zeps(1) write(nout+13,*)'#Length of the searched interval rint=', rint write(nout+13,*)'#Number of the points on the uniform mesh nx=',nx WRITE(NOUT+13,*)'#In columns:' WRITE(NOUT+13,*)'#position, the dipole nonradiative rates' write(nout+13,*)'#for tangentially and radially oriented dipole,' write(nout+13,*)'#and their orientational average, respectively' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+14,FILE='Cnr3o.dat') rewind(NOUT+14) WRITE(NOUT+14,*)'#The 3-term analytic expression for the & GN non-radiative rates as a function of the distance as & a function of the distance normalized to the decay rates in a & host medium' write(nout+14,*)'#Material number =',NMAT write(nout+14,*)'# Homogeneous sphere' * c WRITE(NOUT+14,*) write(nout+14,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+14,*) * WRITE(NOUT+14,*)'#the vacuum (detection) wavelength =', lambda write(nout+14,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+14,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+14,*)'#sphere diel. const.=', zeps(1) write(nout+14,*)'#Length of the searched interval rint=', rint write(nout+14,*)'#Number of the points on the uniform mesh nx=',nx WRITE(NOUT+14,*)'#In columns:' WRITE(NOUT+14,*)'#position, the dipole nonradiative rates' write(nout+14,*)'#for tangentially and radially oriented dipole,' write(nout+14,*)'#and their orientational average, respectively' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+15,FILE='Cnr2o.dat') rewind(NOUT+15) WRITE(NOUT+15,*)'#The 3-term analytic expression for the & GN non-radiative rates as a function of the distance as & a function of the distance normalized to the decay rates in a & host medium' write(nout+15,*)'#Material number =',NMAT write(nout+15,*)'# Homogeneous sphere' * c WRITE(NOUT+15,*) write(nout+15,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+15,*) * WRITE(NOUT+15,*)'#the vacuum (detection) wavelength =', lambda write(nout+15,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+15,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+15,*)'#sphere diel. const.=', zeps(1) write(nout+15,*)'#Length of the searched interval rint=', rint write(nout+15,*)'#Number of the points on the uniform mesh nx=',nx WRITE(NOUT+15,*)'#In columns:' WRITE(NOUT+15,*)'#position, the dipole nonradiative rates' write(nout+15,*)'#for tangentially and radially oriented dipole,' write(nout+15,*)'#and their orientational average, respectively' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+16,FILE='Cnrdo.dat') rewind(NOUT+16) WRITE(NOUT+16,*)'#The 1-term analytic expression for the & GN non-radiative rates as a function of the distance as & a function of the distance normalized to the decay rates in a & host medium' write(nout+16,*)'#Material number =',NMAT write(nout+16,*)'# Homogeneous sphere' * c WRITE(NOUT+16,*) write(nout+16,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+16,*) * WRITE(NOUT+16,*)'#the vacuum (detection) wavelength =', lambda write(nout+16,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+16,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+16,*)'#sphere diel. const.=', zeps(1) write(nout+16,*)'#Length of the searched interval rint=', rint write(nout+16,*)'#Number of the points on the uniform mesh nx=',nx WRITE(NOUT+16,*)'#In columns:' WRITE(NOUT+16,*)'#position, the dipole nonradiative rates' write(nout+16,*)'#for tangentially and radially oriented dipole,' write(nout+16,*)'#and their orientational average, respectively' C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+20,FILE='Coefasmp.dat') rewind(NOUT+20) WRITE(NOUT+20,*)'#The 1-term analytic expression for the & GN non-radiative rates as a function of the distance as & a function of the distance normalized to the decay rates in a & host medium' write(nout+20,*)'#Material number =',NMAT write(nout+20,*)'# Homogeneous sphere' * c WRITE(NOUT+20,*) write(nout+20,*)'#Sphere radius rmuf in relative units=', rmuf c WRITE(NOUT+20,*) * WRITE(NOUT+20,*)'#the vacuum (detection) wavelength =', lambda write(nout+20,*) '# rsnm/wavelength=', rsnm/lambda * write(nout+20,*)'#host dielectric constant=', zeps(2) if (nmat.eq.0) & write(nout+20,*)'#sphere diel. const.=', zeps(1) write(nout+20,*)'#Length of the searched interval rint=', rint write(nout+20,*)'#Number of the points on the uniform mesh nx=',nx WRITE(NOUT+20,*)'#In columns:' WRITE(NOUT+20,*)'#position, the dipole nonradiative rates' write(nout+20,*)'#for tangentially and radially oriented dipole,' write(nout+20,*)'#and their orientational average, respectively' 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 do j=0,7 write(nout+j,*)'#Sphere radius in nm=', rsnm enddo 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 ***************************** ZEPS1 *********************************** * -------------------------------- * Optical constants in the case of a dispersion * 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: * zeps1=CCEPS zeps(2)=zeps0 if (nmat.le.1) then go to 2 !no reading of material data else call readmat(nmat,nfin,rsnm,rmuf,omf,ceps1) end if c lambda=2.d0*pi*rsnm/(omega*RMUF) omega =2.d0*pi*rsnm/(lambda*RMUF) cc xs=RMUF*omega*dble(sqrt(zeps0)) !Equiv. size parameter 2 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) zeps(1)=zeps1 *______________________________________ 8 continue if (imag(zeps1).eq.0.d0) then write(6,*)'Sphere is not absorbing!' write(6,*)'Check the set-up!' pause 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 parameter C >>> * go to 25 CALL GNZBESS(RX(I),LMX,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 go to 40 * 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 * * Assign polarizability: *==== do 150 imod=1,2 !nmd if (imod.eq.1) then !Rayleigh limit write(6,*)'Results for Rayleigh polarization' write(nout,*) write(nout+4,*) else if (imod.eq.2) then write(6,*)'Results for T-asymptotic polarization' write(nout+1,*) write(nout+5,*) else if (imod.eq.3) then write(6,*)'Results for T-polarization' write(nout+2,*) write(nout+6,*) end if zp=zeps1/zeps0 do 100 l=1,lmx if (imod.eq.1) then !Rayleigh limit ztmas=(dble(2*l+1)/l)*imag(zp)/abs(zp+dble(l+1)/l)**2 zal(l,imod)=(zp-cone)/(zp+dble(l+1)/l) ximal(l,imod)=ztmas else if (imod.eq.2) then !Tmat asymptotic zx=rx(2)**2/dble(2*(2*l+3)) !RX(I)=SG(I)*RMF(j) - size parameter zalph=ci*((l+1)*rx(2)**(2*l+1))/(xfac(l)*xfac(l+1)) if (l.eq.1) then zcrmn=(zp-cone)/(zp+ 2.d0 - (zp-cone)*zalph) end if zalph=(zp-cone)*(cone - (zp+cone)*zx)/ & (l*zp+dble(l+1) + zx*(-l*zp**2 - dble(3*(2*l+1))*zp/(2*l-1) + & dble((l+1)*(2*l+3))/(2*l-1)) - (zp-cone)*zalph) zal(l,imod)=dble(l)*zalph ximal(l,imod)=imag(zal(l,imod)) else if (imod.eq.3) then !Tmat exact ztmas=tt2(2,1,l,2,1)/tt2(1,1,l,2,1) zalph=-ci*xfac(l)*xfac(l+1)*ztmas/((l+1)*zk**(2*l+1)) !in gauss; -ci*6.d0*pi*ztmas/zk**3 in SI zal(l,imod)=dble(l)*zalph/rsnm**(2*l+1) ximal(l,imod)=imag(zal(l,imod)) end if 100 continue 150 CONTINUE !over imod write(nout+7,*) write(nout+8,*) write(nout+9,*) write(nout+10,*) write(nout+11,*) write(nout+12,*) write(nout+13,*) write(nout+14,*) write(nout+15,*) write(nout+16,*) OPEN(UNIT=30,FILE='Rnrl612r10l500R.dat') !Cnrl340r100l500f.dat rewind(30) do il=1,nx read(30,*) xd,x1,x2,xr(il) enddo close(30) *>>> scanning the interval RINT: C Prefactor 3/x^3, where x is the size parameter: xkoef=3.d0/rx(2)**3 !implicit veps_h-dependence via size parameter rx(2) c xkoef=3.d0*(lambda/rsnm)**3/(2.d0*pi)**3 ! prefactor 3.d0/(omega*rmuf/c)**3 do 200 il=1,nx !over the interval RINT delo= rint/dble(nx) !=4=800/200 ... too high!! rmff= 1.d0 + dble(il)*delo !distance from the sphere center in rsnm xd=rmff-1.d0 ! d/a cu delo= (rint-1.d0)/dble(nx) ct delo=0.5d0 ct rmff=(rsnm+ dble(il)*delo)/rsnm xq=1.d0/rmff ! a/r_1=1/(r_1/a) zv=zeps(2)/(zeps(1)+zeps(2)) do imod=1,nmd xpar=0.d0 xperp=0.d0 zdpar=czero zdperp=czero zff=czero xli2=0.d0 * Summation over l: do l=1,lmx if (l.ge.2) goto 50 * (B.18') of {GN}: gperp(imod)=abs(cone + 2.d0*xq**3*zal(l,imod))**2 * (B.43') of {GN}: gpar(imod)=abs(cone - xq**3*zal(l,imod))**2 50 continue * (B.16') of {GN} - the sum for Delta(perpendicular): zdperp=zdperp + dble(l+1)**2*zal(l,imod)*xq**(2*l+4) * (B.24') of {GN}: xperp=xperp + dble(l+1)**2*ximal(l,imod)*xq**(2*l+4) * (B.42') of {GN} - the sum for Delta(parallel): zdpar=zdpar + dble(l*(l+1))*zal(l,imod)*xq**(2*l+4) * (B.45') of {GN}: xpar = xpar + dble(l*(l+1))*ximal(l,imod)*xq**(2*l+4) zff=zff+ xq**(2*l)/(dble(l)*(dble(l)+zv)) !F(u,v) xli2=xli2 + xq**(2*l)/(dble(l)**2) !dilograrithm function enddo !over l * Adding prefactor to get the final expressions for the Delta's: c zdperp=xkoef*lambda*zdperp/(4.d0*pi*xc) !Delta(perpendicular) c zdpar =xkoef*lambda*zdpar/(8.d0*pi*xc) !Delta(parallel) zdperp=0.d0 zdpar =0.d0 C Assembling everything to obtain the nonradiative rates gnrpar and gnrperp: gnrperp(imod)=xkoef*xperp/ & (2.d0*abs(cone-zdperp)**2) gnrpar(imod) =xkoef*xpar/ & (4.d0*abs(cone-zdpar)**2) C Assembling everything to obtain the radiative rates gpar and gperp: gperp(imod)=dble(sqrt(zeps0))*gperp(imod)/ & abs(cone-zdperp)**2 gpar(imod) =dble(sqrt(zeps0))*gpar(imod)/ & abs(cone-zdpar)**2 C--------/---------/---------/---------/---------/---------/---------/-- * OUTPUT WRITING: * Nonradiative rates: gnrperp,gnrpar gim=gnrperp(imod)+2.d0*gnrpar(imod) !See Eq. (7) of PRA38_3410 dospar(il+1) = gnrpar(imod) !tangential oscillations dosperp(il+1)= gnrperp(imod) !radial oscillations dos(il+1) = gim/3.d0 !polarization average: there is !one radial and two indep. tangential !orientations possible ct for convergence test: ct dos(il+1)=dos(il+1)/xr(il) c WRITE(NOUT+imod-1,999) rmff, !xd, & dospar(il+1), & dosperp(il+1), & dos(il+1) * Radiative rates: gperp,gpar gim=gperp(imod)+2.d0*gpar(imod) !See Eq. (7) of PRA38_3410 dospar(il+1) = gpar(imod) !tangential oscillations dosperp(il+1)= gperp(imod) !radial oscillations dos(il+1) = gim/3.d0 !polarization average: there is !one radial and two indep. tangential !orientations possible WRITE(NOUT+imod+3,999) rmff, !xd, & dospar(il+1), & dosperp(il+1), & dos(il+1) enddo ! imod C--------/---------/---------/---------/---------/---------/---------/-- c Asymptotic coefficients: za(1)=(zeps(1)-zeps(2))/(4.d0*(zeps(1)+zeps(2))) za(2)=-(zeps(1)-zeps(2))*(zeps(1)+3.d0*zeps(2))/ & (8.d0*(zeps(1)+zeps(2))**2) za(3)=(zeps(1)-zeps(2))*(zeps(1)**2+3.d0*zeps(2)**2)/ & (8.d0*(zeps(1)+zeps(2))**3) za(4)=zeps(2)*zeps(1)**2*(zeps(1)-zeps(2))/(zeps(1)+zeps(2))**4 zb(1)=(zeps(1)-zeps(2))/(4.d0*(zeps(1)+zeps(2))) zb(2)=-(zeps(1)-zeps(2))*(3.d0*zeps(1)+5.d0*zeps(2))/ & (8.d0*(zeps(1)+zeps(2))**2) zb(3)=(zeps(1)-zeps(2))*(3.d0*zeps(1)**2+8.d0*zeps(2)*zeps(1)+ & 9.d0*zeps(2)**2)/(8.d0*(zeps(1)+zeps(2))**3) zb(4)=-zeps(2)**2*zeps(1)*(zeps(1)-zeps(2))/(zeps(1)+zeps(2))**4 xu=1.d0/rmff**2 xf1=xu**3*(1.d0+xu)/(1.d0-xu)**3 xfs=1.d0/(4.d0*xd**3) -5.d0/(8.d0*xd**2) +9.d0/(8.d0*xd) write(nout+20,*) write(nout+20,*) 'xd=', xd write(nout+20,*) 'xf1=', xf1,'xfs1=', xfs xfs=1.d0/(4.d0*xd**3) xf2=xu**3/(1.d0-xu)**2 xfs=1.d0/(4.d0*xd**2) -3.d0/(4.d0*xd) +23.d0/16.d0 write(nout+20,*) 'xf2=', xf2,'xfs2=', xfs xf3=xu**3/(1.d0-xu) xfs=1.d0/(2.d0*xd) -9.d0/4.d0 + 49.d0*xd/8.d0 -209.d0*xd**2/16.d0 write(nout+20,*) 'xf3=', xf3,'xfs3=', xfs C Dual expression for nonradiative rates: zperp=xf1 + xf2*(2.d0-zv) + xf3*(cone-zv)**2 + & zv*(cone-zv)**2*xu**2*log(cone-xu)+ zv**2*(cone-zv)**2*xu**2*zff zpar=xf1 + xf2*(cone-zv) - xf3*zv*(cone-zv) - & zv**2*(cone-zv)*xu**2*log(cone-xu)-zv**3*(cone-zv)*xu**2*zff xperp=imag((zeps(1)-zeps(2))*zperp/(zeps(1)+zeps(2))) xpar=imag((zeps(1)-zeps(2))*zpar/(zeps(1)+zeps(2))) gnrperps=xkoef*xperp/(2.d0*abs(cone-zdperp)**2) gnrpars =xkoef*xpar/(4.d0*abs(cone-zdpar)**2) gim=(gnrperps+2.d0*gnrpars)/3.d0 x1=gnrpars x2=gnrperps x3=gim ct for convergence test: ct gim=gim/xr(il) !comment out the line if not convergence test WRITE(NOUT+7,999) xd,gnrpars,gnrperps,gim C Dual expression for nonradiative rates up to dilogarithm without the remainder: zperp=xf1 + xf2*(2.d0-zv) + xf3*(cone-zv)**2 + & zv*(cone-zv)**2*xu**2*log(cone-xu)+ zv**2*(cone-zv)**2*xu**2*xli2 zpar=xf1 + xf2*(cone-zv) - xf3*zv*(cone-zv) - & zv**2*(cone-zv)*xu**2*log(cone-xu)-zv**3*(cone-zv)*xu**2*xli2 xperp=imag((zeps(1)-zeps(2))*zperp/(zeps(1)+zeps(2))) xpar=imag((zeps(1)-zeps(2))*zpar/(zeps(1)+zeps(2))) gnrperps=xkoef*xperp/(2.d0*abs(cone-zdperp)**2) gnrpars =xkoef*xpar/(4.d0*abs(cone-zdpar)**2) gim=(gnrperps+2.d0*gnrpars)/3.d0 ct for convergence test: ct gim=gim/xr(il) WRITE(NOUT+8,999) xd,gnrpars,gnrperps,gim !absolute values c WRITE(NOUT+8,999) xd,gnrpars/x1,gnrperps/x2,gim/x3 !relative values C Dual expression for nonradiative rates up to logarithm: zperp=xf1 + xf2*(2.d0-zv) + xf3*(cone-zv)**2 + & zv*(cone-zv)**2*xu**2*log(cone-xu) zpar=xf1 + xf2*(cone-zv) - xf3*zv*(cone-zv) - & zv**2*(cone-zv)*xu**2*log(cone-xu) xperp=imag((zeps(1)-zeps(2))*zperp/(zeps(1)+zeps(2))) xpar=imag((zeps(1)-zeps(2))*zpar/(zeps(1)+zeps(2))) gnrperps=xkoef*xperp/(2.d0*abs(cone-zdperp)**2) gnrpars =xkoef*xpar/(4.d0*abs(cone-zdpar)**2) gim=(gnrperps+2.d0*gnrpars)/3.d0 WRITE(NOUT+13,999) xd,gnrpars,gnrperps,gim !absolute values c WRITE(NOUT+13,999) xd,gnrpars/x1,gnrperps/x2,gim/x3 !relative values C Dual 3 terms expression for nonradiative rates: zperp=xf1 + xf2*(2.d0-zv) + xf3*(cone-zv)**2 zpar=xf1 + xf2*(cone-zv) - xf3*zv*(cone-zv) xperp=imag((zeps(1)-zeps(2))*zperp/(zeps(1)+zeps(2))) xpar=imag((zeps(1)-zeps(2))*zpar/(zeps(1)+zeps(2))) gnrperps=xkoef*xperp/(2.d0*abs(cone-zdperp)**2) gnrpars =xkoef*xpar/(4.d0*abs(cone-zdpar)**2) gim=(gnrperps+2.d0*gnrpars)/3.d0 WRITE(NOUT+14,999) xd,gnrpars,gnrperps,gim c WRITE(NOUT+14,999) xd,gnrpars/x1,gnrperps/x2,gim/x3 C Dual 2 terms expression for nonradiative rates: zperp=xf1 + xf2*(2.d0-zv) zpar=xf1 + xf2*(cone-zv) xperp=imag((zeps(1)-zeps(2))*zperp/(zeps(1)+zeps(2))) xpar=imag((zeps(1)-zeps(2))*zpar/(zeps(1)+zeps(2))) gnrperps=xkoef*xperp/(2.d0*abs(cone-zdperp)**2) gnrpars =xkoef*xpar/(4.d0*abs(cone-zdpar)**2) gim=(gnrperps+2.d0*gnrpars)/3.d0 WRITE(NOUT+15,999) xd,gnrpars,gnrperps,gim !absolute values c WRITE(NOUT+15,999) xd,gnrpars/x1,gnrperps/x2,gim/x3 !relative values C Dual first term expression for nonradiative rates: zperp=xf1 zpar=xf1 xperp=imag((zeps(1)-zeps(2))*zperp/(zeps(1)+zeps(2))) xpar=imag((zeps(1)-zeps(2))*zpar/(zeps(1)+zeps(2))) gnrperps=xkoef*xperp/(2.d0*abs(cone-zdperp)**2) gnrpars =xkoef*xpar/(4.d0*abs(cone-zdpar)**2) gim=(gnrperps+2.d0*gnrpars)/3.d0 WRITE(NOUT+16,999) xd,gnrpars,gnrperps,gim !absolute values c WRITE(NOUT+16,999) xd,gnrpars/x1,gnrperps/x2,gim/x3 !relative values C The dominant d**(-3) asymptotic of nonradiative rates: xperp=imag(za(1))/xd**3 xpar =imag(zb(1))/xd**3 gnrperps=xkoef*xperp/(2.d0*abs(cone-zdperp)**2) gnrpars =xkoef*xpar/(4.d0*abs(cone-zdpar)**2) gim=(gnrperps+2.d0*gnrpars)/3.d0 WRITE(NOUT+9,999) xd,gnrpars,gnrperps,gim !absolute values c WRITE(NOUT+9,999) xd,gnrpars/x1,gnrperps/x2,gim/x3 !relative values C Asymptotic nonradiative rates - 2 orders: xperp=imag(za(1))/xd**3 +imag(za(2))/xd**2 xpar=imag(zb(1))/xd**3 +imag(zb(2))/xd**2 gnrperps=xkoef*xperp/(2.d0*abs(cone-zdperp)**2) gnrpars =xkoef*xpar/(4.d0*abs(cone-zdpar)**2) gim=(gnrperps+2.d0*gnrpars)/3.d0 WRITE(NOUT+10,999) xd,gnrpars,gnrperps,gim !absolute values c WRITE(NOUT+10,999) xd,gnrpars/x1,gnrperps/x2,gim/x3 !relative values C Asymptotic nonradiative rates - 3 orders: xperp=imag(za(1))/xd**3 +imag(za(2))/xd**2+imag(za(3))/xd xpar=imag(zb(1))/xd**3 +imag(zb(2))/xd**2+imag(zb(3))/xd gnrperps=xkoef*xperp/(2.d0*abs(cone-zdperp)**2) gnrpars =xkoef*xpar/(4.d0*abs(cone-zdpar)**2) gim=(gnrperps+2.d0*gnrpars)/3.d0 WRITE(NOUT+11,999) xd,gnrpars,gnrperps,gim !absolute values c WRITE(NOUT+11,999) xd,gnrpars/x1,gnrperps/x2,gim/x3 !relative values C Asymptotic nonradiative rates up to dilogarithm: C zff=F(u,v); xli2=dilograrithm function rff1=imag( (zeps(1)-zeps(2))*zv**2*(cone-zv)**2*xu**2*zff & /(zeps(1)+zeps(2)) ) rff2=-imag( (zeps(1)-zeps(2))*zv**3*(cone-zv)*xu**2*zff & /(zeps(1)+zeps(2)) ) xperp=imag(za(1))/xd**3 +imag(za(2))/xd**2+imag(za(3))/xd+ & imag(za(4))*log(2.d0*xd) + rff1 xpar=imag(zb(1))/xd**3 +imag(zb(2))/xd**2+imag(zb(3))/xd+ & imag(zb(4))*log(2.d0*xd) + rff2 gnrperps=xkoef*xperp/(2.d0*abs(cone-zdperp)**2) gnrpars =xkoef*xpar/(4.d0*abs(cone-zdpar)**2) gim=(gnrperps+2.d0*gnrpars)/3.d0 WRITE(NOUT+12,999) xd,gnrpars,gnrperps,gim !absolute values c WRITE(NOUT+12,999) xd,gnrpars/x1,gnrperps/x2,gim/x3 !relative values * Ratios: xrperp= gnrperp(2)/gnrperp(1) xrpar = gnrpar(2) /gnrpar(1) xrperp= gnrperp(3)/gnrperp(1) xrpar = gnrpar(3) /gnrpar(1) WRITE(NOUT+3,998) rmff,xrpar,xrperp 200 CONTINUE !over r 300 CONTINUE !over rsnm 400 CONTINUE !over lambda close(nout) close(nout+1) close(nout+2) close(nout+3) close(nout+4) close(nout+5) close(nout+6) close(nout+7) close(nout+8) close(nout+9) close(nout+10) close(nout+11) close(nout+12) close(nout+13) close(nout+14) close(nout+15) close(nout+16) close(nout+20) write(6,*)'Gersten&Nitzan nonradiative rates in' write(6,*) write(6,*)'Rayleignnr.dat for static Rayleigh polarization' write(6,*) write(6,*)'Tasgnnr.dat for T-asymptotic polarization' write(6,*) write(6,*)'Tmatgnnr.dat for T-matrix polarization' write(6,*)'Ratios.dat for the ratios of T-matrix & to Rayleigh polarization results' write(6,*)'Gersten&Nitzan radiative rates in' write(6,*) write(6,*)'Rlr.dat for static Rayleigh polarization' write(6,*) write(6,*)'Tasr.dat for T-asymptotic polarization' write(6,*) write(6,*)'Tmatr.dat for T-matrix polarization' write(6,*) write(6,*)'Asmpnr.dat for the short-distance asymptotic of & non-radiative rates' write(6,*)'Asmpnrd3.dat for the dominant d**(-3) short-distance & asymptotic of non-radiative rates' 998 format(f9.6,3X,2(e15.8,3x)) 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. 06/2011 Alexander Moroz