xiangyang27.cpp

３D-FDTD防真

/*PML完全匹配层*/
#include <math.h>
#include "shuju8.h"
#include "bianjie.h"

 main()
 {   FILE *out;
     char outfile[3];
	 int m,loop,loop1,i,j,k,start=1000; 
	 float Iz,Uz,Uzz[NN],Izz[NN],fs=0;
	 fuzhi();

   printf("%f,%f,%f,%f,%f\n",dx,dy,dz,c/f);
   printf("%f\n",1/sqrt((1/dx)*(1/dx)+(1/dy)*(1/dy)+(1/dz)*(1/dz)));
   printf("%f\n",1/(dx*dy/dz));
//   printf("%f\n",1.0/(2*Pi*f*dt));
   for(lp=0;lp<1;lp++)
   {// f=f+20000000; 
	   dt=1.0/(f*400);fuzhi();
   printf("%f\n",1/sqrt((1/dx)*(1/dx)+(1/dy)*(1/dy)+(1/dz)*(1/dz)));
 //  printf("%f\n",1/(2*Pi*f*ee*er*dx*dy/dz));
   for(loop=1;loop<=NN;loop++)   /*主时间循环*/
   {
	   m=chuanbo1(loop);
	   m=bianjie();
	   
//	   Uz=-exp(-(loop*dt-100*dt)*(loop*dt-100*dt)/(dt*dt*1000))*sin(2*Pi*f*(loop*dt-100*dt))*1000*dz;
//	   Uz=-(sin(2*Pi*(loop-0.5)*f*dt)+sin(2*Pi*(loop-0.5)*(f+20000000)*dt)+sin(2*Pi*(loop-0.5)*(f+40000000)*dt)+sin(2*Pi*(loop-0.5)*(f+60000000)*dt)+sin(2*Pi*(loop-0.5)*(f+80000000)*dt)+sin(2*Pi*(loop-0.5)*(f+100000000)*dt)+sin(2*Pi*(loop-0.5)*(f+120000000)*dt)+sin(2*Pi*(loop-0.5)*(f+140000000)*dt)+sin(2*Pi*(loop-0.5)*(f+160000000)*dt)+sin(2*Pi*(loop-0.5)*(f+180000000)*dt)+sin(2*Pi*(loop-0.5)*(f+200000000)*dt)+sin(2*Pi*(loop-0.5)*(f+220000000)*dt)+sin(2*Pi*(loop-0.5)*(f+240000000)*dt)+sin(2*Pi*(loop-0.5)*(f+260000000)*dt)+sin(2*Pi*(loop-0.5)*(f+280000000)*dt)+sin(2*Pi*(loop-0.5)*(f+300000000)*dt))*dz*1000;
	   Uz=-sin(2*Pi*f*(loop-0.5)*dt)*1000*dz;
	   Uzz[loop-1]=Uz;
	  Iz=(-dx*(Hxy[Nx/2-DD/2+kx][Ny/2-DD/2+ky-1][NNN+5]-Hxy[Nx/2-DD/2+kx][Ny/2-DD/2+ky][NNN+5])-dy*(Hyz[Nx/2-DD/2+kx][Ny/2-DD/2+ky][NNN+5]-Hyz[Nx/2-DD/2+kx-1][Ny/2-DD/2+ky][NNN+5]))*1000;
	   Izz[loop-1]=Iz;
	  
	   printf("%f,%f\n",Uz,Iz);
	 
	   m=chuanbo2(loop);
		
            
			printf("%d\n",loop);
			printf("%f\n",Ezx[Nx/2][Ny/2][Nz/2]);
			printf("%f\n",Ezy[Nx/2][Ny/2][Nz/2]);

//		    E1[loop]=sqrt(Ezx[Nx/2][Ny/2][Nz/2]*Ezx[Nx/2][Ny/2][Nz/2]+Exy[Nx/2][Ny/2][Nz/2]*Exy[Nx/2][Ny/2][Nz/2]+Eyz[Nx/2][Ny/2][Nz/2]*Eyz[Nx/2][Ny/2][Nz/2]);
//          H1[loop]=sqrt(Hzx[Nx/2][Ny/2][Nz/2]*Hzx[Nx/2][Ny/2][Nz/2]+Hxy[Nx/2][Ny/2][Nz/2]*Hxy[Nx/2][Ny/2][Nz/2]+Hyz[Nx/2][Ny/2][Nz/2]*Hyz[Nx/2][Ny/2][Nz/2]);

   
    flag3++;
/*   if(loop==100||loop==200||loop==300)
   { printf("enter the outfile name:\n");
   scanf("%s",outfile);
   if ((out=fopen(outfile,"w"))==NULL)
   {printf("cannot open outfile\n");
    goto end;
   }
   else  
	   for(j=0;j<Ny;j++)
	      {for(k=0;k<Nz;k++)
	   {fprintf(out,"%f   ",Exy[Nx/2][j][k]);
	   }
	   fprintf(out,";\n");
       }


     fclose(out);

	  printf("enter the outfile name:\n");
   scanf("%s",outfile);
   if ((out=fopen(outfile,"w"))==NULL)
   {printf("cannot open outfile\n");
    goto end;
   }
   else  
	   for(j=0;j<Nx;j++)
	      {for(k=0;k<Nz;k++)
	   {fprintf(out,"%f   ",Eyz[j][Ny/2][k]);
	   }
	   fprintf(out,";\n");
       }


     fclose(out);

   printf("enter the outfile name:\n");
   scanf("%s",outfile);
   if ((out=fopen(outfile,"w"))==NULL)
   {printf("cannot open outfile\n");
    goto end;
   }
   else  
	   for(j=0;j<Nx;j++)
	      {for(k=0;k<Ny;k++)
	   {fprintf(out,"%f   ",Ezx[j][k][Nz/2]);
	   }
	   fprintf(out,";\n");
       }


     fclose(out);

   } 
*/

//	if(loop>=NN)
//		waitui(loop);
		if(loop>=NN-400)
        distill(loop);

 }
  printf("************************************\n");
//  printf("频率为:%f时,第%d 运算",f,lp);
  printf("************************************\n");
   ffft(Uzz,Izz);
   zukang();
   zhubobi(); 
  
   
    for (i=0;i<Nx;i++)
	 for (j=0;j<Ny;j++)
		 for (k=0;k<Nz;k++)
			 Exy[i][j][k]=sqrt(Exy[i][j][k]*Exy[i][j][k]+Eyx[i][j][k]*Eyx[i][j][k]+Ezx[i][j][k]*Ezx[i][j][k]);
        
   
   
   if ((out=fopen("jiexiao1","a"))==NULL)
   {printf("cannot open outfile\n");
    goto end;
   }
   else  

	   for(j=0;j<Ny;j++)
	      {for(k=0;k<Nz;k++)
	   {fprintf(out,"%f   ",Exy[Nx/2-5][j][k]);
	   }
	   fprintf(out,";\n");
       }

    fprintf(out,"\n\n\n");
     fclose(out);


   if ((out=fopen("jiexiao2","a"))==NULL)
   {printf("cannot open outfile\n");
    goto end;
   }
   else  

	   for(j=0;j<Nx;j++)
	      {for(k=0;k<Nz;k++)
	   {fprintf(out,"%f   ",Exy[j][Ny/2-5][k]);
	   }
	   fprintf(out,";\n");
       }

   fprintf(out,"\n\n\n");
     fclose(out);

  
   if ((out=fopen("jiexiao3","a"))==NULL)
   {printf("cannot open outfile\n");
    goto end;
   }
   else  

	   for(j=0;j<Nx;j++)
	      {for(k=0;k<Ny;k++)
	   {fprintf(out,"%f   ",Exy[j][k][Nz/2]);
	   }
	   fprintf(out,";\n");
       }

    fprintf(out,"\n\n\n");
     fclose(out);
  
   waitui(loop,0);
 
   if ((out=fopen("jiexiao4","a"))==NULL)
   {printf("cannot open outfile\n");
    goto end;
   }
   else  

	 {  for(j=0;j<=180;j++)
	      fprintf(out,"%f\n",sqrt(fxhs[j].shi*fxhs[j].shi+fxhs[j].xu*fxhs[j].xu));
     }
     fprintf(out,"\n\n\n");
    fclose(out);

	waitui(loop,90);
  
   if ((out=fopen("jiexiao5","a"))==NULL)
   {printf("cannot open outfile\n");
    goto end;
   }
   else  

	 {  for(j=0;j<=180;j++)
	      fprintf(out,"%f\n",sqrt(fxhs[j].shi*fxhs[j].shi+fxhs[j].xu*fxhs[j].xu));
     }
   fprintf(out,"\n\n\n");
   fclose(out);

   if ((out=fopen("jiexiao100","a"))==NULL)
   {printf("cannot open outfile\n");
    goto end;
   }
   else  

   {     fprintf(out,"%f     %f;\n",fxhs1.sita.shi,fxhs1.sita.xu);
		 fprintf(out,"%f     %f;\n",fxhs1.fai.shi,fxhs1.fai.xu);
   }
      fclose(out);

}

end:  return(0);
 }







void fuzhi()
{int loop1,loop2,loop3;
    for(loop1=0;loop1<Nx;loop1++)
	{ for(loop2=0;loop2<Ny+1;loop2++)
	   {for(loop3=0;loop3<Nz+1;loop3++)
	{      Exy[loop1][loop2][loop3]=0;
	       Exz[loop1][loop2][loop3]=0;
	}
	}
    }   /*给x向电场赋初值*/

    for(loop1=0;loop1<Nx+1;loop1++)
	{ for(loop2=0;loop2<Ny;loop2++)
	   {for(loop3=0;loop3<Nz+1;loop3++)
	{     Eyx[loop1][loop2][loop3]=0;
	      Eyz[loop1][loop2][loop3]=0;
	}
	}
    }   /*给y向电场赋初值*/

   for(loop1=0;loop1<Nx+1;loop1++)
	{ for(loop2=0;loop2<Ny+1;loop2++)
	   {for(loop3=0;loop3<Nz;loop3++)
     
   {       Ezx[loop1][loop2][loop3]=0;
	       Ezy[loop1][loop2][loop3]=0;
   }
       
	}
    }   /*给z向电场赋初值,并加激励源*/

      for(loop1=0;loop1<Nx+1;loop1++)
	{ for(loop2=0;loop2<Ny;loop2++)
	   {for(loop3=0;loop3<Nz;loop3++)
	{Hxy[loop1][loop2][loop3]=0;
     Hxz[loop1][loop2][loop3]=0;

	}
	}
    }   /*给x向磁场赋初值*/

     for(loop1=0;loop1<Nx;loop1++)
   { for(loop2=0;loop2<Ny+1;loop2++)
	   {for(loop3=0;loop3<Nz;loop3++)
	{Hyx[loop1][loop2][loop3]=0;
     Hyz[loop1][loop2][loop3]=0;
	}
	}
    }   /*给y向磁场赋初值*/

      for(loop1=0;loop1<Nx;loop1++)
	{ for(loop2=0;loop2<Ny;loop2++)
	   {for(loop3=0;loop3<Nz+1;loop3++)
	{  Hzx[loop1][loop2][loop3]=0;
       Hzy[loop1][loop2][loop3]=0;
	}
	}
    }   /*给z向磁场赋初值*/
    

	  for(loop1=0;loop1<NN/2;loop1++)
	  {
	    pp1[loop1].shi=0;pp1[loop1].xu=0;
		 pp2[loop1].shi=0;pp2[loop1].xu=0;
		  zk[loop1].shi=0;zk[loop1].xu=0;
	  }

	   fxhs1.sita.shi=0;fxhs1.sita.xu=0;
	   fxhs1.fai.shi=0;fxhs1.fai.xu=0;

	   for(loop1=0;loop1<=180;loop1++)
	   {fxhs[loop1].shi=0,fxhs[loop1].xu=0;}

		for(loop1=0;loop1<=180;loop1++)
		{si[loop1]=sin(loop1/180.0*Pi);co[loop1]=cos(loop1/180.0*Pi);}

		 for(loop1=0;loop1<Nx-2*NNN-2*DDD;loop1++)
			for(loop2=0;loop2<Nx-2*NNN-2*DDD+1;loop2++)
			{
			   Jxu[loop1][loop2].zhengfu=0; Jxu[loop1][loop2].xiangwei=0;
			   Jxd[loop1][loop2].zhengfu=0; Jxd[loop1][loop2].xiangwei=0;
			   Jxl[loop1][loop2].zhengfu=0; Jxl[loop1][loop2].xiangwei=0;
			   Jxr[loop1][loop2].zhengfu=0; Jxr[loop1][loop2].xiangwei=0;
			   Jzl[loop1][loop2].zhengfu=0; Jzl[loop1][loop2].xiangwei=0;
			   Jzr[loop1][loop2].zhengfu=0; Jzr[loop1][loop2].xiangwei=0;
			   Jmxu[loop1][loop2].zhengfu=0; Jmxu[loop1][loop2].xiangwei=0;
			   Jmxd[loop1][loop2].zhengfu=0; Jmxd[loop1][loop2].xiangwei=0;
			   Jmxl[loop1][loop2].zhengfu=0; Jmxl[loop1][loop2].xiangwei=0;
			   Jmxr[loop1][loop2].zhengfu=0; Jmxr[loop1][loop2].xiangwei=0;
			   Jmzl[loop1][loop2].zhengfu=0; Jmzl[loop1][loop2].xiangwei=0;
			   Jmzr[loop1][loop2].zhengfu=0; Jmzr[loop1][loop2].xiangwei=0;