2dtm.cpp

2维Tm极化波的FDTD算法包含吸收边界条件PML

/*FDTD_3.1.c.  2D TM program*/


# include <math.h>
# include <stdlib.h>
# include <stdio.h>

#define IE 60
#define JE 60

main()
{
	 float ga[IE+1][JE+1],dz[IE+1][JE+1],ez[IE+1][JE+1],hx[IE+1][JE+1],hy[IE+1][JE+1];
	 int l,n,i,j,ic,jc,nsteps;
	 float ddx,dt,T,epsz,epsilon,sigma,eaf,pi;
	 float t0,spread,pulse;
	 FILE *fp;

	 /*Initialize */
	 ic = IE/2;
	 jc = JE/2;
	 ddx = .01;        /*cell size*/
	 dt = ddx/6e8;     /*time steps*/
	 epsz = 8.8e-12;
	 pi = 3.14159;
	 
	 t0=20.0;                  /*Center of the incident pulse*/
	 spread=6.0;                /*Width of  the incident pulse*/
	 T=0;
	 nsteps=1;

	  for (j=0;j<=JE;j++)
	 {   
	   printf("%2d",j);
	   for (i=0;i<=IE;i++)
	   {
		   dz[i][j] = 0;
		   ez[i][j] = 0;
           hx[i][j] = 0;
		   hy[i][j] = 0;
	       ga[i][j] = 1.0;
		   printf("%5.2f",ga[i][j]);
	   }
	   printf("\n");
	 }
 
	 
	 
	 
	 while (nsteps>0) 
	 {
		 printf("nsteps --> ");    /*nsteps is the number of times  the main loop has executed*/
		 scanf("%d",&nsteps);      
		 printf("%d \n",nsteps);
		 n=0;

		 for (n=1;n<=nsteps;n++)
		 {
			 T=T+1;               /*T keeps track of the total number of the */
			                     /*main loop has executed*/

		 /* Main FDTD Loop*/
			 /*Calculate the Dz field*/
			 for (j=0;j<JE;j++)
			 {
				 for(i=0;i<IE;i++)
				 {
					 dz[i][j] = dz[i][j]+.5*(hy[i][j]-hy[i-1][j]-hx[i][j]+hx[i][j-1]);
				 }
			 }
			 

			 /*put a Gaussian pulse in the middle */
			 pulse = exp(-.5*(pow((t0-T)/spread,2.0)));
			 dz[ic][jc] = pulse;                       /*put a sinusoidal soft source at cell 5*/
			 
             /*Calculate the Ez field*/
			 for (j=1;j<JE;j++)
			 {
				 for(i=1;i<IE;i++)
				 {
					 ez[i][j] = ga[i][j]*dz[i][j];
				 }
			 }

			 /*Calculate the Hx field*/
			 for (j=0;j<=JE-1;j++)
			 {
				 for(i=0;i<=IE;i++)
				 {
					 hx[i][j] = hx[i][j]+.5*(ez[i][j]-ez[i][j+1]);
				 }
			 }
			 
	         
			 /*Calculate the Hy field*/
			 for (j=0;j<=JE;j++)
			 {
				 for(i=0;i<IE;i++)
				 {
					 hy[i][j] = hy[i][j]+.5*(ez[i+1][j]-ez[i][j]);
				 }
			 }
             
		 }

        /*End of the Main FDTD Loop*/

          /*At the end of the calculation, printf out the Ez fields*/
          for (j=0;j<=jc;j++)
		  {
			  printf("%2d",j);
			  for(i=1;i<ic;i++)
			  {
				  printf("%1.10f\n",ez[i][j]);
			  }

			  printf("\n");
		  }
          printf("T= %5.0f\n",T);

		  /*Write the E field out to a file "Ez"*/
		  fp = fopen("e:\\work\\Ez","w");
		  for (j=0;j<=JE;j++)
		  {
			  for (i=0;i<=IE;i++)
			  {
				  fprintf(fp," %0.10f\n",ez[i][j]);
			  }
			  fprintf(fp,"\n");
		  }

		  fclose(fp);

		  
		  
		 

	 }
	 
}