Hxy[loop1][loop2][loop3]=((2*u-sgmmy*dt)/(2*u+sgmmy*dt))*Hxy[loop1][loop2][loop3]-(2*dt)/(2*u+sgmmy*dt)*((Ezx[loop1][loop2+1][loop3]+Ezy[loop1][loop2+1][loop3]-Ezx[loop1][loop2][loop3]-Ezy[loop1][loop2][loop3])/dy);
          Hxz[loop1][loop2][loop3]=((2*u-sgmmz*dt)/(2*u+sgmmz*dt))*Hxz[loop1][loop2][loop3]+(2*dt)/(2*u+sgmmz*dt)*((Eyz[loop1][loop2][loop3+1]+Eyx[loop1][loop2][loop3+1]-Eyz[loop1][loop2][loop3]-Eyx[loop1][loop2][loop3])/dz);
		}
		   if(loop1>NNN-1&&loop1<Nx-NNN&&loop2>NNN-1&&loop2<Ny-NNN&&loop3>NNN-1&&loop3<Nz-NNN)
		   { Hxy[loop1][loop2][loop3]+= Hxz[loop1][loop2][loop3];
		     Hxz[loop1][loop2][loop3]=0;}
	   } 
	  }
	  }
 
 /**************************************************************************************************************************************************************************************************************************************/
      /*y向磁場(chǎng)*/
	  for(loop1=0;loop1<Nx;loop1++)
	  {for(loop2=0;loop2<Ny+1;loop2++)
	  {for(loop3=0;loop3<Nz;loop3++)
	  {   if(loop1>NNN-1&&loop1<Nx-NNN&&loop2>NNN-1&&loop2<Ny-NNN&&loop3>NNN-1&&loop3<Nz-NNN)
		{sgmx=0;sgmy=0;sgmz=0;
		 sgmmx=0;sgmmy=0;sgmmz=0;
		}
	     else
		 { 
		  if(loop1>=0&&loop1<=NNN-1)
		 {  sgmx=sgmmax*(NNN-loop1-1)/NNN*(NNN-loop1-1)/NNN;
		    sgmmx=sgmx*142122.3034;
		 }
		  else if(loop1>=Nx-NNN&&loop1<Nx)
		 {  sgmx=sgmmax*(NNN+loop1-Nx)/NNN*(NNN+loop1-Nx)/NNN;
		    sgmmx=sgmx*142122.3034;
		 }
		  else
		  {sgmx=0;sgmmx=0;}


           if(loop2>=0&&loop2<=NNN-1)
		 {  sgmy=sgmmax*(NNN-loop2-1)/NNN*(NNN-loop2-1)/NNN;
		    sgmmy=sgmy*142122.3034;
		 }
		  else if(loop2>=Ny-NNN&&loop2<Ny)
		 {  sgmy=sgmmax*(NNN+loop2-Ny)/NNN*(NNN+loop2-Ny)/NNN;
		    sgmmy=sgmy*142122.3034;
		 }
		  else
		  {sgmy=0;sgmmy=0;}


		   if(loop3>=0&&loop3<=NNN-1)
		 {  sgmz=sgmmax*(NNN-loop3-1)/NNN*(NNN-loop3-1)/NNN;
		    sgmmz=sgmz*142122.3034;
		 }
		  else if(loop3>=Nz-NNN&&loop3<Nz)
		 {  sgmz=sgmmax*(NNN+loop3-Nz)/NNN*(NNN+loop3-Nz)/NNN;
		    sgmmz=sgmz*142122.3034;
		 }
		  else
		  {sgmz=0;sgmmz=0;}


		 }
     if(loop1==Nx/2&&loop2==Ny/2&&((loop3>=Nz/2-DD1/2-1&&loop3<Nz/2-1)||(loop3>Nz/2&&loop3<=Nz/2+DD1/2)))
	 {       Hyz[loop1][loop2][loop3]=((2*u-sgmmz*dt)/(2*u+sgmmz*dt))*Hyz[loop1][loop2][loop3]-(2*dt)/(2*u+sgmmz*dt)*((Exy[loop1][loop2][loop3+1]+Exz[loop1][loop2][loop3+1]-Exy[loop1][loop2][loop3]-Exz[loop1][loop2][loop3])/dz);
             Hyx[loop1][loop2][loop3]=((2*u-sgmmx*dt)/(2*u+sgmmx*dt))*Hyx[loop1][loop2][loop3]+(2*dt)/(2*u+sgmmx*dt)*(2.0/log(dz/r0))*((Ezx[loop1+1][loop2][loop3]+Ezy[loop1+1][loop2][loop3])/dx);
	 }
	 else if(loop1==Nx/2-1&&loop2==Ny/2&&((loop3>=Nz/2-DD1/2-1&&loop3<Nz/2-1)||(loop3>Nz/2&&loop3<=Nz/2+DD1/2)))
	 {       Hyz[loop1][loop2][loop3]=((2*u-sgmmz*dt)/(2*u+sgmmz*dt))*Hyz[loop1][loop2][loop3]-(2*dt)/(2*u+sgmmz*dt)*((Exy[loop1][loop2][loop3+1]+Exz[loop1][loop2][loop3+1]-Exy[loop1][loop2][loop3]-Exz[loop1][loop2][loop3])/dz);
             Hyx[loop1][loop2][loop3]=((2*u-sgmmx*dt)/(2*u+sgmmx*dt))*Hyx[loop1][loop2][loop3]-(2*dt)/(2*u+sgmmx*dt)*(2.0/log(dz/r0))*((Ezx[loop1][loop2][loop3]+Ezy[loop1][loop2][loop3])/dx);
	 }
	 else if(((loop1>=Nx/2-DD2/2-1&&loop1<Nx/2-1)||(loop1>Nx/2&&loop1<=Nx/2+DD2/2))&&loop2==Ny/2+movedis&&loop3==Nz/2)
	 {
	         Hyz[loop1][loop2][loop3]=((2*u-sgmmz*dt)/(2*u+sgmmz*dt))*Hyz[loop1][loop2][loop3]+(2*dt)/(2*u+sgmmz*dt)*((Ezx[loop1+1][loop2][loop3]+Ezy[loop1+1][loop2][loop3]-Ezx[loop1][loop2][loop3]-Ezy[loop1][loop2][loop3])/dx);
             Hyx[loop1][loop2][loop3]=((2*u-sgmmx*dt)/(2*u+sgmmx*dt))*Hyx[loop1][loop2][loop3]-(2*dt)/(2*u+sgmmx*dt)*(2.0/log(dz/r1))*((Exy[loop1][loop2][loop3+1]+Exz[loop1][loop2][loop3+1])/dz);
	 }
	 else if(((loop1>=Nx/2-DD2/2-1&&loop1<Nx/2-1)||(loop1>Nx/2&&loop1<=Nx/2+DD2/2))&&loop2==Ny/2+movedis&&loop3==Nz/2-1)
	 {
	         Hyz[loop1][loop2][loop3]=((2*u-sgmmz*dt)/(2*u+sgmmz*dt))*Hyz[loop1][loop2][loop3]+(2*dt)/(2*u+sgmmz*dt)*((Ezx[loop1+1][loop2][loop3]+Ezy[loop1+1][loop2][loop3]-Ezx[loop1][loop2][loop3]-Ezy[loop1][loop2][loop3])/dx);
             Hyx[loop1][loop2][loop3]=((2*u-sgmmx*dt)/(2*u+sgmmx*dt))*Hyx[loop1][loop2][loop3]+(2*dt)/(2*u+sgmmx*dt)*(2.0/log(dz/r1))*((Exy[loop1][loop2][loop3]+Exz[loop1][loop2][loop3])/dz);
	 }
     else
	 {   
          Hyz[loop1][loop2][loop3]=((2*u-sgmmz*dt)/(2*u+sgmmz*dt))*Hyz[loop1][loop2][loop3]-(2*dt)/(2*u+sgmmz*dt)*((Exy[loop1][loop2][loop3+1]+Exz[loop1][loop2][loop3+1]-Exy[loop1][loop2][loop3]-Exz[loop1][loop2][loop3])/dz);
          Hyx[loop1][loop2][loop3]=((2*u-sgmmx*dt)/(2*u+sgmmx*dt))*Hyx[loop1][loop2][loop3]+(2*dt)/(2*u+sgmmx*dt)*((Ezx[loop1+1][loop2][loop3]+Ezy[loop1+1][loop2][loop3]-Ezx[loop1][loop2][loop3]-Ezy[loop1][loop2][loop3])/dx);
	 }
		  if(loop1>NNN-1&&loop1<Nx-NNN&&loop2>NNN-1&&loop2<Ny-NNN&&loop3>NNN-1&&loop3<Nz-NNN)
		  {Hyz[loop1][loop2][loop3]+=Hyx[loop1][loop2][loop3];
		   Hyx[loop1][loop2][loop3]=0;}

	   } 
	  }
	  }
/****************************************************************************************************************************************************************************************************************************************/       
   /*z向磁場(chǎng)*/
	  for(loop1=0;loop1<Nx;loop1++)
	  {for(loop2=0;loop2<Ny;loop2++)
	  {for(loop3=0;loop3<Nz+1;loop3++)
	 {   if(loop1>NNN-1&&loop1<Nx-NNN&&loop2>NNN-1&&loop2<Ny-NNN&&loop3>NNN-1&&loop3<Nz-NNN)
		{sgmx=0;sgmy=0;sgmz=0;
		 sgmmx=0;sgmmy=0;sgmmz=0;
		}
	     else
		 { 
		  if(loop1>=0&&loop1<=NNN-1)
		 {  sgmx=sgmmax*(NNN-loop1-1)/NNN*(NNN-loop1-1)/NNN;
		    sgmmx=sgmx*142122.3034;
		 }
		  else if(loop1>=Nx-NNN&&loop1<Nx)
		 {  sgmx=sgmmax*(NNN+loop1-Nx)/NNN*(NNN+loop1-Nx)/NNN;
		    sgmmx=sgmx*142122.3034;
		 }
		  else
		  {sgmx=0;sgmmx=0;}


           if(loop2>=0&&loop2<=NNN-1)
		 {  sgmy=sgmmax*(NNN-loop2-1)/NNN*(NNN-loop2-1)/NNN;
		    sgmmy=sgmy*142122.3034;
		 }
		  else if(loop2>=Ny-NNN&&loop2<Ny)
		 {  sgmy=sgmmax*(NNN+loop2-Ny)/NNN*(NNN+loop2-Ny)/NNN;
		    sgmmy=sgmy*142122.3034;
		 }
		  else
		  {sgmy=0;sgmmy=0;}


		   if(loop3>=0&&loop3<=NNN-1)
		 {  sgmz=sgmmax*(NNN-loop3-1)/NNN*(NNN-loop3-1)/NNN;
		    sgmmz=sgmz*142122.3034;
		 }
		  else if(loop3>=Nz-NNN&&loop3<Nz)
		 {  sgmz=sgmmax*(NNN+loop3-Nz)/NNN*(NNN+loop3-Nz)/NNN;
		    sgmmz=sgmz*142122.3034;
		 }
		  else
		  {sgmz=0;sgmmz=0;}


		 }
           if(((loop1>=Nx/2-DD2/2-1&&loop1<Nx/2-1)||(loop1>Nx/2&&loop1<=Nx/2+DD2/2))&&loop2==Ny/2+movedis&&loop3==Nz/2)
			{
	               Hzx[loop1][loop2][loop3]=((2*u-sgmmx*dt)/(2*u+sgmmx*dt))*Hzx[loop1][loop2][loop3]-(2*dt)/(2*u+sgmmx*dt)*((Eyz[loop1+1][loop2][loop3]+Eyx[loop1+1][loop2][loop3]-Eyz[loop1][loop2][loop3]-Eyx[loop1][loop2][loop3])/dx);
                   Hzy[loop1][loop2][loop3]=((2*u-sgmmy*dt)/(2*u+sgmmy*dt))*Hzy[loop1][loop2][loop3]+(2*dt)/(2*u+sgmmy*dt)*(2.0/log(dy/r1))*((Exy[loop1][loop2+1][loop3]+Exz[loop1][loop2+1][loop3])/dy);
			}
	        else if(((loop1>=Nx/2-DD2/2-1&&loop1<Nx/2-1)||(loop1>Nx/2&&loop1<=Nx/2+DD2/2))&&loop2==Ny/2-1+movedis&&loop3==Nz/2)
			{
	               Hzx[loop1][loop2][loop3]=((2*u-sgmmx*dt)/(2*u+sgmmx*dt))*Hzx[loop1][loop2][loop3]-(2*dt)/(2*u+sgmmx*dt)*((Eyz[loop1+1][loop2][loop3]+Eyx[loop1+1][loop2][loop3]-Eyz[loop1][loop2][loop3]-Eyx[loop1][loop2][loop3])/dx);
                   Hzy[loop1][loop2][loop3]=((2*u-sgmmy*dt)/(2*u+sgmmy*dt))*Hzy[loop1][loop2][loop3]-(2*dt)/(2*u+sgmmy*dt)*(2.0/log(dy/r1))*((Exy[loop1][loop2][loop3]+Exz[loop1][loop2][loop3])/dy);
			}
            else  if(loop1==Nx/2&&loop2==Ny/2&&((loop3>=Nz/2-DD1/2-1&&loop3<Nz/2)||(loop3>Nz/2&&loop3<=Nz/2+DD1/2+1)))
			{
	               Hzx[loop1][loop2][loop3]=((2*u-sgmmx*dt)/(2*u+sgmmx*dt))*Hzx[loop1][loop2][loop3]-(2*dt)/(2*u+sgmmx*dt)*1.0/(dx*dy-Pi*r0*r0/4)*(Eyz[loop1+1][loop2][loop3]*dy+Eyx[loop1+1][loop2][loop3]*dy-Eyz[loop1][loop2][loop3]*dy/2*log(dy/r0)-Eyx[loop1][loop2][loop3]*dy/2*log(dy/r0));
                   Hzy[loop1][loop2][loop3]=((2*u-sgmmy*dt)/(2*u+sgmmy*dt))*Hzy[loop1][loop2][loop3]+(2*dt)/(2*u+sgmmy*dt)*1.0/(dx*dy-Pi*r0*r0/4)*(Exy[loop1][loop2+1][loop3]*dx+Exz[loop1][loop2+1][loop3]*dx-Exy[loop1][loop2][loop3]*dx/2*log(dx/r0)-Exz[loop1][loop2][loop3]*dx/2*log(dx/r0));	       
			}
		    else  if(loop1==Nx/2&&loop2==Ny/2-1&&((loop3>=Nz/2-DD1/2-1&&loop3<Nz/2)||(loop3>Nz/2&&loop3<=Nz/2+DD1/2+1)))
			{
	                Hzx[loop1][loop2][loop3]=((2*u-sgmmx*dt)/(2*u+sgmmx*dt))*Hzx[loop1][loop2][loop3]-(2*dt)/(2*u+sgmmx*dt)*1.0/(dx*dy-Pi*r0*r0/4)*(Eyz[loop1+1][loop2][loop3]*dy+Eyx[loop1+1][loop2][loop3]*dy-Eyz[loop1][loop2][loop3]*dy/2*log(dy/r0)-Eyx[loop1][loop2][loop3]*dy/2*log(dy/r0));
                    Hzy[loop1][loop2][loop3]=((2*u-sgmmy*dt)/(2*u+sgmmy*dt))*Hzy[loop1][loop2][loop3]+(2*dt)/(2*u+sgmmy*dt)*1.0/(dx*dy-Pi*r0*r0/4)*(Exy[loop1][loop2+1][loop3]*dx/2*log(dx/r0)+Exz[loop1][loop2+1][loop3]*dx/2*log(dx/r0)-Exy[loop1][loop2][loop3]*dx-Exz[loop1][loop2][loop3]*dx);
			}
            else  if(loop1==Nx/2-1&&loop2==Ny/2&&((loop3>=Nz/2-DD1/2-1&&loop3<Nz/2)||(loop3>Nz/2&&loop3<=Nz/2+DD1/2+1)))
		   {
	                Hzx[loop1][loop2][loop3]=((2*u-sgmmx*dt)/(2*u+sgmmx*dt))*Hzx[loop1][loop2][loop3]-(2*dt)/(2*u+sgmmx*dt)*1.0/(dx*dy-Pi*r0*r0/4)*(Eyz[loop1+1][loop2][loop3]*dy/2*log(dy/r0)+Eyx[loop1+1][loop2][loop3]*dy/2*log(dy/r0)-Eyz[loop1][loop2][loop3]*dy-Eyx[loop1][loop2][loop3]*dy);
                    Hzy[loop1][loop2][loop3]=((2*u-sgmmy*dt)/(2*u+sgmmy*dt))*Hzy[loop1][loop2][loop3]+(2*dt)/(2*u+sgmmy*dt)*1.0/(dx*dy-Pi*r0*r0/4)*(Exy[loop1][loop2+1][loop3]*dx+Exz[loop1][loop2+1][loop3]*dx-Exy[loop1][loop2][loop3]*dx/2*log(dx/r0)-Exz[loop1][loop2][loop3]*dx/2*log(dx/r0));
    	   }
            else  if(loop1==Nx/2-1&&loop2==Ny/2-1&&((loop3>=Nz/2-DD1/2-1&&loop3<Nz/2)||(loop3>Nz/2&&loop3<=Nz/2+DD1/2+1)))
		   {
	                Hzx[loop1][loop2][loop3]=((2*u-sgmmx*dt)/(2*u+sgmmx*dt))*Hzx[loop1][loop2][loop3]-(2*dt)/(2*u+sgmmx*dt)*1.0/(dx*dy-Pi*r0*r0/4)*(Eyz[loop1+1][loop2][loop3]*dy/2*log(dy/r0)+Eyx[loop1+1][loop2][loop3]*dy/2*log(dy/r0)-Eyz[loop1][loop2][loop3]*dy-Eyx[loop1][loop2][loop3]*dy);
                    Hzy[loop1][loop2][loop3]=((2*u-sgmmy*dt)/(2*u+sgmmy*dt))*Hzy[loop1][loop2][loop3]+(2*dt)/(2*u+sgmmy*dt)*1.0/(dx*dy-Pi*r0*r0/4)*(Exy[loop1][loop2+1][loop3]*dx/2*log(dx/r0)+Exz[loop1][loop2+1][loop3]*dx/2*log(dx/r0)-Exy[loop1][loop2][loop3]*dx-Exz[loop1][loop2][loop3]*dx);
    	   }

 
           else
		   { 
          Hzx[loop1][loop2][loop3]=((2*u-sgmmx*dt)/(2*u+sgmmx*dt))*Hzx[loop1][loop2][loop3]-(2*dt)/(2*u+sgmmx*dt)*((Eyz[loop1+1][loop2][loop3]+Eyx[loop1+1][loop2][loop3]-Eyz[loop1][loop2][loop3]-Eyx[loop1][loop2][loop3])/dx);
          Hzy[loop1][loop2][loop3]=((2*u-sgmmy*dt)/(2*u+sgmmy*dt))*Hzy[loop1][loop2][loop3]+(2*dt)/(2*u+sgmmy*dt)*((Exy[loop1][loop2+1][loop3]+Exz[loop1][loop2+1][loop3]-Exy[loop1][loop2][loop3]-Exz[loop1][loop2][loop3])/dy);
		   }
		  if(loop1>NNN-1&&loop1<Nx-NNN&&loop2>NNN-1&&loop2<Ny-NNN&&loop3>NNN-1&&loop3<Nz-NNN)
		  {Hzx[loop1][loop2][loop3]+=Hzy[loop1][loop2][loop3];
		   Hzy[loop1][loop2][loop3]=0;}

	   }  
	  }
	  }          /*設(shè)置波動(dòng)源，向前波動(dòng)1個(gè)時(shí)間周期*/


       return(1);    /*此步僅為便于調(diào)試*/
       
}


void waitui(int M)     //huigensi原理
{
	int i,j,k,loop1,loop2;
	float tt1,tt2,Esita,Efai,cansh1,cansh2,cansh3;
	struct ff
	{float shi;float xu;};
   struct ff fx,fy,fz,fmx,fmy,fmz;
     fx.shi=0;fx.xu=0;
	 fy.shi=0;fy.xu=0;
	 fz.shi=0;fz.xu=0;
	 fmx.shi=0;fmx.xu=0;
	 fmy.shi=0;fmy.xu=0;
	 fmz.shi=0;fmz.xu=0;
	 cansh1=2*Pi*f/c*dx;cansh2=2*Pi*f/c*dy;cansh3=2*Pi*f/c*dz;
	 loop2=90;loop1=90;

	  for(i=NNN+DDD;i<Nx-NNN-DDD;i++)                    //上面y向電磁場(chǎng)
	     for(j=NNN+DDD;j<=Ny-NNN-DDD;j++)
		{
		  fx.shi-=Jxu[i-NNN-DDD][j-NNN-DDD].zhengfu*cos(-Jxu[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i+0.5-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		  fx.xu-=Jxu[i-NNN-DDD][j-NNN-DDD].zhengfu*sin(-Jxu[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i+0.5-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		  fmx.shi+=Jmxu[i-NNN-DDD][j-NNN-DDD].zhengfu*cos(-Jmxu[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j+0.5-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		  fmx.xu+=Jmxu[i-NNN-DDD][j-NNN-DDD].zhengfu*sin(-Jmxu[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j+0.5-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		}

	   for(i=NNN+DDD;i<Nx-NNN-DDD;i++)                    //下面y向電磁場(chǎng)
	     for(j=NNN+DDD;j<=Ny-NNN-DDD;j++)
		{
		  fx.shi+=Jxd[i-NNN-DDD][j-NNN-DDD].zhengfu*cos(-Jxd[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i+0.5-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j-Ny/2.0)*si[loop1]*si[loop2]-cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		  fx.xu+=Jxd[i-NNN-DDD][j-NNN-DDD].zhengfu*sin(-Jxd[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i+0.5-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j-Ny/2.0)*si[loop1]*si[loop2]-cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		  fmx.shi-=Jmxd[i-NNN-DDD][j-NNN-DDD].zhengfu*cos(-Jmxd[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j+0.5-Ny/2.0)*si[loop1]*si[loop2]-cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		  fmx.xu-=Jmxd[i-NNN-DDD][j-NNN-DDD].zhengfu*sin(-Jmxd[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j+0.5-Ny/2.0)*si[loop1]*si[loop2]-cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		}

	   for(i=NNN+DDD;i<Nx-NNN-DDD;i++)                    //右面z向電磁場(chǎng)
	     for(k=NNN+DDD;k<=Nz-NNN-DDD;k++)
		{
		  fx.shi+=Jxr[i-NNN-DDD][k-NNN-DDD].zhengfu*cos(-Jxr[i-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(i+0.5-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(Ny/2-NNN-DDD)*si[loop1]*si[loop2]+cansh3*(k-Nz/2.0)*co[loop1]);
		  fx.xu+=Jxr[i-NNN-DDD][k-NNN-DDD].zhengfu*sin(-Jxr[i-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(i+0.5-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(Ny/2-NNN-DDD)*si[loop1]*si[loop2]+cansh3*(k-Nz/2.0)*co[loop1]);
		  fmx.shi-=Jmxr[i-NNN-DDD][k-NNN-DDD].zhengfu*cos(-Jmxr[i-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(i-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(Ny/2-NNN-DDD)*si[loop1]*si[loop2]+cansh3*(k+0.5-Nz/2.0)*co[loop1]);
		  fmx.xu-=Jmxr[i-NNN-DDD][k-NNN-DDD].zhengfu*sin(-Jmxr[i-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(i-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(Ny/2-NNN-DDD)*si[loop1]*si[loop2]+cansh3*(k+0.5-Nz/2.0)*co[loop1]);
		}

	   for(i=NNN+DDD;i<Nx-NNN-DDD;i++)                    //左面y向電磁場(chǎng)
	     for(k=NNN+DDD;k<=Nz-NNN-DDD;k++)
		{
		  fx.shi-=Jxl[i-NNN-DDD][k-NNN-DDD].zhengfu*cos(-Jxl[i-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(i+0.5-Nx/2.0)*si[loop1]*co[loop2]-cansh2*(Ny/2-NNN-DDD)*si[loop1]*si[loop2]+cansh3*(k-Nz/2.0)*co[loop1]);
		  fx.xu-=Jxl[i-NNN-DDD][k-NNN-DDD].zhengfu*sin(-Jxl[i-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(i+0.5-Nx/2.0)*si[loop1]*co[loop2]-cansh2*(Ny/2-NNN-DDD)*si[loop1]*si[loop2]+cansh3*(k-Nz/2.0)*co[loop1]);
		  fmx.shi+=Jmxl[i-NNN-DDD][k-NNN-DDD].zhengfu*cos(-Jmxl[i-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(i-Nx/2.0)*si[loop1]*co[loop2]-cansh2*(Ny/2-NNN-DDD)*si[loop1]*si[loop2]+cansh3*(k+0.5-Nz/2.0)*co[loop1]);
		  fmx.xu+=Jmxl[i-NNN-DDD][k-NNN-DDD].zhengfu*sin(-Jmxl[i-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(i-Nx/2.0)*si[loop1]*co[loop2]-cansh2*(Ny/2-NNN-DDD)*si[loop1]*si[loop2]+cansh3*(k+0.5-Nz/2.0)*co[loop1]);
		}

	   for(i=NNN+DDD;i<=Nx-NNN-DDD;i++)                    //上面x向電磁場(chǎng)
	     for(j=NNN+DDD;j<Ny-NNN-DDD;j++)
		{
		  fy.shi+=Jyu[i-NNN-DDD][j-NNN-DDD].zhengfu*cos(-Jyu[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j+0.5-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		  fy.xu+=Jyu[i-NNN-DDD][j-NNN-DDD].zhengfu*sin(-Jyu[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j+0.5-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		  fmy.shi-=Jmyu[i-NNN-DDD][j-NNN-DDD].zhengfu*cos(-Jmyu[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i+0.5-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		  fmy.xu-=Jmyu[i-NNN-DDD][j-NNN-DDD].zhengfu*sin(-Jmyu[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i+0.5-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		}

	   for(i=NNN+DDD;i<=Nx-NNN-DDD;i++)                    //下面x向電磁場(chǎng)
	     for(j=NNN+DDD;j<Ny-NNN-DDD;j++)
		{
		  fy.shi-=Jyd[i-NNN-DDD][j-NNN-DDD].zhengfu*cos(-Jyd[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j+0.5-Ny/2.0)*si[loop1]*si[loop2]-cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		  fy.xu-=Jyd[i-NNN-DDD][j-NNN-DDD].zhengfu*sin(-Jyd[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j+0.5-Ny/2.0)*si[loop1]*si[loop2]-cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		  fmy.shi+=Jmyd[i-NNN-DDD][j-NNN-DDD].zhengfu*cos(-Jmyd[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i+0.5-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j-Ny/2.0)*si[loop1]*si[loop2]-cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		  fmy.xu+=Jmyd[i-NNN-DDD][j-NNN-DDD].zhengfu*sin(-Jmyd[i-NNN-DDD][j-NNN-DDD].xiangwei+cansh1*(i+0.5-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(j-Ny/2.0)*si[loop1]*si[loop2]-cansh3*(Nz/2-NNN-DDD)*co[loop1]);
		}

	   for(j=NNN+DDD;j<=Ny-NNN-DDD;j++)                    //后面z向電磁場(chǎng)
	     for(k=NNN+DDD;k<Nz-NNN-DDD;k++)
		{
		  fy.shi-=Jyf[j-NNN-DDD][k-NNN-DDD].zhengfu*cos(-Jyb[j-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(Nx/2-NNN-DDD)*si[loop1]*co[loop2]+cansh2*(j+0.5-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(k-Nz/2.0)*co[loop1]);
		  fy.xu-=Jyf[j-NNN-DDD][k-NNN-DDD].zhengfu*sin(-Jyb[j-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(Nx/2-NNN-DDD)*si[loop1]*co[loop2]+cansh2*(j+0.5-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(k-Nz/2.0)*co[loop1]);
		  fmy.shi+=Jmyf[j-NNN-DDD][k-NNN-DDD].zhengfu*cos(-Jmyb[j-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(Nx/2-NNN-DDD)*si[loop1]*co[loop2]+cansh2*(j-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(k+0.5-Nz/2.0)*co[loop1]);
		  fmy.xu+=Jmyf[j-NNN-DDD][k-NNN-DDD].zhengfu*sin(-Jmyb[j-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(Nx/2-NNN-DDD)*si[loop1]*co[loop2]+cansh2*(j-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(k+0.5-Nz/2.0)*co[loop1]);
		}

	   for(j=NNN+DDD;j<=Ny-NNN-DDD;j++)                    //前面z向電磁場(chǎng)
	     for(k=NNN+DDD;k<Nz-NNN-DDD;k++)
		{
		  fy.shi+=Jyb[j-NNN-DDD][k-NNN-DDD].zhengfu*cos(-Jyf[j-NNN-DDD][k-NNN-DDD].xiangwei-cansh1*(Nx/2-NNN-DDD)*si[loop1]*co[loop2]+cansh2*(j+0.5-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(k-Nz/2.0)*co[loop1]);
		  fy.xu+=Jyb[j-NNN-DDD][k-NNN-DDD].zhengfu*sin(-Jyf[j-NNN-DDD][k-NNN-DDD].xiangwei-cansh1*(Nx/2-NNN-DDD)*si[loop1]*co[loop2]+cansh2*(j+0.5-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(k-Nz/2.0)*co[loop1]);
		  fmy.shi-=Jmyb[j-NNN-DDD][k-NNN-DDD].zhengfu*cos(-Jmyf[j-NNN-DDD][k-NNN-DDD].xiangwei-cansh1*(Nx/2-NNN-DDD)*si[loop1]*co[loop2]+cansh2*(j-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(k+0.5-Nz/2.0)*co[loop1]);
		  fmy.xu-=Jmyb[j-NNN-DDD][k-NNN-DDD].zhengfu*sin(-Jmyf[j-NNN-DDD][k-NNN-DDD].xiangwei-cansh1*(Nx/2-NNN-DDD)*si[loop1]*co[loop2]+cansh2*(j-Ny/2.0)*si[loop1]*si[loop2]+cansh3*(k+0.5-Nz/2.0)*co[loop1]);
		}



	   for(i=NNN+DDD;i<Nx-NNN-DDD;i++)                    //右面x向電磁場(chǎng)
	     for(k=NNN+DDD;k<=Nz-NNN-DDD;k++)
		{
		  fz.shi-=Jzr[i-NNN-DDD][k-NNN-DDD].zhengfu*cos(-Jzr[i-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(i-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(Ny/2-NNN-DDD)*si[loop1]*si[loop2]+cansh3*(k+0.5-Nz/2.0)*co[loop1]);
		  fz.xu-=Jzr[i-NNN-DDD][k-NNN-DDD].zhengfu*sin(-Jzr[i-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(i-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(Ny/2-NNN-DDD)*si[loop1]*si[loop2]+cansh3*(k+0.5-Nz/2.0)*co[loop1]);
		  fmz.shi+=Jmzr[i-NNN-DDD][k-NNN-DDD].zhengfu*cos(-Jmzr[i-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(i+0.5-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(Ny/2-NNN-DDD)*si[loop1]*si[loop2]+cansh3*(k-Nz/2.0)*co[loop1]);
		  fmz.xu+=Jmzr[i-NNN-DDD][k-NNN-DDD].zhengfu*sin(-Jmzr[i-NNN-DDD][k-NNN-DDD].xiangwei+cansh1*(i+0.5-Nx/2.0)*si[loop1]*co[loop2]+cansh2*(Ny/2-NNN-DDD)*si[loop1]*si[loop2]+cansh3*(k-Nz/2.0)*co[loop1]);