newtmstair.m

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

clear;
lamda0=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);
%llist=[1500e-9,1510e-9,1520e-9,1530e-9,1540e-9,1550e-9,1560e-9,1570e-9,1580e-9,1590e-9,1600e-9];
%nlist=[1.4,1.4,1.5,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3,3.2];
olist=[0.5,1,2,2.5,3,4,4.5,5,5.5,6,6.5,8,9,10,11,12];
for ii=1:16;
lamda=lamda0;
Neff=1.46;
rela=Neff^2;
order=olist(ii)
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>=0))=e0*rela;
Ep(find((X>=-a/2&X<=-a/4&Y>=X+b/4)|(X>-a/4&X<=a/4&Y>=-X-b/4&Y>=X-b/4)|(X>a/4&X<=a/2&Y>=-X+b/4)))=e0*rela;

Epx(find(Yx>=0))=e0*rela;
Epx(find((Xx>=-a/2&Xx<=-a/4&Yx>=Xx+b/4)|(Xx>-a/4&Xx<=a/4&Yx>=-Xx-b/4&Yx>=Xx-b/4)|(Xx>a/4&Xx<=a/2&Yx>=-Xx+b/4)))=e0*rela;

Epy(find(Yy>=0))=e0*rela;
Epy(find((Xy>=-a/2&Xy<=-a/4&Yy>=Xy+b/4)|(Xy>-a/4&Xy<=a/4&Yy>=-Xy-b/4&Yy>=Xy-b/4)|(Xy>a/4&Xy<=a/2&Yy>=-Xy+b/4)))=e0*rela;
Epy(find(Xy==-a/2-a/(N-1)/2&Yy>-b/4))=e0*rela;
Epy(:,N+1)=Epy(:,1);

%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;
k=W/c*Neff;

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


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



x0=Xx(78,1);
y0=Yx(78,1);

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;

Ezt=Ezi(1:78,:);
Hxt=Hxi(1:78,:);
Hyt=Hyi(1:78,:);

Ezs=Ez(79:N,:);
Hxs=Hx(79:N+1,:);
Hys=Hy(79:N,:);
Xs=X(79:N,:);
Ys=Y(79:N,:);
Xt=X(1:78,:);
Yt=Y(1:78,:);

Eps=Ep(79:N,:);
Ept=Ep(1:78,:);

M=size(Ezs);
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));

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);
TimeSteps=12000;

thita=3*pi/4;
%thita=pi-acos((2*lamda/lamda0-1)/sqrt(2));
number=size(Ezs);
s=number(1)*number(2)*Dx*Dy

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.

	Hxs(2:M,:)=Hxs(2:M,:)-Dt*(Ezs(2:M,:)-Ezs(1:M-1,:))/Dy/mu0;
	Hys(:,2:N)=Hys(:,2:N)+Dt*(Ezs(:,2:N)-Ezs(:,1:N-1))/Dx/mu0;
   
   Hxs(M+1,:)=expboundary*Hxs(M+1,:)-(1-expboundary)*...
      (EzxPMLA(1,:)+EzyPMLA(1,:)-Ezs(M,:))/(SigmaBound*factor*Dy);     %Boundary A
   Hxs(1,:)=Hxs(1,:)-Dt*(Ezs(1,:)-Ezt(78,:)+Ezi(78,:))/Dy/mu0;  %Boundary B
   
   Hys(:,N+1)=Hys(:,2)*exp(-i*a*kx);
   Hys(:,1)=Hys(:,N)*exp(i*a*kx); 

   Hxt(2:78,:)=Hxt(2:78,:)-Dt*(Ezt(2:78,:)-Ezt(1:77,:))/Dy/mu0;
   Hyt(:,2:N)=Hyt(:,2:N)+Dt*(Ezt(:,2:N)-Ezt(:,1:N-1))/Dx/mu0;
   
   Hyt(:,N+1)=Hyt(:,2)*exp(i*a*kx);
   Hyt(:,1)=Hyt(:,N)*exp(-i*a*kx);
    
    Hxt(1,:)=expboundary*Hxt(1,:)-(1-expboundary)*...
      (Ezt(1,:)-EzxPMLB(NPML,:)-EzyPMLB(NPML,:))/(SigmaBound*factor*Dy);  
   
   %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   
   Ezs(:,:)=Ezs(:,:)+Dt*((Hys(:,2:N+1)-Hys(:,1:N))/Dx./Eps-(Hxs(2:M+1,:)-Hxs(1:M,:))/Dy./Eps); 
   Ezt(1:77,:)=Ezt(1:77,:)+Dt*((Hyt(1:77,2:N+1)-Hyt(1:77,1:N))/Dx./Ept(1:77,:)-(Hxt(2:78,:)-Hxt(1:77,:))/Dy./Ept(1:77,:));
   Ezt(78,:)=Ezt(78,:)+Dt*((Hyt(78,2:N+1)-Hyt(78,1:N))/Dx./Ept(78,:)-(Hxs(1,:)+Hxi(79,:)-Hxt(78,:))/Dy./Ept(78,:));
      %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,:)-Hxs(M+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).*...
      (Hxt(1,:)-HxPMLB(NPML,:));    

      
 
  if(mod(m,100)==0)
   m
   Inte=Ezs.*exp(-i*(k*cos(thita)*Xs+k*sin(thita)*Ys))*Dx*Dy;
   temp=(abs(sum(sum(Inte)))/s)^2
   direc(j+1)=temp;
   j=j+1;
            
    %  figure(1);
    % clf reset;
    %  surface(Xt,Yt,real(Ezt));
    %  shading interp;
     % 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
R=((Neff-1)/(Neff+1))^2;  %normal reflectance

%%direction
%dimen=size(direc2);
%aver=sum(direc((dimen(2)-99):dimen(2)))/100
aver=sum(direc(70:120))/51
%load d:\reflect.mat reflect -ASCII;
load d:\db.mat db -ASCII;

%reflect(ii+1)=(aver/s)^2   %reflectance by numerical method
db(ii+1)=10*log10(aver)
%save d:\reflect.mat reflect -ASCII -DOUBLE;
save d:\db.mat db -ASCII ;


end
%save groove250tm direc -ASCII -DOUBLE;
%save d:\reflect.mat reflect -ASCII -DOUBLE;
save d:\db.mat db -ASCII ;