xiangyang27.cpp

3D-FDTD垂直双阵子程序

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

 main()
 {   FILE *out;
     char outfile[3];
	 int m,lp,loop,loop1,i,j,k,start=1000; 
	 float Iz1,Uz1,Uzz1[NN],Izz1[NN],Iz2,Uz2,Uzz2[NN],Izz2[NN],fs=0;
	 fuzhi();

   printf("%f,%f,%f,%f,%f,%f\n",dx,dy,dz,r0,r1,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;
   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;
	   Uz1=-sin(2*Pi*f*(loop-0.5)*dt)*1000*dz;
	   Uzz1[loop-1]=Uz1;
	   Iz1=(dx*(Hxy[Nx/2][Ny/2-1][Nz/2+D/2+DD/4]-Hxy[Nx/2][Ny/2][Nz/2+D/2+DD/4])+dy*(Hyz[Nx/2][Ny/2][Nz/2+D/2+DD/4]-Hyz[Nx/2-1][Ny/2][Nz/2+D/2+DD/4]))*1000;
	   Izz1[loop-1]=Iz1;
	   Uz2=-sin(2*Pi*f*(loop-0.5)*dt-1.5*Pi)*1000*dz;
	   //Uz2=-sin(2*Pi*f*(loop-0.5)*dt)*1000*dz;
	   Uzz2[loop-1]=Uz2;
	   Iz2=(dx*(Hxy[Nx/2][Ny/2-1][Nz/2-D/2-DD/4]-Hxy[Nx/2][Ny/2][Nz/2-D/2-DD/4])+dy*(Hyz[Nx/2][Ny/2][Nz/2-D/2-DD/4]-Hyz[Nx/2-1][Ny/2][Nz/2-D/2-DD/4]))*1000;
	   Izz2[loop-1]=Iz2;
	   printf("%f,%f\n",Uz1,Iz1);
	   printf("%f,%f\n",Uz2,Iz2);
	   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-100)
	{
    for(loop1=0;loop1<=DD/2;loop1++)
	 {
     In1[loop1]=(dx*(Hxy[Nx/2][Ny/2-1][Nz/2+D/2+loop1]-Hxy[Nx/2][Ny/2][Nz/2+D/2+loop1])+dy*(Hyz[Nx/2][Ny/2][Nz/2+D/2+loop1]-Hyz[Nx/2-1][Ny/2][Nz/2+D/2+loop1]))*1000;
	 In2[loop1]=(dx*(Hxy[Nx/2][Ny/2-1][Nz/2-D/2-loop1]-Hxy[Nx/2][Ny/2][Nz/2-D/2-loop1])+dy*(Hyz[Nx/2][Ny/2][Nz/2-D/2-loop1]-Hyz[Nx/2-1][Ny/2][Nz/2-D/2-loop1]))*1000;
	 if(loop==NN-100)
	 {  current1[loop1].zhengfu=In1[loop1];current1[loop1].xiangwei=0;
	    current2[loop1].zhengfu=In2[loop1];current2[loop1].xiangwei=0;
	 }
	 else
	 {
	     if(In1[loop1]>current1[loop1].zhengfu)
	        {current1[loop1].zhengfu=In1[loop1];current1[loop1].xiangwei=(loop-NN+100)/100.0*2*Pi;}
		 if(In2[loop1]>current2[loop1].zhengfu)
	        {current2[loop1].zhengfu=In2[loop1];current2[loop1].xiangwei=(loop-NN+100)/100.0*2*Pi;}
	 }
     }
	}
//	if(loop>=NN)
//		waitui(loop);
		if(loop>=NN-100)
        distill(loop);

 }
  printf("************************************\n");
//  printf("频率为:%f时,第%d 运算",f,lp);
  printf("************************************\n");
   ffft(Uzz1,Izz1);
   zukang();
   zhubobi(); 
  
 }  
    waitui(loop);
    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]);
        
   
   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-5][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   ",Exy[j][Ny/2-5][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   ",Exy[j][k][Nz/2]);
	   }
	   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-2*NNN-2*DDD+1;j++)
	   {        for (k=0;k<Nx-2*NNN-2*DDD;k++)
				   fprintf(out,"%f ",Jyu[j][k].zhengfu);
                   fprintf(out,";\n");
	   }  
	   fprintf(out,"\n\n\n");
       for (j=0;j<Nx-2*NNN-2*DDD+1;j++)
	   {        for (k=0;k<Nx-2*NNN-2*DDD;k++)
				   fprintf(out,"%f ",Jyd[j][k].zhengfu);
                   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-2*NNN-2*DDD+1;j++)
	   {        for (k=0;k<Nx-2*NNN-2*DDD;k++)
				   fprintf(out,"%f ",Jmyu[j][k].zhengfu);
                   fprintf(out,";\n");
	   }  
	   fprintf(out,"\n\n\n");
       for (j=0;j<Nx-2*NNN-2*DDD+1;j++)
	   {        for (k=0;k<Nx-2*NNN-2*DDD;k++)
				   fprintf(out,"%f ",Jmyd[j][k].zhengfu);
                   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-2*NNN-2*DDD;j++)
	   {        for (k=0;k<Nx-2*NNN-2*DDD+1;k++)
				   fprintf(out,"%f ",Jxu[j][k].zhengfu);
                   fprintf(out,";\n");
	   }  
	   fprintf(out,"\n\n\n");
       for (j=0;j<Nx-2*NNN-2*DDD+1;j++)
	   {        for (k=0;k<Nx-2*NNN-2*DDD;k++)
				   fprintf(out,"%f ",Jxd[j][k].zhengfu);
                   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-2*NNN-2*DDD;j++)
	   {        for (k=0;k<Nx-2*NNN-2*DDD+1;k++)
				   fprintf(out,"%f ",Jmxu[j][k].zhengfu);
                   fprintf(out,";\n");
	   }  
	   fprintf(out,"\n\n\n");
       for (j=0;j<Nx-2*NNN-2*DDD+1;j++)
	   {        for (k=0;k<Nx-2*NNN-2*DDD;k++)
				   fprintf(out,"%f ",Jmxd[j][k].zhengfu);
                   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<=DD/2;j++)
	      fprintf(out,"%f      %f\n",current1[j].zhengfu,current1[j].xiangwei);
          fprintf(out,"\n \n");
	    for(j=0;j<=DD/2;j++)
	      fprintf(out,"%f      %f\n",current2[j].zhengfu,current2[j].xiangwei);

	 }

   printf("enter the outfile name:\n");
   scanf("%s",outfile);
   if ((out=fopen(outfile,"w"))==NULL)
   {printf("cannot open outfile\n");
    goto end;
   }
   else