main.c

该程序是用来计算三维光波传播的FDTD程序。 光源为点源。边界有PML。

/* fd3d_4.1c 3D Dipole in free space */


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

#define IE 80          /*IE is the number of cells to the X lable*/
#define JE 80          /*IE is the number of cells to the Y lable*/
#define KE 80          /*IE is the number of cells to the Z lable*/
#define pi 3.1415926
#define epsz 8.85419e-12

main()
{
   static float gax[IE][JE][KE],gay[IE][JE][KE],gaz[IE][JE][KE];
    static float dx[IE][JE][KE],dy[IE][JE][KE],dz[IE][JE][KE];
   static  float hx[IE][JE][KE],hy[IE][JE][KE],hz[IE][JE][KE];
   static float ex[IE][JE][KE],ey[IE][JE][KE],ez[IE][JE][KE];
   static int i,j,n,jc,ic,kc,k,NSTEPS,npml;
   static  float T,ddx,dt,epsilon,omega;
  static  float t0,spread,pulse;
    FILE *fp,*fopen();

    jc=JE/2;
    ic=IE/2;
    kc=KE/2;
    ddx=.01;             /* Set the cell size to 0.01 m */
    dt=ddx/(2*3e8);      /* Calculate the time stem */

    /*Initialize the arrays */
    for(k=0;k<KE;k++){
        for(j=0;j<JE;j++){
         for(i=0;i<IE;i++){
          ex[i][j][k]=0.0;
          ey[i][j][k]=0.0;
          ez[i][j][k]=0.0;
          dx[i][j][k]=0.0;
          dy[i][j][k]=0.0;
          dz[i][j][k]=0.0;
          hx[i][j][k]=0.0;
          hy[i][j][k]=0.0;
          hz[i][j][k]=0.0;
          gax[i][j][k]=1.0;
          gay[i][j][k]=1.0;
          gaz[i][j][k]=1.0;
        }
        }
    }

        /* Specify the dipole */
       for(k=10;k<21;k++){
            gaz[ic][jc][k]=0.;
        }
            gaz[ic][jc][kc]=1.;

     t0=20;              /*Center of the incident pulse*/
     spread=6.0;          /*Width of the incident pulse*/
     T=0;
     NSTEPS=1;

     while (NSTEPS>0) {
       printf("NSTEPS-->");  /*NSTEPS is the number of times the*/
       scanf("%d",&NSTEPS);  /* main loop has executed*/
       printf("%d\n",NSTEPS);

     for (n=1; n<=NSTEPS; n++)
    {
      T=T+1;           /*T keeps track of the total number*/
      /* Start of the main FDTD Loop    */

     /*   Calculate the Dx field  */
      for(k=1;k<KE;k++){
          for(j=1;j<JE;j++){
            for(i=1;i<IE;i++){
             dx[i][j][k] = dx[i][j][k]+0.5*(hz[i][j][k]- hz[i][j-1][k]
               - hy[i][j][k]+ hy[i][j][k-1]);
        }
    }
        }

      /* Calculate the Dy field */
      for(k=1;k<KE;k++){
          for(j=1;j<JE;j++){
            for(i=1;i<IE;i++){
             dy[i][j][k] = dy[i][j][k]+0.5*(hx[i][j][k]- hx[i][j][k-1]
               - hz[i][j][k]+ hz[i-1][j][k]);
            }
          }
      }
            /* Calculate the Dz field */
        for(k=1;k<KE;k++){
          for(j=1;j<JE;j++){
            for(i=1;i<IE;i++){
             dz[i][j][k] = dz[i][j][k]+0.5*(hy[i][j][k]- hy[i-1][j][k]
               - hx[i][j][k]+ hx[i][j-1][k]);
            }
          }
      }

       /*  Source  */

      pulse = exp(-0.5*(pow((t0-T)/spread,2.0)));
      dz[ic][jc][kc]=pulse;

    /* Calculate the E from D field */
      for(k=1;k<KE-1;k++){
          for(j=1;j<JE-1;j++){
            for(i=1;i<IE-1;i++){
               ex[i][j][k]=gax[i][j][k]*dx[i][j][k];
               ey[i][j][k]=gay[i][j][k]*dy[i][j][k];
               ez[i][j][k]=gaz[i][j][k]*dz[i][j][k];
            }
          }
      }

      /* Calculate the Hx field */
       for(k=1;k<KE-1;k++){
          for(j=1;j<JE-1;j++){
            for(i=1;i<IE-1;i++){
               hx[i][j][k] = hx[i][j][k]+0.5*(ey[i][j][k+1]- ey[i][j][k]
               - ez[i][j+1][k]+ ez[i][j][k]);
            }
          }
       }

        /* Calculate the Hy field */
       for(k=1;k<KE-1;k++){
          for(j=1;j<JE-1;j++){
            for(i=1;i<IE-1;i++){
               hy[i][j][k] = hy[i][j][k]+0.5*(ez[i+1][j][k]- ez[i][j][k]
               - ex[i][j][k+1]+ ex[i][j][k]);
            }
          }
       }

        /* Calculate the Hz field */
       for(k=1;k<KE-1;k++){
          for(j=1;j<JE-1;j++){
            for(i=1;i<IE-1;i++){
               hz[i][j][k] = hz[i][j][k]+0.5*(ex[i][j+1][k]- ex[i][j][k]
               - ey[i+1][j][k]+ ey[i][j][k]);
            }
          }
       }
    }
       /*   End of the Main FDTD Loop  */

   /* write the E field out to a field "Ez"  */
   fp=fopen("Ez40.dat","w");
   for(j=0;j<JE;j++){
       for(i=0;i<IE;i++){
        fprintf(fp,"%11.5f",ez[i][j][kc]);
        }
           fprintf(fp,"  \n");
   }
   fclose(fp);
     }
}