fdtd1_3.c

用c语言仿真了FDTD算法,给大家提供一种思路

/*FD1D_1.3.c  1D FDTD simulation */
/*simulation of a pluse hitting a dielectric medium */
#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#define KE 200                /*KE is the number of cells to be used */

main ()
{
  float ex[KE],hy[KE];
  float low_m1,low_m2,high_m1,high_m2;
  int n,k,kc,ke,NSTEPS,kstart;
  float T,ddx,dt,epsz,epsilon,sigma,eaf;
  float cb[KE];
  float m,spread,pulse;
  FILE *fp , *fopen();

  /*initialize*/
  for (k=1; k<KE;k++)
  {
    ex[k]=0;
    hy[k]=0;
    cb[k]=0.5;
  }  
  kc=KE/2;	/*center of the problem space*/
  m=40.0;	/*center of the incident pulse*/
  spread=12;	/*width of the incident pulse*/
  T=0;
  NSTEPS=1;
printf("Dielectric starts at -->");
scanf("%d",&kstart);
printf("Epsilon -->");
scanf("%f",&epsilon);
printf("%d %6.2f  \n",kstart,epsilon);
for (k=kstart;k<=KE;k++)
{
  cb[k]=0.5/epsilon;
}
for (k=1;k<=KE;k++)
{
  printf("%2d %6.3f \n",k,cb[k]);}
  while (NSTEPS>0) {
  printf("NETEPS --> ");	/*NETEPS is the center of times the main loop has excuted*/
  scanf ("%d",&NSTEPS);
  printf("%d\n",NSTEPS);
  n=0;

  for (n=1;n<=NSTEPS;n++)
  {
    
    T=T+1;    /*T keeps track of the totle number of times the main loop is excuted*/
 /*main FDTD loop*/
 /*calculate the EX filed*/
 for (k=1;k<KE;k++)
 {
   ex[k]=ex[k]+cb[k]*(hy[k-1]-hy[k]);}
   /*put a gaussian pulse in the middle*/
   pulse=exp(-0.5*(pow((40-T)/spread,2.0)));
   ex[5]=ex[5]+pulse;
   printf("%5.1f %6.3f %6.3f\n",T,pulse,ex[5]);
   /*Absorbing Boundary Condition */
   ex[0]=low_m2;
   low_m2=low_m1;
   low_m1=ex[1];
   ex[KE-1]=high_m2;
   high_m2=high_m1;
   high_m1=ex[KE-2];
   /*calculate the Hy field */
  for (k=0;k<KE-1;k++)
 {
   hy[k]=hy[k]+0.5*(ex[k]-ex[k+1]);
 }
 /*end of the main FDTD loop*/
 /*at the end of the calculation,print out the ex and hy field*/
 for (k=1;k<KE;k++)
 {
  printf("%3d %6.3f %6.3f\n",k,ex[k],hy[k]); }
  
  /*write the E filed out to a file "EX"*/
  fp=fopen("Ex","w");
  for (k=1;k<=KE;k++)
  {
    fprintf (fp,"%6.3f\n",ex[k]);}
    fclose(fp);
    fp=fopen("Hy","w");
  for (k=1;k<=KE;k++)
  {
    fprintf (fp,"%6.3f\n",hy[k]);}
    fclose(fp);
    printf("T=%5.0f\n",T);
  }
}
}