【DDSCAT——离散偶极近似仿真程序06】模型设计Matlab代码：多六棱柱、多圆柱核六棱柱壳

最新推荐文章于 2021-06-07 18:16:46 发布

XD_Yangf

最新推荐文章于 2021-06-07 18:16:46 发布

阅读量2k

点赞数 1

分类专栏： DDSCAT

本文链接：https://blog.csdn.net/XD_Yangf/article/details/102468583

版权

DDSCAT 专栏收录该内容

16 篇文章 38 订阅

订阅专栏

继续设计文献里面的模型

一、设计多个六棱柱

1）Matlab代码

%%All hex are oriented such that the incident light
%%is along the x axis and strikes the circular face of the hex
R = [20, 20, 20]; %Hex center to side distance
H = [10, 10, 10]; %Hex height
hex = 3; %number of cylinders (including ring)
offX = zeros(1,hex);
offY = zeros(1,hex);
offZ = zeros(1,hex);
%   offX = [0,-40,40]; %User defined offsets
    offY = [18,-18,-18]; %user defined offsets
    offZ = [0,21,-21]; %user defined offsets
%%Calculation variables for output file documantation
X = H;
Y = R*1.5;
Z = R*1.5;
Xd = X;
Yd = Y*2;
Zd = Z*2;
Out = [0 0 0 0 0 0 0];
d = 1; %dipole #
di = 1; %inner dipole #
for off = 1:hex;
    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))
                L1 = (abs(sqrt(3)/2*(j-offY(1,off))+1/2*(i-offZ(1,off))));
                L2 = (abs(sqrt(3)/2*(j-offY(1,off))-1/2*(i-offZ(1,off))));
                L3 = ((i-offZ(1,off))^2 + (j-offY(1,off))^2)^.5;
                if L1 <= R(1,off) && L2 <= R(1,off) && abs(i-offZ(1,off))<= R(1,off);
                    Out = [Out; d k j i 1 1 1];
                    d = d + 1;
                end
            end
        end
    end
end
Out(1,:) = [];
d = d - 1;
di = di - 1;
fid = fopen('output.txt', 'w');
fprintf(fid, ' >TARCEL: Hexagons; AX,AY,AZ= %7.4f %7.4f %7.4f', Xd, Yd, Zd);
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('%dHex%dR%dH',hex,R(1,1),H(1,1));
folderpath = strcat(foldername,'/');
filename = sprintf('%dHex%dR%dH.txt',hex,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文件

先放一个会报错的ddscat.par文件。

下面是在命令行中运行的报错结果。

报错的原因暂时还没有找到，但是可以确定的是，其中X=3.00000D-01中的X表示的是入射光的波长。因为我改变过波长的值，这个值也会跟着改变，继续报错。小范围的波长改动依然会继续报错，但是大范围的波长改动，会导致程序一直迭代下去，不会终止，且error误差值很大。

3）近场图

报这个错误最大的问题是，在ddpostprocess后处理之后，虽然产生了vtr文件，但是在vislt里面无法读取模型的强度值，即无法可视化。

4）正确的ddscat.par文件

在尝试了很多次的改动之后，终于改出了可以正常运行不报错、可以可视化的配置文件。结果其实很简单，就是把介电常数的种类配置为两个。虽然这个纯六棱柱模型是由一种物质构成的，但是我们设置有两种配置文件没有什么影响，只要shape文件里面只出现1，不出现2就行。结果，这样一改，反而结果正确了。

5）成功的近场图

6）模型的3D散点图

之后尽量给大家加上模型的一个三维散点图，除了在电厂图可视化的情况下认识模型，也同时在偶极子构成概念上了解模型。

二、设计多个圆柱核六棱柱壳模型

1）Matlab代码

%%All hex are oriented such that the incident light
%%is along the x axis and strikes the circular face of the hex
R = [20, 20, 20]; %Hex center to side distance
H = [20, 20, 20]; %Hex height
Ri = [10, 10, 10]; %Hex core radius
hex = 3; %number of cylinders (including ring)
offX = zeros(1,hex);
offY = zeros(1,hex);
offZ = zeros(1,hex);
% offX = [0,-40,40];
offY = [18,-18,-18];
offZ = [0,21,-21];
%% calculation parameters used for output documentation file
X = H;
Y = R*1.5;
Z = R*1.5;
Xd = X; %diameter
Yd = Y*2;
Zd = Z*2;
Xs = Ri*2;
Ys = Ri*2;
Zs = Ri*2;
Out = [0 0 0 0 0 0 0];
d = 1; %dipole #
di = 1; %inner dipole #
for off = 1:hex;
    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))
                L1 = (abs(sqrt(3)/2*(j-offY(1,off))+1/2*(i-offZ(1,off))));
                L2 = (abs(sqrt(3)/2*(j-offY(1,off))-1/2*(i-offZ(1,off))));
                L3 = ((i-offZ(1,off))^2 + (j-offY(1,off))^2)^.5;
                if L1 <= R(1,off) && L2 <= R(1,off) && abs(i-offZ(1,off))<= R(1,off);
                    if L3 <= Ri;
                        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: Hexegons with Cores; 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('%dHexCore%dR%dr%dH',hex,R(1,1),Ri(1,1),H(1,1));
folderpath = strcat(foldername,'/');
filename = sprintf('%dHexCore%dR%dr%dH.txt',hex,R(1,1),Ri(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文件

' ========== 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

这里的配置其实和上一个模型的正确配置文件的参数是一模一样的，不过这次是没有出现问题的。

3）近场图