📄 newtegroove.m
字号:
clear;
lamda0=1.55e-6;
%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.6,1.8,2.0,2.2,2.4,2.6,2.8,3,3.2];
%llist=[1520e-9,1510e-9,1520e-9,1530e-9,1540e-9,1550e-9,1560e-9,1570e-9,1580e-9,1590e-9,1600e-9]
olist=[4,4.5,5,5.5,6,6.5,7,8,8.5,9,9.5,10]
for ii=1:12;
Neff=3.2;
rela=Neff^2;
lamda=lamda0;
order=olist(ii)
%N=ceil(151*order/6.5)+40;
N=161;
if(mod(N,2)==0)
N=N+1;
end
a=2*order*lamda0/Neff*sqrt(2); %width of the cell
b=a;
Ep=ones(N,N)*e0;
x=linspace(-a/2,a/2,N);
X=repmat(x,N,1); %2*N+1,N Ep(i,j)'s x cooridnate
y=linspace(-b/2,b/2,N);
Y=repmat(y',1,N); %2*N+1,N Ep(i,j)'s y cooridnate
Epy=ones(N,N+1)*e0;
xy=linspace(-a/2-a/(N-1)/2,a/2+a/(N-1)/2,N+1); %2*N+1,N+1
Xy=repmat(xy,N,1); %Epy(i,j)'s x coordinate
Yy=repmat(y',1,N+1); %Epy(i,j)'s y coordinate
Epx=ones(N+1,N)*e0;
Xx=repmat(x,N+1,1);
yx=linspace(-b/2-b/2/(N-1),b/2+b/2/(N-1),N+1); %2*N+2,N
Yx=repmat(yx',1,N);
Ep(find(Y>=-b/4))=e0*rela;
Ep(find((X>-a/4&X<=0&Y<-X-b/4)|(X>a/4&X<=a/2&Y<-X+b/4)))=e0;
Epx(find(Yx>=-b/4))=e0*rela;
Epx(find((Xx>-a/4&Xx<=0&Yx<-Xx-b/4)|(Xx>a/4&Xx<=a/2&Yx<-Xx+b/4)))=e0;
Epy(find(Yy>=-b/4))=e0*rela;
Epy(find((Xy>-a/4&Xy<=0&Yy<-Xy-b/4)|(Xy>a/4&Xy<=a/2&Yy<-Xy+b/4)))=e0;
Epy(find(Xy==-a/2-a/(N-1)/2&Yy<-Xy-3*b/4))=e0;
%figure(1);
%surf(Xy,Yy,Epy);
Dx=X(N+1)-X(1);
Dy=-(Y(1)-Y(2));
Dt=1/sqrt(1/(Dx*Dx)+1/(Dy*Dy))/c/3; %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(N,N);
Ex=zeros(N+1,N);
Ey=zeros(N,N+1);
%Ey_t=zeros(2*N+1,N+2,M);
%HzReal=zeros(2*N+1,N+1);
Hzi0=zeros(N,N);
Exi0=zeros(N+1,N);
Eyi0=zeros(N,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;
Hzt=Hzi(1:(N+5)/2,:);
Ext=Exi(1:(N+5)/2,:);
Eyt=Eyi(1:(N+5)/2,:);
Hzs=Hz((N+5)/2+1:N,:);
Exs=Ex((N+5)/2+1:N+1,:);
Eys=Ey((N+5)/2+1:N,:);
Xs=X((N+5)/2+1:N,:);
Ys=Y((N+5)/2+1:N,:);
Xt=X(1:(N+5)/2,:);
Yt=Y(1:(N+5)/2,:);
Epsx=Epx((N+5)/2+1:N+1,:);
Epsy=Epy((N+5)/2+1:N,:);
Eptx=Epx(1:(N+5)/2,:);
Epty=Epy(1:(N+5)/2,:);
M=size(Hzs);
M=M(1);
%Parameters about PML:
factor=mu0/e0;
NPML=15;
n=4; %The order of the polynomial that decribes the conductivity profile.
R=1e-7;
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);
thita=3*pi/4;
TimeSteps=15000;
number=size(Hzs);
s=Dx*Dy*number(1)*number(2)
R=((Neff-1)/(Neff+1))^2 %normal reflectance
%thita=pi-acos((2*lamda/lamda0-1)/sqrt(2));
%thita-3*pi/4
%thita=linspace(pi/2+pi/10,pi-pi/10,71);
%thita(1)=pi-acos((14*lamda/7/lamda0-1)/sqrt(2));
%thita(2)=pi-acos((13*lamda/7/lamda0-1)/sqrt(2));
%thita(3)=pi-acos((12*lamda/7/lamda0-1)/sqrt(2));
for m=1:TimeSteps
Exs(2:M,:)=Exs(2:M,:)+Dt*(Hzs(2:M,:)-Hzs(1:M-1,:))/Dy./Epsx(2:M,:);
Eys(:,2:N)=Eys(:,2:N)-Dt*(Hzs(:,2:N)-Hzs(:,1:N-1))/Dx./Epsy(:,2:N);
Exs(M+1,:)=expboundary*Exs(M+1,:)+(1-expboundary)*...
(HzxPMLA(1,:)+HzyPMLA(1,:)-Hzs(M,:))/(SigmaBound*Dy); %Boundary A
Exs(1,:)=Exs(1,:)+Dt*(Hzs(1,:)-Hzt((N+5)/2,:)+Hzi((N+5)/2,:))/Dy./Epsx(1,:); %Boundary B
Eys(:,N+1)=Eys(:,2)*exp(-i*a*kx);
Eys(:,1)=Eys(:,N)*exp(i*a*kx);
Ext(2:(N+5)/2,:)=Ext(2:(N+5)/2,:)+Dt*(Hzt(2:(N+5)/2,:)-Hzt(1:(N+5)/2-1,:))/Dy./Eptx(2:(N+5)/2,:);
Eyt(:,2:N)=Eyt(:,2:N)-Dt*(Hzt(:,2:N)-Hzt(:,1:N-1))/Dx./Epty(:,2:N);
Eyt(:,1)=Eyt(:,N)*exp(-i*a*kx);
Eyt(:,N+1)=Eyt(:,2)*exp(i*a*kx);
Ext(1,:)=expboundary*Ext(1,:)+(1-expboundary)*...
(Hzt(1,:)-HzxPMLB(NPML,:)-HzyPMLB(NPML,:))/(SigmaBound*Dy);
%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
Hzs(:,:)=Hzs(:,:)-Dt*((Eys(:,2:N+1)-Eys(:,1:N))/Dx-(Exs(2:M+1,:)-Exs(1:M,:))/Dy)/mu0;
Hzt(1:(N+5)/2-1,:)=Hzt(1:(N+5)/2-1,:)-Dt*((Eyt(1:(N+5)/2-1,2:N+1)-Eyt(1:(N+5)/2-1,1:N))/Dx-(Ext(2:(N+5)/2,:)-Ext(1:(N+5)/2-1,:))/Dy)/mu0;
Hzt((N+5)/2,:)=Hzt((N+5)/2,:)-Dt*((Eyt((N+5)/2,2:N+1)-Eyt((N+5)/2,1:N))/Dx-(Exs(1,:)+Exi((N+5)/2+1,:)-Ext((N+5)/2,:))/Dy)/mu0;
%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,:)-Exs(M+1,:));
%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).*...
(Ext(1,:)-ExPMLB(NPML,:));
%Hz(find(abs(Y-k*X-0.25)<Dx))=Hz(find(abs(Y-k*X-0.25)<Dx))+exp(i*(-W*m*Dt));
%H components.
if(mod(m,100)==0)
%surface(Xt,Yt,real(Hzt));
% shading interp;
%pause(0.1);
%for l=1:71
Inte1=Hzs.*exp(-i*(k*cos(thita)*Xs+k*sin(thita)*Ys))*Dx*Dy;
temp1=(abs(sum(sum(Inte1)))/s)^2
direc1(j+1)=temp1;
%Inte2=Hzs.*exp(-i*(k*cos(thita(2))*Xs+k*sin(thita(2))*Ys))*Dx*Dy;
%temp2=(abs(sum(sum(Inte2)))/s)^2
%direc2(j+1)=temp2;
%Inte3=Hzs.*exp(-i*(k*cos(thita(3))*Xs+k*sin(thita(3))*Ys))*Dx*Dy;
%temp3=(abs(sum(sum(Inte3)))/s)^2
%direc3(j+1)=temp3;
%end
j=j+1;
m
end
Hzi=Hzi0*exp(-i*W*m*Dt);
Exi=Exi0*exp(-i*W*m*Dt);
Eyi=Eyi0*exp(-i*W*m*Dt);
%Hzix=Hzix0*exp(-i*W*m*Dt);
%Hziy=Hziy0*exp(-i*W*m*Dt);
%Exiz=Exiz0*exp(-i*W*m*Dt);
%Eyiz=Eyiz0*exp(-i*W*m*Dt);
end
%%direction
aver1=sum(direc1(120:150))/31;
%load d:\reflect.mat reflect -ASCII;
load d:\db.mat db -ASCII;
db(ii+1)=10*log10(aver1)
%save d:\reflect.mat reflect -ASCII -DOUBLE;
save d:\db.mat db -ASCII;
%save d:\direct.dat direc -ASCII -DOUBLE
end
%save d:\reflect.mat reflect -ASCII -DOUBLE;
save d:\db.mat db -ASCII;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -