fdtd.m

一个参照葛德彪书编写的3D FDTD源文件

    Ex(Ismin+1:Ismax-1,Ismin,Ismax-1)=Ab1x(Ismin+1:Ismax-1,3,Ismax-1)-(Ex(Ismin+1:Ismax-1,Ismin+1,Ismax-1)-Ab1x(Ismin+1:Ismax-1,1,Ismax-1))/3;
    Ez(Ismin+1,Ismin,Ismin+1:Ismax-1)=Ab1z(Ismin+1,3,Ismin+1:Ismax-1)-(Ez(Ismin+1,Ismin+1,Ismin+1:Ismax-1)-Ab1z(Ismin+1,1,Ismin+1:Ismax-1))/3;
    Ex(Ismin+1,Ismin,Ismin+1:Ismax-1)=Ab1x(Ismin+1,3,Ismin+1:Ismax-1)-(Ex(Ismin+1,Ismin+1,Ismin+1:Ismax-1)-Ab1x(Ismin+1,1,Ismin+1:Ismax-1))/3;
    Ez(Ismax-1,Ismin,Ismin+1:Ismax-1)=Ab1z(Ismax-1,3,Ismin+1:Ismax-1)-(Ez(Ismax-1,Ismin+1,Ismin+1:Ismax-1)-Ab1z(Ismax-1,1,Ismin+1:Ismax-1))/3;
    Ex(Ismax-1,Ismin,Ismin+1:Ismax-1)=Ab1x(Ismax-1,3,Ismin+1:Ismax-1)-(Ex(Ismax-1,Ismin+1,Ismin+1:Ismax-1)-Ab1x(Ismax-1,1,Ismin+1:Ismax-1))/3;
    Ab1z(:,2,:)=Ab1z(:,1,:);         Ab1x(:,2,:)=Ab1x(:,1,:);%存Pn-1
    Ab1z(:,4,:)=Ab1z(:,3,:);         Ab1x(:,4,:)=Ab1x(:,3,:);%存Qn-1
    Ab1z(:,1,:)=Ez(:,Ismin,:);       Ab1x(:,1,:)=Ex(:,Ismin,:);%存Pn
    Ab1z(:,3,:)=Ez(:,Ismin+1,:);     Ab1x(:,3,:)=Ex(:,Ismin+1,:);%存Qn
    %---------------面2(y=Ismax,吸收Ez、Ex)-------------------------
    for i=Ismin+2:Ismax-2
        for j=Ismin+2:Ismax-2
            %--Ez--
            T1=Ez(i,Ismax-1,j)+Ab2z(i,2,j);
            T2=Ab2z(i,3,j)+Ab2z(i,1,j);
            T3=Ab2z(i+1,1,j)+Ab2z(i-1,1,j)+Ab2z(i,1,j+1)+Ab2z(i,1,j-1) + Ab2z(i+1,3,j)+Ab2z(i-1,3,j)+Ab2z(i,3,j+1)+Ab2z(i,3,j-1);
            Ez(i,Ismax,j)=-Ab2z(i,4,j)-T1/3+T2+T3/12;
            %--Ex--
            T1=Ex(i,Ismax-1,j)+Ab2x(i,2,j);
            T2=Ab2x(i,3,j)+Ab2x(i,1,j);
            T3=Ab2x(i+1,1,j)+Ab2x(i-1,1,j)+Ab2x(i,1,j+1)+Ab2x(i,1,j-1) + Ab2x(i+1,3,j)+Ab2x(i-1,3,j)+Ab2x(i,3,j+1)+Ab2x(i,3,j-1);
            Ex(i,Ismax,j)=-Ab2x(i,4,j)-T1/3+T2+T3/12;
        end
    end
    Ez(Ismin+1:Ismax-1,Ismax,Ismin+1)=Ab2z(Ismin+1:Ismax-1,3,Ismin+1)-(Ez(Ismin+1:Ismax-1,Ismax-1,Ismin+1)-Ab2z(Ismin+1:Ismax-1,1,Ismin+1))/3;
    Ex(Ismin+1:Ismax-1,Ismax,Ismin+1)=Ab2x(Ismin+1:Ismax-1,3,Ismin+1)-(Ex(Ismin+1:Ismax-1,Ismax-1,Ismin+1)-Ab2x(Ismin+1:Ismax-1,1,Ismin+1))/3;
    Ez(Ismin+1:Ismax-1,Ismax,Ismax-1)=Ab2z(Ismin+1:Ismax-1,3,Ismax-1)-(Ez(Ismin+1:Ismax-1,Ismax-1,Ismax-1)-Ab2z(Ismin+1:Ismax-1,1,Ismax-1))/3;
    Ex(Ismin+1:Ismax-1,Ismax,Ismax-1)=Ab2x(Ismin+1:Ismax-1,3,Ismax-1)-(Ex(Ismin+1:Ismax-1,Ismax-1,Ismax-1)-Ab2x(Ismin+1:Ismax-1,1,Ismax-1))/3;
    Ez(Ismin+1,Ismax,Ismin+1:Ismax-1)=Ab2z(Ismin+1,3,Ismin+1:Ismax-1)-(Ez(Ismin+1,Ismax-1,Ismin+1:Ismax-1)-Ab2z(Ismin+1,1,Ismin+1:Ismax-1))/3;
    Ex(Ismin+1,Ismax,Ismin+1:Ismax-1)=Ab2x(Ismin+1,3,Ismin+1:Ismax-1)-(Ex(Ismin+1,Ismax-1,Ismin+1:Ismax-1)-Ab2x(Ismin+1,1,Ismin+1:Ismax-1))/3;
    Ez(Ismax-1,Ismax,Ismin+1:Ismax-1)=Ab2z(Ismax-1,3,Ismin+1:Ismax-1)-(Ez(Ismax-1,Ismax-1,Ismin+1:Ismax-1)-Ab2z(Ismax-1,1,Ismin+1:Ismax-1))/3;
    Ex(Ismax-1,Ismax,Ismin+1:Ismax-1)=Ab2x(Ismax-1,3,Ismin+1:Ismax-1)-(Ex(Ismax-1,Ismax-1,Ismin+1:Ismax-1)-Ab2x(Ismax-1,1,Ismin+1:Ismax-1))/3;
    Ab2z(:,2,:)=Ab2z(:,1,:);         Ab2x(:,2,:)=Ab2x(:,1,:);%存Pn-1
    Ab2z(:,4,:)=Ab2z(:,3,:);         Ab2x(:,4,:)=Ab2x(:,3,:);%存Qn-1
    Ab2z(:,1,:)=Ez(:,Ismax,:);       Ab2x(:,1,:)=Ex(:,Ismax,:);%存Pn
    Ab2z(:,3,:)=Ez(:,Ismax-1,:);     Ab2x(:,3,:)=Ex(:,Ismax-1,:);%存Qn
    %---------------面3(z=Ismin,吸收Ey、Ex)----------------------
    for i=Ismin+2:Ismax-2
        for j=Ismin+2:Ismax-2
            %--Ey--
            T1=Ey(i,j,Ismin+1)+Ab3y(i,j,2);
            T2=Ab3y(i,j,3)+Ab3y(i,j,1);
            T3=Ab3y(i+1,j,1)+Ab3y(i-1,j,1)+Ab3y(i,j+1,1)+Ab3y(i,j-1,1) + Ab3y(i+1,j,3)+Ab3y(i-1,j,3)+Ab3y(i,j+1,3)+Ab3y(i,j-1,3);
            Ey(i,j,Ismin)=-Ab3y(i,j,4)-T1/3+T2+T3/12;
            %--Ex--
            T1=Ex(i,j,Ismin+1)+Ab3x(i,j,2);
            T2=Ab3x(i,j,3)+Ab3x(i,j,1);
            T3=Ab3x(i+1,j,1)+Ab3x(i-1,j,1)+Ab3x(i,j+1,1)+Ab3x(i,j-1,1) + Ab3x(i+1,j,3)+Ab3x(i-1,j,3)+Ab3x(i,j+1,3)+Ab3x(i,j-1,3);
            Ex(i,j,Ismin)=-Ab3x(i,j,4)-T1/3+T2+T3/12;
        end
    end
    Ey(Ismin+1:Ismax-1,Ismin+1,Ismin)=Ab3y(Ismin+1:Ismax-1,Ismin+1,3)-(Ey(Ismin+1:Ismax-1,Ismin+1,Ismin+1)-Ab3y(Ismin+1:Ismax-1,Ismin+1,1))/3;
    Ex(Ismin+1:Ismax-1,Ismin+1,Ismin)=Ab3x(Ismin+1:Ismax-1,Ismin+1,3)-(Ex(Ismin+1:Ismax-1,Ismin+1,Ismin+1)-Ab3x(Ismin+1:Ismax-1,Ismin+1,1))/3;
    Ey(Ismin+1:Ismax-1,Ismax-1,Ismin)=Ab3y(Ismin+1:Ismax-1,Ismax-1,3)-(Ey(Ismin+1:Ismax-1,Ismax-1,Ismin+1)-Ab3y(Ismin+1:Ismax-1,Ismax-1,1))/3;
    Ex(Ismin+1:Ismax-1,Ismax-1,Ismin)=Ab3x(Ismin+1:Ismax-1,Ismax-1,3)-(Ex(Ismin+1:Ismax-1,Ismax-1,Ismin+1)-Ab3x(Ismin+1:Ismax-1,Ismax-1,1))/3;
    Ey(Ismin+1,Ismin+1:Ismax-1,Ismin)=Ab3y(Ismin+1,Ismin+1:Ismax-1,3)-(Ey(Ismin+1,Ismin+1:Ismax-1,Ismin+1)-Ab3y(Ismin+1,Ismin+1:Ismax-1,1))/3;
    Ex(Ismin+1,Ismin+1:Ismax-1,Ismin)=Ab3x(Ismin+1,Ismin+1:Ismax-1,3)-(Ex(Ismin+1,Ismin+1:Ismax-1,Ismin+1)-Ab3x(Ismin+1,Ismin+1:Ismax-1,1))/3;
    Ey(Ismax-1,Ismin+1:Ismax-1,Ismin)=Ab3y(Ismax-1,Ismin+1:Ismax-1,3)-(Ey(Ismax-1,Ismin+1:Ismax-1,Ismin+1)-Ab3y(Ismax-1,Ismin+1:Ismax-1,1))/3;
    Ex(Ismax-1,Ismin+1:Ismax-1,Ismin)=Ab3x(Ismax-1,Ismin+1:Ismax-1,3)-(Ex(Ismax-1,Ismin+1:Ismax-1,Ismin+1)-Ab3x(Ismax-1,Ismin+1:Ismax-1,1))/3;
    Ab3y(:,:,2)=Ab3y(:,:,1);         Ab3x(:,:,2)=Ab3x(:,:,1);%存Pn-1
    Ab3y(:,:,4)=Ab3y(:,:,3);         Ab3x(:,:,4)=Ab3x(:,:,3);%存Qn-1
    Ab3y(:,:,1)=Ey(:,:,Ismin);       Ab3x(:,:,1)=Ex(:,:,Ismin);%存Pn
    Ab3y(:,:,3)=Ey(:,:,Ismin+1);     Ab3x(:,:,3)=Ex(:,:,Ismin+1);%存Qn
    %---------------面4(z=Ismax,吸收Ey、Ex)----------------------
    for i=Ismin+2:Ismax-2
        for j=Ismin+2:Ismax-2
            %--Ey--
            T1=Ey(i,j,Ismax-1)+Ab4y(i,j,2);
            T2=Ab4y(i,j,3)+Ab4y(i,j,1);
            T3=Ab4y(i+1,j,1)+Ab4y(i-1,j,1)+Ab4y(i,j+1,1)+Ab4y(i,j-1,1) + Ab4y(i+1,j,3)+Ab4y(i-1,j,3)+Ab4y(i,j+1,3)+Ab4y(i,j-1,3);
            Ey(i,j,Ismax)=-Ab4y(i,j,4)-T1/3+T2+T3/12;
            %--Ex--
            T1=Ex(i,j,Ismax-1)+Ab4x(i,j,2);
            T2=Ab4x(i,j,3)+Ab4x(i,j,1);
            T3=Ab4x(i+1,j,1)+Ab4x(i-1,j,1)+Ab4x(i,j+1,1)+Ab4x(i,j-1,1) + Ab4x(i+1,j,3)+Ab4x(i-1,j,3)+Ab4x(i,j+1,3)+Ab4x(i,j-1,3);
            Ex(i,j,Ismax)=-Ab4x(i,j,4)-T1/3+T2+T3/12;
        end
    end
    Ey(Ismin+1:Ismax-1,Ismin+1,Ismax)=Ab4y(Ismin+1:Ismax-1,Ismin+1,3)-(Ey(Ismin+1:Ismax-1,Ismin+1,Ismax-1)-Ab4y(Ismin+1:Ismax-1,Ismin+1,1))/3;
    Ex(Ismin+1:Ismax-1,Ismin+1,Ismax)=Ab4x(Ismin+1:Ismax-1,Ismin+1,3)-(Ex(Ismin+1:Ismax-1,Ismin+1,Ismax-1)-Ab4x(Ismin+1:Ismax-1,Ismin+1,1))/3;
    Ey(Ismin+1:Ismax-1,Ismax-1,Ismax)=Ab4y(Ismin+1:Ismax-1,Ismax-1,3)-(Ey(Ismin+1:Ismax-1,Ismax-1,Ismax-1)-Ab4y(Ismin+1:Ismax-1,Ismax-1,1))/3;
    Ex(Ismin+1:Ismax-1,Ismax-1,Ismax)=Ab4x(Ismin+1:Ismax-1,Ismax-1,3)-(Ex(Ismin+1:Ismax-1,Ismax-1,Ismax-1)-Ab4x(Ismin+1:Ismax-1,Ismax-1,1))/3;
    Ey(Ismin+1,Ismin+1:Ismax-1,Ismax)=Ab4y(Ismin+1,Ismin+1:Ismax-1,3)-(Ey(Ismin+1,Ismin+1:Ismax-1,Ismax-1)-Ab4y(Ismin+1,Ismin+1:Ismax-1,1))/3;
    Ex(Ismin+1,Ismin+1:Ismax-1,Ismax)=Ab4x(Ismin+1,Ismin+1:Ismax-1,3)-(Ex(Ismin+1,Ismin+1:Ismax-1,Ismax-1)-Ab4x(Ismin+1,Ismin+1:Ismax-1,1))/3;
    Ey(Ismax-1,Ismin+1:Ismax-1,Ismax)=Ab4y(Ismax-1,Ismin+1:Ismax-1,3)-(Ey(Ismax-1,Ismin+1:Ismax-1,Ismax-1)-Ab4y(Ismax-1,Ismin+1:Ismax-1,1))/3;
    Ex(Ismax-1,Ismin+1:Ismax-1,Ismax)=Ab4x(Ismax-1,Ismin+1:Ismax-1,3)-(Ex(Ismax-1,Ismin+1:Ismax-1,Ismax-1)-Ab4x(Ismax-1,Ismin+1:Ismax-1,1))/3;
    Ab4y(:,:,2)=Ab4y(:,:,1);         Ab4x(:,:,2)=Ab4x(:,:,1);%存Pn-1
    Ab4y(:,:,4)=Ab4y(:,:,3);         Ab4x(:,:,4)=Ab4x(:,:,3);%存Qn-1
    Ab4y(:,:,1)=Ey(:,:,Ismax);       Ab4x(:,:,1)=Ex(:,:,Ismax);%存Pn
    Ab4y(:,:,3)=Ey(:,:,Ismax-1);     Ab4x(:,:,3)=Ex(:,:,Ismax-1);%存Qn
    %---------------面5(x=Ismin,吸收Ey、Ez)----------------------
    for i=Ismin+2:Ismax-2
        for j=Ismin+2:Ismax-2
            %--Ey--
            T1=Ey(Ismin+1,i,j)+Ab5y(2,i,j);
            T2=Ab5y(3,i,j)+Ab5y(1,i,j);
            T3=Ab5y(1,i+1,j)+Ab5y(1,i-1,j)+Ab5y(1,i,j+1)+Ab5y(1,i,j-1) + Ab5y(3,i+1,j)+Ab5y(3,i-1,j)+Ab5y(3,i,j+1)+Ab5y(3,i,j-1);
            Ey(Ismin,i,j)=-Ab5y(4,i,j)-T1/3+T2+T3/12;
            %--Ez--
            T1=Ez(Ismin+1,i,j)+Ab5z(2,i,j);
            T2=Ab5z(3,i,j)+Ab5z(1,i,j);
            T3=Ab5z(1,i+1,j)+Ab5z(1,i-1,j)+Ab5z(1,i,j+1)+Ab5z(1,i,j-1) + Ab5z(3,i+1,j)+Ab5z(3,i-1,j)+Ab5z(3,i,j+1)+Ab5z(3,i,j-1);
            Ez(Ismin,i,j)=-Ab5z(4,i,j)-T1/3+T2+T3/12;
        end
    end
    Ey(Ismin,Ismin+1:Ismax-1,Ismin+1)=Ab5y(3,Ismin+1:Ismax-1,Ismin+1)-(Ey(Ismin+1,Ismin+1:Ismax-1,Ismin+1)-Ab5y(1,Ismin+1:Ismax-1,Ismin+1))/3;
    Ez(Ismin,Ismin+1:Ismax-1,Ismin+1)=Ab5z(3,Ismin+1:Ismax-1,Ismin+1)-(Ez(Ismin+1,Ismin+1:Ismax-1,Ismin+1)-Ab5z(1,Ismin+1:Ismax-1,Ismin+1))/3;
    Ey(Ismin,Ismin+1:Ismax-1,Ismax-1)=Ab5y(3,Ismin+1:Ismax-1,Ismax-1)-(Ey(Ismin+1,Ismin+1:Ismax-1,Ismax-1)-Ab5y(1,Ismin+1:Ismax-1,Ismax-1))/3;
    Ez(Ismin,Ismin+1:Ismax-1,Ismax-1)=Ab5z(3,Ismin+1:Ismax-1,Ismax-1)-(Ez(Ismin+1,Ismin+1:Ismax-1,Ismax-1)-Ab5z(1,Ismin+1:Ismax-1,Ismax-1))/3;
    Ey(Ismin,Ismin+1,Ismin+1:Ismax-1)=Ab5y(3,Ismin+1,Ismin+1:Ismax-1)-(Ey(Ismin+1,Ismin+1,Ismin+1:Ismax-1)-Ab5y(1,Ismin+1,Ismin+1:Ismax-1))/3;
    Ez(Ismin,Ismin+1,Ismin+1:Ismax-1)=Ab5z(3,Ismin+1,Ismin+1:Ismax-1)-(Ez(Ismin+1,Ismin+1,Ismin+1:Ismax-1)-Ab5z(1,Ismin+1,Ismin+1:Ismax-1))/3;
    Ey(Ismin,Ismax-1,Ismin+1:Ismax-1)=Ab5y(3,Ismax-1,Ismin+1:Ismax-1)-(Ey(Ismin+1,Ismax-1,Ismin+1:Ismax-1)-Ab5y(1,Ismax-1,Ismin+1:Ismax-1))/3;
    Ez(Ismin,Ismax-1,Ismin+1:Ismax-1)=Ab5z(3,Ismax-1,Ismin+1:Ismax-1)-(Ez(Ismin+1,Ismax-1,Ismin+1:Ismax-1)-Ab5z(1,Ismax-1,Ismin+1:Ismax-1))/3;
    Ab5y(2,:,:)=Ab5y(1,:,:);         Ab5z(2,:,:)=Ab5z(1,:,:);%存Pn-1
    Ab5y(4,:,:)=Ab5y(3,:,:);         Ab5z(4,:,:)=Ab5z(3,:,:);%存Qn-1
    Ab5y(1,:,:)=Ey(Ismin,:,:);       Ab5z(1,:,:)=Ez(Ismin,:,:);%存Pn
    Ab5y(3,:,:)=Ey(Ismin+1,:,:);     Ab5z(3,:,:)=Ez(Ismin+1,:,:);%存Qn
    %---------------面6(x=Ismax,吸收Ey、Ez)----------------------
    for i=Ismin+2:Ismax-2
        for j=Ismin+2:Ismax-2
            %--Ey--
            T1=Ey(Ismax-1,i,j)+Ab6y(2,i,j);
            T2=Ab6y(3,i,j)+Ab6y(1,i,j);
            T3=Ab6y(1,i+1,j)+Ab6y(1,i-1,j)+Ab6y(1,i,j+1)+Ab6y(1,i,j-1) + Ab6y(3,i+1,j)+Ab6y(3,i-1,j)+Ab6y(3,i,j+1)+Ab6y(3,i,j-1);
            Ey(Ismax,i,j)=-Ab6y(4,i,j)-T1/3+T2+T3/12;
            %--Ez--
            T1=Ez(Ismax-1,i,j)+Ab6z(2,i,j);
            T2=Ab6z(3,i,j)+Ab6z(1,i,j);
            T3=Ab6z(1,i+1,j)+Ab6z(1,i-1,j)+Ab6z(1,i,j+1)+Ab6z(1,i,j-1) + Ab6z(3,i+1,j)+Ab6z(3,i-1,j)+Ab6z(3,i,j+1)+Ab6z(3,i,j-1);
            Ez(Ismax,i,j)=-Ab6z(4,i,j)-T1/3+T2+T3/12;
        end
    end
    Ey(Ismax,Ismin+1:Ismax-1,Ismin+1)=Ab6y(3,Ismin+1:Ismax-1,Ismin+1)-(Ey(Ismax-1,Ismin+1:Ismax-1,Ismin+1)-Ab6y(1,Ismin+1:Ismax-1,Ismin+1))/3;
    Ez(Ismax,Ismin+1:Ismax-1,Ismin+1)=Ab6z(3,Ismin+1:Ismax-1,Ismin+1)-(Ez(Ismax-1,Ismin+1:Ismax-1,Ismin+1)-Ab6z(1,Ismin+1:Ismax-1,Ismin+1))/3;
    Ey(Ismax,Ismin+1:Ismax-1,Ismax-1)=Ab6y(3,Ismin+1:Ismax-1,Ismax-1)-(Ey(Ismax-1,Ismin+1:Ismax-1,Ismax-1)-Ab6y(1,Ismin+1:Ismax-1,Ismax-1))/3;
    Ez(Ismax,Ismin+1:Ismax-1,Ismax-1)=Ab6z(3,Ismin+1:Ismax-1,Ismax-1)-(Ez(Ismax-1,Ismin+1:Ismax-1,Ismax-1)-Ab6z(1,Ismin+1:Ismax-1,Ismax-1))/3;
    Ey(Ismax,Ismin+1,Ismin+1:Ismax-1)=Ab6y(3,Ismin+1,Ismin+1:Ismax-1)-(Ey(Ismax-1,Ismin+1,Ismin+1:Ismax-1)-Ab6y(1,Ismin+1,Ismin+1:Ismax-1))/3;
    Ez(Ismax,Ismin+1,Ismin+1:Ismax-1)=Ab6z(3,Ismin+1,Ismin+1:Ismax-1)-(Ez(Ismax-1,Ismin+1,Ismin+1:Ismax-1)-Ab6z(1,Ismin+1,Ismin+1:Ismax-1))/3;
    Ey(Ismax,Ismax-1,Ismin+1:Ismax-1)=Ab6y(3,Ismax-1,Ismin+1:Ismax-1)-(Ey(Ismax-1,Ismax-1,Ismin+1:Ismax-1)-Ab6y(1,Ismax-1,Ismin+1:Ismax-1))/3;
    Ez(Ismax,Ismax-1,Ismin+1:Ismax-1)=Ab6z(3,Ismax-1,Ismin+1:Ismax-1)-(Ez(Ismax-1,Ismax-1,Ismin+1:Ismax-1)-Ab6z(1,Ismax-1,Ismin+1:Ismax-1))/3;
    Ab6y(2,:,:)=Ab6y(1,:,:);         Ab6z(2,:,:)=Ab6z(1,:,:);%存Pn-1
    Ab6y(4,:,:)=Ab6y(3,:,:);         Ab6z(4,:,:)=Ab6z(3,:,:);%存Qn-1
    Ab6y(1,:,:)=Ey(Ismax,:,:);       Ab6z(1,:,:)=Ez(Ismax,:,:);%存Pn
    Ab6y(3,:,:)=Ey(Ismax-1,:,:);     Ab6z(3,:,:)=Ez(Ismax-1,:,:);%存Qn
    %------------------------------------------------6个吸收面 end--------------------------------------------
    %------------------------------------------------end of E FDTD--------------------------------------------
    %
    %-------------------------------------------------start H FDTD--------------------------------------------
    for i=Ismin+1:Ismax-1
        for j=Ismin+1:Ismax-1
            for k=Ismin+1:Ismax-1
                T1=Ey(i,j,k+1)-Ey(i,j,k)-Ez(i,j+1,k)+Ez(i,j,k);
                Hx(i,j,k)=CP(obj(i+1,j,k),obj(i,j,k))*Hx(i,j,k)+CQ(obj(i+1,j,k),obj(i,j,k))*T1;
                T1=Ez(i+1,j,k)-Ez(i,j,k)-Ex(i,j,k+1)+Ex(i,j,k);
                Hy(i,j,k)=CP(obj(i,j+1,k),obj(i,j,k))*Hy(i,j,k)+CQ(obj(i,j+1,k),obj(i,j,k))*T1;
                T1=Ex(i,j+1,k)-Ex(i,j,k)-Ey(i+1,j,k)+Ey(i,j,k);
                Hz(i,j,k)=CP(obj(i,j,k+1),obj(i,j,k))*Hz(i,j,k)+CQ(obj(i,j,k+1),obj(i,j,k))*T1;
            end
        end
    end
    %----------start to add Hin wave---------
    %--面 1---
    Hx(Izmin+1:Izmax-1,Izmin,Izmin+1:Izmax-1)=Hx(Izmin+1:Izmax-1,Izmin,Izmin+1:Izmax-1)+FH*Ein1(Izmin+1:Izmax-1,1,Izmin+1:Izmax-1);
    %--面 2---
    Hx(Izmin+1:Izmax-1,Izmax,Izmin+1:Izmax-1)=Hx(Izmin+1:Izmax-1,Izmax,Izmin+1:Izmax-1)-FH*Ein2(Izmin+1:Izmax-1,1,Izmin+1:Izmax-1);
    %--面 5---
    Hy(Izmin,Izmin+1:Izmax-1,Izmin+1:Izmax-1)=Hy(Izmin,Izmin+1:Izmax-1,Izmin+1:Izmax-1)-FH*Ein5(1,Izmin+1:Izmax-1,Izmin+1:Izmax-1);
    %--面 6---
    Hy(Izmax,Izmin+1:Izmax-1,Izmin+1:Izmax-1)=Hy(Izmax,Izmin+1:Izmax-1,Izmin+1:Izmax-1)+FH*Ein6(1,Izmin+1:Izmax-1,Izmin+1:Izmax-1);
    %---------end of adding Hin wave---------
    %-------------------------------------------------end of H FDTD-------------------------------------------
    %
    %---------------------------start Far-zone trans&out--------------------------
    %-----------------------------end of Far-zone trans&out-----------------------
    
end
%======================Time loop end==========================