2d_te.m

来自「二维波导中计算反射系数」· M 代码 · 共 198 行

clear all;
clc;

dx=1.5e-2;
dy=dx;
dt=1.81e-11;       

%NX=52;
%NY=7;
for i=1:53
    for j=1:8
        ex(i,j)=0.0;
        D(i,j)=0.0;
        ey(i,j)=0.0;
        hz(i,j)=0.0;
    end
end
%initial field points

eps0=8.85e-12;
mu0=1.2566e-006;%H/m
c=3e8;%m/s
sigma=0.0; 

%free space parameter
m=1;
d=8;  % thickness of PML layer 
R0=1e-8;  %Page167
sigmax=log(R0).*(m+1).*eps0.*c;
sigmax=-0.8.*sigmax./(2.*8.*dx);
for i=1:9
    sigex(i)=sigmax.*(((d+1-0.5-i)/(d-0.5)).^m);
    sighz(i)=sigex(i);
    sigey(i)=sigmax.*(((d+1-i)/d).^m);
end
    sigex(9)=0.0;
    sighz(9)=0.0;
    epsr=1.0;

    for i=1:51                   
    epsrx(i)=1.0;
    epsry(i)=1.0;
end


%epsry(30)=2.5;   
%epsry(40)=2.5;
%for i=31:39
%    epsrx(i)=4.0;
%    epsry(i)=4.0;
%end

%Add the substrate


f=3e9;          %define the source: Gauss Pulse
tao=2./3e9;
tao=tao./dt;
for i=1:400
    if i<60
        source(i)=exp(-4.*pi.*((i-0.8.*tao).^2)./(tao.*tao));
    else source(i)=0.0;
    end
end 
%plot((1:400),source);

%main loop
for n=1:400
    %free space

    for i=10:43
        for j=3:6 %先不算ＰＥＣ　boundary
            ex(i,j) =  ex(i,j)+(hz(i,j)-hz(i,j-1)).*dt.*c./(dy*epsrx(i));
        end
    end
    for i=11:43
        for j=2:6
            ey(i,j) =  ey(i,j)-(hz(i,j)-hz(i-1,j)).*dt.*c./(dx*epsry(i));
        end
    end
    for i=10:43
        for j=2:6
            hz(i,j) =  hz(i,j)+((ex(i,j+1)-ex(i,j))-(ey(i+1,j)-ey(i,j))).*dt.*c/(dx);
        end
    end
    %pec boudary
    for i=2:51 %10:43
            ex(i,2)=0.0;
            ex(i,7)=0.0;
    end

    %pml boundary
    for i=2:10
        temp=sigey(i-1)*dt/(2*eps0);
        for j=2:6
            ey(i,j)=((1-temp)./(1+temp)).*ey(i,j)-(hz(i,j)-hz(i-1,j)).*dt.*c./(dx.*(1+temp));
        end
    end
    for i=44:52
        temp=sigey(52+1-i)*dt/(2*eps0);
        for j=2:6           
            ey(i,j)=((1-temp)./(1+temp)).*ey(i,j)-(hz(i,j)-hz(i-1,j)).*dt.*c./(dx.*(1+temp));
        end
    end
    %ey
    for i=2:9
        temp=sigex(i-1)*dt/(2*eps0);
        for j=3:6 %未算pec边界处ex值
            temp_D=D(i,j);
            D(i,j)=D(i,j)+(hz(i,j)-hz(i,j-1)).*dt.*c./(dy);  
            ex(i,j)=ex(i,j)+((1+temp).*D(i,j)-(1-temp).*temp_D);
        end
    end
    for i=44:51
        temp=sigex(52-i)*dt/(2*eps0);
        for j=3:6 %2:7
            temp_D=D(i,j);
            D(i,j)=D(i,j)+(hz(i,j)-hz(i,j-1)).*dt.*c./(dy);
            ex(i,j)=ex(i,j)+((1+temp).*D(i,j)-(1-temp).*temp_D);
        end
    end
    %ex
    
    for i=2:9 
        temp=sighz(i-1)*dt/(2*eps0);
        for j=2:6          
            hz(i,j)=((1-temp)./(1+temp)).*hz(i,j)+((ex(i,j+1)-ex(i,j))-(ey(i+1,j)-ey(i,j))).*dt.*c/(dx.*(1+temp));
        end
    end
    for i=44:51
        temp=sighz(52-i)*dt/(2*eps0);
        for j=2:6          
            hz(i,j)=((1-temp)./(1+temp)).*hz(i,j)+((ex(i,j+1)-ex(i,j))-(ey(i+1,j)-ey(i,j))).*dt.*c/(dx.*(1+temp));
        end
    end
    %hz
    
%    for i=2
%        for j=2:6
%            ey(i,j)=0.0;
%        end
%    end
%    for i=52
%       for j=2:6
%           ey(i,j)=0.0;
%       end
%   end
      

    %Add source
    if n<100
        for j=2:6
            ey(15,j)=source(n);
        end
    end
    %draw
    
   

    plot((1:53),ey(:,4));   %display
    xlim([1 50]);
    ylim([-1 1]);
    M(n)=getframe;

    
    %note
  notein(n)=ey(20,4);
%  notetotle(n)=ey(20,4);
end

 figure(2);plot((1:400),notein);
%figure(3);plot((1:400),notetotle);


%notere=notetotle-notein;
%figure(4);plot((1:400),notere);

%nf=2^15;
%df=1./(dt.*nf);
%ff=1:df:df*2^15;
%fin=fft(notein,nf);
%fin=abs(fin);
%plot(ff(1:3000),fin(1:3000));

%fre=fft(notere,nf);
%fre=abs(fre);
%plot(ff(1:3000),fre(1:3000));

%fco=fre./fin;
%figure(5);
%plot(ff(1:3000),fco(1:3000));

%figure(1);plot((1:500),note10);
%figure(2);plot((1:500),note20);
%figure(3);plot((1:500),note30);
%figure(4);plot((1:500),note40);