xiangyang24.cpp

全参数模板在FDTD数值算法下．微带天线方向系数、半空间方向函数、输入阻抗、驻波比、固定方位角的极化轴比的计算方法。

		   }
		   else*/
		   {  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=7.6347*(NNN-loop3-1)/NNN*(NNN-loop3-1)/NNN;
		    sgmmz=sgmz*142122.3034;
		 }
		  else if(loop3>=Nz-NNN&&loop3<Nz)
		 {  sgmz=7.6347*(NNN+loop3-Nz)/NNN*(NNN+loop3-Nz)/NNN;
		    sgmmz=sgmz*142122.3034;
		 }
		  else
		  {sgmz=0;sgmmz=0;}


		 }

          Ezx[loop1][loop2][loop3]=((2*e-sgmx*dt)/(2*e+sgmx*dt))*Ezx[loop1][loop2][loop3]+(2*dt)/(2*e+sgmx*dt)*((Hyz[loop1][loop2][loop3]+Hyx[loop1][loop2][loop3]-Hyz[loop1-1][loop2][loop3]-Hyx[loop1-1][loop2][loop3])/dx);
		  Ezy[loop1][loop2][loop3]=((2*e-sgmy*dt)/(2*e+sgmy*dt))*Ezy[loop1][loop2][loop3]-(2*dt)/(2*e+sgmy*dt)*((Hxy[loop1][loop2][loop3]+Hxz[loop1][loop2][loop3]-Hxy[loop1][loop2-1][loop3]-Hxz[loop1][loop2-1][loop3])/dy);

		   if(loop1>NNN-1&&loop1<Nx-NNN&&loop2>NNN-1&&loop2<Ny-NNN&&loop3>NNN-1&&loop3<Nz-NNN)
		   {Ezx[loop1][loop2][loop3]+= Ezy[loop1][loop2][loop3];
		    Ezy[loop1][loop2][loop3]=0;}
		   e=ee;
	     
	  }		 
	  }
	  }
	  return(1);
	  }


/***********************************************************************************************************************************************************************************************************************************/
   /*x向磁场*/
int chuanbo2(int M)
{int loop1,loop2,loop3;
 for(loop1=0;loop1<Nx+1;loop1++)
	{for(loop2=0;loop2<Ny;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=7.6347*(NNN-loop3-1)/NNN*(NNN-loop3-1)/NNN;
		    sgmmz=sgmz*142122.3034;
		 }
		  else if(loop3>=Nz-NNN&&loop3<Nz)
		 {  sgmz=7.6347*(NNN+loop3-Nz)/NNN*(NNN+loop3-Nz)/NNN;
		    sgmmz=sgmz*142122.3034;
		 }
		  else
		  {sgmz=0;sgmmz=0;}


		 }
/*	if(loop1==Nx/2-DD/2+kx&&loop2==Ny/2-DD/2+ky&&loop3>NNN+5&&loop3<NNN+5+5)  
	{
	      Hxy[loop1][loop2][loop3]=((2*u-sgmmy*dt)/(2*u+sgmmy*dt))*Hxy[loop1][loop2][loop3]-(2*dt)/(2*u+sgmmy*dt)*(2.0/log(dy/r0))*((Ezx[loop1][loop2+1][loop3]+Ezy[loop1][loop2+1][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); 

	}
    else if(loop1==Nx/2-DD/2+kx&&loop2==Ny/2-DD/2+ky-1&&loop3>NNN+5&&loop3<NNN+5+5)  
	{
	      Hxy[loop1][loop2][loop3]=((2*u-sgmmy*dt)/(2*u+sgmmy*dt))*Hxy[loop1][loop2][loop3]+(2*dt)/(2*u+sgmmy*dt)*(2.0/log(dy/r0))*((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); 
	}

    else*/
	{   
          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向磁场*/
	  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=7.6347*(NNN-loop3-1)/NNN*(NNN-loop3-1)/NNN;
		    sgmmz=sgmz*142122.3034;
		 }
		  else if(loop3>=Nz-NNN&&loop3<Nz)
		 {  sgmz=7.6347*(NNN+loop3-Nz)/NNN*(NNN+loop3-Nz)/NNN;
		    sgmmz=sgmz*142122.3034;
		 }
		  else
		  {sgmz=0;sgmmz=0;}


		 }
  /*   if(loop1==Nx/2-DD/2+kx&&loop2==Ny/2-DD/2+ky&&loop3>NNN+5&&loop3<NNN+5+5) 
	 {     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(dx/r0))*((Ezx[loop1+1][loop2][loop3]+Ezy[loop1+1][loop2][loop3])/dx);
	 }

	 else if(loop1==Nx/2-DD/2+kx-1&&loop2==Ny/2-DD/2+ky&&loop3>NNN+5&&loop3<NNN+5+5) 
	 {       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(dx/r0))*((Ezx[loop1][loop2][loop3]+Ezy[loop1][loop2][loop3])/dx);
	 }

     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向磁场*/
	  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=7.6347*(NNN-loop3-1)/NNN*(NNN-loop3-1)/NNN;
		    sgmmz=sgmz*142122.3034;
		 }
		  else if(loop3>=Nz-NNN&&loop3<Nz)
		 {  sgmz=7.6347*(NNN+loop3-Nz)/NNN*(NNN+loop3-Nz)/NNN;
		    sgmmz=sgmz*142122.3034;
		 }
		  else
		  {sgmz=0;sgmmz=0;}


		 }
/*		   if(loop1==Nx/2-DD/2+kx&&loop2==Ny/2-DD/2+ky&&loop3>=NNN+6&&loop3<NNN+5+5)
		   {
	               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-DD/2+kx&&loop2==Ny/2-DD/2+ky-1&&loop3>=NNN+6&&loop3<NNN+5+5)
		   {
	                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-DD/2+kx-1&&loop2==Ny/2-DD/2+ky&&loop3>=NNN+6&&loop3<NNN+5+5)
		   {
	                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-DD/2+kx-1&&loop2==Ny/2-DD/2+ky-1&&loop3>=NNN+6&&loop3<NNN+5+5)
		   {
	                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;}

	   }  
	  }
	  }          /*设置波动源，向前波动1个时间周期*/


       return(1);    /*此步仅为便于调试*/
       
}

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;
	 loop1=90;loop2=M;
	  for(loop1=0;loop1<180;loop1++)
	  {
	  for(i=NNN+DDD;i<Nx-NNN-DDD;i++)                    //上面y向电磁场
	     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向电磁场
	     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向电磁场