fdbpm_3d_paper.m

有限差分法

function FDBPM_3D_paper
clc
clear
%计算常量
r=0.5;
lamda=1.55e-6;
k0=2*pi/lamda;
dx=0.4;
dy=0.4;
dz=0.2;
thickness=8;
height=2;
width=12;
xwidth=60;
ywidth=60;
M=xwidth/dx;
N=ywidth/dy;
L=500;
er =ones(M, N); 
ed =-9.89*ones(M, N);
exy=zeros(M,N);
er( 1:M, N/2+thickness/2+1: N ) = 10.6;
ed( 1:M, N/2+thickness/2+1: N ) = -0.03;
er( M/2-width/2:M/2+width/2, N/2-thickness/2-1: N/2+thickness/2 ) = 11.7;
ed( M/2-width/2:M/2+width/2, N/2-thickness/2-1: N/2+thickness/2) = 0.07;
er (M/2-width/2:M/2+width/2, N/2-thickness/2-height: N/2-thickness/2 ) = 4.37;
ed (M/2-width/2:M/2+width/2, N/2-thickness/2-height: N/2-thickness/2 ) = -1.49;
exy( M/2-width/2:M/2+width/2, N/2-thickness/2-height: N/2-thickness/2 ) = 0;
%mesh( er(:,:) )
%mesh( ed(:,:) )
%axis([1 M 1 N -10 12])
%view(0,90)
Power1 = zeros(1,L);
Power2 = zeros(1,L);

d=dx/dy;         %dx/dy
for m=1:M
    for n=1:N
        P(m,n,1)=1*exp( -(m-M/2)^2/25)*exp( -(n-N/2)^2/25);
    end
end
      for n=2:N-1
        for m=2:M-1
            c=2i*dy^2/dz-2+(1-r)*er(m,n)*ed(m,n)*(dy^2);
        E1(m,n)=d*(P(m,n+1,1)+c*P(m,n,1)+P(m,n-1,1)); 
       end
     end

 Ey2=zeros(M,N);
%边界条件
for m=1:M
        if    P(m,2,1)~=0
              klxy1(m)=i*(log(P(m,1,1)/P(m,2,1)))/dy;
              klxy(m)=abs(real(klxy1(m)))-i*abs(imag(klxy1(m)));
        else  klxy(m)=0;
        end
        if    P(m,N-1,1)~=0
              krxy1(m)=i*(log(P(m,N,1)/P(m,N-1,1)))/dy;
              krxy(m)=abs(real(krxy1(m)))-i*abs(imag(krxy1(m)));
        else  krxy(m)=0;
        end
        klyy(m)=0;
        kryy(m)=0;
 end
    
    for n=1:N
        if    P(2,n,1)~=0
              klxx1(n)=i*(log(P(1,n,1)/P(2,n,1)))/dx;
              klxx(n)=abs(real(klxx1(n)))-i*abs(imag(klxx1(n)));
        else  klxx(n)=0;
        end
        if    P(M-1,n,1)~=0
              krxx1(n)=i*(log(P(M,n,1)/P(M-1,n,1)))/dx;
              krxx(n)=abs(real(krxx1(n)))-i*abs(imag(krxx1(n)));
        else  krxx(n)=0;
        end
        klyx(n)=0;
        kryx(n)=0;    
    end 
  
%求解矩阵,迭代计算
%z->z+dz/2
for  l=1:L
  l
%l->l+0.5,即z->z+dz/2
for n=2:N-1
     a=2i*dx^2/dz+2-r*er(2,n)*ed(2,n)*dx^2;
         b11=a-exp(-i*klxx(n)*dx);
         x1(2,n)=E1(2,n)/b11;      
        for   p=3:(M-2)
           a=2i*dx^2/dz+2-r*er(p,n)*ed(p,n)*dx^2;
            c11(p)=-1/b11;
            b11=a+c11(p); 
            x1(p,n)=(E1(p,n)+x1(p-1,n))/b11;
        end  
        a=2i*dx^2/dz+2-r*er(M-1,n)*ed(M-1,n)*dx^2;
         c11(M-1)=-1/b11;
         b11=a-exp(-i*krxx(n)*dx)+c11(M-1);
         x1(M-1,n)=(E1(M-1,n)+x1(M-2,n))/b11;
          Ex1(M-1,n)=x1(M-1,n);
        for v=(M-2):-1:2
            Ex1(v,n)=x1(v,n)-c11(v+1)*Ex1(v+1,n);
        end   
            Ex1(1,n)=Ex1(2,n)*exp(-i*klxx(n)*dx);
            Ex1(M,n)=Ex1(M-1,n)*exp(-i*krxx(n)*dx);      
   end        
%计算Ex1边界 
for m=2:M-1
    Ex1(m,1)=Ex1(m,2)*exp(-i*klxy(m)*dy);
    Ex1(m,N)=Ex1(m,N-1)*exp(-i*krxy(m)*dy);
end
    Ex1(1,1)=Ex1(2,1)*exp(-i*klxx(1)*dx);
    Ex1(1,N)=Ex1(2,N)*exp(-i*klxx(N)*dx);
    Ex1(M,1)=Ex1(M,2)*exp(-i*klxy(M)*dy);
    Ex1(M,N)=Ex1(M,N-1)*exp(-i*krxy(M)*dy);

%计算头半步Ey1
      for n=2:N-1
        for m=2:M-1
            c=2i*dy^2/dz-2+(1-r)*er(m,n)*ed(m,n)*(dy^2);
            b1=er(m,n)*exy(m,n)*(dx^2);
        F1(m,n)=d*(Ey2(m,n+1)+c*Ey2(m,n)+Ey2(m,n-1))-b1*Ex1(m,n);
       end
     end
for n=2:N-1
        a=2i*dx^2/dz+2-r*er(2,n)*ed(2,n)*dx^2;
         b12=a-exp(-i*klyx(n)*dx);
         y1(2,n)=F1(2,n)/b12;      
        for   p=3:(M-2)
            a=2i*dx^2/dz+2-r*er(p,n)*ed(p,n)*dx^2;
            c12(p)=-1/b12;
            b12=a+c12(p); 
            y1(p,n)=(F1(p,n)+y1(p-1,n))/b12;
        end  
        a=2i*dx^2/dz+2-r*er(M-1,n)*ed(M-1,n)*dx^2;
         c12(M-1)=-1/b12;
         b12=a-exp(-i*kryx(n)*dx)+c12(M-1);
         y1(M-1,n)=(F1(M-1,n)+y1(M-2,n))/b12;
          Ey1(M-1,n)=y1(M-1,n);
        for v=(M-2):-1:2
            Ey1(v,n)=y1(v,n)-c12(v+1)*Ey1(v+1,n);
        end   
            Ey1(1,n)=Ey1(2,n)*exp(-i*klyx(n)*dx);
            Ey1(M,n)=Ey1(M-1,n)*exp(-i*kryx(n)*dx);      
   end 
%计算Ey1边界 
for m=2:M-1
    Ey1(m,1)=Ey1(m,2)*exp(-i*klyy(m)*dy);
    Ey1(m,N)=Ey1(m,N-1)*exp(-i*kryy(m)*dy);
end
    Ey1(1,1)=Ey1(2,1)*exp(-i*klyx(1)*dx);
    Ey1(1,N)=Ey1(2,N)*exp(-i*klyx(N)*dx);
    Ey1(M,1)=Ey1(M,2)*exp(-i*klyy(M)*dy);
    Ey1(M,N)=Ey1(M,N-1)*exp(-i*kryy(M)*dy);
%计算H(:,:)     
for n=2:N-1
    for m=2:M-1
         c=2i*(dy^2)/dz-2+(1-r)*er(m,n)*ed(m,n)*(dy^2);
         b2=er(m,n)*exy(m,n)*(dy^2);
        H(m,n)=(Ey1(m+1,n)+c*Ey1(m,n)+Ey1(m-1,n))/d-b2*Ex1(m,n);
    end
end   

%x,y的边界   
%x的边界
    for m=1:M
        if    Ex1(m,2)~=0
              klxy1(m)=i*(log(Ex1(m,1)/Ex1(m,2)))/dy; 
              klxy(m)=abs(real(klxy1(m)))-i*abs(imag(klxy1(m)));
        else  klxy(m)=0;
        end
        if    Ex1(m,N-1)~=0
              krxy1(m)=i*(log(Ex1(m,N)/Ex1(m,N-1)))/dy;
              krxy(m)=abs(real(krxy1(m)))-i*abs(imag(krxy1(m)));
        else  krxy(m)=0;
        end
    end
  for n=1:N
        if    Ex1(2,n)~=0
              klxx1(n)=i*(log(Ex1(1,n)/Ex1(2,n)))/dx;
              klxx(n)=abs(real(klxx1(n)))-i*abs(imag(klxx1(n)));
        else  klxx(n)=0;
        end
        if    Ex1(M-1,n)~=0
              krxx1(n)=i*(log(Ex1(M,n)/Ex1(M-1,n)))/dx;
              krxx(n)=abs(real(krxx1(n)))-i*abs(imag(krxx1(n)));
        else  krxx(n)=0;
        end
    end
%y的边界     
   for m=1:M
        if    Ey1(m,2)~=0
              klyy1(m)=i*(log(Ey1(m,1)/Ey1(m,2)))/dy;
              klyy(m)=abs(real(klyy1(m)))-i*abs(imag(klyy1(m)));
        else  klyy(m)=0;
        end
        if    Ey1(m,N-1)~=0
              kryy1(m)=i*(log(Ey1(m,N)/Ey1(m,N-1)))/dy;
              kryy(m)=abs(real(kryy1(m)))-i*abs(imag(kryy1(m)));
        else  kryy(m)=0;
        end
    end
    
  for n=1:N
        if    Ey1(2,n)~=0
              klyx1(n)=i*(log(Ey1(1,n)/Ey1(2,n)))/dx;
              klyx(n)=abs(real(klyx1(n)))-i*abs(imag(klyx1(n)));
        else  klyx(n)=0;
        end
        if    Ey1(M-1,n)~=0
              kryx1(n)=i*(log(Ey1(M,n)/Ey1(M-1,n)))/dx;
              kryx(n)=abs(real(kryx1(n)))-i*abs(imag(kryx1(n)));
        else  kryx(n)=0;
        end
    end
%再计算后半步长
%l+0.5->l+1
%计算Ey2
for m=2:M-1
        a=2i*dx^2/dz+2-r*er(m,2)*ed(m,2)*dx^2;
         b22=a-exp(-i*klyy(m)*dy);
         y2(m,2)=H(m,2)/b22; 
        for   p=3:(N-2)
            a=2i*dx^2/dz+2-r*er(m,p)*ed(m,p)*dx^2;
            c22(p)=-1/b22;
            b22=a+c22(p); 
            y2(m,p)=(H(m,p)+y2(m,p-1))/b22;
        end  
        a=2i*dx^2/dz+2-r*er(m,N-1)*ed(m,N-1)*dx^2;
         c22(N-1)=-1/b22;
         b22=a-exp(-i*kryy(m)*dy)+c22(N-1);
         y2(m,N-1)= (H(m,N-1)+y2(m,N-2))/b22;
          Ey2(m,N-1)=y2(m,N-1);
        for v=(N-2):-1:2
            Ey2(m,v)=y2(m,v)-c22(v+1)*Ey2(m,v+1);
        end   
            Ey2(m,1)=Ey2(m,2)*exp(-i*klyy(m)*dy);
            Ey2(m,N)=Ey2(m,N-1)*exp(-i*kryy(m)*dy);      
   end 
%计算Ey2边界
for n=2:N-1
    Ey2(1,n)=Ey2(2,n)*exp(-i*klyx(n)*dx);
    Ey2(M,n)=Ey2(M-1,n)*exp(-i*kryx(n)*dx);       
end  
    Ey2(1,1)=Ey2(2,1)*exp(-i*klyx(1)*dx);
    Ey2(1,N)=Ey2(2,N)*exp(-i*klyx(N)*dx);
    Ey2(M,1)=Ey2(M,2)*exp(-i*klyy(M)*dy);
    Ey2(M,N)=Ey2(M,N-1)*exp(-i*kryy(M)*dy);

%计算Ex2
    for n=2:N-1
    for m=2:M-1
         c=2i*(dy^2)/dz-2+(1-r)*er(m,n)*ed(m,n)*(dy^2);
         b2=er(m,n)*exy(m,n)*(dy^2);
        G(m,n)=(Ex1(m+1,n)+c*Ex1(m,n)+Ex1(m-1,n))/d+b2*Ey2(m,n);
    end
end 
for m=2:M-1
     a=2i*dx^2/dz+2-r*er(m,2)*ed(m,2)*dx^2;
         b21=a-exp(-i*klxy(m)*dy);
         x2(m,2)=G(m,2)/b21;      
        for   p=3:(N-2)
           a=2i*dx^2/dz+2-r*er(m,p)*ed(m,p)*dx^2;
            c21(p)=-1/b21;
            b21=a+c21(p); 
            x2(m,p)=(G(m,p)+x2(m,p-1))/b21;
        end  
        a=2i*dx^2/dz+2-r*er(m,N-1)*ed(m,N-1)*dx^2;
         c21(N-1)=-1/b21;
         b21=a-exp(-i*krxy(m)*dy)+c21(N-1);
         x2(m,N-1)=(G(m,N-1)+x2(m,N-2))/b21;
          Ex2(m,N-1)=x2(m,N-1);
        for v=(M-2):-1:2
            Ex2(m,v)=x2(m,v)-c21(v+1)*Ex2(m,v+1);
        end   
            Ex2(m,1)=Ex2(m,2)*exp(-i*klxy(m)*dy);
            Ex2(m,N)=Ex2(m,N-1)*exp(-i*krxy(m)*dy);      
   end   
%计算Ex2边界
for n=2:N-1
    Ex2(1,n)=Ex2(2,n)*exp(-i*klxx(n)*dx);
    Ex2(M,n)=Ex2(M-1,n)*exp(-i*krxx(n)*dx);      
end  
    Ex2(1,1)=Ex2(2,1)*exp(-i*klxx(1)*dx);
    Ex2(1,N)=Ex2(2,N)*exp(-i*klxx(N)*dx);
    Ex2(M,1)=Ex2(M,2)*exp(-i*klxy(M)*dy);
    Ex2(M,N)=Ex2(M,N-1)*exp(-i*krxy(M)*dy);
  
%计算E1(:,:)     
for m=2:M-1
    for n=2:N-1
        c=2i*dy^2/dz-2+(1-r)*er(m,n)*ed(m,n)*(dy^2);
           b1=er(m,n)*exy(m,n)*(dx^2);
        E1(m,n)=(Ex2(m,n+1)+c*Ex2(m,n)+Ex2(m,n-1))*d+b1*Ey2(m,n);
    end
end

%x,y的边界   
%x的边界
    for m=1:M
        if    Ex2(m,2)~=0
              klxy2(m)=i*(log(Ex2(m,1)/Ex2(m,2)))/dy;
              klxy(m)=abs(real(klxy2(m)))-i*abs(imag(klxy2(m)));
        else  klxy(m)=0;
        end
        if    Ex2(m,N-1)~=0
              krxy2(m)=i*(log(Ex2(m,N)/Ex2(m,N-1)))/dy;
              krxy(m)=abs(real(krxy2(m)))-i*abs(imag(krxy2(m)));
        else  krxy(m)=0;
        end
    end
  for n=1:N
        if    Ex2(2,n)~=0
              klxx2(n)=i*(log(Ex2(1,n)/Ex2(2,n)))/dx;
              klxx(n)=abs(real(klxx2(n)))-i*abs(imag(klxx2(n)));
        else  klxx(n)=0;
        end
        if    Ex2(M-1,n)~=0
              krxx2(n)=i*(log(Ex2(M,n)/Ex2(M-1,n)))/dx;
              krxx(n)=abs(real(krxx2(n)))-i*abs(imag(krxx2(n)));
        else  krxx(n)=0;
        end
    end
%y的边界     
   for m=1:M
        if    Ey2(m,2)~=0
              klyy2(m)=i*(log(Ey2(m,1)/Ey2(m,2)))/dy;
              klyy(m)=abs(real(klyy2(m)))-i*abs(imag(klyy2(m)));
        else  klyy(m)=0;
        end
        if    Ey2(m,N-1)~=0
              kryy2(m)=i*(log(Ey2(m,N)/Ey2(m,N-1)))/dy;
              kryy(m)=abs(real(kryy2(m)))-i*abs(imag(kryy2(m)));
        else  kryy(m)=0;
        end
    end
    
  for n=1:N
        if    Ey2(2,n)~=0
              klyx2(n)=i*(log(Ey2(1,n)/Ey2(2,n)))/dx;
              klyx(n)=abs(real(klyx2(n)))-i*abs(imag(klyx2(n)));
        else  klyx(n)=0;
        end
        if    Ey2(M-1,n)~=0
              kryx2(n)=i*(log(Ey2(M,n)/Ey2(M-1,n)))/dx;
              kryx(n)=abs(real(kryx2(n)))-i*abs(imag(kryx2(n)));
        else  kryx(n)=0;
        end
    end
    
    for ii = 1:N
        for jj = 1:M
            Power1(l) = Power1(l) + abs(Ex2(ii, jj))^2;
        end
    end
    
    for ii = 1:N
        for jj = 1:M
            Power2(l) = Power2(l) + abs(Ey2(ii, jj))^2;
        end
    end
    subplot(2,1,1);
    mesh( abs(Ex2(:, :)) )
    axis([1 M 1 N 0 1])
    view(45,45)
    subplot(2,1,2);
    mesh( abs(Ey2(:, :)) )
    axis([1 M 1 N 0 1])
    view(45,45)
    pause(1e-200)
   
end  
%save paper_L=0.2_150X150;

figure
 xaxis=dz*6.6/k0:dz*6.6/k0:L*dz*6.6/k0;
plot(xaxis,Power1/max(Power1),'-b')
grid on
hold on

plot(xaxis,Power2/max(Power1),'--g')
Power = Power1/max(Power1) + Power2/max(Power1);
plot(xaxis,Power,'-r')
xlabel('传输距离','FontSize',12,'FontWeight','bold');
ylabel('归一化功率','FontSize',12,'FontWeight','bold');