%------ 均匀三线互联系统程序 -------%
% 模型图为6-5(P168)
% 单根线(导带为PEC,且无限薄)时的结构参数
% d=20um,w=80um,s=60um,L=4um,
% 导带上方8倍d,导带左右方8倍s
clear
clc
close all
% 记录程序运行时间
tic
% 脉冲最高频率(为了增大时间步长,将高斯脉冲在时域压缩)及相关参数
f_max=100*10^9;
c=3*10^8;
lamda_min=c/f_max;
% 介质参量
mju=4*pi*10^(-7);
epslong_0=8.45*10^(-12);
epslong_r=3.2;
% 脉冲周期参数
tao=1/(2*f_max);
t0=3*tao;
% 模型尺寸设置
dim_w=80*10^(-6);
dim_s=60*10^(-6);
dim_d=20*10^(-6);
dim_L=4*10^(-3);
% 网格尺寸
deltay=dim_d/4;
deltax=dim_s/3;
deltaz=dim_L/20;
% 空间步数
Nx=(3*dim_w+2*dim_s+8*dim_s+8*dim_s)/deltax;
Ny=(dim_d+dim_d*8)/deltay;
Nz=dim_L/deltaz;
% 介质参量的矩阵存储
epslong_x=zeros(Nx,Ny+1,Nz+1);
epslong_y=zeros(Nx+1,Ny,Nz+1);
epslong_z=zeros(Nx+1,Ny+1,Nz);
Ny_fenjiemian=dim_d/deltay;
epslong_x(1:Nx,1:Ny_fenjiemian,1:Nz+1)=epslong_r*epslong_0;
epslong_x(1:Nx,Ny_fenjiemian+1,1:Nz+1)=(epslong_r+1)/2*epslong_0;
epslong_x(1:Nx,Ny_fenjiemian+2:Ny+1,1:Nz+1)=epslong_0;
epslong_y(1:Nx+1,1:Ny_fenjiemian,1:Nz+1)=epslong_r*epslong_0;
epslong_y(1:Nx+1,Ny_fenjiemian+1:Ny,1:Nz+1)=epslong_0;
epslong_z(1:Nx+1,1:Ny_fenjiemian,1:Nz)=epslong_r*epslong_0;
epslong_z(1:Nx+1,Ny_fenjiemian+1,1:Nz)=(epslong_r+1)/2*epslong_0;
epslong_z(1:Nx+1,Ny_fenjiemian+2:Ny+1,1:Nz)=epslong_0;
% 初始化所有场分量
Hx(1:Nx+1,1:Ny,1:Nz)=0;
Hy(1:Nx,1:Ny+1,1:Nz)=0;
Hz(1:Nx,1:Ny,1:Nz+1)=0;
Ex(1:Nx,1:Ny+1,1:Nz+1)=0;
Ey(1:Nx+1,1:Ny,1:Nz+1)=0;
Ez(1:Nx+1,1:Ny+1,1:Nz)=0;
% 时间步长
deltat=1/c/sqrt(1/deltax^2+1/deltay^2+1/deltaz^2)/1.1;
% 时间步数
Nt=5400;
% 吸收边界条件存储数组
Ex_qian=[];
Ey_qian=[];
Ez_qian=[];
% 观测面
Nz_qian_observe=4;
Nz_hou_observe=16;
% 记录观测点电压电流
v=0;
I=0;
V_t_11=zeros(Nt,1);
V_t_12=zeros(Nt,1);
V_t_13=zeros(Nt,1);
V_t_14=zeros(Nt,1);
V_t_15=zeros(Nt,1);
V_t_16=zeros(Nt,1);
I_t_11=zeros(Nt,1);
I_t_12=zeros(Nt,1);
I_t_13=zeros(Nt,1);
V_t_21=zeros(Nt,1);
V_t_22=zeros(Nt,1);
V_t_23=zeros(Nt,1);
V_t_24=zeros(Nt,1);
V_t_25=zeros(Nt,1);
V_t_26=zeros(Nt,1);
I_t_21=zeros(Nt,1);
I_t_22=zeros(Nt,1);
I_t_23=zeros(Nt,1);
% 脉冲完全送进端口时间步数
Nt_pulse=2100;
% 高斯脉冲
J_zz=[];
j_zz=0;
for t=1:Nt
j_zz=exp(-((t*deltat-t0)^2)/(tao^2));
J_zz=[J_zz,j_zz];
end
J_zz;
Nt*deltat*10^9
%高斯电流密度作图
figure(1),
plot(1*deltat*10^(9):1*deltat*10^(9):Nt*deltat*10^(9),J_zz);
title('时域高斯脉冲信号');
xlabel('时间/(ns)');
ylabel('电流幅度/(1)');
grid;
% 开始循环迭代
% 第一次激励
for t=1:Nt
Hx(1:Nx+1,1:Ny,1:Nz)=Hx(1:Nx+1,1:Ny,1:Nz)+deltat/mju*((Ez(1:Nx+1,1:Ny,1:Nz)-Ez(1:Nx+1,2:Ny+1,1:Nz))/deltay+...
(Ey(1:Nx+1,1:Ny,2:Nz+1)-Ey(1:Nx+1,1:Ny,1:Nz))/deltaz);
Hy(1:Nx,1:Ny+1,1:Nz)=Hy(1:Nx,1:Ny+1,1:Nz)+deltat/mju*((Ez(2:Nx+1,1:Ny+1,1:Nz)-Ez(1:Nx,1:Ny+1,1:Nz))/deltax+...
(Ex(1:Nx,1:Ny+1,1:Nz)-Ex(1:Nx,1:Ny+1,2:Nz+1))/deltaz);
Hz(1:Nx,1:Ny,1:Nz+1)=Hz(1:Nx,1:Ny,1:Nz+1)+deltat/mju*((Ex(1:Nx,2:Ny+1,1:Nz+1)-Ex(1:Nx,1:Ny,1:Nz+1))/deltay+...
(Ey(1:Nx,1:Ny,1:Nz+1)-Ey(2:Nx+1,1:Ny,1:Nz+1))/deltax);
% 吸收边界条件存储
Ex_qian(1:Nx,1:Ny+1,1:Nz+1)=Ex(1:Nx,1:Ny+1,1:Nz+1);
Ey_qian(1:Nx+1,1:Ny,1:Nz+1)=Ey(1:Nx+1,1:Ny,1:Nz+1);
Ez_qian(1:Nx+1,1:Ny+1,1:Nz)=Ez(1:Nx+1,1:Ny+1,1:Nz);
Ex(1:Nx,2:Ny,2:Nz)=Ex(1:Nx,2:Ny,2:Nz)+(deltat/epslong_x(1:Nx,2:Ny,2:Nz)).*((Hz(1:Nx,2:Ny,2:Nz)-Hz(1:Nx,1:Ny-1,2:Nz))/deltay+...
(Hy(1:Nx,2:Ny,1:Nz-1)-Hy(1:Nx,2:Ny,2:Nz))/deltaz);
% z=1
Ex(1:Nx,2:Ny_fenjiemian,1)=Ex_qian(1:Nx,2:Ny_fenjiemian,2)+(c/sqrt(epslong_r)*deltat-deltaz)/(c/sqrt(epslong_r)*deltat+deltaz)*...
(Ex(1:Nx,2:Ny_fenjiemian,2)-Ex(1:Nx,2:Ny_fenjiemian,1));
Ex(1:Nx,Ny_fenjiemian+1,1)=Ex_qian(1:Nx,Ny_fenjiemian+1,2)+(c/sqrt(epslong_r/2+1/2)*deltat-deltaz)/(c/sqrt(epslong_r/2+1/2)*deltat+deltaz)*...
(Ex(1:Nx,Ny_fenjiemian+1,2)-Ex(1:Nx,Ny_fenjiemian+1,1));
Ex(1:Nx,Ny_fenjiemian+2:Ny,1)=Ex_qian(1:Nx,Ny_fenjiemian+2:Ny,2)+(c*deltat-deltaz)/(c*deltat+deltaz)*...
(Ex(1:Nx,Ny_fenjiemian+2:Ny,2)-Ex(1:Nx,Ny_fenjiemian+2:Ny,1));
% z=Nz+1
Ex(1:Nx,2:Ny_fenjiemian,Nz+1)=Ex_qian(1:Nx,2:Ny_fenjiemian,Nz)+(c/sqrt(epslong_r)*deltat-deltaz)/(c/sqrt(epslong_r)*deltat+deltaz)*...
(Ex(1:Nx,2:Ny_fenjiemian,Nz)-Ex(1:Nx,2:Ny_fenjiemian,Nz+1));
Ex(1:Nx,Ny_fenjiemian+1,Nz+1)=Ex_qian(1:Nx,Ny_fenjiemian+1,Nz)+(c/sqrt(epslong_r/2+1/2)*deltat-deltaz)/(c/sqrt(epslong_r/2+1/2)*deltat+deltaz)*...
(Ex(1:Nx,Ny_fenjiemian+1,Nz)-Ex(1:Nx,Ny_fenjiemian+1,Nz+1));
Ex(1:Nx,Ny_fenjiemian+2:Ny,Nz+1)=Ex_qian(1:Nx,Ny_fenjiemian+2:Ny,Nz)+(c*deltat-deltaz)/(c*deltat+deltaz)*...
(Ex(1:Nx,Ny_fenjiemian+2:Ny,Nz)-Ex(1:Nx,Ny_fenjiemian+2:Ny,Nz+1));
% y=1
Ex(1:Nx,1,1:Nz+1)=0;
% y=Ny+1
Ex(1:Nx,Ny+1,1:Nz+1)=0;
% 导带
Ex(25:28,Ny_fenjiemian+1,1:Nz+1)=0;
Ex(32:35,Ny_fenjiemian+1,1:Nz+1)=0;
Ex(39:42,Ny_fenjiemian+1,1:Nz+1)=0;
Ey(2:Nx,1:Ny,2:Nz)=Ey(2:Nx,1:Ny,2:Nz)+(deltat/epslong_y(2:Nx,1:Ny,2:Nz)).*((Hx(2:Nx,1:Ny,2:Nz)-Hx(2:Nx,1:Ny,1:Nz-1))/deltaz+...
(Hz(1:Nx-1,1:Ny,2:Nz)-Hz(2:Nx,1:Ny,2:Nz))/deltax);
if t<=Nt_pulse
% 端口1激励
Ey(39:43,1:Ny_fenjiemian,2)=-J_zz(t);
end
% z=1
Ey(2:Nx,1:Ny_fenjiemian,1)=Ey_qian(2:Nx,1:Ny_fenjiemian,2)+...
(c/sqrt(epslong_r)*deltat-deltaz)/(c/sqrt(epslong_r)*deltat+deltaz)*(Ey(2:Nx,1:Ny_fenjiemian,2)-Ey(2:Nx,1:Ny_fenjiemian,1));
Ey(2:Nx,Ny_fenjiemian+1:Ny,1)=Ey_qian(2:Nx,Ny_fenjiemian+1:Ny,2)+...
(c*deltat-deltaz)/(c*deltat+deltaz)*(Ey(2:Nx,Ny_fenjiemian+1:Ny,2)-Ey(2:Nx,Ny_fenjiemian+1:Ny,1));
% z=Nz+1
Ey(2:Nx,1:Ny_fenjiemian,Nz+1)=Ey_qian(2:Nx,1:Ny_fenjiemian,Nz)+(c/sqrt(epslong_r)*deltat-deltaz)/(c/sqrt(epslong_r)*deltat+deltaz)*...
(Ey(2:Nx,1:Ny_fenjiemian,Nz)-Ey(2:Nx,1:Ny_fenjiemian,Nz+1));
Ey(2:Nx,Ny_fenjiemian+1:Ny,Nz+1)=Ey_qian(2:Nx,Ny_fenjiemian+1:Ny,Nz)+(c*deltat-deltaz)/(c*deltat+deltaz)*...
(Ey(2:Nx,Ny_fenjiemian+1:Ny,Nz)-Ey(2:Nx,Ny_fenjiemian+1:Ny,Nz+1));
% x=1
Ey(1,1:Ny,1:Nz+1)=0;
% x=Nx+1
Ey(Nx+1,1:Ny,1:Nz+1)=0;
Ez(2:Nx,2:Ny,1:Nz)=Ez(2:Nx,2:Ny,1:Nz)+(deltat/epslong_z(2:Nx,2:Ny,1:Nz)).*((Hy(2:Nx,2:Ny,1:Nz)-Hy(1:Nx-1,2:Ny,1:Nz))/deltax+...
(Hx(2:Nx,1:Ny-1,1:Nz)-Hx(2:Nx,2:Ny,1:Nz))/deltay);
% x=1
Ez(1,1:Ny,1:Nz)=0;
% x=Nx+1
Ez(Nx+1,1:Ny,1:Nz)=0;
% y=1
Ez(1:Nx+1,1,1:Nz)=0;
% y=Ny+1
Ez(1:Nx+1,Ny+1,1:Nz)=0;
% 导带
Ez(25:29,Ny_fenjiemian+1,1:Nz)=0;
Ez(32:36,Ny_fenjiemian+1,1:Nz)=0;
Ez(39:43,Ny_fenjiemian+1,1:Nz)=0;
% 记录观测点信息
% 端口1电压
v=0;
for y=Ny_fenjiemian:-1:1
v=v-Ey(41,y,Nz_qian_observe)*deltay;
end
V_t_11(t)=v;
% 端口2电压
v=0;
for y=Ny_fenjiemian:-1:1
v=v-Ey(34,y,Nz_qian_observe)*deltay;
end
V_t_12(t)=v;
% 端口3电压
v=0;
for y=Ny_fenjiemian:-1:1
v=v-Ey(27,y,Nz_qian_observe)*deltay;
end
V_t_13(t)=v;
% 端口4电压
v=0;
for y=Ny_fenjiemian:-1:1
v=v-Ey(41,y,Nz_hou_observe)*deltay;
end
V_t_14(t)=v;
% 端口5电压
v=0;
for y=Ny_fenjiemian:-1:1
v=v-Ey(34,y,Nz_hou_observe)*deltay;
end
V_t_15(t)=v;
% 端口6电压
v=0;
for y=Ny_fenjiemian:-1:1
v=v-Ey(27,y,Nz_hou_observe)*deltay;
end
V_t_16(t)=v;
% 端口1电流
I=0;
for x=38:44
I=I+(Hx(x,Ny_fenjiemian-1,Nz_qian_observe)-Hx(x,Ny_fenjiemian+2,Nz_qian_observe))*deltax;
end
for y=Ny_fenjiemian-1:Ny_fenjiemian+3
I=I+(Hy(43,y,Nz_qian_observe)-Hy(38,y,Nz_qian_observe))*deltay;
end
I_t_11(t)=I;
% 端口2电流
I=0;
for x=31:37
I=I+(Hx(x,Ny_fenjiemian-1,Nz_qian_observe)-Hx(x,Ny_fenjiemian+2,Nz_qian_observe))*deltax;
end
for y=Ny_fenjiemian