fdtd1_5.c

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

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

#define KE 200                /*KE is the number of cells to be used */
#define pi 3.1415926535
main ()
{
  float ex[KE],hy[KE];
  float low_m1,low_m2,high_m1,high_m2;
  int   n,k,NSTEPS,kstart;
  float T,ddx,dt,epsz,epsilon,sigma,eaf;
  float cb[KE],ca[KE];
  float pulse;
  float freq_in;
  FILE *fp , *fopen();

  /*initialize*/
  for (k=1; k<KE;k++)
  {
    ex[k]=0;
    hy[k]=0;
    ca[k]=1;
    cb[k]=0.5;
  }  
  
        T=0;
	NSTEPS=1;
	ddx=0.01;     /*set the cell size to 1 cm*/
	dt=ddx/(2*3e8);   /*calculate the time step*/
  	epsz=8.85419e-12;
printf("Dielectric starts at -->");
scanf("%d",&kstart);
printf("Epsilon -->");
scanf("%f",&epsilon);
printf("Conductivity -->");
scanf("%f",&sigma);
printf("%d %6.2f %6.2f \n",kstart,epsilon,sigma);
eaf=dt*sigma/(2*epsz*epsilon);
printf("%6.4f\n",eaf);
for (k=kstart;k<=KE;k++)
{
ca[k]=(1.-eaf)/(1.+eaf);  
cb[k]=0.5/(epsilon*(1+eaf));
}
for (k=1;k<=KE;k++)
{
  printf("%2d %6.3f %6.3f \n",k,ca[k],cb[k]);}
	/*These parameters specify the input pulse*/
	printf("Input freq_in (MHz) -->");
	scanf("%f",&freq_in);
	freq_in=freq_in*1e6;
	printf("%8.0f\n",freq_in);

  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;
	ddx=0.01;     /*set the cell size to 1 cm*/
	dt=ddx/(2*3e8);   /*calculate the time step*/
  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]=ca[k]*ex[k]+cb[k]*(hy[k-1]-hy[k]);}
   /*put a sinusoidal source at the cell 5 p0*/
   pulse=sin(2*pi*freq_in*dt*T);
   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);
  }
}
}