FMz_fy1=FMz_fy1+Mz_fy1*phex;FMz_fy2=FMz_fy2+Mz_fy2*phex;
        FMx_fy1=FMx_fy1+Mx_fy1*phex;FMx_fy2=FMx_fy2+Mx_fy2*phex;
        FMy_fz1=FMy_fz1+My_fz1*phex;FMy_fz2=FMy_fz2+My_fz2*phex;
        FMx_fz1=FMx_fz1+Mx_fz1*phex;FMx_fz2=FMx_fz2+Mx_fz2*phex;
        
        phex=exp(-i*2*pi*(n-1/2-TT*(SS-1))*f*dt);
        FJz_fx1=FJz_fx1+Jz_fx1*phex;FJz_fx2=FJz_fx2+Jz_fx2*phex;
        FJy_fx1=FJy_fx1+Jy_fx1*phex;FJy_fx2=FJy_fx2+Jy_fx2*phex;
        FJz_fy1=FJz_fy1+Jz_fy1*phex;FJz_fy2=FJz_fy2+Jz_fy2*phex;
        FJx_fy1=FJx_fy1+Jx_fy1*phex;FJx_fy2=FJx_fy2+Jx_fy2*phex;
        FJy_fz1=FJy_fz1+Jy_fz1*phex;FJy_fz2=FJy_fz2+Jy_fz2*phex;
        FJx_fz1=FJx_fz1+Jx_fz1*phex;FJx_fz2=FJx_fz2+Jx_fz2*phex;
    end
      set(gca,'FontName','TimesNewRoman');
      set(gca,'FontSize',14);
      title(['時間步為:',num2str(n),'步']);  

    pause(.01);  
end
%**********************************************************
nx=[(SX+1/2-midx):1:(XX-1/2-midx)]*dx;     %積分變量設(shè)置
ny=[(SY+1/2-midy):1:(YY-1/2-midy)]*dx;
nz=[(SZ+1/2-midz):1:(ZZ-1/2-midz)]*dx;
%*************************
xf1_x=repmat((SX-midx)*dx,[YY-SY ZZ-SZ]);      
xf2_x=repmat((XX-midx)*dx,[YY-SY ZZ-SZ]);
yf1_y=repmat((SY-midy)*dx,[XX-SX ZZ-SZ]);
yf2_y=repmat((XX-midy)*dx,[XX-SX ZZ-SZ]);
zf1_z=repmat((SZ-midz)*dx,[XX-SX YY-SY]);
zf2_z=repmat((ZZ-midz)*dx,[XX-SX YY-SY]);

xf1_y=repmat(ny',[1 ZZ-SZ]);xf2_y=repmat(ny',[1 ZZ-SZ]);
yf1_x=repmat(nx',[1 ZZ-SZ]);yf2_x=repmat(nx',[1 ZZ-SZ]);
zf1_x=repmat(nx',[1 YY-SY]);zf2_x=repmat(nx',[1 YY-SY]);

xf1_z=repmat(nz,[YY-SY 1]);xf2_z=repmat(nz,[YY-SY 1]);
yf1_z=repmat(nz,[XX-SX 1]);yf2_z=repmat(nz,[XX-SX 1]);
zf1_y=repmat(ny,[XX-SX 1]);zf2_y=repmat(ny,[XX-SX 1]);
%************************
n=40;
phi = (0:1:n)/n*2*pi;
theta = (0:1:n)/n*pi;
for l=1:length(phi),
    for m=1:length(theta),    
        index1=i*k*sin(theta(m))*cos(phi(l));
        index2=i*k*sin(theta(m))*sin(phi(l));
        index3=i*k*cos(theta(m));
       
       xf1_fz(l,m)= sum(sum(FJz_fx1.*exp(index1*xf1_x+index2*xf1_y+index3*xf1_z)*dx*dx));        %x向第一個面
       xf1_fy(l,m)= sum(sum(FJy_fx1.*exp(index1*xf1_x+index2*xf1_y+index3*xf1_z)*dx*dx)); 
       xf1_fmz(l,m)=sum(sum(FMz_fx1.*exp(index1*xf1_x+index2*xf1_y+index3*xf1_z)*dx*dx));
       xf1_fmy(l,m)=sum(sum(FMy_fx1.*exp(index1*xf1_x+index2*xf1_y+index3*xf1_z)*dx*dx));
       
       xf2_fz(l,m)= sum(sum(FJz_fx2.*exp(index1*xf2_x+index2*xf2_y+index3*xf2_z)*dx*dx));       %x向第二個面
       xf2_fy(l,m)= sum(sum(FJy_fx2.*exp(index1*xf2_x+index2*xf2_y+index3*xf2_z)*dx*dx)); 
       xf2_fmz(l,m)=sum(sum(FMz_fx2.*exp(index1*xf2_x+index2*xf2_y+index3*xf2_z)*dx*dx));
       xf2_fmy(l,m)=sum(sum(FMy_fx2.*exp(index1*xf2_x+index2*xf2_y+index3*xf2_z)*dx*dx));
       
       yf1_fx(l,m)= sum(sum(FJx_fy1.*exp(index1*yf1_x+index2*yf1_y+index3*yf1_z)*dx*dx));
       yf1_fz(l,m)= sum(sum(FJz_fy1.*exp(index1*yf1_x+index2*yf1_y+index3*yf1_z)*dx*dx));
       yf1_fmx(l,m)=sum(sum(FMx_fy1.*exp(index1*yf1_x+index2*yf1_y+index3*yf1_z)*dx*dx));
       yf1_fmz(l,m)=sum(sum(FMz_fy1.*exp(index1*yf1_x+index2*yf1_y+index3*yf1_z)*dx*dx));
       
       yf2_fx(l,m)= sum(sum(FJx_fy2.*exp(index1*yf2_x+index2*yf2_y+index3*yf2_z)*dx*dx));
       yf2_fz(l,m)= sum(sum(FJz_fy2.*exp(index1*yf2_x+index2*yf2_y+index3*yf2_z)*dx*dx));
       yf2_fmx(l,m)=sum(sum(FMx_fy2.*exp(index1*yf2_x+index2*yf2_y+index3*yf2_z)*dx*dx));
       yf2_fmz(l,m)=sum(sum(FMz_fy2.*exp(index1*yf2_x+index2*yf2_y+index3*yf2_z)*dx*dx));
       
       zf1_fx(l,m)= sum(sum(FJx_fz1.*exp(index1*zf1_x+index2*zf1_y+index3*zf1_z)*dx*dx));
       zf1_fy(l,m)= sum(sum(FJy_fz1.*exp(index1*zf1_x+index2*zf1_y+index3*zf1_z)*dx*dx));
       zf1_fmx(l,m)=sum(sum(FMx_fz1.*exp(index1*zf1_x+index2*zf1_y+index3*zf1_z)*dx*dx));
       zf1_fmy(l,m)=sum(sum(FMy_fz1.*exp(index1*zf1_x+index2*zf1_y+index3*zf1_z)*dx*dx));
       
       zf2_fx(l,m)= sum(sum(FJx_fz2.*exp(index1*zf2_x+index2*zf2_y+index3*zf2_z)*dx*dx));
       zf2_fy(l,m)= sum(sum(FJy_fz2.*exp(index1*zf2_x+index2*zf2_y+index3*zf2_z)*dx*dx));
       zf2_fmx(l,m)=sum(sum(FMx_fz2.*exp(index1*zf2_x+index2*zf2_y+index3*zf2_z)*dx*dx));
       zf2_fmy(l,m)=sum(sum(FMy_fz2.*exp(index1*zf2_x+index2*zf2_y+index3*zf2_z)*dx*dx));
    end

end 

    fx=yf1_fx+yf2_fx+zf1_fx+zf2_fx;                                      %各遠區(qū)剖分點電場
    fmx=yf1_fmx+yf2_fmx+zf1_fmx+zf2_fmx;
    fy=xf1_fy+xf2_fy+zf2_fy+zf1_fy;
    fmy=xf1_fmy+xf2_fmy+zf2_fmy+zf1_fmy;
    fz=xf1_fz+xf2_fz+yf1_fz+yf2_fz;
    fmz=xf1_fmz+xf2_fmz+yf1_fmz+yf2_fmz;
    
  for l=1:length(phi),
        for m=1:length(theta),            
            index1=cos(theta(m))*cos(phi(l));
            index2=cos(theta(m))*sin(phi(l));
            index3=sin(theta(m));
            
            f_theta(l,m)=fx(l,m)*index1+fy(l,m)*index2-fz(l,m)*index3;
            f_m_theta(l,m)=fmx(l,m)*index1+fmy(l,m)*index2-fmz(l,m)*index3;
            
            f_phi(l,m)=-fx(l,m)*sin(phi(l))+fy(l,m)*cos(phi(l));
            f_m_phi(l,m)=-fmx(l,m)*sin(phi(l))+fmy(l,m)*cos(phi(l));
        end
    end
    
   r=abs(Z*f_theta+f_m_phi).*abs(Z*f_theta+f_m_phi) +abs(-Z*f_phi+f_m_theta).*abs(-Z*f_phi+f_m_theta);
   r=r/max(max(r));
  %***********************************************************************
  figure(1);
  r_h1=squeeze(r(1,:)); r_h2=squeeze(r(fix(length(phi)/2),:));
  polar([theta  -theta],[r_h1 r_h2],'r');
  
  figure(2);
  r_e=squeeze(r(:,fix(length(theta)/2)));
  polar(phi,r_e','r');
  
  figure(3);
for l=1:length(phi),
    for m=1:length(theta),    
        x(l,m)=sin(theta(m))*cos(phi(l))*r(l,m);
        y(l,m)=sin(theta(m))*sin(phi(l))*r(l,m);
        z(l,m)=cos(theta(m))*r(l,m);
    end
end
surf(x,y,z);
%*****************************************************方向系數(shù)計算
d_theta=1/n*pi;
d_phi=1/n*2*pi;
solidangle=0;
for l=1:length(phi),
     for m=1:length(theta),     
          solidangle=solidangle+r(l,m)*sin(theta(m))*d_theta*d_phi;
      end
  end
D=4*pi/solidangle;
toc