fd2d_32.cpp

二维的光子晶体能带计算

/* fd2d_3.2.c  2D TM  program with the PML*/
# include <math.h>
# include <stdlib.h>
# include <stdio.h>

# define IE  60
# define JE  60

void main()
{ 
	float ga[IE][JE],dz[IE][JE],ez[IE][JE],hx[IE][JE],hy[IE][JE];
	int   l,n,i,j,ic,jc,nsteps,npml;
    float ddx,dt,T,epsz,pi,epsilon,sigma,eaf;
	float xn,xxn,xnum,xd,curl_e;
	float t0,spread,pulse;  
	float gi2[IE],gi3[IE];
	float gj2[JE],gj3[JE];
	float fi1[IE],fi2[IE],fi3[IE];
	float fj1[JE],fj2[JE],fj3[JE];
	float ihx[IE][JE],ihy[IE][JE];
	FILE  *fp;
   
    ic=IE/2-5;
	jc=JE/2-5;
	pi=3.14159;
    ddx=0.01;       // 设置网格尺寸为1cm
	dt=ddx/(2*3e8); // 计算时间步长 
	epsz=8.8e-12; 


    /* 初始化 */
	for(j=0; j<JE; j++)
	{ printf("%2d  ",j);
	  for (i=0; i<IE; i++)
	  {dz[i][j]=0.0;
	   ez[i][j]=0.0;
	   hx[i][j]=0.0;
	   hy[i][j]=0.0;
	   ihx[i][j]=0.0;
	   ihy[i][j]=0.0;
	   ga[i][j]=1.0;
	   printf("%5.2f ",ga[i][j]);
	  }
	  printf(" \n");
	  
	}
   /*Calculate the PML parameters */
	for (i=0; i<IE; i++){
		gi2[i] = 1.0;
		gi3[i] = 1.0;
		fi1[i] = 0.0;
		fi2[i] = 1.0;
		fi3[i] = 1.0;
	}
	for (j=0; j<JE; j++){
		gj2[j] = 1.0;
		gj3[j] = 1.0;
		fj1[j] = 0.0;
		fj2[j] = 1.0;
		fj3[j] = 1.0;
	}

   printf( "Number of PML cells -->");
   scanf("%d",&npml);
   printf("Number = %d\n", npml);

   for (i=0; i<=npml; i++){
	   xnum = npml - i;
	   xd  = npml;
	   xxn = xnum/xd;
	   xn  = 0.33*pow(xxn,3.0);
	   printf("%d  %7.4f  %7.4f \n",i,xxn,xn);
	   gi2[i] = 1.0/(1.0+xn);
	   gi2[IE-1-i] =  1.0/(1.0+xn);
	   gi3[i] = (1.0-xn)/(1.0+xn);
	   gi3[IE-i-1] = (1.0-xn)/(1.0+xn);
	xxn = (xnum-.5)/xd;
	xn = 0.25*pow(xxn,3.0);
	   fi1[i] = xn;
	   fi1[IE-2-i] = xn;
	   fi2[i] = 1.0/(1.0+xn);
	   fi2[IE-2-i] =  1.0/(1.0+xn);
	   fi3[i] = (1.0-xn)/(1.0+xn);
	   fi3[IE-2-i] = (1.0-xn)/(1.0+xn);
	   
   }
   for(i=0; i<IE; i++)
  { printf("%d  %7.4f  %7.4f  %7.4f  \n", i,fi1[i],fi2[i],fi3[i]);}

      for (j=0; j<=npml; j++){
	   xnum = npml - j;
	   xd  = npml;
	   xxn = xnum/xd;
	   xn  = 0.33*pow(xxn,3.0);
	   printf("%d  %7.4f  %7.4f \n",j,xxn,xn);
	   gj2[j] = 1.0/(1.0+xn);
	   gj2[JE-1-j] =  1.0/(1.0+xn);
	   gj3[j] = (1.0-xn)/(1.0+xn);
	   gj3[JE-j-1] = (1.0-xn)/(1.0+xn);
	xxn = (xnum-.5)/xd;
	xn = 0.25*pow(xxn,3.0);
	   fj1[j] = xn;
	   fj1[JE-2-j] = xn;
	   fj2[j] = 1.0/(1.0+xn);
	   fj2[JE-2-j] =  1.0/(1.0+xn);
	   fj3[j] = (1.0-xn)/(1.0+xn);
	   fj3[JE-2-j] = (1.0-xn)/(1.0+xn);
   }
/*************这些参数说明输入脉冲******************/
	t0=40.0;
	spread=15;
	T=0;
	nsteps=1; 
 
 
	while( nsteps > 0)
	{
		printf( "nsteps--> ");
		scanf("%d",&nsteps);
		printf("%d \n",nsteps);

    

   /*********FDTD主程序********/
		for( n=1; n<= nsteps; n++)
		{
			T=T+1;

			/**** Calculate the Dz field****/

			for (j=1; j<JE; j++) { 
			  for (i=1; i<IE; i++){
				  dz[i][j]=gi3[i]*gj3[j]*dz[i][j]+gi2[i]*gj2[j]*0.5*(hy[i][j]-hy[i-1][j]-hx[i][j]+hx[i][j-1]);
			  }
			}

            /**** Put a Sinusoidal pulse in the middle****/
			//pulse =  exp(-.5*(pow((t0-T)/spread,2.0)));
			pulse = sin(2*pi*1500*1e6*dt*T);
			dz[ic][jc]= pulse;
		   
			/******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];
			  }
			}

			/*Set the Ez edges to 0, as part of the PML */
			for (j=0; j<JE-1; j++){
			ez[0][j] = 0.0;
			ez[IE-1][j] = 0.0;}
			for (i=0; i<IE-1; i++){
			ez[i][0] = 0.0;
			ez[i][JE-1] = 0.0;}

            /******Calculate the Hx field*******/
			for (j=0; j<JE-1; j++) { 
			  for (i=0; i<IE; i++){
                curl_e = ez[i][j] -ez[i][j+1];
				ihx[i][j] = ihx[i][j] +fi1[i]*curl_e;
				hx[i][j]=fj3[j]*hx[i][j]+fj2[j]*0.5*(curl_e + ihx[i][j]);
			  }
			}   
              /******Calculate the Hy field*******/
			for (j=0; j<=JE-1; j++) { 
			  for (i=0; i<IE-1; i++){
                curl_e = ez[i+1][j] -ez[i][j];
				ihy[i][j] = ihy[i][j] + fj1[j]*curl_e;
				hy[i][j]=fi3[i]*hy[i][j]+fi2[i]*0.5*(curl_e + ihy[i][j]);
			  }
			}   

		}

   /*********FDTD主程序结束********/


	   /***********输出电场和磁场值************/

	    for (j=1; j<JE; j++ ){
		  printf( "%2d  ",j);
		  for (i=1; i<=IE; i++){
			  printf("%4.1f ",ez[i][j]);
		  }
		  printf("\n");

		}

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

       /**********把电场值写到Ez文件************/
     	fp=fopen( "Ez.dat","w");
		for (j=1; j<JE; j++) { 
			  for (i=1; i<IE; i++){
				fprintf(fp," %6.3f",ez[i][j]);
			  }
            fprintf(fp," \n");
		} 
        fclose(fp);

	}
 }