doubleg.m

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

clear;
lamda=1.55e-6;
N=151;
if(mod(N,2)==0)
    N=N+1;
end
order=7;
nlist=[1.3,1.4,1.5,1.6,1.8,2.0,2.2,2.4,2.6,2.8,3,3.2];
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);

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

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);
Hzi=zeros(N,N);
Exi=zeros(N+1,N);
Hzix0=zeros(N+1,N);
Eyi=zeros(N,N+1);
Hziy0=zeros(N,N+1);
%Exiz0=zeros(N,N);
%Eyiz0=zeros(N,N);

flag=find(Y==b/2);
Hzi(flag)=exp(i*(kx*X(flag)+ky*Y(flag)));
Hz(flag)=Hzi(flag);

flag=find(Yx>=(b/2-b/(N-1))&Yx~=b/2+b/(N-1)/2);
Exi(flag)=-ky/W/(e0*rela)*exp(i*(kx*Xx(flag)+ky*Yx(flag)));
Ex(flag)=Exi(flag);

flag=find(Yy==b/2);
Eyi(flag)=kx/W/(e0*rela)*exp(i*(kx*Xy(flag)+ky*Yy(flag)));
Ey(flag)=Eyi(flag);



%Exiz=Exiz0;
%Eyiz=Eyiz0;

%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=25000;


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.

	Ex(2:N,:)=Ex(2:N,:)+Dt*(Hz(2:N,:)-Hz(1:N-1,:))/Dy./Epx(2:N,:);
   %i*Dt*ky*Hzix(2:N,:)./Epx(2:N,:);
	Ey(:,2:N)=Ey(:,2:N)-Dt*(Hz(:,2:N)-Hz(:,1:N-1))/Dx./Epy(:,2:N);
   %-i*kx*Dt*Hziy(:,2:N)./Epy(:,2:N);
   
         
   Ex(N+1,:)=expboundary*Ex(N+1,:)+(1-expboundary)*...
      (HzxPMLA(1,:)+HzyPMLA(1,:)-Hz(N,:))/(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:N+1,:)-Ex(1:N,:))/Dy)/mu0;
   %-i*Dt/mu0*(kx*(Eyi(:,2:N+1)+Eyi(:,1:N))/2-ky*(Exi(2:N+1,:)+Exi(1:N,:))/2);
  
        %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(N+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).*...
      (Ex(1,:)-ExPMLB(NPML,:));    
   
   if(mod(m,100)==0)
      m
      %surface(X,Y,real(Hz));
      %shading interp;
      %pause(0.1);
      flag=find(Ep==e0*rela&Y<(b/2-b/(N-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,100)+1)=temp;
      j=j+1;



   end

 
   
flag=find(Y==b/2);
Hz(flag)=Hzi(flag)*exp(-i*W*m*Dt);

flag=find(Yx>=(b/2-b/(N-1))&Yx~=b/2+b/(N-1)/2);
Ex(flag)=Exi(flag)*exp(-i*W*m*Dt);

flag=find(Yy==b/2);
Ey(flag)=Eyi(flag)*exp(-i*W*m*Dt);

end


thita=-pi/4;
flag=find(Ep==e0*rela&Y<(b/2-b/(N-1)));
aa=zeros(N,N);
aa(flag)=1;

s=abs(sum(sum(aa)))*Dx*Dy;
R=((Neff-1)/(Neff+1))^2;  %normal reflectance

%%direction
aver=sum(direc)/100;

load d:\reflect.dat reflect -ASCII;
load d:\db.dat db -ASCII;

reflect(ii+1)=(aver/s)^2   %reflectance by numerical method
db(ii+1)=10*log10(reflect(ii+1))

save d:\reflect.dat reflect -ASCII -DOUBLE;
save d:\db.dat db -ASCII -DOUBLE;


end
save d:\reflect.dat reflect -ASCII -DOUBLE;
save d:\db.dat db -ASCII -DOUBLE;