%%%%%%%%三維計算%%%%%%%%%%%%
clc;
clear;
%%%%%%%%%%%%%%%%%
IE=40;JE=40;KE=40;
ic=IE/2;
jc=JE/2;
kc=KE/2;
ddx=0.01;
dt=ddx/6e8;
epsz=8.8e-12;
%%%%%%%%%%%%%
ex=zeros(IE,JE,KE);
ey=zeros(IE,JE,KE);
ez=zeros(IE,JE,KE);
dx=zeros(IE,JE,KE);
dy=zeros(IE,JE,KE);
dz=zeros(IE,JE,KE);
hx=zeros(IE,JE,KE);
hy=zeros(IE,JE,KE);
hz=zeros(IE,JE,KE);
gax=ones(IE,JE,KE);
gay=ones(IE,JE,KE);
gaz=ones(IE,JE,KE);
%%%%%%%%%%%%%
for k=12:30
    gaz(ic,jc,k)=0;
end
gaz(ic,jc,kc)=0;
%%%%%%%%
t0=20;
spread=6;
%%%%%%%
for T=1:300
    %%%%%%calculate the Dx field%%%%%%%
    for k=2:KE
        for j=2:JE
            for i=2:IE
                dx(i,j,k)=dx(i,j,k)+0.5*(hz(i,j,k)-hz(i,j,k-1)...
                    -hy(i,j,k)+hy(i-1,j,k));
            end
        end
    end
    %%%%%%calculate the Dy field%%%%%%%
    for k=2:KE
        for j=2:JE
            for i=2:IE
                dy(i,j,k)=dy(i,j,k)+0.5*(hx(i,j,k)-hx(i,j,k-1)...
                    -hz(i,j,k)+hz(i-1,j,k));
            end
        end
    end
    %%%%%%%%calculate the Dz field%%%%%%%%
    for k=2:KE
        for j=2:JE
            for i=2:IE
                dz(i,j,k)=dz(i,j,k)+0.5*(hy(i,j,k)-hy(i-1,j,k)...
                    -hx(i,j,k)+hx(i,j-1,k));
            end
        end
    end
    %%%%%%%%source%%%%%%%
    pulse=exp(-0.00005*((t0-T)/spread)^2);
    dz(ic,jc,kc)=pulse;
    %%%%%%%%%%%%calculate the E from D field%%%%%
    for k=2:KE-1
        for j=2:JE-1
            for i=2:IE-1
                ex(i,j,k)=gax(i,j,k)*dx(i,j,k);
                ey(i,j,k)=gay(i,j,k)*dy(i,j,k);
                ez(i,j,k)=gaz(i,j,k)*dz(i,j,k);
            end
        end
    end
    %%%%%%%%%%%calculate the Hx field%%%%%%%%
    for k=2:KE-1
        for j=2:JE-1
            for i=2:IE-1
                hx(i,j,k)=hx(i,j,k)+0.5*(ey(i,j,k+1)-ey(i,j,k)...
                    -ez(i,j+1,k)+ez(i,j,k));
            end
        end
    end
    %%%%%%%%%calculate the Hy field%%%%%%%%%%
    for k=2:KE-1
        for j=2:JE-1
            for i=2:IE-1
                hy(i,j,k)=hy(i,j,k)+0.5*(ez(i+1,j,k)-ez(i,j,k)...
                    -ex(i,j,k+1)+ex(i,j,k));
            end
        end
    end
    %%%%%%%%calculate the Hz field%%%%%%%%%
    for k=2:KE
        for j=2:JE-1
            for i=2:IE-1
                hz(i,j,k)=hz(i,j,k)+0.5*(ex(i,j+1,k)-ex(i,j,k)...
                    -ey(i+1,j,k)+ey(i,j,k));
            end
        end
    end
    %%%%%%%%%
 end