clear;
%**********************************************************
pi=3.1415926;
mu=4*pi*1e-7;yb=8.851*1e-12;c=1/sqrt(mu*yb);f=1e8;
k=2*pi*f*sqrt(mu*yb);Z=sqrt(mu/yb);
wl=c/f;dx=wl/40.0;dt=dx/c/2.0;TT=1.0/f/dt;
%**********************************************************
NX=60;	NY=60;NZ=60;AA=2;SS=5;TIME=TT*SS;           %外推邊界距離吸收邊界的距離為AA,時間循環為ss個周期
%**********************************************************     
hx=zeros(NX+1,NY,NZ);hy=zeros(NX,NY+1,NZ);hz=zeros(NX,NY,NZ+1);
ex=zeros(NX,NY+1,NZ+1);ey=zeros(NX+1,NY,NZ+1);ez=zeros(NX+1,NY+1,NZ);
%**********************************************************
waste=5/1000;
hx_1=zeros(NX+1,NY,NZ);hy_1=zeros(NX,NY+1,NZ);hz_1=zeros(NX,NY,NZ+1);                     %三個時刻的場（邊界問題）
hx_2=zeros(NX+1,NY,NZ);hy_2=zeros(NX,NY+1,NZ);hz_2=zeros(NX,NY,NZ+1);
hx_3=zeros(NX+1,NY,NZ);hy_3=zeros(NX,NY+1,NZ);hz_3=zeros(NX,NY,NZ+1);
ex_1=zeros(NX,NY+1,NZ+1);ey_1=zeros(NX+1,NY,NZ+1);ez_1=zeros(NX+1,NY+1,NZ);
ex_2=zeros(NX,NY+1,NZ+1);ey_2=zeros(NX+1,NY,NZ+1);ez_2=zeros(NX+1,NY+1,NZ);
ex_3=zeros(NX,NY+1,NZ+1);ey_3=zeros(NX+1,NY,NZ+1);ez_3=zeros(NX+1,NY+1,NZ);
%**********************************************************
midx=NX/2+1;midy=NY/2+1;midz=NZ/2+1;mx=fix(midx);my=fix(midy);mz=fix(midz);
SX=1+AA;SY=1+AA;SZ=1+AA;XX=NX+1-AA;YY=NY+1-AA;ZZ=NZ+1-AA;                  %各個方向外推邊界距離吸收邊界的距離
FMz_fx1=0;FMy_fx1=0;FMz_fy1=0;FMx_fy1=0;FMy_fz1=0;FMx_fz1=0;
FMz_fx2=0;FMy_fx2=0;FMz_fy2=0;FMx_fy2=0;FMy_fz2=0;FMx_fz2=0;
FJz_fx1=0;FJy_fx1=0;FJz_fy1=0;FJx_fy1=0;FJy_fz1=0;FJx_fz1=0;
FJz_fx2=0;FJy_fx2=0;FJz_fy2=0;FJx_fy2=0;FJy_fz2=0;FJx_fz2=0;
%**********************************************************
r0=1/3*dx;c1=dt/mu/(dx*dx-pi*r0*r0/4);c2=dx/2*log(dx/r0);ll=9;
ia=mx;ja=my;ka=mz;star=fix(mz-ll);endn=fix(mz+ll);
%**********************************************************
tic
figure(1);clf;
set(gcf,'Renderer','OpenGL');
for n=1:TIME
      %-------------------電磁場循環---------------------------
    hx_3=hx_2;hx_2=hx_1;hx_1=hx;hy_3=hy_2;hy_2=hy_1;hy_1=hy;hz_3=hz_2;hz_2=hz_1;hz_1=hz;
    hx=hx+dt/mu/dx*(ey(:,:,2:NZ+1)-ey(:,:,1:NZ)+ez(:,1:NY,:)-ez(:,2:NY+1,:));
    hy=hy+dt/mu/dx*(ez(2:NX+1,:,:)-ez(1:NX,:,:)+ex(:,:,1:NZ)-ex(:,:,2:NZ+1));
    hz=hz+dt/mu/dx*(ex(:,2:NY+1,:)-ex(:,1:NY,:)+ey(1:NX,:,:)-ey(2:NX+1,:,:));
    %*************************
%     hy(ia,ja,star:endn)=hy_1(ia,ja,star:endn)-dt/mu/dx*(ex(ia,ja,star+1:endn+1)-ex(ia,ja,star:endn))+dt/mu/dx*2/log(dx/r0)*ez(ia+1,ja,star:endn);
%     hx(ia,ja,star:endn)=hx_1(ia,ja,star:endn)+dt/mu/dx*(ey(ia,ja,star+1:endn+1)-ey(ia,ja,star:endn))-dt/mu/dx*2/log(dx/r0)*ez(ia,ja+1,star:endn);
%     hy(ia-1,ja,star:endn)=hy_1(ia-1,ja,star:endn)-dt/mu/dx*(ex(ia-1,ja,star+1:endn+1)-ex(ia-1,ja,star:endn))-dt/mu/dx*2/log(dx/r0)*ez(ia-1,ja,star:endn);
%     hx(ia,ja-1,star:endn)=hx_1(ia,ja-1,star:endn)+dt/mu/dx*(ey(ia,ja-1,star+1:endn+1)-ey(ia,ja-1,star:endn))+dt/mu/dx*2/log(dx/r0)*ez(ia,ja-1,star:endn);
%     %*************************
%     hy(ia,ja,ka)=hy_1(ia,ja,ka)-dt/mu/dx*(ex(ia,ja,ka+1)-ex(ia,ja,ka))+dt/mu/dx*2/log(dx/r0)*(ez(ia+1,ja,ka)-ez_1(ia,ja,ka));
%     hx(ia,ja,ka)=hx_1(ia,ja,ka)+dt/mu/dx*(ey(ia,ja,ka+1)-ey(ia,ja,ka))-dt/mu/dx*2/log(dx/r0)*(ez(ia,ja+1,ka)-ez_1(ia,ja,ka));
%     hy(ia-1,ja,ka)=hy_1(ia-1,ja,ka)-dt/mu/dx*(ex(ia-1,ja,ka+1)-ex(ia-1,ja,ka))-dt/mu/dx*2/log(dx/r0)*(ez(ia-1,ja,ka)-ez_1(ia,ja,ka));
%     hx(ia,ja-1,ka)=hx_1(ia,ja-1,ka)+dt/mu/dx*(ey(ia,ja-1,ka+1)-ey(ia,ja-1,ka))+dt/mu/dx*2/log(dx/r0)*(ez(ia,ja-1,ka)-ez_1(ia,ja,ka));
%    %*************************  
%     hz(ia,ja,star:endn+1)=hz_1(ia,ja,star:endn+1)+c1*(c2*ey(ia,ja,star:endn+1)-c2*ex(ia,ja,star:endn+1)+ex(ia,ja+1,star:endn+1)*dx-dx*ey(ia+1,ja,star:endn+1));
%     hz(ia,ja-1,star:endn+1)=hz_1(ia,ja-1,star:endn+1)+c1*(c2*ey(ia,ja-1,star:endn+1)+c2*ex(ia,ja,star:endn+1)-ex(ia,ja-1,star:endn+1)*dx-dx*ey(ia+1,ja-1,star:endn+1));
%     hz(ia-1,ja,star:endn+1)=hz_1(ia-1,ja,star:endn+1)+c1*(-c2*ey(ia,ja,star:endn+1)-c2*ex(ia-1,ja,star:endn+1)+ex(ia-1,ja+1,star:endn+1)*dx+dx*ey(ia-1,ja,star:endn+1));
%     hz(ia-1,ja-1,star:endn+1)=hz_1(ia-1,ja-1,star:endn+1)+c1*(-c2*ey(ia,ja-1,star:endn+1)+c2*ex(ia-1,ja,star:endn+1)-ex(ia-1,ja-1,star:endn+1)*dx+dx*ey(ia-1,ja-1,star:endn+1));
      
    ex_3=ex_2;ex_2=ex_1;ex_1=ex;ey_3=ey_2;ey_2=ey_1;ey_1=ey;ez_3=ez_2;ez_2=ez_1;ez_1=ez;
    ex(:,2:NY,2:NZ)=ex(:,2:NY,2:NZ)+dt/yb/dx*(hz(:,2:NY,2:NZ)-hz(:,1:NY-1,2:NZ)+hy(:,2:NY,1:NZ-1)-hy(:,2:NY,2:NZ));
    ey(2:NX,:,2:NZ)=ey(2:NX,:,2:NZ)+dt/yb/dx*(hx(2:NX,:,2:NZ)-hx(2:NX,:,1:NZ-1)+hz(1:NX-1,:,2:NZ)-hz(2:NX,:,2:NZ));
    ez(2:NX,2:NY,:)=ez(2:NX,2:NY,:)+dt/yb/dx*(hy(2:NX,2:NY,:)-hy(1:NX-1,2:NY,:)+hx(2:NX,1:NY-1,:)-hx(2:NX,2:NY,:));
    ez(ia,ja,ka)=ez(ia,ja,ka)+sin(2*pi*f*n*dt);
    %************吸收邊界條件*************
    ey(1,:,:)=+(3*ey_1(1+1,:,:)-3*ey_2(1+2,:,:)+ey_3(1+3,:,:));
    ey(NX+1,:,:)=+(3*ey_1(NX+1-1,:,:)-3*ey_2(NX+1-2,:,:)+ey_3(NX+1-3,:,:));
    ez(1,:,:)=+(3*ez_1(1+1,:,:)-3*ez_2(1+2,:,:)+ez_3(1+3,:,:));
    ez(NX+1,:,:)=+(3*ez_1(NX+1-1,:,:)-3*ez_2(NX+1-2,:,:)+ez_3(NX+1-3,:,:));
    ex(:,1,:)=+(3*ex_1(:,1+1,:)-3*ex_2(:,1+2,:)+ex_3(:,1+3,:));
    ex(:,NY+1,:)=+(3*ex_1(:,NY+1-1,:)-3*ex_2(:,NY+1-2,:)+ex_3(:,NY+1-3,:));
    ez(:,1,:)=+(3*ez_1(:,1+1,:)-3*ez_2(:,1+2,:)+ez_3(:,1+3,:));
    ez(:,NY+1,:)=+(3*ez_1(:,NY+1-1,:)-3*ez_2(:,NY+1-2,:)+ez_3(:,NY+1-3,:));
    ex(:,:,1)=+(3*ex_1(:,:,1+1)-3*ex_2(:,:,1+2)+ex_3(:,:,1+3));
    ex(:,:,NZ+1)=+(3*ex_1(:,:,NZ+1-1)-3*ex_2(:,:,NZ+1-2)+ex_3(:,:,NZ+1-3));
    ey(:,:,1)=+(3*ey_1(:,:,1+1)-3*ey_2(:,:,1+2)+ey_3(:,:,1+3));
    ey(:,:,NZ+1)=+(3*ey_1(:,:,NZ+1-1)-3*ey_2(:,:,NZ+1-2)+ey_3(:,:,NZ+1-3));
    %***********存儲**************
    if n>TT*(SS-1)        
        ey_fx1=(ey(SX,SY:YY-1,SZ:ZZ-1)+ey(SX,SY:YY-1,SZ+1:ZZ))/2;
        ez_fx1=(ez(SX,SY:YY-1,SZ:ZZ-1)+ez(SX,SY+1:YY,SZ:ZZ-1))/2;
        ex_fy1=(ex(SX:XX-1,SY,SZ:ZZ-1)+ex(SX:XX-1,SY,SZ+1:ZZ))/2;
        ez_fy1=(ez(SX:XX-1,SY,SZ:ZZ-1)+ez(SX+1:XX,SY,SZ:ZZ-1))/2;
        ex_fz1=(ex(SX:XX-1,SY:YY-1,SZ)+ex(SX:XX-1,SY+1:YY,SZ))/2;
        ey_fz1=(ey(SX:XX-1,SY:YY-1,SZ)+ey(SX+1:XX,SY:YY-1,SZ))/2;
        
        ey_fx2=(ey(XX,SY:YY-1,SZ:ZZ-1)+ey(XX,SY:YY-1,SZ+1:ZZ))/2;
        ez_fx2=(ez(XX,SY:YY-1,SZ:ZZ-1)+ez(XX,SY+1:YY,SZ:ZZ-1))/2;
        ex_fy2=(ex(SX:XX-1,YY,SZ:ZZ-1)+ex(SX:XX-1,YY,SZ+1:ZZ))/2;
        ez_fy2=(ez(SX:XX-1,YY,SZ:ZZ-1)+ez(SX+1:XX,YY,SZ:ZZ-1))/2;
        ex_fz2=(ex(SX:XX-1,SY:YY-1,ZZ)+ex(SX:XX-1,SY+1:YY,ZZ))/2;
        ey_fz2=(ey(SX:XX-1,SY:YY-1,ZZ)+ey(SX+1:XX,SY:YY-1,ZZ))/2;
        
        hy_fx1=(hy(SX,SY:YY-1,SZ:ZZ-1)+hy(SX-1,SY:YY-1,SZ:ZZ-1)+hy(SX,SY+1:YY,SZ:ZZ-1)+hy(SX-1,SY+1:YY,SZ:ZZ-1))/4;
        hz_fx1=(hz(SX,SY:YY-1,SZ:ZZ-1)+hz(SX-1,SY:YY-1,SZ:ZZ-1)+hz(SX,SY:YY-1,SZ+1:ZZ)+hz(SX-1,SY:YY-1,SZ+1:ZZ))/4;
        hx_fy1=(hx(SX:XX-1,SY,SZ:ZZ-1)+hx(SX:XX-1,SY-1,SZ:ZZ-1)+hx(SX+1:XX,SY,SZ:ZZ-1)+hx(SX+1:XX,SY-1,SZ:ZZ-1))/4;
        hz_fy1=(hz(SX:XX-1,SY,SZ:ZZ-1)+hz(SX:XX-1,SY-1,SZ:ZZ-1)+hz(SX:XX-1,SY,SZ+1:ZZ)+hz(SX:XX-1,SY-1,SZ+1:ZZ))/4;
        hx_fz1=(hx(SX:XX-1,SY:YY-1,SZ)+hx(SX:XX-1,SY:YY-1,SZ-1)+hx(SX+1:XX,SY:YY-1,SZ)+hx(SX+1:XX,SY:YY-1,SZ-1))/4;
        hy_fz1=(hy(SX:XX-1,SY:YY-1,SZ)+hy(SX:XX-1,SY:YY-1,SZ-1)+hy(SX:XX-1,SY+1:YY,SZ)+hy(SX:XX-1,SY+1:YY,SZ-1))/4;
        
        hy_fx2=(hy(XX,SY:YY-1,SZ:ZZ-1)+hy(XX-1,SY:YY-1,SZ:ZZ-1)+hy(XX,SY+1:YY,SZ:ZZ-1)+hy(XX-1,SY+1:YY,SZ:ZZ-1))/4;
        hz_fx2=(hz(XX,SY:YY-1,SZ:ZZ-1)+hz(XX-1,SY:YY-1,SZ:ZZ-1)+hz(XX,SY:YY-1,SZ+1:ZZ)+hz(XX-1,SY:YY-1,SZ+1:ZZ))/4;
        hx_fy2=(hx(SX:XX-1,YY,SZ:ZZ-1)+hx(SX:XX-1,YY-1,SZ:ZZ-1)+hx(SX+1:XX,YY,SZ:ZZ-1)+hx(SX+1:XX,YY-1,SZ:ZZ-1))/4;
        hz_fy2=(hz(SX:XX-1,YY,SZ:ZZ-1)+hz(SX:XX-1,YY-1,SZ:ZZ-1)+hz(SX:XX-1,YY,SZ+1:ZZ)+hz(SX:XX-1,YY-1,SZ+1:ZZ))/4;
        hx_fz2=(hx(SX:XX-1,SY:YY-1,ZZ)+hx(SX:XX-1,SY:YY-1,ZZ-1)+hx(SX+1:XX,SY:YY-1,ZZ)+hx(SX+1:XX,SY:YY-1,ZZ-1))/4;
        hy_fz2=(hy(SX:XX-1,SY:YY-1,ZZ)+hy(SX:XX-1,SY:YY-1,ZZ-1)+hy(SX:XX-1,SY+1:YY,ZZ)+hy(SX:XX-1,SY+1:YY,ZZ-1))/4;
        
        ey_fx1=squeeze(ey_fx1);ey_fx2=squeeze(ey_fx2);
        ez_fx1=squeeze(ez_fx1);ez_fx2=squeeze(ez_fx2);
        ex_fy1=squeeze(ex_fy1);ex_fy2=squeeze(ex_fy2);
        ez_fy1=squeeze(ez_fy1);ez_fy2=squeeze(ez_fy2);
        ex_fz1=squeeze(ex_fz1);ex_fz2=squeeze(ex_fz2);
        ey_fz1=squeeze(ey_fz1);ey_fz2=squeeze(ey_fz2);
        
        hy_fx1=squeeze(hy_fx1);hy_fx2=squeeze(hy_fx2);
        hz_fx1=squeeze(hz_fx1);hz_fx2=squeeze(hz_fx2);
        hx_fy1=squeeze(hx_fy1);hx_fy2=squeeze(hx_fy2);
        hz_fy1=squeeze(hz_fy1);hz_fy2=squeeze(hz_fy2);
        hx_fz1=squeeze(hx_fz1);hx_fz2=squeeze(hx_fz2);
        hy_fz1=squeeze(hy_fz1);hy_fz2=squeeze(hy_fz2);
        %*************************
        Mz_fx1=--ey_fx1;Mz_fx2=-+ey_fx2;
        My_fx1=-+ez_fx1;My_fx2=--ez_fx2;
        Mz_fy1=-+ex_fy1;Mz_fy2=--ex_fy2;
        Mx_fy1=--ez_fy1;Mx_fy2=-+ez_fy2;
        My_fz1=--ex_fz1;My_fz2=-+ex_fz2;
        Mx_fz1=-+ey_fz1;Mx_fz2=--ey_fz2;
        
        Jz_fx1=-+hy_fx1;Jz_fx2=--hy_fx2;
        Jy_fx1=--hz_fx1;Jy_fx2=-+hz_fx2;
        Jz_fy1=--hx_fy1;Jz_fy2=-+hx_fy2;
        Jx_fy1=-+hz_fy1;Jx_fy2=--hz_fy2;
        Jy_fz1=-+hx_fz1;Jy_fz2=--hx_fz2;
        Jx_fz1=--hy_fz1;Jx_fz2=-+hy_fz2;
        %*************************
        phex=exp(-i*2*pi*(n-TT*(SS-1))*f*dt);
        FMz_fx1=FMz_fx1+Mz_fx1*phex;FMz_fx2=FMz_fx2+Mz_fx2*phex;            %傅立葉變換提起幅值和相位  
        FMy_fx1=FMy_fx1+My_fx1*phex;FMy_fx2=FMy_fx2+My_fx2*phex;