program sphere C--------/---------/---------/---------/---------/---------/---------/-- C This routines calculates the single sphere scattering properties C (including coated spheres) C C The comments lines below are intended to provide a self-consistent C description of the program parameters C C Requires the following supplementary routines: C C 1) the Bessel function package bessel.zip available at C http://www.wave-scattering.com/bessel.zip C C 2) material data routines sordalc.f available at C http://www.wave-scattering.com/sordalc.f C and znsrefind.f available at C http://www.wave-scattering.com/znsrefind.f C C 3) Complex arctan routine zartan.f available at C http://www.wave-scattering.com/zartan.f C C 4) If NMAT.gt.1 you will need to supply the required C material data. C C The gold data file "Audat.dat" (the option NMAT=3) can be downloaded at C http://www.wave-scattering.com/Audat.dat.txt C For security reasons, *.dat files cannot be posted; therefore C the additional *.txt extension. You have to rename the data file C back to Audat.dat after the download. C C A makefile "makesssc" for UNIX/LINUX compilers, which C for some computer security reasons, which do not allow to distribute C *.dat files and files without an extenmsion, can be downloaded as C http://www.wave-scattering.com/makesssc.txt C (After a download it has to be renamed to makesssc) C Check your local compiler options and, if necessary, C amend the makefile appropriately. C Save all the required files within the same directory and issue C the command C C make -f makesssc C C to compile the source code. C C C Outputs, among other, the following cross sections (cs): C C Extinction versus wavelength in ssext.dat C Scattering (elastic) cs versus wavelength in sssccs.dat C Absorption versus wavelength in ssabs.dat C Albedo versus wavelength in ssalbedo.dat C Dipole extinction in ssdipolext.dat C Quadrupole extinction in ssquadrext.dat C C There is also a possibility to output elements of the scattering matrix C which are stored in the array CX. The lines are commented out C (see the lines preceded by ccx string, but you can activate them if you want. 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,ILCS,i,j,l,ij1,ikl,ieps,istep integer NOUT,NSTEP,NFIN,NMAT real*8 TOL COMPLEX*16 ZEPS0,CCEPS,CSEPS,ZARTAN character*1 ync logical ynperfcon,yntest,ynbrug external ZARTAN c Parameters: C ::: number of the output unit PARAMETER (NOUT=35) c number of spherical harmonics used PARAMETER (lmx=50) * If sphere is coated, ync=y, otherwise ync=n parameter (ync='n') * * ynperfcon=.true. if core is a perfect conductor, otherwise * ynperfcon=.false. * PARAMETER (ynperfcon=.false.) * * yntest=.true. only for debugging, otherwise * ynperfcon=.false. * PARAMETER (yntest=.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/sphere radius' c (If lcs.ne.1, program is singular for rff=0. and 1.) - use homogeneous c sphere instead!!! c PARAMETER (rff=0.95d0) c background dielectric constant PARAMETER (ZEPS0=1.D0**2) c sphere (core) dielectric constant (depending whether lcs=1 or lcs=2) PARAMETER (CCEPS=(1.45D0,0.d0)**2) c PARAMETER (CCEPS=(-10.84D0,0.762d0)) C >>> SPHERE (OUTER SHELL SCATTERER) PERMITTIVITY <<< * n(silica)=1.45 <---> EPS(1)=2.1025D0 * n(ZnS)=2. <---> EPS(1)=4.D0 PARAMETER (CSEPS=(1.8d0,0.d0)**2) * * 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 ZnS c NMAT=9 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=66 * PARAMETER (NFIN=53) * C ::: relative error allowed for the TCS. If the convergence * within TOL is not reached, program issues warning PARAMETER (TOL=1.d-6) * c If ynbrug=.true., performs Bruggeman approximation for ZEPS1. Otherwise c ynbrug=false. parameter (ynbrug=.false.) * ****************************************************************** c Declarations: REAL*8 RMF(lcs),rff(lcs),RMUF,ff,filfrac !ff for Bruggemann; filfrac for ZnS real*8 xs,acs,tcs,tsc,xcs1,xcs2,xcs3,pi real*8 rsnm,lambda real*8 enw,xstep real*8 omf(NFIN),omxf,reepsz,plasma,omxp real*8 delo,omega0,omega complex*16 ceps1(NFIN),ZEPS1,Z1,Z2 COMPLEX*16 ci,cx(2,lmx),zeps(lcs+1),cqeps(2) COMPLEX*16 RX(2),SG(2),zs COMPLEX*16 KMT(2,lmx),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) *--------------------------------------------------------------- c Data: DATA PI/3.141592653589793d0/ DATA ci/(0.d0,1.d0)/ 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-------------------------------------------------------------------- * 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: C C C RMUF = 1 .... CLOSED PACKED SC (LAT=0) C --> f_{max}= 0.52359878 C C RMUF = 1./SQRT(2.) .... CLOSED PACKED FCC (LAT=1) C RMUF = 0.70710678118654746 --> f_{max}= 0.74048049 C C RMUF = SQRT(3.)/2. .... CLOSED PACKED BC (LAT=2) C RMUF = 0.8660254037844386 --> f_{max}= 0.68017476 c c FCC LATTICE c f=0.05 (--) c RMUF=0.2879411911484862d0 c f=0.06 c RMUF=0.3059831741945870d0 c f=0.07 c RMUF=0.3221166265075572d0 c f=0.08 c RMUF=0.3367780601921260d0 c f=0.09 c RMUF=0.3502632974822208d0 c f=0.1 (0f) c RMUF=0.3627831678597810d0 c f=0.12 c RMUF=0.3855146420814098d0 c f=0.15 (--) c RMUF=0.4152830592077074d0 c f=0.2 (f) c RMUF=0.4570781497340833d0 c f=0.25 (--) c RMUF=0.4923725109213483d0 c f=0.3 (--) c RMUF=0.5232238679605294d0 c f=0.35 (--) c RMUF=0.5508116833525640d0 c f=0.36 (vATL) c RMUF=0.5560083268891277d0 c f=0.4 (ff) c RMUF=0.5758823822969722d0 c f=0.45 (--) c RMUF=0.5989418136982620d0 c f=0.5 (--) c RMUF=0.6203504908994001d0 c f=0.55 (--) c RMUF=0.6403754763690468d0 c f=0.58 (2w) c RMUF=0.6518131631368212d0 c f=0.6 (3f) c RMUF=0.6592207650508868d0 c f=0.62 (3w) c RMUF=0.6664655293794289d0 c f=0.64 (4f) c RMUF=0.6735561203842518d0 c f=0.65 (--) c RMUF=0.6770461107710686d0 c f=0.66 (4w) c RMUF=0.6805004874579641d0 c f=0.68 (5f) c RMUF=0.6873059437985093d0 c f=0.7 (6f) c RMUF=0.6939792343839249d0 c f=0.72 (7f) c RMUF=0.7005265949644416d0 c close packed : c RMUF=1.d0/DSQRT(2.D0) c c DIAMOND LATTICE c c diamond: (close packed) c PARAMETER (RMUF=0.4330127018922193d0) c RMUF=1.d0 rmf(lcs)=rmuf *oooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooooo * Scanning over frequency interval: *------------------------------------------------- write(6,*)'Read sphere radius rsnm in nm' read(5,*) rsnm cc rsnm=63.3d0 cc write(6,*)'Read the sphere (core) dielectric constant' cc read(5,*) zeps1 * n(silica)=1.45 <---> ZEPS(1)=2.1025D0 * n(ZnS)=2. <---> ZEPS(1)=4.D0 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 * 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 c rff(1)=0.75d0 rmf(ikl)=rff(ikl)*rmuf 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 in the lth-sphere layer diel. const. for l=',ikl read (5,*) zeps(ikl) end if END IF enddo !IKL.GE.2 * end if !lcs.ge.2 * write(6,*)'Read initial (maximal) lambda (in the vacuum) in nm' read(5,*) lambda cc lambda=500d0 * * 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*rsnm*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 rsnm is the sphere radius (in nm) and lambda is the wavelengths (in nm) c in the vacuum: omega=2.d0*pi*rsnm/(lambda*rmuf) 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,*)'Size parameter x=2*pi*rsnm*n_0/lambda=', xs * write(6,*)'Scan down to wavelength (in nm)' read(5,*) enw c enw=500 write(6,*)'Scanning step in nm' read(5,*) xstep c xstep=5 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/(enw*rmuf) enw=enw-omega0 delo=enw/dble(nstep) end if * -------------------------------- * output initial statements 11 OPEN(UNIT=NOUT,FILE='sssccs.dat') rewind(NOUT) WRITE(NOUT,*)'#Scattering cs for a single (coated) sphere' WRITE(NOUT,*)'#(cross sections normalized per sphere effective & surface S=pi*rsnm**2)' write(nout,*) write(nout,*)'#Sphere radius in nm=', rsnm write(nout,*)'#Sphere radius rmuf in relative units=', rmuf write(nout,*) if (.not.ynperfcon) then write(nout,*)'#Material number =',NMAT else write(nout,*)'#Perfectly conducting sphere (core)' end if * if (ynbrug) write(6,*)'Bruggeman approximation performed' if (ynbrug) write(nout,*)'#Bruggeman approximation performed' * if (ync .eq.'n') then write(nout,*)'#Homogeneous sphere' else if (ync.eq.'y') then write(nout,*)'#coated sphere' end if write(nout,*) write(nout,*)'#host dielectric constant=', zeps(lcs+1) if (nmat.ge.1) write(nout,*)'#Dispersive layer number=', ilcs if ((lcs.eq.1).and.(nmat.eq.0)) & write(nout,*)'#sphere diel. const.=', zeps(1) if (lcs.ge.2) then write(nout,*)'#core radius/sphere radius =',rff(1) if ((ilcs.ne.1).or.(nmat.eq.0)) & write(nout,*)'#sphere core diel. const.=', zeps(1) if ((ilcs.ne.lcs).or.(nmat.eq.0)) & write(nout,*)'#coating diel. const.=',zeps(lcs) C--------/---------/---------/---------/---------/---------/---------/-- end if WRITE(NOUT,*)'#r_s/lambda_0, sg_sc in columns' write(nout,*) OPEN(UNIT=NOUT+1,FILE='ssext.dat') rewind(NOUT+1) WRITE(NOUT+1,*)'#Extinction for a single (coated) sphere' WRITE(NOUT+1,*)'#(cross sections normalized per sphere effective & surface S=pi*rsnm**2)' write(nout+1,*) write(nout+1,*)'#Sphere radius in nm=', rsnm write(nout+1,*) if (.not.ynperfcon) then write(nout+1,*)'#Material number =',NMAT else write(nout+1,*)'#Perfectly conducting sphere (core)' end if if (ynbrug) write(nout+1,*)'#Bruggeman approximation performed' if (ync .eq.'n') then write(nout+1,*)'#Homogeneous sphere' else if (ync.eq.'y') then write(nout+1,*)'#coated sphere' end if write(nout+1,*) write(nout+1,*)'#host dielectric constant=', zeps(lcs+1) if (nmat.ge.1) write(nout+1,*)'#Dispersive layer number=', ilcs if ((lcs.eq.1).and.(nmat.eq.0)) & write(nout+1,*)'#sphere diel. const.=', zeps(1) if (lcs.ge.2) then write(nout+1,*)'#core radius/sphere radius =',rff(1) if ((ilcs.ne.1).or.(nmat.eq.0)) & write(nout+1,*)'#sphere core diel. const.=', zeps(1) if ((ilcs.ne.lcs).or.(nmat.eq.0)) & write(nout+1,*)'#coating diel. const.=',zeps(lcs) C--------/---------/---------/---------/---------/---------/---------/-- end if WRITE(NOUT+1,*)'#r_s/lambda_0, sg_tot in columns' write(nout+1,*) OPEN(UNIT=NOUT+2,FILE='ssabs.dat') rewind(NOUT+2) WRITE(NOUT+2,*)'#Absorption for a single (coated) sphere' WRITE(NOUT+2,*)'#(cross sections normalized per sphere effective & surface S=pi*rsnm**2)' write(nout+2,*) write(nout+2,*)'#Sphere radius in nm=', rsnm write(nout+2,*) if (.not.ynperfcon) then write(nout+2,*)'#Material number =',NMAT else write(nout+2,*)'#Perfectly conducting sphere (core)' end if if (ynbrug) write(nout+2,*)'#Bruggeman approximation performed' if (ync .eq.'n') then write(nout+2,*)'#Homogeneous sphere' else if (ync.eq.'y') then write(nout+2,*)'#coated sphere' end if write(nout+2,*) write(nout+2,*)'# host dielectric constant=', zeps(lcs+1) if (nmat.ge.1) write(nout+2,*)'#Dispersive layer number=', ilcs if ((lcs.eq.1).and.(nmat.eq.0)) & write(nout+2,*)'#sphere diel. const.=', zeps(1) if (lcs.ge.2) then write(nout+2,*)'#core radius/sphere radius =',rff(1) if ((ilcs.ne.1).or.(nmat.eq.0)) & write(nout+2,*)'#sphere core diel. const.=', zeps(1) if ((ilcs.ne.lcs).or.(nmat.eq.0)) & write(nout+2,*)'#coating diel. const.=',zeps(lcs) C--------/---------/---------/---------/---------/---------/---------/-- end if WRITE(NOUT+2,*)'#r_s/lambda_0 and sg_abs in columns' write(nout+2,*) OPEN(UNIT=NOUT+5,FILE='ssalbedo.dat') rewind(NOUT+5) WRITE(NOUT+5,*)'#Albedo for a single (coated) sphere' write(nout+5,*)'#Sphere radius in nm=', rsnm write(nout+5,*) if (.not.ynperfcon) then write(nout+5,*)'#Material number =',NMAT else write(nout+5,*)'#Perfectly conducting sphere (core)' end if * if (ynbrug) write(nout+5,*)'#Bruggeman approximation performed' * if (ync .eq.'n') then write(nout+5,*)'#Homogeneous sphere' else if (ync.eq.'y') then write(nout+5,*)'#coated sphere' end if write(nout+5,*) write(nout+5,*)'#host dielectric constant=', zeps(lcs+1) if (nmat.ge.1) write(nout+5,*)'#Dispersive layer number=', ilcs if ((lcs.eq.1).and.(nmat.eq.0)) & write(nout+5,*)'#sphere diel. const.=', zeps(1) if (lcs.ge.2) then write(nout+5,*)'#core radius/sphere radius =',rff(1) if ((ilcs.ne.1).or.(nmat.eq.0)) & write(nout+5,*)'#sphere core diel. const.=', zeps(1) if ((ilcs.ne.lcs).or.(nmat.eq.0)) & write(nout+5,*)'#coating diel. const.=',zeps(lcs) C--------/---------/---------/---------/---------/---------/---------/-- end if WRITE(NOUT+5,*)'#r_s/lambda_0, albedo, tcs-(acs+tsc) in columns' write(nout+5,*) OPEN(UNIT=NOUT+6,FILE='ssdipolext.dat') rewind(NOUT+6) WRITE(NOUT+6,*)'#Dipole ext. for a single (coated) sphere' write(nout+6,*)'#Sphere radius in nm=', rsnm write(nout+6,*) if (.not.ynperfcon) then write(nout+6,*)'#Material number =',NMAT else write(nout+6,*)'#Perfectly conducting sphere (core)' end if if (ynbrug) write(nout+6,*)'#Bruggeman approximation performed' if (ync .eq.'n') then write(nout+6,*)'#Homogeneous sphere' else if (ync.eq.'y') then write(nout+6,*)'#coated sphere' end if write(nout+6,*) write(nout+6,*)'#host dielectric constant=', zeps(lcs+1) if (nmat.ge.1) write(nout+6,*)'#Dispersive layer number=', ilcs if ((lcs.eq.1).and.(nmat.eq.0)) & write(nout+6,*)'#sphere diel. const.=', zeps(1) if (lcs.ge.2) then write(nout+6,*)'#core radius/sphere radius =',rff(1) if ((ilcs.ne.1).or.(nmat.eq.0)) & write(nout+6,*)'#sphere core diel. const.=', zeps(1) if ((ilcs.ne.lcs).or.(nmat.eq.0)) & write(nout+6,*)'#coating diel. const.=',zeps(lcs) C--------/---------/---------/---------/---------/---------/---------/-- end if WRITE(NOUT+6,*)'#r_s/lambda_0 and dipole extinction columns' write(nout+6,*) C >>> OPEN(UNIT=NOUT+9,FILE='ssdipolscs.dat') rewind(NOUT+9) WRITE(NOUT+9,*)'#Dipole scs for a single (coated) sphere' write(nout+9,*)'#Sphere radius in nm=', rsnm write(NOUT+9,*) if (.not.ynperfcon) then write(NOUT+9,*)'#Material number =',NMAT else write(NOUT+9,*)'#Perfectly conducting sphere (core)' end if if (ynbrug) write(NOUT+9,*)'#Bruggeman approximation performed' if (ync .eq.'n') then write(NOUT+9,*)'#Homogeneous sphere' else if (ync.eq.'y') then write(NOUT+9,*)'#coated sphere' end if write(NOUT+9,*) write(NOUT+9,*)'#host dielectric constant=', zeps(lcs+1) if (nmat.ge.1) write(NOUT+9,*)'#Dispersive layer number=', ilcs if ((lcs.eq.1).and.(nmat.eq.0)) & write(NOUT+9,*)'#sphere diel. const.=', zeps(1) if (lcs.ge.2) then write(NOUT+9,*)'#core radius/sphere radius =',rff(1) if ((ilcs.ne.1).or.(nmat.eq.0)) & write(NOUT+9,*)'#sphere core diel. const.=', zeps(1) if ((ilcs.ne.lcs).or.(nmat.eq.0)) & write(NOUT+9,*)'#coating diel. const.=',zeps(lcs) C--------/---------/---------/---------/---------/---------/---------/-- end if WRITE(NOUT+9,*)'#r_s/lambda_0 and dipole scs columns' write(NOUT+9,*) C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+7,FILE='ssquadrext.dat') rewind(NOUT+7) WRITE(NOUT+7,*)'#Quadrupole ext. for a single (coated) sphere' write(nout+7,*)'#Sphere radius in nm=', rsnm write(nout+7,*) if (.not.ynperfcon) then write(nout+7,*)'#Material number =',NMAT else write(nout+7,*)'#Perfectly conducting sphere (core)' end if if (ynbrug) write(nout+7,*)'#Bruggeman approximation performed' if (ync .eq.'n') then write(nout+7,*)'#Homogeneous sphere' else if (ync.eq.'y') then write(nout+7,*)'#coated sphere' end if write(nout+7,*) write(nout+7,*)'#host dielectric constant=', zeps(lcs+1) if (nmat.ge.1) write(nout+7,*)'#Dispersive layer number=', ilcs if ((lcs.eq.1).and.(nmat.eq.0)) & write(nout+7,*)'#sphere diel. const.=', zeps(1) if (lcs.ge.2) then write(nout+7,*)'#core radius/sphere radius =',rff(1) if ((ilcs.ne.1).or.(nmat.eq.0)) & write(nout+7,*)'#sphere core diel. const.=', zeps(1) if ((ilcs.ne.lcs).or.(nmat.eq.0)) & write(nout+7,*)'#coating diel. const.=',zeps(lcs) C--------/---------/---------/---------/---------/---------/---------/-- end if WRITE(NOUT+7,*)'#r_s/lambda_0, quadrupole extinction in columns' write(nout+7,*) C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+10,FILE='ssquadrscs.dat') rewind(NOUT+10) WRITE(NOUT+10,*)'#Quadrupole scs for a single (coated) sphere' write(nout+10,*)'#Sphere radius in nm=', rsnm write(NOUT+10,*) if (.not.ynperfcon) then write(NOUT+10,*)'#Material number =',NMAT else write(NOUT+10,*)'#Perfectly conducting sphere (core)' end if if (ynbrug) write(NOUT+10,*)'#Bruggeman approximation performed' if (ync .eq.'n') then write(NOUT+10,*)'#Homogeneous sphere' else if (ync.eq.'y') then write(NOUT+10,*)'#coated sphere' end if write(NOUT+10,*) write(NOUT+10,*)'#host dielectric constant=', zeps(lcs+1) if (nmat.ge.1) write(NOUT+10,*)'#Dispersive layer number=', ilcs if ((lcs.eq.1).and.(nmat.eq.0)) & write(NOUT+10,*)'#sphere diel. const.=', zeps(1) if (lcs.ge.2) then write(NOUT+10,*)'#core radius/sphere radius =',rff(1) if ((ilcs.ne.1).or.(nmat.eq.0)) & write(NOUT+10,*)'#sphere core diel. const.=', zeps(1) if ((ilcs.ne.lcs).or.(nmat.eq.0)) & write(NOUT+10,*)'#coating diel. const.=',zeps(lcs) C--------/---------/---------/---------/---------/---------/---------/-- end if WRITE(NOUT+10,*)'#r_s/lambda_0, quadrupole scs in columns' write(NOUT+10,*) C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+8,FILE='ssoctupext.dat') rewind(NOUT+8) WRITE(NOUT+8,*)'#Octupole ext. for a single (coated) sphere' write(nout+8,*)'#Sphere radius in nm=', rsnm write(NOUT+8,*) if (.not.ynperfcon) then write(NOUT+8,*)'#Material number =',NMAT else write(NOUT+8,*)'#Perfectly conducting sphere (core)' end if if (ynbrug) write(NOUT+8,*)'#Bruggeman approximation performed' if (ync .eq.'n') then write(NOUT+8,*)'#Homogeneous sphere' else if (ync.eq.'y') then write(NOUT+8,*)'#coated sphere' end if write(NOUT+8,*) write(NOUT+8,*)'#host dielectric constant=', zeps(lcs+1) if (nmat.ge.1) write(NOUT+8,*)'#Dispersive layer number=', ilcs if ((lcs.eq.1).and.(nmat.eq.0)) & write(NOUT+8,*)'#sphere diel. const.=', zeps(1) if (lcs.ge.2) then write(NOUT+8,*)'#core radius/sphere radius =',rff(1) if ((ilcs.ne.1).or.(nmat.eq.0)) & write(NOUT+8,*)'#sphere core diel. const.=', zeps(1) if ((ilcs.ne.lcs).or.(nmat.eq.0)) & write(NOUT+8,*)'#coating diel. const.=',zeps(lcs) C--------/---------/---------/---------/---------/---------/---------/-- end if WRITE(NOUT+8,*)'#r_s/lambda_0, octupole extinction in columns' write(NOUT+8,*) C--------/---------/---------/---------/---------/---------/---------/-- OPEN(UNIT=NOUT+11,FILE='ssoctupscs.dat') rewind(NOUT+11) WRITE(NOUT+11,*)'#Octupole scs for a single (coated) sphere' write(nout+11,*)'#Sphere radius in nm=', rsnm write(NOUT+11,*) if (.not.ynperfcon) then write(NOUT+11,*)'#Material number =',NMAT else write(NOUT+11,*)'#Perfectly conducting sphere (core)' end if if (ynbrug) write(NOUT+11,*)'#Bruggeman approximation performed' if (ync .eq.'n') then write(NOUT+11,*)'#Homogeneous sphere' else if (ync.eq.'y') then write(NOUT+11,*)'#coated sphere' end if write(NOUT+11,*) write(NOUT+11,*)'#host dielectric constant=', zeps(lcs+1) if (nmat.ge.1) write(NOUT+11,*)'#Dispersive layer number=', ilcs if ((lcs.eq.1).and.(nmat.eq.0)) & write(NOUT+11,*)'#sphere diel. const.=', zeps(1) if (lcs.ge.2) then write(NOUT+11,*)'#core radius/sphere radius =',rff(1) if ((ilcs.ne.1).or.(nmat.eq.0)) & write(NOUT+11,*)'#sphere core diel. const.=', zeps(1) if ((ilcs.ne.lcs).or.(nmat.eq.0)) & write(NOUT+11,*)'#coating diel. const.=',zeps(lcs) C--------/---------/---------/---------/---------/---------/---------/-- end if WRITE(NOUT+11,*)'#r_s/lambda_0, octupole scs in columns' write(NOUT+11,*) if (yntest) then OPEN(UNIT=NOUT+15,FILE='tr2diag.dat') rewind(NOUT+15) end if ***************************** 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/(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') !Aluminium 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) end if ! material constant reading ********************* * begin main scanning loop: 2 do 200 istep=1,nstep+1 omega=omega0 + dble(istep-1)*delo lambda=2.d0*pi*rsnm/(omega*rmuf) cc xs=RMUF*omega*dble(sqrt(zeps0)) !Equiv. size parameter if ((nmat.eq.0).or.(ynperfcon)) 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. reepsz=2.d0*pi*RSNM/(323.83d0*rmuf) IF (NMAT.EQ.1) THEN !Material decision IF - Drude metal plasma=reepsz omxp=plasma/omega zeps1=1.d0-omxp**2/(1.d0+ci*plasma/(144.d0*omega)) go to 5 * ELSE IF (nmat.eq.4) then !Material decision IF - ZnS * filfrac=0.62d0 ! filfrac of ZnS in ZnS core call znsrefind(LAMBDA,FILFRAC,zeps1) go to 5 * ELSE IF (NMAT.EQ.2) THEN !Material decision IF - Ag c >>> real material data: !silver * lambda_z=323.83d0 * lambda_p=164.d0 * When real material data are used, * reepsz differs from plasma!!! The plasma wavelength is * calculated below: plasma=reepsz*7.2d0/3.8291d0 * security trap - remainder (not optimized!) omxf=omega/reepsz if (omxf.gt.omf(1)) then write(6,*)'Calculation of has to stop with' write(6,*)' OMF(1)' write(6,*)' OMXF=', omxf stop end if if (omxf.lt.omf(nfin)) then omxp=plasma/omega zeps1=1.d0-omxp**2/(1.d0+ci*plasma/(144.d0*omega)) * damping coefficient for silver is plasma/144 where plasma is different from * the Re eps zero crossing at 3.8291 eV according to Palik!!! go to 5 else if (omxf.eq.omf(1)) then zeps1=ceps1(1) go to 5 else do ieps=2,nfin * data file ordered with the increased wavelength * omxf increases in the loop and is oriented opposite to the data file if (omxf.gt.omf(ieps)) then ! linear interpolation zeps1=ceps1(ieps)+(omxf-omf(ieps))*(ceps1(ieps-1)-ceps1(ieps)) 1 /(omf(ieps-1)-omf(ieps)) go to 5 end if enddo end if * ELSE IF ((NMAT.GE.3).or.((nmat.ge.5).and.(nmat.le.7))) then !Material decision IF !Au,Cu,Al,Pt c >>> * data file ordered with the decreased wavelength * omega increases in the loop and is oriented along the data file * if ( (omega.lt.omf(1)).or.(omega.gt.omf(nfin)) ) then cc write(6,*)'Material data not available for this wavelength' cc stop * call sordalc(NMAT,lambda,ZEPS1) go to 5 * end if * if (omega.eq.omf(nfin)) then zeps1=ceps1(nfin) go to 5 else do ieps=1,nfin-1 if (omega.lt.omf(ieps+1)) then ! linear interpolation zeps1=ceps1(ieps)+(omega-omf(ieps))*(ceps1(ieps+1)-ceps1(ieps)) 1 /(omf(ieps+1)-omf(ieps)) go to 5 end if enddo end if END IF ! END of Material decision IF * The end of reading real data according to Palik's book *_____________________________________ * activate Bruggeman: 5 if (ynbrug) then ff=0.8d0 z1 = (3.d0*ff-1.d0)*zeps1+(2.d0 - 3.d0*ff)*zeps0 z2 = sqrt(z1*z1 + 8.d0*zeps1*zeps0) * if (IMAG(z2).GE.0.0) then zeps1= (z1 + z2)/4.d0 else zeps1= (z1 - z2)/4.d0 end if end if 7 zeps(ilcs)=zeps1 *______________________________________ 8 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)) SG(2)=omega*CQEPS(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)) SG(1)=omega*CQEPS(1) * CQEPS(2)=SQRT(ZEPS(j+1)) SG(2)=omega*CQEPS(2) * DO 25 I=1,2 * RX(I)=SG(I)*RMF(j) c 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)= UL(1,L) tt1(1,2,l,1)= VL(1,L) tt1(2,1,l,1)= DRUL(1,L) tt1(2,2,l,1)= DRVL(1,L) * tt2(1,1,l,1)= sg(2)*DRVL(2,L) tt2(1,2,l,1)= - sg(2)*VL(2,L) tt2(2,1,l,1)= - sg(2)*DRUL(2,L) tt2(2,2,l,1)= sg(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=-tan \eta_l below * KMT(1,L)=bm(l)/am(l) != -tan \eta_l KMT(2,L)=be(l)/ae(l) != -tan \eta_l * 40 CONTINUE * write(6,*)'kmt(1,1)=', kmt(1,1) write(6,*)'kmt(2,1)=', kmt(2,1) * write(6,*)'kmt(1,2)=', kmt(1,2) write(6,*)'kmt(2,2)=', kmt(2,2) * C******************************************************************** * Calculation of the cross sections tsc=0.d0 tcs=0.d0 acs=0.d0 xs=omega*rmuf*dble(sqrt(zeps0)) * write(6,*)'Size parameter x=2*pi*rsnm*n_0/lambda=', xs DO 60 j=1,lmx * >>> extracting phase-shifts from KMT *>>> magnetic part: zs=-kmt(1,j) alpha(j)=zartan(zs) zs=ci*sin(alpha(j))*exp(ci*alpha(j)) cc write(6,*)'j, i*sin(eta_m)*exp(i*eta_m)=', j, zs if ((yntest).and.(j.le.5)) then write(nout+15,*) 'j, sin^2\eta_M', j, (sin(alpha(j)))**2 end if * * Under normal circumstances, Im. part of a phase shift \geq 0 !!!) * if(imag(alpha(j)).lt.0.d0) then write(6,*)'imag(alpha(j))=',imag(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) zs=ci*sin(beta(j))*exp(ci*beta(j)) cc write(6,*)'j, i*sin(eta_e)*exp(i*eta_e)=', j, zs if ((yntest).and.(j.le.5)) then write(nout+15,*) 'j, sin^2\eta_E', j, (sin(beta(j)))**2 end if * * Under normal circumstances, Im. part of a phase shift \geq 0 !!!) * if(imag(beta(j)).lt.0.d0) then write(6,*)'imag(alpha(j))=',imag(beta(j)) ,' is negative' write(6,*)'omega, j=', omega, j pause end if c write(6,*)'beta(j)=', beta(j) * total elastic scattering cross section: * [for example, formula (2.137) of Newton] xcs1=dble(2*j+1)*(2.d0+exp(-4.d0*imag(alpha(j)))- 1 2.d0*exp(-2.d0*imag(alpha(j)))*cos(2.d0*dble(alpha(j)))+ 2 exp(-4.d0*imag(beta(j)))- 4 2.d0*exp(-2.d0*imag(beta(j)))*cos(2.d0*dble(beta(j)))) tsc=tsc+xcs1 * absorption cross section: * [for example, formula (2.138) of Newton] xcs2=dble(2*j+1)*(2.d0 - 1 exp(-4.d0*imag(alpha(j))) - exp(-4.d0*imag(beta(j)))) acs=acs+xcs2 * total scattering/extinction cross section: * [for example, formula (2.136) of Newton] xcs3=dble(2*j+1)*(2.d0 - 1 exp(-2.d0*imag(alpha(j)))*cos(2.d0*dble(alpha(j)))- 2 exp(-2.d0*imag(beta(j)))*cos(2.d0*dble(beta(j)))) *<<< if (j.eq.1) then WRITE(NOUT+6,*)rsnm/lambda, xcs3/xs**2 WRITE(NOUT+9,*)rsnm/lambda, xcs1/(2.d0*xs**2) end if if (j.eq.2) then WRITE(NOUT+7,*)rsnm/lambda, xcs3/xs**2 WRITE(NOUT+10,*)rsnm/lambda, xcs1/(2.d0*xs**2) end if if (j.eq.3) then WRITE(NOUT+8,*)rsnm/lambda, xcs3/xs**2 WRITE(NOUT+11,*)rsnm/lambda, xcs1/(2.d0*xs**2) end if *<<< tcs=tcs+xcs3 c write(6,*)'j, tsc, acs, tcs =', tsc, acs, tcs ccx DO 50 i=1,2 ccx cx(i,j)=(dcmplx(1.d0,0.d0) - ci*kmt(i,j))/ ccx 1 (dcmplx(1.d0,0.d0) + ci*kmt(i,j)) ccx 50 CONTINUE 60 CONTINUE xcs1=xcs1/tsc if (acs.ne.0.) then xcs2=xcs2/acs else xcs2=0.d0 end if xcs3=xcs3/tcs if (xcs1.gt.tol ) then write(6,*)'Warning - convergence worse than', xcs1 go to 200 else if (xcs2.gt.tol ) then write(6,*)'Warning - convergence worse than', xcs2 go to 200 else if (xcs3.gt.tol ) then write(6,*)'Warning - convergence worse than', xcs3 go to 200 end if c write(6,*)'bare tcs=', tcs c write(6,*)'Convergence within', xcs *-------------------------------------------- * Scattering coefficients (cross sections normalized * per sphere effective surface S=pi*rsnm**2. * According, for example, formulae (2.136-8) of Newton: * Elastic scattering coefficient(2.137) of Newton: tsc=tsc/(2.d0*xs**2) * Absorption coefficient (2.138) of Newton: acs=acs/(2.d0*xs**2) * QEXT=(2*PI/k**2) \sum_{AL} \sin^2\eta_{AL} * Total scattering coefficient (2.136) of Newton: tcs=tcs/xs**2 *-------------------------------------------- * Note factor 2 here for the vector waves instead of 4 * for the scalar ones * write(nout,*) rsnm/lambda,tsc write(nout+1,*) rsnm/lambda,tcs write(nout+2,*) rsnm/lambda,acs write(nout+5,*) rsnm/lambda,tsc/tcs,tcs -(tsc+acs) if(nstep.gt.10) go to 200 c write(6,*) 'istep=', istep write(6,*) 'lambda=',lambda write(6,*) write(6,*)'Scattering coefficient=' write(6,*)'tsc=', tsc write(6,*) write(6,*)'Absorption coefficient=' write(6,*)'acs=', acs write(6,*) write(6,*)'Extinction coefficient=' write(6,*)'tcs=', tcs write(6,*) write(6,*)'Scattering identity tcs=tsc+acs' write(6,*)'tcs-(tsc+acs)=', tcs -(tsc+acs) write(6,*) 200 continue close(nout) close(nout+1) close(nout+2) close(nout+5) close(nout+6) close(nout+7) close(nout+8) close(nout+9) close(nout+10) close(nout+11) close(nout+15) if (ync .eq.'n') then write(6,*)'Homogeneous sphere' else if (ync.eq.'y') then write(6,*)'coated sphere' end if * write(6,*)'Sphere parameters:' write(6,*) write(6,*)'radius =', rsnm if (ync.eq.'n') write(6,*)'sphere diel. constant=', zeps(1) write(6,*)'background dielectric constant ZEPS0=', zeps0 if (ync.eq.'y') write(6,*)'core diel. constant=', zeps(1) if (ync.eq.'y') write(6,*)'coating diel. constant=', zeps(lcs) if (ync.eq.'y') write(6,*)'core radius/sphere radius =', & rff(1) write(6,*) write(6,*)'Extinction versus wavelength in ssext.dat' write(6,*)'Scattering cs versus wavelength in sssccs.dat' write(6,*)'Absorption versus wavelength in ssabs.dat' write(6,*)'Albedo versus wavelength in ssalbedo.dat' write(6,*)'Dipole extinction in ssdipolext.dat' write(6,*)'Quadrupole extinction in ssquadrext.dat' write(6,*)'Octupole extinction in ssoctupext.dat' write(6,*)'Dipole scs in ssdipolscs.dat' write(6,*)'Quadrupole scs in ssquadrscs.dat' write(6,*)'Octupole scs in ssoctupscs.dat' *--------/---------/---------/---------/---------/---------/---------/-- end C (C) Copr. 1/1999 Alexander Moroz