📄 tmgroovewhq.m
字号:
clear;
lamda=1.55e-6;
N=151;
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=2*order*lamda/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);
W=2*pi*c/lamda;
Ez=zeros(N,N);
Hx=zeros(N+1,N);
Hy=zeros(N,N+1);
%Ey_t=zeros(18*N+1,N+2,M);
%HzReal=zeros(18*N+1,N+1);
Ezi=zeros(N,N);
Hxi=zeros(N+1,N);
Hyi=zeros(N,N+1);
Ezi0=zeros(N,N);
Hxi0=zeros(N+1,N);
Hyi0=zeros(N,N+1);
k=W/c*Neff;
kx=k*cos(-pi/4);
ky=k*sin(-pi/4);
x0=Xx(N+1,1);
y0=Yx(N+1,1);
flag=find(Ep==e0*rela);
Ezi0(flag)=exp(i*(kx*X(flag)+ky*Y(flag))).*exp(-i*W*((X(flag)-x0)*cos(-pi/4)+(Y(flag)-y0)*sin(-pi/4))*Neff/c);
flag=find(Epx==e0*rela);
Hxi0(flag)=ky/W/mu0*exp(i*(kx*Xx(flag)+ky*Yx(flag))).*exp(-i*W*((Xx(flag)-x0)*cos(-pi/4)+(Yx(flag)-y0)*sin(-pi/4))*Neff/c);
flag=find(Epy==e0*rela);
Hyi0(flag)=-kx/W/mu0*exp(i*(kx*Xy(flag)+ky*Yy(flag))).*exp(-i*W*((Xy(flag)-x0)*cos(-pi/4)+(Yy(flag)-y0)*sin(-pi/4))*Neff/c);
Ezi=Ezi0;
Hxi=Hxi0;
Hyi=Hyi0;
%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));
EzxPMLA=zeros(NPML,N);
EzyPMLA=zeros(NPML,N);
HxPMLA=zeros(NPML,N);
HyPMLA=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
EzxPMLB=zeros(NPML,N);
EzyPMLB=zeros(NPML,N);
HxPMLB=zeros(NPML,N);
HyPMLB=zeros(NPML,N+1); %Zone B
Sigmay_zB=flipud(Sigmay_zA);
Sigmay_xB=flipud(Sigmay_xA);
Sigmay_yB=flipud(Sigmay_yA); %Zone B
j=0;
expboundary=exp(-SigmaBound*Dt/e0);
flag=find(Ep==e0*rela);
TimeSteps=20000;
thita=3*pi/4;
for m=1:TimeSteps
%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.
Hx(2:N,:)=Hx(2:N,:)-Dt*(Ez(2:N,:)-Ez(1:N-1,:))/Dy/mu0+...
Dt*i*W*Hxi(2:N,:)-Dt*(Ezi(2:N,:)-Ezi(1:N-1,:))/Dy/mu0;
Hy(:,2:N)=Hy(:,2:N)+Dt*(Ez(:,2:N)-Ez(:,1:N-1))/Dx/mu0+...
Dt*i*W*Hyi(:,2:N)+Dt*(Ezi(:,2:N)-Ezi(:,1:N-1))/Dx/mu0;
Hx(N+1,:)=expboundary*Hx(N+1,:)-(1-expboundary)*...
(EzxPMLA(1,:)+EzyPMLA(1,:)-Ez(N,:))/(SigmaBound*factor*Dy); %Boundary A
Hx(1,:)=expboundary*Hx(1,:)-(1-expboundary)*...
(Ez(1,:)-EzxPMLB(NPML,:)-EzyPMLB(NPML,:))/(SigmaBound*factor*Dy); %Boundary B
Hy(:,N+1)=Hy(:,2)*exp(i*a*kx);
Hy(:,1)=Hy(:,N)*exp(-i*a*kx);
%The following part is for ZONE A. H components.%%%%%%%%%%????????????????????e0?u0?
HxPMLA(1:NPML-1,:)=exp(-Sigmay_xA(1:NPML-1,:)*Dt/e0).*HxPMLA(1:NPML-1,:)-...
(1-exp(-Sigmay_xA(1:NPML-1,:)*Dt/e0))./(Sigmay_xA(1:NPML-1,:)*factor*Dy).*...
(EzxPMLA(2:NPML,:)+EzyPMLA(2:NPML,:)-EzxPMLA(1:NPML-1,:)-EzyPMLA(1:NPML-1,:));
HyPMLA(:,2:N)=HyPMLA(:,2:N)+...
Dt*(EzxPMLA(:,2:N)+EzyPMLA(:,2:N)-EzxPMLA(:,1:N-1)-EzyPMLA(:,1:N-1))/(mu0*Dx);
HyPMLA(:,1)=HyPMLA(:,N)*exp(-i*a*kx); %Boundary Left
HyPMLA(:,N+1)=HyPMLA(:,2)*exp(i*a*kx); %Boundary Right
HxPMLA(NPML,:)=HxPMLA(NPML-1,:); %Boundary Upper;
%The following part is for ZONE B.
%H components.
HxPMLB(2:NPML,:)=exp(-Sigmay_xB(2:NPML,:)*Dt/e0).*HxPMLB(2:NPML,:)-...
(1-exp(-Sigmay_xB(2:NPML,:)*Dt/e0))./(Sigmay_xB(2:NPML,:)*factor*Dy).*...
(EzxPMLB(2:NPML,:)+EzyPMLB(2:NPML,:)-EzxPMLB(1:NPML-1,:)-EzyPMLB(1:NPML-1,:));
HyPMLB(:,2:N)=HyPMLB(:,2:N)+...
Dt*(EzxPMLB(:,2:N)+EzyPMLB(:,2:N)-EzxPMLB(:,1:N-1)-EzyPMLB(:,1:N-1))/(mu0*Dx);
HyPMLB(:,1)=HyPMLB(:,N)*exp(-i*a*kx); %Boundary Left
HyPMLB(:,N+1)=HyPMLB(:,2)*exp(i*a*kx); %Boundary Right
HxPMLB(1,:)=HxPMLB(2,:); %Boundary Bottom;
%E components
Ez(:,:)=Ez(:,:)+Dt*((Hy(:,2:N+1)-Hy(:,1:N))/Dx./Ep-(Hx(2:N+1,:)-Hx(1:N,:))/Dy./Ep)+...
i*Dt*W*Ezi(:,:)+Dt*((Hyi(:,2:N+1)-Hyi(:,1:N))/Dx./Ep-(Hxi(2:N+1,:)-Hxi(1:N,:))/Dy./Ep);
%The following part is for ZONE A.
%Ez component.
EzxPMLA(:,:)=EzxPMLA(:,:)+Dt*(HyPMLA(:,2:N+1)-HyPMLA(:,1:N))/(e0*Dx);
EzyPMLA(2:NPML,:)=exp(-Sigmay_zA(2:NPML,:)*Dt/e0).*EzyPMLA(2:NPML,:)-...
(1-exp(-Sigmay_zA(2:NPML,:)*Dt/e0))./(Sigmay_zA(2:NPML,:)*Dy).*...
(HxPMLA(2:NPML,:)-HxPMLA(1:NPML-1,:));
EzyPMLA(1,:)=exp(-Sigmay_zA(1,:)*Dt/e0).*EzyPMLA(1,:)-...
(1-exp(-Sigmay_zA(1,:)*Dt/e0))./(Sigmay_zA(1,:)*Dy).*...
(HxPMLA(1,:)-Hx(N+1,:));
%The following part is for ZONE B.
%Ez component.
EzxPMLB=EzxPMLB+Dt*(HyPMLB(:,2:N+1)-HyPMLB(:,1:N))/(e0*Dx);
EzyPMLB(1:NPML-1,:)=exp(-Sigmay_zB(1:NPML-1,:)*Dt/e0).*EzyPMLB(1:NPML-1,:)-...
(1-exp(-Sigmay_zB(1:NPML-1,:)*Dt/e0))./(Sigmay_zB(1:NPML-1,:)*Dy).*...
(HxPMLB(2:NPML,:)-HxPMLB(1:NPML-1,:));
EzyPMLB(NPML,:)=exp(-Sigmay_zB(NPML,:)*Dt/e0).*EzyPMLB(NPML,:)-...
(1-exp(-Sigmay_zB(NPML,:)*Dt/e0))./(Sigmay_zB(NPML,:)*Dy).*...
(Hx(1,:)-HxPMLB(NPML,:));
%HzAll=[HzxPMLA+HzyPMLA;Hz;HzxPMLB+HzyPMLB];
% if mod(m,100)==0
%if (m>500&mod(m,100)==0)
% m
% flag=find(Ep==e0&Y<0.2);
% Inte=Hz(find(flag)).*exp(i*(kx*X(flag)+ky*Y(flag)))*Dx*Dy;
% abs(sum(Inte))
%end
if(mod(m,100)==0)
m
Inte=Ez(flag).*exp(-i*(k*cos(thita)*X(flag)+k*sin(thita)*Y(flag)))*Dx*Dy;
temp=abs(sum(Inte))
direc(mod(j,100)+1)=temp;
j=j+1;
%for kk=1:21
%thita=pi/4+pi/20*(kk-1);
%Inte=Ez(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;
% figure(1);
% clf reset;
% surface(X,Y,real(Ez));
% shading interp;
% pause(0.1);
%Ez(203*30+150)
%abs(Hzi(203*30+150))
%pause(0.1);
end
Ezi=Ezi0*exp(-i*W*m*Dt);
Hxi=Hxi0*exp(-i*W*m*Dt);
Hyi=Hyi0*exp(-i*W*m*Dt);
end
thita1=-pi/4;
Inte=Ezi(flag).*exp(-i*(k*cos(thita1)*X(flag)+k*sin(thita1)*Y(flag)))*Dx*Dy;
In=abs(sum(Inte));
R=((Neff-1)/(Neff+1))^2; %normal reflectance
%%direction
%dimen=size(direc);
aver=sum(direc)/100
%load d:\reflect.mat reflect -ASCII;
load d:\db.dat db -ASCII;
reflect(ii+1)=(aver/In)^2 %reflectance by numerical method
db(ii+1)=10*log10(reflect(ii+1))
%save d:\reflect.mat reflect -ASCII -DOUBLE;
save d:\db.dat 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 d:\reflect.mat reflect -ASCII -DOUBLE;
save d:\db.dat db -ASCII -DOUBLE;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -