【DDSCAT——离散偶极近似仿真程序05】模型设计Matlab代码：多核壳球 、多核壳圆柱
简洁明了，后面及部分给出的都是模型设计的Matlab代码，这些代码都来自于一篇文章：Chin C W S. Localized Surface Plasmon Resonance with the use of Silver and Titanium Oxide Nanostructures[J]. 2011.
一、设计多个核壳球形结构
1）Matlab代码
Matlab代码生成shape文件，shape.dat会被自动放在当前路径下的一个文件夹里。
%**************************************************************************
% 希望此代码能帮助各位学习DDSCAT的网友更近一步
% 作者：XD_Yangf
% 来源：CSDN
% 原文：Chin C W S. Localized Surface Plasmon Resonance with the use of Silver and Titanium Oxide Nanostructures[J]. 2011.
% 版权声明：本文为博主原创文章，转载请附上博文链接！
% Date:2019-05-15
%**************************************************************************
clear;
clc;
%**************************************************************************
Ls = 15; %sphere radius for calculations
Li = 8; %inner sphere radius
spheres = 2; %number of spheres
offX = zeros(1,spheres);
offY = zeros(1,spheres);
offZ = zeros(1,spheres);
offX = [0,0,40,0,0,0];
offY = [-20,20,0,40,80,120];
% offZ = [0,20,-20,0,20];
%% calculation parameters for output file documentation
X = Ls; %sphere radius
Y = Ls;
Z = Ls;
Xd = X*2; %diameter
Yd = Y*2;
Zd = Z*2;
Xs = Li*2;
Ys = Li*2;
Zs = Li*2;
Out = [0 0 0 0 0 0 0];
d = 1; %dipole #
di = 1; %inner dipole #
off = 1;
for off = 1:spheres
for i = (-Z+offZ(1,off)):(Z+offZ(1,off))
for j = (-Y+offY(1,off)):(Y+offY(1,off))
for k = (-X+offX(1,off)):(X+offX(1,off))
L = ((i-offZ(1,off))^2 + (j-offY(1,off))^2 + (k-offX(1,off))^2)^.5;
if L <= Ls
if L <= Li
Out = [Out; d k j i 1 1 1];
d = d + 1;
di = di + 1;
else
Out = [Out; d k j i 2 2 2];
d = d + 1;
end
end
end
end
end
end
Out(1,:) = [];
d = d - 1;
di = di - 1;
fid = fopen('output.txt', 'w');
fprintf(fid, ' >TARCEL: concentric spheres; AX,AY,AZ= %7.4f %7.4f %7.4f BX,BY,BZ= %7.4f %7.4f %7.4f', Xd, Yd, Zd, Xs, Ys, Zs);
fprintf(fid, '\r\n %8.0f %8.0f = NAT, NIN,', d, di);
fprintf(fid, '\r\n 1.000000 0.000000 0.000000 = A_1 vector\r\n 0.000000 1.000000 0.000000 = A_2 vector');
fprintf(fid, '\r\n 1.000000 1.000000 1.000000 = lattice spacings (d_x,d_y,d_z)/d');
fprintf(fid, '\r\n -0.50000 -0.50000 -0.50000 = lattice offset x0(1-3) = (x_TF,y_TF,z_TF)/d for dipole 0 0 0');
fprintf(fid, '\r\n JA IX IY IZ ICOMP(x,y,z)\r\n');
fprintf(fid,'%7d %3d %3d %3d %1d %1d %1d\r\n', Out');
fclose(fid);
foldername = sprintf('%dSpheres%dR%dr',spheres,Ls,Li);
folderpath = strcat(foldername,'/');
filename = sprintf('%dSpheres%dR%dr.txt',spheres,Ls,Li);
mkdir(foldername);
copyfile('output.txt',filename);
copyfile('output.txt','shape.dat');
movefile(filename, folderpath);
movefile('shape.dat',folderpath); % 将shape.dat文件移动到 %dSpheres%dR%dr/ 文件夹中
% 后面的也是一些文件复制移动的操作，大家可以注释掉，需要用的话得去掉注释自行修改代码，否则会报错
% copyfile('Ag.txt',folderpath);
% copyfile('Ti.txt',folderpath);
% copyfile('ddscat.exe',folderpath);
% copyfile('ddscat2.par','ddscat.par');
% movefile('ddscat.par',folderpath);
% disp(['Sucess, file saved to: ', foldername]);2）ddscat.par文件
在这里跑一下这个双核壳球模型在300nm入射光下的消光计算，以便跟大家对比。
' ========== Parameter file for v7.3 ==================='
'**** Preliminaries ****'
'NOTORQ' = CMTORQ*6 (DOTORQ, NOTORQ) -- either do or skip torque calculations
'PBCGS2' = CMDSOL*6 (PBCGS2, PBCGST, GPBICG, QMRCCG, PETRKP) -- CCG method
'GPFAFT' = CMETHD*6 (GPFAFT, FFTMKL) -- FFT method
'GKDLDR' = CALPHA*6 (GKDLDR, LATTDR, FLTRCD) -- DDA method
'NOTBIN' = CBINFLAG (NOTBIN, ORIBIN, ALLBIN)
'**** Initial Memory Allocation ****'
150 150 150 = dimensioning allowance for target generation
'**** Target Geometry and Composition ****'
'FROM_FILE' = CSHAPE*9 shape directive
no SHPAR parameters needed
2 = NCOMP = number of dielectric materials
'../diel/SiO2' = file with refractive index 1
'../diel/Au_evap' = file with refractive index 2
'**** Additional Nearfield calculation? ****'
1 = NRFLD (=0 to skip nearfield calc., =1 to calculate nearfield E)
0.0 0.0 0.0 0.0 0.0 0.0 (fract. extens. of calc. vol. in -x,+x,-y,+y,-z,+z)
'**** Error Tolerance ****'
1.00e-5 = TOL = MAX ALLOWED (NORM OF |G>=AC|E>-ACA|X>)/(NORM OF AC|E>)
'**** Maximum number of iterations ****'
2000 = MXITER
'**** Integration cutoff parameter for PBC calculations ****'
1.00e-2 = GAMMA (1e-2 is normal, 3e-3 for greater accuracy)
'**** Angular resolution for calculation of <cos>, etc. ****'
0.5 = ETASCA (number of angles is proportional to [(3+x)/ETASCA]^2 )
'**** Vacuum wavelengths (micron) ****'
0.300 0.300 1 'LIN' = wavelengths (first,last,how many,how=LIN,INV,LOG)
'**** Refractive index of ambient medium'
1.0000 = NAMBIENT
'**** Effective Radii (micron) **** '
0.03 0.03 1 'LIN' = eff. radii (first, last, how many, how=LIN,INV,LOG)
'**** Define Incident Polarizations ****'
(0,0) (1.,0.) (0.,0.) = Polarization state e01 (k along x axis)
2 = IORTH (=1 to do only pol. state e01; =2 to also do orth. pol. state)
'**** Specify which output files to write ****'
1 = IWRKSC (=0 to suppress, =1 to write ".sca" file for each target orient.
'**** Specify Target Rotations ****'
0. 0. 1 = BETAMI, BETAMX, NBETA (beta=rotation around a1)
0. 0. 1 = THETMI, THETMX, NTHETA (theta=angle between a1 and k)
0. 0. 1 = PHIMIN, PHIMAX, NPHI (phi=rotation angle of a1 around k)
'**** Specify first IWAV, IRAD, IORI (normally 0 0 0) ****'
0 0 0 = first IWAV, first IRAD, first IORI (0 0 0 to begin fresh)
'**** Select Elements of S_ij Matrix to Print ****'
6 = NSMELTS = number of elements of S_ij to print (not more than 9)
11 12 21 22 31 41 = indices ij of elements to print
'**** Specify Scattered Directions ****'
'LFRAME' = CMDFRM (LFRAME, TFRAME for Lab Frame or Target Frame)
1 = NPLANES = number of scattering planes
0. 0. 180. 5 = phi, theta_min, theta_max (deg) for plane A
90. 0. 180. 5 = phi, theta_min, theta_max (deg) for plane B
对于这个配置文件，我对核壳的选择分别是SiO2和Au，Au_evap折射率文件是DDSCAT自带的，SiO2是我在网上找的，然后自己写到文件中的。下面是SiO2的折射率文件。
来源网站：refractiveindex.info
SiO2, (from refractiveindex.info)
1 2 3 0 0 = columns for wave, Re(n), Im(n), eps1, eps2
wave(um) Re(n) Im(n) eps1 eps2
0.9840 1.45 0.000 -40.27 2.79
0.8920 1.45 0.000 -32.04 1.93
0.8211 1.45 0.000 -25.81 1.63
0.7560 1.45 0.000 -20.61 1.27
0.7045 1.46 0.000 -16.82 1.07
0.6595 1.46 0.000 -13.65 1.04
0.6168 1.46 0.000 -10.66 1.37
0.5821 1.46 0.000 -8.11 1.66
0.5486 1.46 0.000 -5.84 2.11
0.5209 1.46 0.000 -3.95 2.58
0.4959 1.46 0.000 -2.28 3.81
0.4714 1.46 0.000 -1.70 4.84
0.4509 1.47 0.000 -1.76 5.28
0.4305 1.47 0.000 -1.69 5.65
0.4133 1.47 0.000 -1.70 5.72
0.3974 1.47 0.000 -1.65 5.74
0.3815 1.47 0.000 -1.60 5.64
0.3679 1.47 0.000 -1.40 5.61
0.3542 1.48 0.000 -1.23 5.60
0.3425 1.48 0.000 -1.31 5.54
0.3315 1.48 0.000 -1.36 5.57
0.3204 1.48 0.000 -1.23 5.85
0.3107 1.49 0.000 -1.24 5.79
0.3009 1.49 0.000 -1.23 5.78
0.2924 1.49 0.000 -1.31 5.60
0.2844 1.49 0.000 -1.33 5.49
0.2761 1.50 0.000 -1.37 5.28
0.2689 1.50 0.000 -1.35 4.98
3）近场图
为了计算快点，我就选了300 nm处的计算（这里有个规律，就是入射光波长越长，计算时间越长），并将其在300 nm处的近场图呈现了。
正视近场图。
切面近场图。
二、设计多个核壳圆柱形结构
1）Matlab代码
%**************************************************************************
% 希望此代码能帮助各位学习DDSCAT的网友更近一步
% 作者：XD_Yangf
% 来源：CSDN
% 原文：Chin C W S. Localized Surface Plasmon Resonance with the use of Silver and Titanium Oxide Nanostructures[J]. 2011.
% 版权声明：本文为博主原创文章，转载请附上博文链接！
% Date:2019-05-15
%**************************************************************************
clear;
clc;
%**************************************************************************
%%All cylinders are oriented such that the incident light
%%is along the x axis and strikes the circular face of the cylinder
R = [10, 10, 10]; %outer cylinder radius
r = [5, 5, 10]; %inner cylinder radius
H = [20, 20, 10]; %outer cylinder height
h = [10, 10, 10]; %inner cylinder height
hoff = [0, 0, 0]; %inner cylinder offset along x axis (0 means centered)
cylinders = 2; %number of cylinders
offX = zeros(1,cylinders);
offY = zeros(1,cylinders);
offZ = zeros(1,cylinders);
% offX = [0,-40,40];
offY = [-12,12];
% offZ = [0,20,-20,0,20];
X = H;
Y = R;
Z = R;
Rc = R; %cylinder radius for calculations
Xi = H/2 + hoff; %inner cylinder centered at half outer height (with offset)
Xo = Xi - h/2; %inner cylinder lower bound
Xmax = Xi + h/2;
Xd = X; %diameter
Yd = Y*2;
Zd = Z*2;
ri = r; %inner cylinder radius
Xs = ri*2;
Ys = ri*2;
Zs = ri*2;
Out = [0 0 0 0 0 0 0];
d = 1; %dipole #
di = 1; %inner dipole #
off = 1;
for off = 1:cylinders
for i = (-Z(1,off)+offZ(1,off)):(Z(1,off)+offZ(1,off))
for j = (-Y(1,off)+offY(1,off)):(Y(1,off)+offY(1,off))
for k = (0+offX(1,off)):(X(1,off)+offX(1,off))
L = ((i-offZ(1,off))^2 + (j-offY(1,off))^2)^.5;
if L <= Rc(1,off)
if L <= ri(1,off)
if k >= Xo(1,off) && k <= Xmax(1,off)
Out = [Out; d k j i 1 1 1];
d = d + 1;
di = di + 1;
else
Out = [Out; d k j i 2 2 2];
d = d + 1;
end
else
Out = [Out; d k j i 2 2 2];
d = d + 1;
end
end
end
end
end
end
Out(1,:) = [];
d = d - 1;
di = di - 1;
fid = fopen('output.txt', 'w');
fprintf(fid, ' >TARCEL: concentric cylinders; AX,AY,AZ= %7.4f %7.4f %7.4f BX,BY,BZ= %7.4f %7.4f %7.4f', Xd, Yd, Zd, Xs, Ys, Zs);
fprintf(fid, '\r\n %8.0f %8.0f = NAT, NIN,', d, di);
fprintf(fid, '\r\n 1.000000 0.000000 0.000000 = A_1 vector\r\n 0.000000 1.000000 0.000000 = A_2 vector')
fprintf(fid, '\r\n 1.000000 1.000000 1.000000 = lattice spacings (d_x,d_y,d_z)/d')
fprintf(fid, '\r\n -0.50000 -0.50000 -0.50000 = lattice offset x0(1-3) = (x_TF,y_TF,z_TF)/d for dipole 0 0 0')
fprintf(fid, '\r\n JA IX IY IZ ICOMP(x,y,z)\r\n');
fprintf(fid,'%7d %3d %3d %3d %1d %1d %1d\r\n', Out');
fclose(fid);
foldername = sprintf('%dCylinder%dR%dr%dHh',cylinders,R(1,1),r(1,1),H(1,1));
folderpath = strcat(foldername,'/');
filename = sprintf('%dCylinder%dR%dr%dHh.txt',cylinders,R(1,1),r(1,1),H(1,1));
mkdir(foldername);
copyfile('output.txt',filename);
copyfile('output.txt','shape.dat');
movefile(filename, folderpath);
movefile('shape.dat',folderpath); % 将shape.dat文件移动到 %dSpheres%dR%dr/ 文件夹中
% 后面的也是一些文件复制移动的操作，大家可以注释掉，需要用的话得去掉注释自行修改代码，否则会报错
% copyfile('Ag.txt',folderpath);
% copyfile('Ti.txt',folderpath);
% copyfile('ddscat.exe',folderpath);
% copyfile('ddscat2.par','ddscat.par');
% movefile('ddscat.par',folderpath);
% disp(['Sucess, file saved to: ', foldername]);2）ddscat.par文件
配置同上
3）近场图
正视近场图。
切面近场图。
这个圆柱有点扁，不好看出来。
后面再改改这一部分中模型的参数。
