📄 stair.m
字号:
clear;
lamda=1.55e-6;
N=101;
if(mod(N,2)==0)
N=N+1;
end
order=7;
mu0=4*pi*1.0e-7; %Epsilon Zero, if using Gauss Unit, it equals to 1.
e0=8.85*1e-12; %Mu Zero, if using Gauss Unit, it equals to 1.
c=1/sqrt(mu0*e0);
nlist=[1.3,1.4,1.5,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3,3.2];
for ii=1:12;
Neff=nlist(ii);
rela=Neff^2;
a=order*lamda/Neff*sqrt(2); %width of the cell
b=2*a;
Ep=ones(2*N+1,N)*e0;
x=linspace(-a/2+a/(N-1)/2,a/2+a/(N-1)/2,N);
X=repmat(x,2*N+1,1); %2*N+1,N Ep(i,j)'s x cooridnate
y=linspace(b/2/(2*N),b+b/(2*N)/2,2*N+1);
Y=repmat(y',1,N); %2*N+1,N Ep(i,j)'s y cooridnate
Epy=ones(2*N+1,N+1)*e0;
xy=linspace(-a/2,a/2+a/(N-1),N+1); %2*N+1,N+1
Xy=repmat(xy,2*N+1,1); %Epy(i,j)'s x coordinate
Yy=repmat(y',1,N+1); %Epy(i,j)'s y coordinate
Epx=ones(2*N+2,N)*e0;
Xx=repmat(x,2*N+2,1);
yx=linspace(0,b+b/2/N,2*N+2); %2*N+2,N
Yx=repmat(yx',1,N);
Ep(find(Y>=b/2+b/4/N))=e0*rela;
Ep(find((X>=-a/2+a/2/(N-1)&X<=a/2/(N-1)&Y>=X+b/4/N+b/2)|(X>a/2/(N-1)&X<=a/2+a/2/(N-1)&Y>=-X+b/4/N+b/2)))=e0*rela;
Epx(find(Yx>=b/2+b/4/N))=e0*rela;
Epx(find((Xx>=-a/2+a/2/(N-1)&Xx<=a/2/(N-1)&Yx>=Xx+b/4/N+b/2)|(Xx>a/2/(N-1)&Xx<=a/2+a/2/(N-1)&Yx>=-Xx+b/4/N+b/2)))=e0*rela;
Epy(find(Yy>=b/2+b/4/N))=e0*rela;
Epy(find((Xy>=-a/2+a/2/(N-1)&Xy<=a/2/(N-1)&Yy>=Xy+b/4/N+b/2)|(Xy>a/2/(N-1)&Xy<=a/2+a/2/(N-1)&Yy>=-Xy+b/4/N+b/2)))=e0*rela;
Epy(find(Xy==-a/2&Yy>b/4/N+b/4))=e0*rela;
Epy(:,N+1)=Epy(:,1);
X=X-a/2/(N-1);
Xx=Xx-a/2/(N-1);
Xy=Xy-a/2/(N-1);
Y=Y-b/2-b/4/N;
Yx=Yx-b/2-b/4/N;
Yy=Yy-b/2-b/4/N;
%figure(1);
%surf(Xy,Yy,Epy);
Dx=X(2*N+2)-X(1);
Dy=-(Y(1)-Y(2));
Dt=1/sqrt(1/(Dx*Dx)+1/(Dy*Dy))/c; %Time interval
%Dt=1.667592769157701e-011;
W=2*pi*c/lamda;
kx=W/c*Neff*sin(pi/4);
ky=-W/c*Neff*cos(pi/4);
Hz=zeros(2*N+1,N);
Ex=zeros(2*N+2,N);
Ey=zeros(2*N+1,N+1);
%Ey_t=zeros(2*N+1,N+2,M);
%HzReal=zeros(2*N+1,N+1);
Hzi0=zeros(2*N+1,N);
Exi0=zeros(2*N+2,N);
Eyi0=zeros(2*N+1,N+1);
flag=find(Ep==e0*rela);
Hzi0(flag)=exp(i*(kx*X(flag)+ky*Y(flag)));
flag=find(Epx==e0*rela);
Exi0(flag)=-ky/W/(e0*rela)*exp(i*(kx*Xx(flag)+ky*Yx(flag)));
flag=find(Epy==e0*rela);
Eyi0(flag)=kx/W/(e0*rela)*exp(i*(kx*Xy(flag)+ky*Yy(flag)));
Hzi=Hzi0;
Exi=Exi0;
Eyi=Eyi0;
%Parameters about PML:
factor=mu0/e0;
NPML=22;
n=4; %The order of the polynomial that decribes the conductivity profile.
R=1e-8;
Delta=NPML*Dy;
SigmaMax=-(n+1)*e0*c*log(R)/(Delta*2);
NUM=NPML*2:-1:1;
Sigmay=SigmaMax*((NUM*Dy/2+Dy/2).^(n+1)-(NUM*Dy/2-Dy/2).^(n+1))/(Delta^n*Dy*(n+1));
SigmaBound=SigmaMax*(Dy/2).^(n+1)/(Delta^n*Dy*(n+1));
Sigmax=Sigmay;
Sigmay_z1=fliplr(repmat(Sigmax(2:2:NPML*2),NPML,1));
Sigmay_x1=fliplr(repmat(Sigmax(1:2:NPML*2-1),NPML,1));
Sigmay_y1=fliplr(repmat(Sigmax(2:2:NPML*2),NPML,1));
HzxPMLA=zeros(NPML,N);
HzyPMLA=zeros(NPML,N);
ExPMLA=zeros(NPML,N);
EyPMLA=zeros(NPML,N+1); %Zone A
Sigmay_zA=repmat(Sigmay_z1(1,:)',1,N);
Sigmay_xA=repmat(Sigmay_x1(1,:)',1,N);
Sigmay_yA=repmat(Sigmay_y1(1,:)',1,N+1); %Zone A
HzxPMLB=zeros(NPML,N);
HzyPMLB=zeros(NPML,N);
ExPMLB=zeros(NPML,N);
EyPMLB=zeros(NPML,N+1); %Zone B
Sigmay_zB=flipud(Sigmay_zA);
Sigmay_xB=flipud(Sigmay_xA);
Sigmay_yB=flipud(Sigmay_yA); %Zone B
k=W/c*Neff;
j=0;
expboundary=exp(-SigmaBound*Dt/e0);
flag=find(Ep==e0*rela);
TimeSteps=40000;
thita=3*pi/4;
for m=1:TimeSteps
Ex(2:2*N+1,:)=Ex(2:2*N+1,:)+Dt*(Hz(2:2*N+1,:)-Hz(1:2*N,:))/Dy./Epx(2:2*N+1,:)+...
Dt*i*W*Exi(2:2*N+1,:)+i*Dt*ky*Hzi(2:2*N+1,:)./Epx(2:2*N+1,:);
Ey(:,2:N)=Ey(:,2:N)-Dt*(Hz(:,2:N)-Hz(:,1:N-1))/Dx./Epy(:,2:N)+...
Dt*i*W*Eyi(:,2:N)-i*kx*Dt*Hzi(:,2:N)./Epy(:,2:N);
Ex(2*N+2,:)=expboundary*Ex(2*N+2,:)+(1-expboundary)*...
(HzxPMLA(1,:)+HzyPMLA(1,:)-Hz(2*N+1,:))/(SigmaBound*Dy); %Boundary A
Ex(1,:)=expboundary*Ex(1,:)+(1-expboundary)*...
(Hz(1,:)-HzxPMLB(NPML,:)-HzyPMLB(NPML,:))/(SigmaBound*Dy); %Boundary B
Ey(:,N+1)=Ey(:,2)*exp(i*a*kx);
Ey(:,1)=Ey(:,N)*exp(-i*a*kx);
%The following part is for ZONE A. E components.%%%%%%%%%%????????????????????e0?u0?
ExPMLA(1:NPML-1,:)=exp(-Sigmay_xA(1:NPML-1,:)*Dt/e0).*ExPMLA(1:NPML-1,:)+...
(1-exp(-Sigmay_xA(1:NPML-1,:)*Dt/e0))./(Sigmay_xA(1:NPML-1,:)*Dy).*...
(HzxPMLA(2:NPML,:)+HzyPMLA(2:NPML,:)-HzxPMLA(1:NPML-1,:)-HzyPMLA(1:NPML-1,:));
EyPMLA(:,2:N)=EyPMLA(:,2:N)-...
Dt*(HzxPMLA(:,2:N)+HzyPMLA(:,2:N)-HzxPMLA(:,1:N-1)-HzyPMLA(:,1:N-1))/(e0*Dx);
EyPMLA(:,N+1)=EyPMLA(:,2)*exp(i*a*kx);
EyPMLA(:,1)=EyPMLA(:,N)*exp(-i*a*kx);%Boundary Right
ExPMLA(NPML,:)=ExPMLA(NPML-1,:); %Boundary Upper;
%The following part is for ZONE B.
%E components.
ExPMLB(2:NPML,:)=exp(-Sigmay_xB(2:NPML,:)*Dt/e0).*ExPMLB(2:NPML,:)+...
(1-exp(-Sigmay_xB(2:NPML,:)*Dt/e0))./(Sigmay_xB(2:NPML,:)*Dy).*...
(HzxPMLB(2:NPML,:)+HzyPMLB(2:NPML,:)-HzxPMLB(1:NPML-1,:)-HzyPMLB(1:NPML-1,:));
EyPMLB(:,2:N)=EyPMLB(:,2:N)-...
Dt*(HzxPMLB(:,2:N)+HzyPMLB(:,2:N)-HzxPMLB(:,1:N-1)-HzyPMLB(:,1:N-1))/(e0*Dx);
EyPMLB(:,N+1)=EyPMLB(:,2)*exp(i*a*kx);
EyPMLB(:,1)=EyPMLB(:,N+1)*exp(-i*a*kx);%Boundary Right
ExPMLB(1,:)=ExPMLB(2,:); %Boundary Bottom;
%H components
Hz(:,:)=Hz(:,:)-Dt*((Ey(:,2:N+1)-Ey(:,1:N))/Dx-(Ex(2:2*N+2,:)-Ex(1:2*N+1,:))/Dy)/mu0+...
i*Dt*W*Hzi(:,:)-i*Dt/mu0*(kx*Eyi(:,1:N)-ky*Exi(1:2*N+1,:));
%The following part is for ZONE A.
%Hz component.
HzxPMLA=HzxPMLA-Dt*(EyPMLA(:,2:N+1)-EyPMLA(:,1:N))/(mu0*Dx);
HzyPMLA(2:NPML,:)=exp(-Sigmay_zA(2:NPML,:)*Dt/e0).*HzyPMLA(2:NPML,:)+...
(1-exp(-Sigmay_zA(2:NPML,:)*Dt/e0))./(Sigmay_zA(2:NPML,:)*factor*Dy).*...
(ExPMLA(2:NPML,:)-ExPMLA(1:NPML-1,:));
HzyPMLA(1,:)=exp(-Sigmay_zA(1,:)*Dt/mu0).*HzyPMLA(1,:)+...
(1-exp(-Sigmay_zA(1,:)*Dt/mu0))./(Sigmay_zA(1,:)*factor*Dy).*...
(ExPMLA(1,:)-Ex(2*N+2,:));
%The following part is for ZONE B.
%Hz component.
HzxPMLB=HzxPMLB-Dt*(EyPMLB(:,2:N+1)-EyPMLB(:,1:N))/(mu0*Dx);
HzyPMLB(1:NPML-1,:)=exp(-Sigmay_zB(1:NPML-1,:)*Dt/e0).*HzyPMLB(1:NPML-1,:)+...
(1-exp(-Sigmay_zB(1:NPML-1,:)*Dt/e0))./(Sigmay_zB(1:NPML-1,:)*factor*Dy).*...
(ExPMLB(2:NPML,:)-ExPMLB(1:NPML-1,:));
HzyPMLB(NPML,:)=exp(-Sigmay_zB(NPML,:)*Dt/e0).*HzyPMLB(NPML,:)+...
(1-exp(-Sigmay_zB(NPML,:)*Dt/e0))./(Sigmay_zB(NPML,:)*factor*Dy).*...
(Ex(1,:)-ExPMLB(NPML,:));
if(mod(m,100)==0)
m
%for kk=1:21
%thita=pi/4+pi/20*(kk-1);
%Inte=Hz(flag).*exp(-i*(k*cos(thita)*X(flag)+k*sin(thita)*Y(flag)))*Dx*Dy;
%abs(sum(Inte))
%direc(j,kk)=abs(sum(Inte));
%end
%j=j+1;
Inte=Hz(flag).*exp(-i*(k*cos(thita)*X(flag)+k*sin(thita)*Y(flag)))*Dx*Dy;
temp=abs(sum(Inte))
direc(mod(j,120)+1)=temp;
j=j+1;
%Hz(203*30+150)
%abs(Hzi(203*30+150))
%pause(0.1);
end
%if mod(m,10)==0
%m
% Epss=Ep;
% Epss=Epss(:,1:N+1-mod(fix(m/step),loop));
% Epss=[temp Epss];
% figure(2);
% clf reset;
% surface(X,Y,real(Hz));
% shading interp;
%surface(X,Y,real(Hz));
%pause(0.2);
% end
Hzi=Hzi0*exp(-i*W*m*Dt);
Exi=Exi0*exp(-i*W*m*Dt);
Eyi=Eyi0*exp(-i*W*m*Dt);
end
thita=-pi/4;
Inte=Hzi(flag).*exp(-i*(k*cos(thita)*X(flag)+k*sin(thita)*Y(flag)))*Dx*Dy;
In=abs(sum(Inte))
R=((Neff-1)/(Neff+1))^2 %normal reflectance
%%direction
aver=sum(direc)/120
load e:\reflect.mat reflect -ASCII;
load e:\db.mat db -ASCII;
reflect(ii+1)=(aver/In)^2 %reflectance by numerical method
db(ii+1)=10*log10(reflect(ii+1))
save e:\reflect.mat reflect -ASCII -DOUBLE;
save e:\db.mat db -ASCII -DOUBLE;
end
%xx=1/4:1/20:5/4;
%hold on
%for j=(dimen(1)-10):dimen(1)
% plot(xx,direc(j,:));
%end
%save groove250tm direc -ASCII -DOUBLE;
save e:\reflect.mat reflect -ASCII -DOUBLE;
save e:\db.mat db -ASCII -DOUBLE;⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -