1_4_test.cpp

FDTD 1-d C++ sourece

/* FDTD_1_4.cpp*/
/*simulation of a sinusoidal wave 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 */
 # define pi 3.1415926

 void main ()
 {
  float ex[KE],hy[KE],cb[KE];
  float ex_low_m1,ex_low_m2,ex_high_m1,ex_high_m2,ex_high_m3,ex_high_m4;
  int n,k,kc,ke,kstart,NSTEPS;
  int freq_in;
  float T,ddx,dt,espz,epsilon,sigma,eaf;

  float t0,spread,pulse;
  FILE *fp;

  /* Initialize */
  for (k = 0; k <KE; k++)
  {
     ex[k]=0;
     hy[k]=0;
     cb[k]=.5;                  /*Initialize to free space*/
  }

  printf ("Dielectric starts at -->");
  scanf ("%d",&kstart);
  printf ("Epsilon -->");
  scanf ("%f",&epsilon);
  printf ("%d %6.2f \n",kstart,epsilon);
  
  
  ddx = .01;                     /*Set the cell size to 1 cm*/
  dt = ddx/(2*3e8);              /*Calculate the time step*/


  /* These parameters specify the input pulse */
  printf ( "Input freq (MHz) --> ");
  scanf ( "%f",&freq_in);
  freq_in = freq_in*1e6;
  printf ( "%8.0f  \n", freq_in);


  for (k = kstart; k<KE; k++)
  {
     cb[k] = .5/epsilon;
  }


//  for (k = 1; k <= KE; k ++)
//  {
//	  printf ("%2d %4.2f\n", k, cb[k]);
 // }


  kc=KE/2;                    /* Center of the problem space */
  t0=40.0;                    /* Center of the incident pulse */
  spread=12;                  /* Width of the incident pulse */
  T=0;
  NSTEPS=1;
  ex_low_m1=0;
  ex_low_m2=0;
  ex_high_m1=0;
  ex_high_m2=0;
  ex_high_m3 = 0;
  ex_high_m4 = 0;

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

   for (n=1;n<=NSTEPS;n++)
   {
    T=T+1;                        /* T keeps track of the total number */
				 /* of times the main loop is executed */
    /* Main FDTD Loop */

    /* Calculate the Ex field */
    for (k=0;k<KE;k++)
    {
     ex[k]=ex[k]+cb[k]*(hy[k-1]-hy[k]);
    }

    /* Put a sinusoidal source at cell 5 p0 */

    pulse = sin (2*pi*freq_in*dt*T);
    ex[5] = ex[5] + pulse;
    printf ("%5.1f %6.2f  %6.2f\n",T,pulse,ex[5]);

    /*Absorbing Boundary Conditions*/
    ex[0]      =       ex_low_m2;
    ex_low_m2  =       ex_low_m1;
    ex_low_m1  =       ex[1];

    ex[KE-1]   =       ex_high_m4;
    ex_high_m4 =       ex_high_m3;
	ex_high_m3 =       ex_high_m2;
	ex_high_m2 =       ex_high_m1;
    ex_high_m1 =       ex[KE-2];

    /* Calculatee the Hy field */
    for (k=0;k<KE-1;k++)
    {
     hy[k]=hy[k]+.5*(ex[k]-ex[k+1]);
    }
   }

    /* End of the Main FDTD Loop */

	/* At the end of the calculation,print out
	  the Ex and Hy fields */
    for (k=0;k<KE;k++)
    {
     printf ("%3d  %6.2f  %6.2f \n",k,ex[k],hy[k]);
    }

    /* Write the E field out ot a file "Ex" */
    fp=fopen( "Ex","w");
    for (k=0;k<KE;k++)
    {
     fprintf (fp," %6.2f \n",ex[k]);
    }
    fclose(fp);

    /* Write the H field out ot a file "Hy" */
    fp=fopen( "Hy","w");
    for (k=0;k<KE;k++)
    {
     fprintf (fp," %6.2f \n",hy[k]);
    }
    fclose(fp);

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