groovetmwl.m

二微时域有限差分的matlab模拟 二微时域有限差分的matlab模拟

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>=b/4/N+b/4)|(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>=b/4/N+b/4)|(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>=b/4/N+b/4)|(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(find(Xy==a/2+a/(N-1)&Yy>=b/4/N+b/4))=e0*rela;

Y=Y-b/2-b/4/N;
Yx=Yx-b/2-b/4/N;
Yy=Yy-b/2-b/4/N;
X=X-a/2/(N-1);
Xx=Xx-a/2/(N-1);
Xy=Xy-a/2/(N-1);

%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/10; %Time interval
%Dt=1.667592769157701e-011;

W=2*pi*c/lamda;

Ez=zeros(2*N+1,N);
Hx=zeros(2*N+2,N);
Hy=zeros(2*N+1,N+1);
%Ey_t=zeros(18*N+1,N+2,M);
%HzReal=zeros(18*N+1,N+1);
Ezi=zeros(2*N+1,N);
Hxi=zeros(2*N+2,N);
Hyi=zeros(2*N+1,N+1);

Ezi0=zeros(2*N+1,N);
Hxi0=zeros(2*N+2,N);
Hyi0=zeros(2*N+1,N+1);


k=W/c*Neff;

kx=k*cos(-pi/4);
ky=k*sin(-pi/4);

flag=find(Ep==e0*rela);
Ezi0(flag)=exp(i*(kx*X(flag)+ky*Y(flag)));
flag=find(Epx==e0*rela);
Hxi0(flag)=ky/W/mu0*exp(i*(kx*Xx(flag)+ky*Yx(flag)));
flag=find(Epy==e0*rela);
Hyi0(flag)=-kx/W/mu0*exp(i*(kx*Xy(flag)+ky*Yy(flag)));


Ezi=Ezi0;
Hxi=Hxi0;
Hyi=Hyi0;

%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));

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=40000;

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:2*N+1,:)=Hx(2:2*N+1,:)-Dt*(Ez(2:2*N+1,:)-Ez(1:2*N,:))/Dy/mu0+...
   	Dt*i*W*Hxi(2:2*N+1,:)-Dt*i*ky*Ezi(2:2*N+1,:)/mu0;
	Hy(:,2:N)=Hy(:,2:N)+Dt*(Ez(:,2:N)-Ez(:,1:N-1))/Dx/mu0+...
   	Dt*i*W*Hyi(:,2:N)+i*Dt*kx*Ezi(:,2:N)/mu0;
   
   Hx(2*N+2,:)=expboundary*Hx(2*N+2,:)-(1-expboundary)*...
      (EzxPMLA(1,:)+EzyPMLA(1,:)-Ez(2*N+1,:))/(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:2*N+2,:)-Hx(1:2*N+1,:))/Dy./Ep)+...
      i*Dt*W*Ezi(:,:)+i*kx*Dt*Hyi(:,2:N+1)./Ep-i*ky*Dt*Hxi(2:2*N+2,:)./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(2*N+2,:));    
   
      %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,120)+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)/120

load d:\reflect.mat reflect -ASCII;
load d:\db.mat 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.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 d:\reflect.mat reflect -ASCII -DOUBLE;
save d:\db.mat db -ASCII -DOUBLE;