ģ

模拟二维空间中TM波的传播过程

%模拟TM波传播的MATLAB源程序
%function main
fprintf('please wait patiently ,some files are being created ...\n');
v=3*10^8;
delta_t=2.5*10^(-11);
delta_s=2*v*delta_t;
steps_t=100;
cicle_N=2*steps_t+1;
sourse_location_x=51;
sourse_location_y=51;
grids_M=100;
grids_N=100;
Ez_M=grids_M+1;
Ez_N=grids_N+1;
Hx_N=grids_N+1;
Hy_M=grids_M+1;
Hy_N=grids_N;

%初始化
temp=zeros(Ez_M,Ez_N);
temp(sourse_location_x,sourse_location_y)=sourse(0);
Ez_data_getfile;
fid=fopen(Ez_data_file{1},'w');
fprintf(fid,'%e',temp);
fclose(fid);
%value_Ez(1)=temp;           %cicle=1

Hx_data_getfile;
fid=fopen(Hx_data_file{1},'w');
fprintf(fid,'%e',zeros(Hx_M,Hx_N));
fclose(fid);
%value_Hx(1)=zeros(Hx_M,Hx_N);     %cicle=2

Hy_data_getfile;
fid=fopen(Hy_data_file{1},'w');
fprintf(fid,'%e',zeros(Hy_M,Hy_N));
fclose(fid);
%value_Hy(1)=zeros(Hy_M,Hy_N);     %cicle=2

for cicle_n=3:cicle_N
    
%*****************计算Ez********************
%下面将所有的电场分量都更新
if mod(cicle_n,2)==1
    temp_Ez_this=zeros(Ez_M,Ez_N);         %初始化temp_Ez_this;
    fid=fopen(Ez_data_file{(cicle_n+1)/2-1},'r');
    temp_Ez_last=fscanf(fid,'%e',[Ez_M,Ez_N]);
    fclose(fid);
    fid=fopen(Hx_data_file{(cicle_n-1)/2},'r');
    temp_Hx_last=fscanf(fid,'%e',[Hx_M,Hx_N]);
    fclose(fid);
    fid=fopen(Hy_data_file{(cicle_n-1)/2},'r');
    temp_Hy_last=fscanf(fid,'%e',[Hy_M,Hy_N]);
    fclose(fid);
    for i=1:Ez_M
        for j=1:Ez_N
            %第一行的处理
            if i==1
                if j==1
                    if cicle_n==3
                        temp_Ez_this(1,1)=0;
                    else
                        f_rad=sqrt(1-1/(50*sqrt(2)));
                        alfa_11=pi/4;
                        fid=fopen(Ez_data_file{(cicle_n+1)/2-2},'r');
                        temp_Ez_last2=fscanf(fid,'%e',[Ez_M,Ez_N]);
                        fclose(fid);
%                        temp_Ez_last2=value_Ez((cicle_n+1)/2-2);
temp_Ez_this(1,1)=f_rad*((1-sin(alfa_11))*(1-cos(alfa_11))*temp_Ez_last2(1,1)+(1-sin(alfa_11))*cos(alfa_11)*tenp_Ez_last2(1,2)+sin(alfa_11)*(1-cos(alfa_11))*temp_Ez_last2(2,1)+sin(alfa_11)*cos(alfa_11)*temp_Ez_last2(2,2));
                    end
% 第一行，第二列到第Ez_N-1列的处理
                elseif j>=2 & j<=(Ez_N-1)
% 先求出第二行第j列的值
temp_Ez_this(2,j)=temp_Ez_last(2,j)+delta_t*v/delta_s*(temp_Hy_last(2,j)-temp_Hy_last(2,(j-1))+temp_Hx_last(2,j)-temp_Hx_last(1,j));
                    if cicle_n==3
A1=temp_Ez_this(2,j);
A2=temp_Ez_this(2,j)-2*temp_Ez_last(2,j)-2*temp_Ez_last(1,j);
A3=temp_Ez_last(2,(j+1))-2*temp_Ez_last(2,j)+temp_Ez_last(2,(j-1))+temp_Ez_last(1,(j+1))-2*temp_Ez_last(1,j)+temp_Ez_this(1,(j-1));
                    else
                        fid=fopen(Ez_data_file{(cicle_n+1)/2-2},'r');
                        temp_Ez_last2=fscanf(fid,'%e',[Ez_M,Ez_N]);
                        fclose(fid);
                        temp_Ez_last2=value_Ez((cicle_n+1)/2-2);
                        A1=temp_Ez_this(2,j)-temp_Ez_last2(2,j)+temp_Ez_last2(1,j);
A2=temp_Ez_this(2,j)-2*temp_Ez_last(2,j)-2*temp_Ez_last(1,j)+temp_Ez_last2(2,j)+temp_Ez_last2(1,j);
A3=temp_Ez_last(2,(j+1))-2*temp_Ez_last(2,j)+temp_Ez_last(2,(j-1))+temp_Ez_last(1,(j+1))-2*temp_Ez_last(1,j)+temp_Ez_this(1,(j-1));
                    end
                    p0=-1;
                    p2=1/2;
temp_Ez_this(1,j)=1/(1*p0+1)*A1-2*p0/(2*p0+1)*A2-p2/(4*p0+2)*A3;
%第一行第Ez_N列的处理
                elseif j==Ez_N
                    if cicle_n==3
                        temp_Ez_this(1,Ez_N)=0;
                    else
                        f_rad=sqrt(1-1/(50*sqrt(2)));
                        alfa_11=pi/4;
                        fid=fopen(Ez_data_file{(cicle_n+1)/2-2},'r');
                        temp_Ez_last2=fscanf(fid,'%e',[Ez_M,Ez_N]);
                        fclose(fid);
temp_Ez_this(1,Ez_N)=f_rad*((1-sin(alfa_11))*(1-cos(alfa_11))*temp_Ez_last2(1,Ez_N)+(1-sin(alfa_11))*cos(alfa_11)*tenp_Ez_last2(1,(Ez_N-1))+sin(alfa_11)*(1-cos(alfa_11))*temp_Ez_last2(2,Ez_N)+sin(alfa_11)*cos(alfa_11)*temp_Ez_last2(2,(Ez_N-1)));
                    end
                end
%第二行到第Ez_M-1行的处理
            elseif i>=2&&i<=(Ez_M-1)
%第二行到第Ez_M-1行，第一列的处理
                if j==1
%先求出第i行，第二列的值
temp_Ez_this(i,2)=temp_Ez_last(i,2)+delta_t*v/delta_s*(temp_Hy_last(i,2)-temp_Hy_last(i,1)+temp_Hx_last(i,2)-temp_Hx_last((i-1),2));
                   if cicle_n==3
                       A1=temp_Ez_this(i,2);
A2=temp_Ez_this(i,2)-2*temp_Ez_last(i,2)-2*temp_Ez_last(i,1);
A3=temp_Ez_last((i-1),2)-2*temp_Ez_last(i,2)+temp_Ez_last((i+1),2)+temp_Ez_last((i-1),1)-2*temp_Ez_last(i,1)+temp_Ez_this((i+1),1);
                   else
                       fid=fopen(Ez_data_file{(cicle_n+1)/2-2},'r');
                        temp_Ez_last2=fscanf(fid,'%e',[Ez_M,Ez_N]);
                        fclose(fid);
                        A1=temp_Ez_this(i,2)-temp_Ez_last2(i,2)+temp_Ez_last2(i,1);
A2=temp_Ez_this(i,2)-2*temp_Ez_last(i,2)-2*temp_Ez_last(i,1)+temp_Ez_last2(i,2)+temp_Ez_last2(i,1);
A3=temp_Ez_last((i-1),2)-2*temp_Ez_last(i,2)+temp_Ez_last((i+1),2)+temp_Ez_last((i-1),1)-2*temp_Ez_last(i,1)+temp_Ez_this((i+1),1);
                    end
                    p0=-1;
                    p2=1/2;
temp_Ez_this(i,1)=1/(2*p0+1)*A1-2*p0/(2*p0+1)*A2-p2/(4*p0+2)*A3;
%第三行到第Ez_M-2行，第三列到第Ez_N-2列的处理
                elseif j>=3&&j<=Ez_N-2
%如果是激励源的位置
                    if i==sourse_location_x&&j==sourse_location_y
                        temp_Ez_this(i,j)=sourse((cicle_n-1)/2);
                    else
temp_Ez_this(i,j)=temp_Ez_last(i,j)+delta_t*v/delta_s*(temp_Hy_last(i,j)-temp_Hy_last(i,(j-1))+temp_Hx_last(i,j)-temp_Hx_last((i-1),j));
                    end
%第二行到第Ez_M-1行，第Ez_N列的处理
                elseif j==Ez_N
%先求出第二行到第Ez_M-1行，第Ez_N-1列的值
temp_Ez_this(i,Ez_N-1)=temp_Ez_last(i,Ez_N-1)+delta_t*v/delta_s*(temp_Hy_last(i,Ez_N-1)-temp_Hy_last(i,Ez_N-2)+temp_Hx_last(i,Ez_N-1)-temp_Hx_last((i-1),Ez_N-1));
                    if cicle_n==3
                        A1=temp_Ez_this(i,Ez_N-1)
A2=temp_Ez_this(i,Ez_N-1)-2*temp_Ez_last(i,Ez_N-1)-2*temp_Ez_last(i,Ez_N);
A3=temp_Ez_last((i-1),Ez_N-1)-2*temp_Ez_last(i,Ez_N-1)+temp_Ez_last((i+1),Ez_N-1)+temp_Ez_last((i-1),Ez_N)-2*temp_Ez_last(i,Ez_N)+temp_Ez_this((i+1),Ez_N);
                    else
                        fid=fopen(Ez_data_file{(cicle_n+1)/2-2},'r');
                        temp_Ez_last2=fscanf(fid,'%e',[Ez_M,Ez_N]);
                        fclose(fid);
                        A1=temp_Ez_this(i,Ez_N-1)-temp_Ez_last2(i,Ez_N-1)+temp_Ez_last2(i,Ez_N);
A2=temp_Ez_this(i,Ez_N-1)-2*temp_Ez_last(i,Ez_N-1)-2*temp_Ez_last(i,Ez_N-1)+temp_Ez_last2(i,Ez_N-1)+temp_Ez_last2(i,Ez_N);
A3=temp_Ez_last((i-1),Ez_N-1)-2*temp_Ez_last(i,Ez_N-1)+temp_Ez_last((i+1),Ez_N-1)+temp_Ez_last((i-1),Ez_N)-2*temp_Ez_last(i,Ez_N)+temp_Ez_this((i+1),Ez_N);
                    end
                    p0=-1;
                    p2=1/2;
temp_Ez_this(i,Ez_N)=1/(2*p0+1)*A1-2*p0/(2*p0+1)*A2-p2/(4*p0+2)*A3;
                end
%第Ez_M行的处理
            elseif i==Ez_M
%第Ez_M行第一列的处理
                if j==1
                    if cicle_n==3
                        temp_Ez_this(Ez_M,1)=0;
                    else
                        f_rad=sqrt(1-1/(50*sqrt(2)));
                        alfa_11=pi/4;
                        fid=fopen(Ez_data_file{(cicle_n+1)/2-2},'r');
                        temp_Ez_last2=fscanf(fid,'%e',[Ez_M,Ez_N]);
                        fclose(fid);
temp_Ez_this(Ez_M,1)=f_rad*((1-sin(alfa_11))*(1-cos(alfa_11))*temp_Ez_last2(Ez_M,1)+(1-sin(alfa_11))*cos(alfa_11)*tenp_Ez_last2(Ez_M,2)+sin(alfa_11)*(1-cos(alfa_11))*temp_Ez_last2((Ez_M-1),1)+sin(alfa_11)*cos(alfa_11)*temp_Ez_last2((Ez_M-1),2));
                    end
%第Ez_M行第二列到第Ez_N-1列的处理
                elseif j>=2&&j<=(Ez_N-1)
                    if cicle_n==3
                        A1=temp_Ez_this(Ez_M-1,j);
A2=temp_Ez_this(Ez_M-1,j)-2*temp_Ez_last(Ez_M-1,j)-2*temp_Ez_last(Ez_M,j);
A3=temp_Ez_last(Ez_M-1,j+1)-2*temp_Ez_last(Ez_M-1,j)+temp_Ez_last(Ez_M-1,j-1)+temp_Ez_last(Ez_M,j+1)-2*temp_Ez_last(Ez_M,j)+temp_Ez_this(Ez_M,j-1);
                    else
                        fid=fopen(Ez_data_file{(cicle_n+1)/2-2},'r');
                        temp_Ez_last2=fscanf(fid,'%e',[Ez_M,Ez_N]);
                        fclose(fid);
                        A1=temp_Ez_this(Ez_M-1,j)-temp_Ez_last2(Ez_M-1,j)+temp_Ez_last2(Ez_M,j);
A2=temp_Ez_this(Ez_M-1,j)-2*temp_Ez_last(Ez_M-1,j)-2*temp_Ez_last(Ez_M,j)+temp_Ez_last2(Ez_M-1,j)+temp_Ez_last2(Ez_M,j);
A3=temp_Ez_last(Ez_M-1,j+1)-2*temp_Ez_last(Ez_M-1,j)+temp_Ez_last(Ez_M-1,j-1)+temp_Ez_last(Ez_M+1,j)-2*temp_Ez_last(Ez_M,j)+temp_Ez_this(Ez_M,j-1);
                    end
                    p0=-1;
                    p2=1/2;
temp_Ez_this(Ez_M,j)=1/(2*p0+1)*A1-2*p0/(2*p0+1)*A2-p2/(4*p0+2)*A3;
%第Ez_M行第Ez_N列的处理
                elseif j==Ez_N
                    if cicle_n==3
                        temp_Ez_this(Ez_M,Ez_N)=0;
                    else
                        f_rad=sqrt(1-1/(50*sqrt(2)));
                        alfa_11=pi/4;
                        fid=fopen(Ez_data_file{(cicle_n+1)/2-2},'r');
                        temp_Ez_last2=fscanf(fid,'%e',[Ez_M,Ez_N]);
                        fclose(fid);
temp_Ez_this(Ez_M,Ez_N)=f_rad*((1-sin(alfa_11))*(1-cos(alfa_11))*temp_Ez_last2(Ez_M,Ez_N)+(1-sin(alfa_11))*cos(alfa_11)*tenp_Ez_last2(Ez_M,(Ez_N-1))+sin(alfa_11)*(1-cos(alfa_11))*temp_Ez_last2((Ez_M-1),Ez_N)+sin(alfa_11)*cos(alfa_11)*temp_Ez_last2((Ez_M-1),(Ez_N-1)));
                    end
                end
            end
        end
    end
    fid=fopen(Ez_data_file{(cicle_n+1)/2},'w');
    fprintf(fid,'%e',temp_Ez_this);
    fclose(fid);
%      value_Ez((cicle_n+1)/2)=temp_Ez_this;

%*********************************** 计算 Hx 和 Hy****************************
%下面将所有的磁场分量都更新
   else
       %初始化 Hx
       temp_Hx_this=zeros(Hx_M,Hx_N);
       fid=fopen(Hx_data_file{cicle_n/2-1},'r');
       temp_Hx_last=fscanf(fid,'%e',[Hx_M,Hx_N]);
       fclose(fid);
       %初始化 Hy
       temp_Hy_this=zeros(Hy_M,Hy_N);
       fid=fopen(Hy_data_file{cicle_n/2-1},'r');
       temp_Hy_last=fscanf(fid,'%e',[Hy_M,Hy_N]);
       fclose(fid);
       fid=fopen(Ez_data_file{cicle_n/2},'r');
       temp_Ez_last=fscanf(fid,'%e',[Ez_M,Ez_N]);
       fclose(fid);
    
       for i=1:Hx_M
           for j=1:Hx_M

temp_Hx_this(i,j)=temp_Hx_last(i,j)+delta_t*v/delta_s*(temp_Ez_last(i+1,j)-temp_Ez_last(i,j));
           end
       end
       fid=fopen(Hx_data_file{cicle_n/2},'w');
       fprintf(fid,'%e',temp_Hx_this);
       fclose(fid);
       for i=1:Hy_M
           for j=1:Hy_N
temp_Hy_this(i,j)=temp_Hy_last(i,j)+delta_t*v/delta_s*(temp_Ez_last(i,j+1)-temp_Ez_last(i,j));
           end
       end
       fid=fopen(Hy_data_file{cicle_n/2},'w');
       fprintf(fid,'%e',temp_Hy_this);
       fclose(fid);
   end
end
fprintf('done!\n');