Mode_sover_2D_me.zip_Micro–Nano_couple166_光场分布_模场分布_纳米光场

clear all; %% 物理常数及结构参数 ybxn0=8.854187817e-12; %真空介电常数 miu0=4*pi*1e-7; %真空磁导率 n1=3.34; %衬底折射率 n2=3.44; %脊折射率 n3=1.00; %空气折射率 a=2e-8; %相邻网格间隔20nm lambda=1.55e-6; %波长 k=2*pi/lambda; %波矢 L=399; %网格列数 W=409; %网格行数 size1=4*2+6*(L-2);%计算第1/w行Hx对应的系数矩阵非零元素个数 size2=5*2+7*(L-2);%计算第其他行Hx时，每一行对应的系数矩阵非零元素个数 size3=4*2+5*(L-2);%计算第1/w行Hy对应的系数矩阵非零元素个数 size4=6*2+7*(L-2);%计算第其他行Hy时，每一行对应的系数矩阵非零元素个数 num=2*(size1+size3)+(W-2)*(size2+size4);%系数矩阵非零元素总个数 %% ebsn=zeros(W+1,L+1); %波导的介电常数矩阵 for i=1:L+1 for j=1:200 ebsn(j,i)=n1^2; end end for i=1:L+1 for j=201:210 ebsn(j,i)=n2^2; end end for i=1:L+1 for j=211:W+1 ebsn(j,i)=n3^2; end end for i=151:250 for j=211:265 ebsn(j,i)=n2^2; end end %% A_width=zeros(1,num); % 稀疏矩阵的行下标向量（宽） B_height=zeros(1,num); % 稀疏矩阵的列下标向量（高） S_value=zeros(1,num); % 稀疏矩阵的值向量 %% 确定稀疏矩阵行列下标向量以及值 %% Hx for i=1:W for j=1:L %% 下边界 if i==1 if j==1 A_width(1:2)=1:2; A_width(3)=L+1; A_width(4)=2; B_height(1:3)=1; B_height(4)=W*L+1; end if j>1&&j<L num_1=5+(j-2)*6; A_width(num_1:num_1+2)=j-1:j+1; A_width(num_1+3:num_1+5)=[L*i+j,j-1,j+1]; B_height(num_1:num_1+3)=j; B_height(num_1+4:num_1+5)=W*L+j; end if j==L A_width(size1-3:size1-2)=j-1:j; A_width(size1-1:size1)=[L+j,j-1]; B_height(size1-3:size1-1)=j; B_height(size1)=W*L+L*(i-1)+j; end ebsn1=ebsn(i+1,j); ebsn2=ebsn(i,j); ebsn3=ebsn(i,j+1); ebsn4=ebsn(i+1,j+1); xxw=1/a^2; xxe=1/a^2; xxn=1/a^2*(ebsn3/(ebsn3+ebsn4)+ebsn2/(ebsn1+ebsn2)); xxs=1/a^2*(ebsn4/(ebsn3+ebsn4)+ebsn1/(ebsn1+ebsn2)); xxp=(ebsn1*ebsn2/(ebsn1+ebsn2)+ebsn3*ebsn4/(ebsn3+ebsn4))*k^2-4/a^2; xye=1/(2*a^2)*((ebsn1-ebsn2)/(ebsn1+ebsn2)+(ebsn4-ebsn3)/(ebsn3+ebsn4)); xyw=-xye; switch j case 1 %左下角点 S_value(1:4)=[xxp xxe xxn xye]; case L %右下角点 S_value(size1-3:size1)=[xxw xxp xxn xyw]; otherwise S_value(6*j-7:6*j-2)=[xxw xxp xxe xxn xyw xye]; end end %% 中间 if i>1&&i<W if j==1 A_width(size1+size2*(i-2)+1)=L*(i-2)+j; A_width(size1+size2*(i-2)+2:size1+size2*(i-2)+3)=L*(i-1)+j:L*(i-1)+j+1; A_width(size1+size2*(i-2)+4:size1+size2*(i-2)+5)=[i*L+j,L*(i-1)+j+1]; B_height(size1+size2*(i-2)+1:size1+size2*(i-2)+4)=L*(i-1)+j; B_height(size1+size2*(i-2)+5)=W*L+L*(i-1)+j; end if j>1&&j<L num_1=size1+size2*(i-2)+7*(j-2)+6; A_width(num_1)=L*(i-2)+j; A_width(num_1+1:num_1+3)=L*(i-1)+j-1:L*(i-1)+j+1; A_width(num_1+4:num_1+6)=[L*i+j,L*(i-1)+j-1,L*(i-1)+j+1]; B_height(num_1:num_1+4)=L*(i-1)+j; B_height(num_1+5:num_1+6)=W*L+L*(i-1)+j; end if j==L num_1=size1+size2*(i-2)+7*(j-2)+6; A_width(num_1)=L*(i-2)+j; A_width(num_1+1:num_1+2)=L*(i-1)+j-1:L*(i-1)+j; A_width(num_1+3:num_1+4)=[L*i+j,L*(i-1)+j-1]; B_height(num_1:num_1+3)=L*(i-1)+j; B_height(num_1+4)=W*L+L*(i-1)+j; end ebsn1=ebsn(i+1,j); ebsn2=ebsn(i,j); ebsn3=ebsn(i,j+1); ebsn4=ebsn(i+1,j+1); xxw=1/a^2; xxe=1/a^2; xxn=1/a^2*(ebsn3/(ebsn3+ebsn4)+ebsn2/(ebsn1+ebsn2)); xxs=1/a^2*(ebsn4/(ebsn3+ebsn4)+ebsn1/(ebsn1+ebsn2)); xxp=(ebsn1*ebsn2/(ebsn1+ebsn2)+ebsn3*ebsn4/(ebsn3+ebsn4))*k^2-4/a^2; xye=1/(2*a^2)*((ebsn1-ebsn2)/(ebsn1+ebsn2)+(ebsn4-ebsn3)/(ebsn3+ebsn4)); xyw=-xye; aa=size1+size2*(i-2); if(j==1) %左边点 S_value(aa+1:aa+5)=[xxs xxp xxe xxn xye]; else if(j==L) %右边点 S_value(aa+size2-4:aa+size2)=[xxs xxw xxp xxn xyw]; else S_value(aa+7*j-8:aa+7*j-2)=[xxs xxw xxp xxe xxn xyw xye]; end end end %% 上边界 if i==W if j==1 num_1=size1+size2*(i-2)+1; A_width(num_1)=L*(i-2)+j; A_width(num_1+1:num_1+2)=L*(i-1)+j:L*(i-1)+j+1; A_width(num_1+3)=L*(i-1)+j+1; B_height(num_1:num_1+2)=L*(i-1)+j; B_height(num_1+3)=W*L+L*(i-1)+j; end if j>1&&j<L num_1=size1+size2*(i-2)+6*(j-2)+5; A_width(num_1)=L*(i-2)+j; A_width(num_1+1:num_1+3)=L*(i-1)+j-1:L*(i-1)+j+1; A_width(num_1+4:num_1+5)=[L*(i-1)+j-1,L*(i-1)+j+1]; B_height(num_1:num_1+3)=L*(i-1)+j; B_height(num_1+4:num_1+5)=W*L+L*(i-1)+j; end if j==L num_1=size1+size2*(i-2)+6*(j-2)+5; A_width(num_1)=L*(i-2)+j; A_width(num_1+1:num_1+2)=L*(i-1)+j-1:L*(i-1)+j; A_width(num_1+3)=L*(i-1)+j-1; B_height(num_1:num_1+2)=L*(i-1)+j; B_height(num_1+3)=W*L+L*(i-1)+j; end ebsn1=ebsn(i+1,j); ebsn2=ebsn(i,j); ebsn3=ebsn(i,j+1); ebsn4=ebsn(i+1,j+1); xxw=1/a^2; xxe=1/a^2; xxn=1/a^2*(ebsn3/(ebsn3+ebsn4)+ebsn2/(ebsn1+ebsn2)); xxs=1/a^2*(ebsn4/(ebsn3+ebsn4)+ebsn1/(ebsn1+ebsn2)); xxp=(ebsn1*ebsn2/(ebsn1+ebsn2)+ebsn3*ebsn4/(ebsn3+ebsn4))*k^2-4/a^2; xye=1/(2*a^2)*((ebsn1-ebsn2)/(ebsn1+ebsn2)+(ebsn4-ebsn3)/(ebsn3+ebsn4)); xyw=-xye; aa1=size1+(i-2)*size2; if(j==1) %左上角点 S_value(aa1+1:aa1+4)=[xxs xxp xxe xye]; else if(j==L) %右上角点 S_value(aa1+size1-3:aa1+size1)=[xxs xxw xxp xyw]; else S_value(aa1+6*j-7:aa1+6*j-2)=[xxs xxw xxp xxe xyw xye]; end end end end end %% Hy for i=1:W for j=1:L %% 上边界 if i==1 if j==1 num_1=2*size1+(W-2)*size2+1; A_width(num_1:num_1+3)=[W*L+(i-1)*L+j,W*L+(i-1)*L+j+1,W*L+i*L+j,W*L+i*L+j]; B_height(num_1:num_1+2)=W*L+(i-1)*L+j; B_height(num_1+3)=(i-1)*L+j; end if j>1&&j<L num_1=2*size1+(W-2)*size2+5*(j-2)+5; A_width(num_1:num_1+4)=[W*L+(i-1)*L+j-1:W*L+(i-1)*L+j+1,W*L+i*L+j,W*L+i*L+j]; B_height(num_1:num_1+3)=W*L+(i-1)*L+j; B_height(num_1+4)=(i-1)*L+j; end if j==L num_1=2*size1+(W-2)*size2+size3-3; A_width(num_1:num_1+3)=[W*L+(i-1)*L+j-1:W*L+(i-1)*L+j,W*L+i*L+j,W*L+i*L+j]; B_height(num_1:num_1+2)=W*L+(i-1)*L+j; B_height(num_1+3)=(i-1)*L+j; end ebsn1=ebsn(i+1,j); ebsn2=ebsn(i,j); ebsn3=ebsn(i,j+1); ebsn4=ebsn(i+1,j+1); yyw=1/a^2*(eb