doubles.m

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

clear;
lamda=1.55e-6;
N=201;
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);
Neff=1.5
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>=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;

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);
Hzix0=zeros(N+1,N);
Eyi0=zeros(N,N+1);
Hziy0=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)));

flag=find(Epx==e0*rela);
Hzix0(flag)=exp(i*(kx*Xx(flag)+ky*Yx(flag)));
flag=find(Epy==e0*rela);
Hziy0(flag)=exp(i*(kx*Xy(flag)+ky*Yy(flag)));

Hzi=Hzi0;
Exi=Exi0;
Eyi=Eyi0;
Hzix=Hzix0;
Hziy=Hziy0;
%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=1;

expboundary=exp(-SigmaBound*Dt/e0);
thita=3*pi/4;
TimeSteps=30000;
flag=find(Ep==e0*rela);

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,:)+...
   	Dt*i*W*Exi(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)+...
   	Dt*i*W*Eyi(:,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*W*Hzi(:,:)-i*Dt/mu0*(kx*Eyi(:,1:N)-ky*Exi(1:N,:));
  
        %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,10)==0)
      m
      surface(X,Y,real(Hz));
      shading interp;
      pause(0.1);

      Inte=Hz(flag).*exp(-i*(k*cos(thita)*X(flag)+k*sin(thita)*Y(flag)))*Dx*Dy;
      temp=abs(sum(Inte))
      direc(j)=temp;
      j=j+1;


%Hz(203*30+150)
%abs(Hzi(203*30+150))
%pause(0.1);
   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);
end

thita=-pi/4;
Inte=Hzi(flag).*exp(-i*(k*cos(thita)*X(flag)+k*sin(thita)*Y(flag)))*Dx*Dy;
In=abs(sum(Inte))
R=((Neff-1)/(Neff+1))^2  %normal reflectance

%%direction
dimen=size(direc);
aver=sum(direc(dimen(2)-100:dimen(2)))/101;
reflect=(aver/In)^2    %reflectance by numerical method
db=10*log10(reflect)

%save e:\double200.mat direc -ASCII -DOUBLE;