C (C) Copr. 06/2005 Alexander Moroz C Wave-scattering.com C http://www.wave-scattering.com/ C wavescattering@yahoo.com C C Last update: February 2005 C C This is EFFE2P, a computer code to calculate the effective C medium parameters within the context of the extended C Maxwell - Garnett theory (a quasi-static extention of the c Maxwell Garnett theory) in the case of a binary composite C (two different types of spheres embedded in a matrix) C C This code has been employed in the article C C V. Yannopapas and A. Moroz, Negative refractive index metamaterials C from inherently non-magnetic materials for deep infrared to terahertz C frequency ranges, published in C J. Phys.: Condens. Matter. 17, 3717-3734 (2005) C (http://www.iop.org/EJ/abstract/0953-8984/17/25/002). C C where you can also find an additional teoretical background for C the extended Maxwell - Garnett theory employed here. C C The program default options correspond to the cases considered in the above C article. (The generated results are shown in Fig. 4a of the article.) C C In order to modify the default option, read through the comments C in the program and description of the parameters and variables. C C OUTPUT: C Returns file effmed.dat containing the effective medium parameters C versus wavelength. C ===== C C TESTS: C C 1) THE RESULTS OBTAINED BY DEFAULT SETUP OF PARAMETERS C SHOULD COINCIDE WITH THOSE DISPLAYED IN FIG. 4a of the article. C In particular, the first result lines of effmed.dat with C the program default options should read as follows: C C 13209.000 -50.1126 80.2440 1.3153 0.00080 5.4064 9.7573 C 13206.646 -51.6593 69.7463 1.3162 0.00080 4.8056 9.5469 C 13204.292 -51.5292 60.5496 1.3171 0.00081 4.2896 9.2908 C 13201.939 -50.3046 52.6734 1.3180 0.00081 3.8505 9.0095 C 13199.587 -48.4108 46.0136 1.3189 0.00082 3.4787 8.7172 C 13197.236 -46.1410 40.4189 1.3198 0.00082 3.1645 8.4229 C 13194.886 -43.6896 35.7316 1.3208 0.00083 2.8993 8.1326 C 13192.536 -41.1812 31.8070 1.3217 0.00083 2.6756 7.8495 C 13190.188 -38.6937 28.5202 1.3226 0.00084 2.4876 7.5756 C 13187.840 -36.2743 25.7673 1.3236 0.00084 2.3301 7.3118 C C======================================================================== C C REQUIREMENTS & NOT DISTRIBUTED DEPENDENT ROUTINES: C C You should supply your own routine C C GNZBESS(RX(I),LMAX,jl,drjl,nl,drnl) C C which generates the Bessel functions JL(0:LMAX),NL(0:LMAX) C of the first and second order and their derivatives C DJL(0:LMAX),DNL(0:LMAX) for a complex argument RX(I) in the C subroutine SPHEREF. (SPHEREF routine calculates the single C sphere scattering properties (including coated spheres) and C returns the elements of the diagonal T matrix.) 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 ===================================================================== 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 (that's me). I'd appreciate a reprint or e-copy C of any paper where my code has been used. C C - You can not sell it or use it in any way to get money out of it C (save for the research grants you might get thanks to it) C C - I'm making this code available to the best of my knowledge. C It has been tested as throughly as possible. C It work for Windows and Linux platforms. However, it's C your responsibility to check for its accuracy in the environment C and size/shape ranges you use. C C - Disclaimer: The usual stuff. I claim no responsibility for any C misuse or damage arising from its use. If this program erases your C hard disk, makes your computer cry out for Linux, C makes you lose your boy/girlfriend, melts your TV set, spends out C your lifesaving or forces you to get to decaf coffee, it's your C problem. For your tranquility, no such side effects has been C reported ... yet. 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 I would highly appreciate informing me of any problems encountered C with this code. C C Improvements, comments, additions: all welcome at C C wavescattering@yahoo.com C C ===================================================================== PROGRAM effe2p C Variables declared but never referenced: C YNTEST declared at line 12 file effe.f C Variables set but never used: C LMX set at line 20 file effe.f C TOL set at line 90 file effe.f C Variable ZEPS is used before its value has been defined C--------/---------/---------/---------/---------/---------/---------/-- C PROGRAM TO CALCULATE EFFECTIVE REFRACTIVE INDEX USING THE C RENORMALIZED MAXWELL-GARNETT FORMULA FOR A (COATED) SPHERE C C f77 -g -C effe.f -o rneffe C--------/---------/---------/---------/---------/---------/---------/-- implicit none integer NOUT,LMX,LCS,ILCS,NMAT,NFIN,NSTEP,NPC real*8 TOL logical ynperfcon,ynperfconv,yntest,ynbrug COMPLEX*16 ZEPS0,CCEPS,CSEPS COMPLEX*16 ZMU0,CCMU,CSMUS character*1 ync,yncv C ::: output unit number PARAMETER (NOUT=60) c number of spherical harmonics used PARAMETER (lmx=1) c number of different phases parameter (npc=2) * * If particle is coated, ync=y, otherwise ync=n (lower case only!) parameter (ync='n') * * ynperfcon=.true. if core is a perfect conductor, otherwise * ynperfcon=.false. * PARAMETER (ynperfcon=.false.) * c number of coatings parameter (lcs=1) c The coating layer to which material data are read in parameter (ilcs=1) c if coated, the ratio 'core radius/particle radius' c (If lcs.ne.1, program is singular for rff=0. and 1.) - use homogeneous c particle instead!!! c PARAMETER (rff=0.95d0) c background dielectric constant PARAMETER (ZEPS0=1.D0**2) c c background permeability PARAMETER (ZMU0=1.D0**2) c c Temporarily, before the coated part is finished, c read in particle core dielectric function to CSEPS c particle (core) dielectric constant (depending whether lcs=1 or lcs=2) c If LCS.EQ.1, CCEPS is particle dielectric constant c PARAMETER (CCEPS=(1.45D0,0.d0)**2) !SiO2 PARAMETER (CCEPS=(-50.84D0,0.762d0)) !JAP89_5774 ellipsoid for ld=633 c PARAMETER (CCEPS=(-2.03D0,0.602d0)) !JAP89_5774 sphere for ld=354 C >>> SPHERE (OUTER SHELL SCATTERER) PERMITTIVITY <<< * n(silica)=1.45 <---> EPS(1)=2.1025D0 * n(ZnS)=2. <---> EPS(1)=4.D0 PARAMETER (CSEPS=(1.45d0,0.d0)**2) c PARAMETER (CSEPS=(-50.84D0,0.762d0)) !=(-10.84D0,0.762d0)) JAP89_5774 c If LCS.EQ.1, CCMU is particle permeability PARAMETER (CCMU=(1.D0,0.d0)) C >>> SPHERE (OUTER SHELL SCATTERER) PERMEABILIITY <<< PARAMETER (CSMUS=(1.D0,0.d0)) * 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=1) * c Temporarily option for reading of the real data for the dielectric constant c The number of the entries in a material data file to be read below c AGC.DAT NFIN=73 ! from Palik c Audat.dat NFIN=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=66 * PARAMETER (NFIN=66) * C ::: relative error allowed for the TCS. If the convergence * within TOL is not reached, program issues warning PARAMETER (TOL=7.d-2) * c If ynbrug=.true., performs Bruggeman approximation for ZEPS1. Otherwise c ynbrug=false. parameter (ynbrug=.false.) ****************************************************************** c Declarations: integer ikl,ieps,istep,ipc REAL*8 rmuff,lambda,ff(npc),ffb,fft,pi,filfrac,enw,xstep, !ffb,filfrac is for Bruggemann & RMF(lcs),rff(lcs),RMUF(npc),rsnm(npc) real*8 omf(NFIN),omxf,reepsz,plasma,omxp real*8 delo,omega0,omega,omegal,p1,p2 COMPLEX*16 ci,cone,xs,zeff,zneff,zmuff,z1,z2,zx,zy, & zalphae(npc),zalpham(npc),ca(npc),cb(npc),cc(npc),cd(npc), & cam(npc),cbm(npc),ccm(npc),cdm(npc) complex*16 ceps1(NFIN),ZEPS1,zeps(lcs+1),zmu(lcs+1) COMPLEX*16 TMT(4,3,3) * COMMON /TOSPHERECR/ rmuff COMMON /TOSPHERECH/ yncv COMMON /TOSPHERECL/ ynperfconv * * From here to spheref C C--------/---------/---------/---------/---------/---------/---------/-- c Data: DATA PI/3.141592653589793d0/ DATA ci/(0.d0,1.d0)/,cone/(1.d0,0.d0)/ *--------------------------------------------------------------- CCCCCCCCCCCCCCCCCC Assignement of common variables CCCCCCCCCCC * ynperfconv=ynperfcon yncv=ync ***************************************************************** c write(6,*)'The (coated) sphere radius in nanometers' c read(5,*) rsnm rsnm(1)=446d0 !almost close-packed LiTaO_{3} spheres p1=0.5d0 !probability of sphere 1 rsnm(2)=358d0/sqrt(zeps0) !n-type Ge spheres p2=1-p1 !probability of sphere 2 * sphere radius in the units, where * the side ALPHA of a conventional unit cell of a cubic lattice * is equal to 2 (then OMEGA=PI*ALPHA/LAMBDA). Various examples * of the sphere radius for different sphere filling fractions * and different lattices are listed below: RMUF(1)=rsnm(1)*2.D0/1264d0 ff(1)=p1*2.d0*pi*RMUF(1)**3/3.d0 RMUF(2)=rsnm(2)*2.D0/1264d0 ff(2)=p2*2.d0*pi*RMUF(2)**3/3.d0 fft=ff(1)+ff(2) !total filling fraction of spheres *oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo * Scanning over frequency interval: *------------------------------------------------- * ct write(6,*)'READ INITIAL (MAXIMAL) LAMBDA (in the vacuum) in nm' ct read(5,*) lambda cv lambda=633.d0 !JAP89_5776 for gold ellipsod c lambda=354.d0 !JAP89_5776 for silver sphere cc lambda=500d0 lambda=13209d0 !299792.458d0/26.7d0/0.85d0 !lambda_T=11228,1819 !lambda_L= 6392,1632 * * * size parameter is customarily defined as the ratio of * circumference of particle to the wavelength in the host medium * in which the particle is embedded * x=kr=\sg a=2*pi*r/lambda * xs=2.d0*pi*rs*dble(sqrt(zeps0))/lambda * convert lambda to the lambda in vacuum: c lambda=lambda*dble(sqrt(zeps0)) c omega=2.d0*pi*rs/(lambda*rmuf)=xs/(rmuf*dble(sqrt(zeps0))), c where rs is the particle radius (in nm) and lambda c is the wavelengths (in nm) c in the vacuum: omega=2.d0*pi*rsnm(1)/(lambda*rmuf(1)) omega0=omega c write(6,*)'Read omega' c read(5,*) omega c WRITE(6,*)'omega=', omega c if (1.eq.1) nstep=0 ! temporarily c delo=0.d0 c omega0=omega c go to 11 * * Option for omega input: c write(6,*)'Read omega =' c read(5,*) omega * c xs=RMUF*omega*dble(sqrt(zeps0)) * c write(6,*)'Equiv. size parameter x=2*pi*rs*n_0/lambda=',xs * ct write(6,*)'Scan down to wavelength (in nm)' ct read(5,*) enw c enw=500 enw=11224d0 !299792.458d0/46.9d0/0.95d0 ct write(6,*)'Scanning step in nm' ct read(5,*) xstep xstep=2 * C ::: number of steps on frequency interval: if (lambda.eq.enw) then NSTEP=0 delo=0.d0 else NSTEP=(lambda-enw)/xstep C ::: width of the searched frequency interval ENW=2.d0*pi*rsnm(1)/(enw*rmuf(1)) enw=enw-omega0 delo=enw/dble(nstep) end if * ZEPS(1)=cceps if((lcs.gt.1).and. & ((nmat.eq.0).or.((nmat.ne.0).and.(ilcs.ne.lcs)))) & zeps(lcs)=cseps zeps(lcs+1)=zeps0 ZMU(1)=ccmu if((lcs.gt.1).and. & ((nmat.eq.0).or.((nmat.ne.0).and.(ilcs.ne.lcs)))) & zmu(lcs)=csmus zmu(lcs+1)=zmu0 * -------------------------------- * output initial statements OPEN(UNIT=NOUT,FILE='effmed.dat') rewind(NOUT) WRITE(NOUT,*)'#EFFECTIVE MEDIUM PARAMETERS' write(nout,*) write(nout,*)'The effective dielectric constant obtained in the & non-quasi-static case is complex even in the case of & non-absorbing spheres embedded in a non-absorbing medium. & Therefore, except for the quasi-static case, the imaginary & of the effective dielectric constant of the extended & MG theories cannot be used to evaluate the absorption & or the heating rate.' write(nout,*) write(nout,*)'#Spheres filling fraction=', fft write(nout,*)'#Sphere 1 radius in nm=', rsnm(1) write(nout,*)'#Sphere 1 radius rmuf in relative units=', rmuf(1) write(nout,*)'#Sphere 2 radius in nm=', rsnm(2) write(nout,*)'#Sphere 2 radius rmuf in relative units=', rmuf(2) if (.not.ynperfcon) then write(nout,*)'#Material number =',NMAT else write(nout,*)'#Perfectly conducting sphere (core)' end if if (ynbrug) write(nout,*)'#Bruggeman approximation performed' if (ync .eq.'n') then write(nout,*)'#Homogeneous particle' else if (ync.eq.'y') then write(nout,*)'#Coated particle with (including core)', & LCS,'shells' end if if (nmat.ge.1) write(nout,*)'#Dispersive layer number=', ilcs write(nout,*) write(nout,*)'#host dielectric constant=', zeps(lcs+1) 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) C--------/---------/---------/---------/---------/---------/---------/-- end if WRITE(NOUT,*) WRITE(NOUT,*)'#In columns: lambda_0 and the respective Re and Im' WRITE(NOUT,*)'#parts of eff. diel. constant, eff. mu, and n_eff' write(nout,*) C--------/---------/---------/---------/---------/---------/---------/-- C--------/---------/---------/---------/---------/---------/---------/-- * Checking set up: * if (NMAT.GT.1) write(6,*)'Real material data & are to be provided' * if ((ync.eq.'y'.and.lcs.eq.1).or.(ync.eq.'n'.and.lcs.ne.1)) then write(6,*)'Check compatibility of YNC and LCS' stop end if C-------------------------------------------------------------------- * if (lcs.ge.2) then !coated particle * cc write(6,*)'Read equal-volume-core radius in nm' cc read(5,*) rff(1) cc rff(1)=204.9d0 cc rff(1)=rff(1)/rs write(6,*)'Coated sphere shell 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=1,lcs-1 * write(6,*)'Read r(l) for l=',ikl read (5,*) rff(ikl) rff(ikl)=rff(ikl)/rsnm(ipc) c rff(1)=0.75d0 rmf(ikl)=rff(ikl)*rmuff IF ((IKL.GE.2).and.(IKL.LT.LCS)) THEN if ((nmat.eq.0).or.((nmat.ne.0).and.(ilcs.ne.ikl))) then write(6,*)'Read the lth-sphere layer diel. const. for l=',ikl read (5,*) zeps(ikl) end if * write(6,*)'Read the lth-sphere layer magnetic permeability & for l=',ikl read (5,*) zmu(ikl) END IF enddo !IKL.GE.2 * end if !lcs.ge.2 ***************************** ZEPS1 *********************************** * -------------------------------- * Sphere optical constants in the case of a dispersion * READING IN MATERIAL DATA: * Reading real silver data according to Palik's book * requires reading data files OMF and CEPS1 of dimension NFIN * OMF is reepsz/omega and CEPS1 contains the sphere EPS * 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(1)/(omf(ieps)*rmuf(1)) 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(1)/(omf(ieps)*rmuf(1)) 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(1)/(omf(ieps)*rmuf(1)) 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(1)/(omf(ieps)*rmuf(1)) enddo close(30) end if ! material constant reading ********************* * -------------------------------- * begin main scanning loop: do 200 istep=1,nstep+1 omega=omega0 + dble(istep-1)*delo C--------/---------/---------/---------/---------/---------/---------/-- do 100 ipc=1,npc xs=RMUF(ipc)*omega*dble(sqrt(zeps0)) !Equiv. size parameter rmuff=RMUF(ipc) rmf(lcs)=rmuff 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*particle radius in nm/(lambda_z in nm*rmuf) * where lambda_z is the wavelength for which Re eps_s=0. C reepsz=2.d0*pi*rsnm(1)/(324.269d0*rmuf(1)) !setting to Ag plasma frequency C plasma=reepsz * omega_p=omega_T of 2*pi*26.7 THZ corresponds to * in local frequency units IF (NMAT.EQ.1) THEN !Material decision IF - Drude metal * c = 2.99792458d8 m/s * THz is equivalent to $\ld=2.99792458\times 10^{-4} m$, or to $299,792458$ $\mu$m. plasma=2.d0*pi*rsnm(1)*26.7d0/(299792.458d0*rmuf(1)) !omega_T omegal=2.d0*pi*rsnm(1)*46.9d0/(299792.458d0*rmuf(1)) !omega_L if (ipc.eq.1) then c omegal=2.d0*pi*rsnm(1)*46.9d0/(299792.458d0*rmuf(1)) zeps1=13.4d0*(1.d0 + (omegal**2-plasma**2)/(plasma**2-omega**2)) else if (ipc.eq.2) then c plasma=omegal omxp=plasma*1.d0/omega !=plasma/omega zeps1=15.8d0*(1.d0 - omxp**2) end if 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 !setting to Au plasma frequency * 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 END IF ! END of Material decision IF *The end of reading real data according to Palik's book *_____________________________________ * activate Bruggeman: 5 if (ynbrug) then ffb=0.8d0 z1 = (3.d0*ffb-1.d0)*zeps1+(2.d0 - 3.d0*ffb)*zeps0 z2 = sqrt(z1*z1 + 8.d0*zeps1*zeps0) * if (IMAG(z2).GE.0.0) then zeps1= (z1 + z2)/4.d0 else zeps1= (z1 - z2)/4.d0 end if end if zeps(ilcs)=zeps1 *______________________________________ 8 lambda=2.d0*pi*rsnm(1)/(omega*RMUF(1)) * call spheref(1,lcs,lambda,rsnm(ipc),rmuf(ipc),zeps,zmu,tmt) C--------/---------/---------/---------/---------/---------/---------/-- C >>> lmax,lcs,lambda,rsnm(ipc),rmf(ipc),zeps,zmu C <<< tmt C =============== C This routine calculates the single sphere scattering properties C (including coated spheres) C C Outputs the elements of the diagonal T matrix C C | TMT(1,*) | 0 | C TMT = | ----------+-------------| C | 0 | TMT(2,*) | C C TMT(1,*) terms corresponds to TEE scattering matrices C TMT(2,*) terms corresponds to TMM scattering matrices C C TMT's equal to i*sin(eta)*exp(i*eta), where eta is a phase-shift C--------/---------/---------/---------/---------/---------/---------/-- c dielectric polarization factor: zalphae(ipc)=-ci*3.d0*tmt(1,1,1)/(2.d0*xs**3) ca(ipc)= zeps0*(cone +2.d0*zalphae(ipc)*fft) cb(ipc)= cone - zalphae(ipc)*fft cc(ipc)= zeps0*(cone + 2.d0*zalphae(ipc)*fft) cd(ipc)= 2.d0*cone -2.d0*zalphae(ipc)*fft c magnetic polarization factor: zalpham(ipc)=-ci*3.d0*tmt(2,1,1)/(2.d0*xs**3) zmuff=zmu0*(cone + 2.d0*zalpham(ipc)*ff(ipc))/ & (cone - zalpham(ipc)*ff(ipc)) write(6,*)'MG zmuff=', zmuff cam(ipc)= zmu0*(cone +2.d0*zalpham(ipc)*fft) cbm(ipc)= cone - zalpham(ipc)*fft ccm(ipc)= zmu0*(cone + 2.d0*zalpham(ipc)*fft) cdm(ipc)= 2.d0*cone -2.d0*zalpham(ipc)*fft 100 continue ! ipc-loop over different spheres c Luo two-phase Maxwell-Garnett formula zx=-(p1*(ca(1)*cd(2)-cb(1)*cc(2))+p2*(ca(2)*cd(1)-cb(2)*cc(1)))/ & (p1*cb(1)*cd(2)+p2*cb(2)*cd(1)) zy=-(p1*ca(1)*cc(2)+p2*ca(2)*cc(1))/ & (p1*cb(1)*cd(2)+p2*cb(2)*cd(1)) C--------/---------/---------/---------/---------/---------/---------/-- zeff=(-zx + sqrt(zx**2-4.d0*zy))/2.d0 if (imag(zeff).lt.0) zeff=(-zx - sqrt(zx**2-4.d0*zy))/2.d0 zx=-(p1*(cam(1)*cdm(2)-cbm(1)*ccm(2))+p2*(cam(2)*cdm(1) & -cbm(2)*ccm(1)))/(p1*cbm(1)*cdm(2)+p2*cbm(2)*cdm(1)) zy=-(p1*cam(1)*ccm(2)+p2*cam(2)*ccm(1))/ & (p1*cbm(1)*cdm(2)+p2*cbm(2)*cdm(1)) C--------/---------/---------/---------/---------/---------/---------/-- zmuff=(-zx + sqrt(zx**2-4.d0*zy))/2.d0 if (imag(zmuff).lt.0) zmuff=(-zx - sqrt(zx**2-4.d0*zy))/2.d0 zneff=sqrt(zeff*zmuff) if (imag(zneff).lt.0) zneff=-sqrt(zeff*zmuff) C write(6,*) write(6,*) write(6,*)'LAMBDA=',lambda write(6,*) write(6,*)'QUASISTATIC VALUES:' write(6,*)'effective dielectric constant e_eff=', zeff write(6,*)'effective permeability mu_eff=', zmuff write(6,*)'effective complex refractive index n_eff=',zneff C write(nout,1105) lambda, zeff, zmuff, zneff C 200 continue * close(nout) write(6,*) write(6,*)'Effective medium parameters versus wavelength in & effmed.dat' write(6,*) write(6,*)'The effective dielectric constant obtained in the & non-quasi-static case is complex even in the case of & non-absorbing spheres embedded in a non-absorbing medium. & Therefore, except for the quasi-static case, the imaginary & of the effective dielectric constant of the extended & MG theories cannot be used to evaluate the absorption & or the heating rate.' C--------/---------/---------/---------/---------/---------/---------/-- 1105 FORMAT(F12.3,4X,F10.4,2X,F10.4,4X,F10.4,2X,F8.5,4X,F10.4,2X,F10.4) END C******************************************************************** subroutine spheref(lmax,lcs,lambda,rsnm,rmf,zeps,zmu,tmt) C--------/---------/---------/---------/---------/---------/---------/-- C >>> lmax,lcs,lambda,rsnm,rmf,zeps,zmu C <<< tmt C C Outputs the elements of the diagonal T matrix C C | TMT(1,*) | 0 | C TMT = | ----------+-------------| C | 0 | TMT(2,*) | C C TMT(1,*) terms corresponds to TEE scattering matrices C TMT(2,*) terms corresponds to TMM scattering matrices C C TMT's equal to i*sin(eta)*exp(i*eta), where eta is a phase-shift C =============== C This routine calculates the single sphere scattering properties C (including coated spheres) C C C k_l length in units (2*PI/A=) PI: xkl= 0.8660254037844386d0 C C KMT ... elements of the K matrix C ALPHA, BETA ... arrays of the electric and magnetic phase shifts C C Partial wave expansion is used, which is badly convergent C for large size parameters $x> 100$. In numerical applications, the C series is to be cut off after C LMAX (LMX parameter here) \approx x+4x^{1/3}+2$. C In the case of LMX=50 this means that x <= 35. C If one wants to observe ripples, the cutoff for a given x has to C be even larger C--------------------------------------------------------------------- implicit none integer LMX,LCS,LMAX,LMAXD1,LMTD real*8 TOL character*1 ync logical ynperfcon c Parameters: c number of spherical harmonics used PARAMETER (lmx=1,LMAXD1=LMX+1,LMTD=LMAXD1*LMAXD1-1) c maximal number of sphere coatings cc parameter (lcs=3) c C ::: relative error allowed for the TCS. If the convergence * within TOL is not reached, program issues warning PARAMETER (TOL=1.d-6) * C--------/---------/---------/---------/---------/---------/---------/-- c Declarations: integer i,j,l,ikl,ij,ij1,m real*8 ZEPS0 real*8 rmf(lcs),rmuf,xs,rsnm,pi real*8 omega,lambda * COMPLEX*16 ci,czero,zeps(lcs+1),zmu(lcs+1),cqeps(2),cqmu(2) COMPLEX*16 RX(2),SG(2),zs,zz1,zz2 COMPLEX*16 KMT(2,lmx),TMT(4,LMTD,LMTD),alpha(lmx),beta(lmx) * COMPLEX*16 cm(lcs,lmx),dm(lcs,lmx),ce(lcs,lmx),de(lcs,lmx) COMPLEX*16 tt1(2,2,lmx,2),tt2(2,2,lmx,2) COMPLEX*16 AM(lmx),AE(lmx),BM(lmx),BE(lmx) * COMPLEX*16 JL(0:lmx),NL(0:lmx) COMPLEX*16 DRJL(0:lmx),DRNL(0:lmx) COMPLEX*16 UL(2,0:lmx),VL(2,0:lmx) COMPLEX*16 DRUL(2,0:lmx),DRVL(2,0:lmx) COMPLEX*16 ZARTAN * external ZARTAN * * COMMON /TOSPHERECR/ rmuf !rmuf(ipc) COMMON /TOSPHERECH/ ync COMMON /TOSPHERECL/ ynperfcon * * From the main here *--------------------------------------------------------------- c Data: DATA PI/3.141592653589793d0/ DATA ci/(0.d0,1.d0)/,czero/(0.D0,0.D0)/ C--------/---------/---------/---------/---------/---------/---------/-- * Initialization: do ij=1,lmtd do ikl=1,lmtd tmt(3,ikl,ij)=czero tmt(4,ikl,ij)=czero if (ikl.ne.ij) then tmt(1,ikl,ij)=czero tmt(2,ikl,ij)=czero end if enddo enddo C--------/---------/---------/---------/---------/---------/---------/-- * Checking set up: if ((ync.eq.'y'.and.lcs.eq.1).or.(ync.eq.'n'.and.lcs.ne.1)) then write(6,*)'Check compatibility of YNC and LCS' stop end if C-------------------------------------------------------------------- * Reading in the input data: c close packed : c RMUF=1.d0/DSQRT(2.D0) * * size parameter is customarily defined as the ratio of * circumference of sphere to the wavelength in the host medium * in which the sphere is embedded * x=kr=\sg a=2*pi*r/lambda * xs=2.d0*pi*rs*dble(sqrt(zeps0))/lambda * convert lambda to the lambda in vacuum: c lambda=lambda*dble(sqrt(zeps0)) c omega=2.d0*pi*rsnm/(lambda*rmuf)=xs/(rmuf*dble(sqrt(zeps0))), c where rs is the sphere radius (in nm) and lambda is the wavelengths (in nm) c in the vacuum: * * Option for omega input: c write(6,*)'Read omega =' c read(5,*) omega * zeps0=zeps(lcs+1) omega=2.d0*pi*rsnm/(lambda*rmuf) !=2.d0*pi*rsnm(ipc)/(lambda*rmf(ipc)) c xs=2.d0*pi*rsnm*sqrt(zeps0)/lambda c xs=RMUF*omega*dble(sqrt(zeps0)) * c write(6,*)'Size parameter x=2*pi*rs*n_0/lambda=', xs * * -------------------------------- *______________________________________ * if (.not.ynperfcon) then ij1=1 do l=1,lmx AM(l)=dcmplx(1.d0,0.d0) AE(l)=dcmplx(1.d0,0.d0) BM(l)=dcmplx(0.d0,0.d0) BE(l)=dcmplx(0.d0,0.d0) enddo else if (ynperfcon) then CQEPS(2)=SQRT(ZEPS(2)) CQMU(2)=SQRT(ZMU(2)) SG(2)=omega*SQRT(ZEPS(2)*ZMU(2)) RX(1)=SG(2)*RMF(1) * * call gnzbess(RX(1),LMX,jl,drjl,nl,drnl) * DO 10 L=1,lmx C >>> (AS 10.1.22): UL(1,L)=RMF(1)*JL(L) VL(1,L)=RMF(1)*NL(L) DRJL(L)=SG(2)*DRJL(L) DRNL(L)=SG(2)*DRNL(L) DRUL(1,L)=JL(L)+RMF(1)*DRJL(L) DRVL(1,L)=NL(L)+RMF(1)*DRNL(L) AM(l)= NL(L) ! cm(1,l) BM(l)=-JL(L) ! dm(1,l) AE(l)= DRVL(1,L) ! ce(1,l) BE(l)=-DRUL(1,L) ! de(1,l) * cf. Jackson 1962, p. 571, Eqs. (16.147); * B/A should yield -tan(phase shift) 10 continue if (lcs.eq.1) go to 30 ij1=2 end if C******************************************************************** c Execution: * Calculation of the phase shifts * DO 28 j=ij1,lcs CQEPS(1)=SQRT(ZEPS(j)) CQMU(1)=SQRT(ZMU(j)) SG(1)=omega*SQRT(ZEPS(j)*ZMU(j)) * CQEPS(2)=SQRT(ZEPS(j+1)) CQMU(2)=SQRT(ZMU(j+1)) SG(2)=omega*SQRT(ZEPS(j+1)*ZMU(j+1)) * DO 25 I=1,2 * RX(I)=SG(I)*RMF(j) c WRITE(6,*)'i, rx(i)=', i, rx(i) C >>> * call gnzbess(RX(I),lmx,jl,drjl,nl,drnl) * c write(6,*)'jl=', jl DO 15 L=1,lmx C >>> (AS 10.1.22): UL(I,L)=RMF(j)*JL(L) VL(I,L)=RMF(j)*NL(L) DRJL(L)=SG(I)*DRJL(L) DRNL(L)=SG(I)*DRNL(L) DRUL(I,L)=JL(L)+RMF(j)*DRJL(L) DRVL(I,L)=NL(L)+RMF(j)*DRNL(L) 15 continue 25 CONTINUE * c write(6,*)'ul=', ul * C >>> END OF THE LOOP TO ASSIGN VALUES OF BESSEL FUNCTIONS C JL and NL start to oscillate after RX.GT. approx 2.5 C******************************************************************** C Transfer matrix for a layered (coated) sphere C******************************************************************** * do l=1,lmx * * magnetic part * tt1(1,1,l,1)= cqmu(1)*UL(1,L) tt1(1,2,l,1)= cqmu(1)*VL(1,L) tt1(2,1,l,1)= DRUL(1,L)/cqmu(1) tt1(2,2,l,1)= DRVL(1,L)/cqmu(1) * tt2(1,1,l,1)= sg(2)*DRVL(2,L)/cqmu(2) tt2(1,2,l,1)= - sg(2)*cqmu(2)*VL(2,L) tt2(2,1,l,1)= - sg(2)*DRUL(2,L)/cqmu(2) tt2(2,2,l,1)= sg(2)*cqmu(2)*UL(2,L) * * electric part * tt1(1,1,l,2)=cqeps(1)*UL(1,L) tt1(1,2,l,2)=cqeps(1)*VL(1,L) tt1(2,1,l,2)=DRUL(1,L)/cqeps(1) tt1(2,2,l,2)= DRVL(1,L)/cqeps(1) * tt2(1,1,l,2)= sg(2)*DRVL(2,L)/cqeps(2) tt2(1,2,l,2)= -sg(2)*cqeps(2)*VL(2,L) tt2(2,1,l,2)= -sg(2)*DRUL(2,L)/cqeps(2) tt2(2,2,l,2)= sg(2)*cqeps(2)*UL(2,L) * * m-part * cm(j,l)=AM(l)*(tt2(1,1,l,1)*tt1(1,1,l,1) 1 +tt2(1,2,l,1)*tt1(2,1,l,1))+BM(l)*( 2 tt2(1,1,l,1)*tt1(1,2,l,1)+tt2(1,2,l,1)*tt1(2,2,l,1)) * dm(j,l)=AM(l)*(tt2(2,1,l,1)*tt1(1,1,l,1) 1 +tt2(2,2,l,1)*tt1(2,1,l,1))+BM(l)*( 2 tt2(2,1,l,1)*tt1(1,2,l,1)+tt2(2,2,l,1)*tt1(2,2,l,1)) * * e-part * ce(j,l)=AE(l)*(tt2(1,1,l,2)*tt1(1,1,l,2) 1 +tt2(1,2,l,2)*tt1(2,1,l,2))+BE(l)*( 2 tt2(1,1,l,2)*tt1(1,2,l,2)+tt2(1,2,l,2)*tt1(2,2,l,2)) * de(j,l)=AE(l)*(tt2(2,1,l,2)*tt1(1,1,l,2) 1 +tt2(2,2,l,2)*tt1(2,1,l,2))+BE(l)*( 2 tt2(2,1,l,2)*tt1(1,2,l,2)+tt2(2,2,l,2)*tt1(2,2,l,2)) * AM(l)=cm(j,l) BM(l)=dm(j,l) AE(l)=ce(j,l) BE(l)=de(j,l) c write(6,*) AM(l), BM(l) c write(6,*) AE(l), BE(l) * enddo * 28 CONTINUE c write(6,*)'am=', am c write(6,*)'bm=', bm C--------/---------/---------/---------/---------/---------/---------/-- C ASSIGNING VALUES TO ELEMENTS OF THE K-MATRIX C >>> 30 CONTINUE * DO 40 L=1,lmx * In the following, one needs only phase shifts, so that division * by SG(2) is omitted: * KMT(1,L)=bm(l)/am(l) KMT(2,L)=be(l)/ae(l) * 40 CONTINUE c write(6,*)'kmt=', kmt C******************************************************************** DO 60 j=1,lmx * >>> Extracting phase-shifts from KMT * >>> magnetic part: zs=-kmt(1,j) alpha(j)=zartan(zs) * * Under normal circumstances, Im. part of a phase shift \geq 0 !!!) * if(dimag(alpha(j)).lt.0.d0) then write(6,*)'dimag(alpha(j))=',dimag(alpha(j)) ,' is negative' write(6,*)'omega, j=', omega, j pause end if c write(6,*)'j, alpha(j)=', j, alpha(j) *>>> electric part: zs=-kmt(2,j) beta(j)=zartan(zs) * * Under normal circumstances, Im. part of a phase shift \geq 0 !!!) * if(dimag(beta(j)).lt.0.d0) then write(6,*)'dimag(alpha(j))=',dimag(beta(j)) ,' is negative' write(6,*)'omega, j=', omega, j pause end if c write(6,*)'beta(j)=', beta(j) 60 CONTINUE do 100 l=1,lmax zz1 = ci*sin(beta(l))*exp(ci*beta(l)) zz2 = ci*sin(alpha(l))*exp(ci*alpha(l)) do 100 m=-l,l ij=l*(l+1)+m tmt(1,ij,ij)= zz1 tmt(2,ij,ij)= zz2 100 continue * return end C (C) Copr. 11/2003 Alexander Moroz C******************************************************************** function zartan(zs) *--------/---------/---------/---------/---------/---------/---------/-- * For a complex argument using that * atan (z)= -ln((1+iz)/(1-iz))*i/2.d0 * Hower, a direct application of this formula often * results in a nonzero imaginary part of $\arctan z$ even * for a purely real $z=x$. * Therefore, in order to avoid this, in a numerical * implementation, the complex logarithm is written * explicitly as * * ln\left(\fr{1+iz}{1-iz}\right)= (ar1,ar2)= log( * (1.d0 + abs(z)**2 -2.d0*imag(z))/1.d0+abs(z)**2 +2.d0*imag(z)))/2.d0 * + * ci*atan(2.d0*dble(z)/(1-abs(z)**2). * * For a real z=x, |x|<1, log here is purely imaginary and * equals to i*atan(2x/(1-x**2))\equiv 2*i*atan x *--------/---------/---------/---------/---------/---------/---------/-- complex*16 zartan,zs real*8 xxs,xas,ar1,ar2,pi DATA PI/3.141592653589793d0/ if (dimag(zs).eq.0.) then zartan=dcmplx(atan(dble(zs)),0.d0) return end if xas=abs(zs)**2 ar1=log((1.d0+xas-2.d0*dimag(zs))/(1.d0+xas+2.d0*dimag(zs)))/2.d0 xxs=dble(zs) * special case: if(xas.eq.1.) then if(xxs.ge.0.) then ar2=pi/2.d0 else if (xxs.lt.0.) then ar2=-pi/2.d0 end if zartan =dcmplx(ar2,- ar1)/2.d0 end if * remaining cases: ar2=2.d0*xxs/(1.d0-xas) if(xas.lt.1.d0) then ! 1st and 4th quadrant zartan=dcmplx(atan(ar2),- ar1)/2.d0 else if (xas.gt.1. .and. xxs.ge.0.) then ! 2nd quadrant zartan=dcmplx(pi+atan(ar2),- ar1)/2.d0 else if(xas.gt.1. .and. xxs.lt.0.) then ! 3rd quadrant zartan=dcmplx(-pi+atan(ar2),- ar1)/2.d0 end if return end C (C) Copr. 1/1999 Alexander Moroz C******************************************************************** 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 * C (C) Copr. 04/2002 Alexander Moroz