# include <math.h>
# include <stdlib.h>
# include <stdio.h>
# define IE 60
# define JE 60
int 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;
float dx, dy, dt, T, epsz, pi, epsilon, sigma, mu, eaf;
float t0, spread, pulse;
FILE *fp, *fopen();
ic = IE/2;
jc = JE/2;
dx = .01; /* cell size */
dy = .01; /* cell size */
dt = dx/6e8; /* Time steps */
epsz = 8.8e-12;
pi = 3.14159;
mu = 1.0e-7*4.0*pi;
sigma = 1.0e-8;
for( j = 0; j < JE; j++) {
printf("%2d", j);
for(i = 0; i < IE; i++) {
dz[i][j] = 0.;
ez[i][j] = 0.;
hx[i][j] = 0.;
hy[i][j] = 0.;
ga[i][j] = 1.0;
printf("%5.2f", ga[i][j]);
}
printf("\n");
}
t0 = 20.0;
spread = 6.0;
T = 0;
nsteps = 1;
while(nsteps > 0){
printf("nsteps -> ");
scanf("%d", &nsteps);
printf("%d \n", nsteps);
for(n = 1; n <= nsteps; n++) {
T = T+1;
/* - Start of Main FDTD loop - */
/* - Put a Gaussian pulse in the middle - */
pulse = sin(2*pi*1500*1e6*dt*T);
ez[ic][jc] = pulse;
/* - Calculate Ez field - */
for(j = 1; j < JE; j++){
for(i = 1; i < IE; i++) {
ez[i][j] = (1.0 - sigma*dt/(2.0*epsz))/(1.0+sigma*dt)/(2.0*epsz)*ez[i][j]+(dt/epsz)
/(1.0+sigma*dt/(2.0*epsz*(1/dx))*(hy[i][j] - hy[i - 1][j]) - (dt/epsz)/((1.0+sigma*dt)/(2.0*epsz))*(1/dy)
*(hx[i][j] - hx[i][j - 1]));
}
}
/* - Calculate Hx field - */
for(j = 0; j < JE - 1; j++){
for(i = 0; i < IE - 1; i++){
hx[i][j] = hx[i][j] - (dt/mu*1/dy)*(ez[i][j+1] - ez[i][j]);
}
}
/* Calculate Hy field */
for(j = 0; j < JE - 1; j++){
for(i = 0; i <= IE - 1; i++) {
hy[i][j] = hy[i][j]+(dt/mu*1/dx)*(ez[i+1][j] - ez[i][j]);
}
}
}
/* - End of main FDTD loop - */
for(j = 1; j < jc; j++) {
printf("%2d", j);
for(i = 1; i < ic; i++) {
printf("%5.2f", ez[2*i][2*j]);
}
printf("\n");
}
printf("T = %5.0f \n", T);
/* Write E field out to file "Ez" */
fp = fopen("Ez", "w");
for(j = 0; j < JE; j++) {
for(i = 0; i < IE; i++) {
fprintf(fp,"%d\t%d\t%6.3f",i, j, ez[i][j]);
}
fprintf(fp, "\n");
}
fclose(fp);
}
}