tm_point.m

fdtd模拟点源在自由空间的传播

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%TM波的传播
%今日完成 tm波的传播

clear all; clc;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 真空常数设定
epsz=1/(4*pi*9*10^9);          %真空介电常数
mu=4*pi*10^(-7);               %真空磁导率
Z=sqrt(mu/epsz);               %真空阻抗
epsilon=1;                     %真空相对介电常数
sigma=0;                       %真空电导率

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 设定截止频率
v=3e8;                         %波速
f=0.3e9;                       %截止频率
lamda=v/f;                     %波长
k=2*pi/lamda;                  %波数

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 设定网格
IE=150;                        %x向网格数
JE=150;                        %y向网格数
ISteps=2000;                   %迭代次数
ddx=lamda/20;                  %网格尺寸
dt=ddx/(2*v);                  %时间步大小

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 设定总场
ia=30;
ib=IE-ia;                      %x向总场边界位置
ja=30;                         
jb=JE-ja;                      %y向总场边界位置

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PML层数
npml=8;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  设定脉冲源参数
spread=8;                      %脉冲宽度
t0=25;                         %脉冲高度
%x0=IE/2;
%y0=JE/2;                       %脉冲位置
x0=50;
y0=50;
pulse=0;                       %记录脉冲大小


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 计算中的变量内存分配
dz=zeros(IE,JE);               %z向电场密度，归一化后的
ez=zeros(IE,JE);               %z向电场强度
iz=zeros(IE,JE);               %迭代中间变量
hx=zeros(IE,JE);               %x向磁场变量
hy=zeros(IE,JE);               %y向磁场变量
ihx=zeros(IE,JE);              %x向hx计算的中间变量
ihy=zeros(IE,JE);              %y向hy计算的中间变量

%% PML初始化设置
gi2=ones(1,IE);                  %x向系数
gi3=ones(1,IE);
fi1=zeros(1,IE);
fi2=ones(1,IE);
fi3=ones(1,IE);

gj2=ones(1,JE);                  %y向系数
gj3=ones(1,JE);
fj1=zeros(1,JE);
fj2=ones(1,JE);
fj3=ones(1,JE);

%% PML层中阻抗渐变设置
for i=1:npml+1;
    xnum=npml-i+1;
    xxn=xnum/npml;
    xn=0.333*(xxn^3);
    gi2(i)=1/(1+xn);
    gi2(IE-i+1)=1/(1+xn);
    gi3(i)=(1-xn)/(1+xn);
    gi3(IE-i+1)=(1-xn)/(1+xn);
    
    xxn=(xnum-0.5)/npml;
    xn=0.25*(xxn^3);
    fi1(i)=xn;
    fi1(IE-i)=xn;
    fi2(i)=1/(1+xn);
    fi2(IE-i)=1/(1+xn);
    fi3(i)=(1-xn)/(1+xn);
    fi3(IE-i)=(1-xn)/(1+xn);
end;

for j=1:npml+1;
    xnum=npml-j+1;   
    xxn=xnum/npml;      
    xn=0.33*(xxn^3);   
    gj2(j)=1/(1+xn);   
    gj2(JE-j+1)=1/(1+xn);   
    gj3(j)=(1-xn)/(1+xn);   
    gj3(JE-j+1)=(1-xn)/(1+xn);   
       
    xxn=(xnum-0.5)/npml;   
    xn=0.25*(xxn^3);   
    fj1(j)=xn;   
    fj1(JE-j)=xn;   
    fj2(j)=1/(1+xn);   
    fj2(JE-j)=1/(1+xn);   
    fj3(j)=(1-xn)/(1+xn);   
    fj3(JE-j)=(1-xn)/(1+xn);  
end;



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  FDTD计算循环
for T=1:ISteps;
    disp(T);
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% 计算Dz
     for i=2:IE;
        for j=2:JE;
            dz(i,j)=gi3(i)*gj3(j)*dz(i,j)+...
                gi2(i)*gj2(j)*0.5*(hy(i,j)-hy(i-1,j)-hx(i,j)+hx(i,j-1));
        end;
    end;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% 添加信号源
    %pulse=sin(2*pi*600e6*dt*n);  %使用正弦波
    pulse=exp(-0.5*((t0-T)/spread)^2);
    dz(x0,y0)=dz(x0,y0)+pulse;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% 计算Ez     后期回波很大，可能与此有关系，在边界层的地方两者是不等的。以后要修改，添加ga，gb
    ez=dz;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%  设置Ez边界为0，作为PML辅助
    ez(1,:)=0;
    ez(IE,:)=0;
    ez(:,1)=0;
    ez(:,JE)=0;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% 计算hx
    for i=1:IE;
        for j=1:JE-1;
            curl_e=ez(i,j)-ez(i,j+1);
            ihx(i,j)=ihx(i,j)+fi1(i)*curl_e;
            hx(i,j)=fj3(j)*hx(i,j)+fj2(j)*0.5*(curl_e+ihx(i,j));
        end;
    end;             
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %% 计算hy
    for i=1:IE-1;
        for j=1:JE;
            curl_e=ez(i+1,j)-ez(i,j);
            ihy(i,j)=ihy(i,j)+fj1(j)*curl_e;
            hy(i,j)=fj3(i)*hy(i,j)+fi2(i)*0.5*(curl_e+ihy(i,j));
        end;
    end;
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    %%  显示
    subplot(1,2,1);
    mesh(ez);
    axis([1 IE 1 JE -1 1]);
    subplot(1,2,2);
    plot(ez(:,50));
    %axis([1 200]);
    %pause(0.01);
    drawnow;
end;