A question about Sullivan's book for 2d PML

2d pml

I am learning FDTD through Sullivan's book---Electromagnetic Simulation Using the FDTD Method, and I can not get the right result from the code listed in the book for fd2d_3.2.c. Also I downloaded this code from his website, it doesn't work too. Can any body help me to figure it out!

Thanks a lot!

Below is the source code from the book:

*******************************************************

/* Fd2d_3.2.c. 2D TM program with the PML */
# include <math.h>
# include <stdlib.h>
# include <stdio.h>

#define IE 140
#define JE 140

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,npml;
float ddx,dt,T,epsz,pi,epsilon,sigma,eaf;
float xn,xxn,xnum,xd,curl_e;
float t0,spread,pulse;
float gi2[IE],gi3[IE];
float gj2[JE],gj3[IE];
float fi1[IE],fi2[IE],fi3[JE];
float fj1[JE],fj2[JE],fj3[JE];
float ihx[IE][JE],ihy[IE][JE];
FILE *fp, *fopen();

ic = IE/2-20;
jc = JE/2-20;
ddx = .01; /* Cell size */
dt =ddx/6e8; /* Time steps */
epsz = 8.8e-12;
pi=3.14159;

/* Initialize the arrays */
for ( j=0; j < JE; j++ ) {
printf( "%2d ",j);
for ( i=0; i < IE; i++ ) {
dz[j]= 0.0 ;
hx[j]= 0.0 ;
hy[j]= 0.0 ;
ihx[j]= 0.0 ;
ihy[j]= 0.0 ;
ga[j]= 1.0 ;
printf( "%5.2f ",ga[j]);
}
printf( " \n");
}

/* Calculate the PML parameters */

for (i=0;i< IE; i++) {
gi2 = 1.0;
gi3 = 1.0;
fi1 = 0.0;
fi2 = 1.0;
fi3 = 1.0;
}
for (j=0;j< IE; j++) {
gj2[j] = 1.0;
gj3[j] = 1.0;
fj1[j] = 0.0;
fj2[j] = 1.0;
fj3[j] = 1.0;
}

printf( "Number of PML cells --> ");
scanf("%d", &npml);

for (i=0;i<= npml; i++) {
xnum = npml - i;
xd = npml;
xxn = xnum/xd;
xn = 0.25*pow(xxn,3.0);
printf(" %d %7.4f %7.4f \n",i,xxn,xn);
gi2 = 1.0/(1.0+xn);
gi2[IE-1-i] = 1.0/(1.0+xn);
gi3 = (1.0 - xn)/(1.0 + xn);
gi3[IE-i-1] = (1.0 - xn)/(1.0 + xn);
xxn = (xnum-.5)/xd;
xn = 0.25*pow(xxn,3.0);
fi1 = xn;
fi1[IE-2-i] = xn;
fi2 = 1.0/(1.0+xn);
fi2[IE-2-i] = 1.0/(1.0+xn);
fi3 = (1.0 - xn)/(1.0 + xn);
fi3[IE-2-i] = (1.0 - xn)/(1.0 + xn);
}

for (j=0;j<= npml; j++) {
xnum = npml - j;
xd = npml;
xxn = xnum/xd;
xn = 0.25*pow(xxn,3.0);
printf(" %d %7.4f %7.4f \n",i,xxn,xn);
gj2[j] = 1.0/(1.0+xn);
gj2[JE-1-j] = 1.0/(1.0+xn);
gj3[j] = (1.0 - xn)/(1.0 + xn);
gj3[JE-j-1] = (1.0 - xn)/(1.0 + xn);
xxn = (xnum-.5)/xd;
xn = 0.25*pow(xxn,3.0);
fj1[j] = xn;
fj1[JE-2-j] = xn;
fj2[j] = 1.0/(1.0+xn);
fj2[JE-2-j] = 1.0/(1.0+xn);
fj3[j] = (1.0 - xn)/(1.0 + xn);
fj3[JE-2-j] = (1.0 - xn)/(1.0 + xn);
}

printf("gi + fi \n");
for (i=0; i< IE; i++) {
printf( "%2d %5.2f %5.2f \n",
i,gi2,gi3),
printf( " %5.2f %5.2f %5.2f \n ",
fi1,fi2,fi3);
}

printf("gj + fj \n");
for (j=0; j< JE; j++) {
printf( "%2d %5.2f %5.2f \n",
j,gj2[j],gj3[j]),
printf( " %5.2f %5.2f %5.2f \n ",
fj1[j],fj2[j],fj3[j]);
}

t0 = 40.0;
spread = 12.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 the Main FDTD loop ---- */

/* Calculate the Dz field */
for ( j=1; j < IE; j++ ) {
for ( i=1; i < IE; i++ ) {
dz[j] = gi3*gj3[j]*dz[j]
+ gi2*gj2[j]*.5*( hy[j] - hy[i-1][j]
- hx[j] + hx[j-1]) ;
}
}

/* Sinusoidal Source */

/* pulse = sin(2*pi*1500*1e6*dt*T);; */
pulse = exp(-.5*pow( (T-t0)/spread,2.));
dz[ic][jc] = pulse;

/* Calculate the Ez field */
/* Leave the Ez edges to 0, as part of the PML */
for ( j=1; j < JE-1; j++ ) {
for ( i=1; i < IE-1; i++ ) {
ez[j] = ga[j]*dz[j] ;
}
}

printf("%3f %6.2f \n ",T,ez[ic][jc]);

/* Calculate the Hx field */
for ( j=0; j < JE-1; j++ ) {
for ( i=0; i < IE; i++ ) {
curl_e = ez[j] - ez[j+1] ;
ihx[j] = ihx[j] + fi1*curl_e ;
hx[j] = fj3[j]*hx[j]
+ fj2[j]*.5*(curl_e + ihx[j]);
}
}

/* Calculate the Hy field */
for ( j=0; j <= JE-1; j++ ) {
for ( i=0; i < IE-1; i++ ) {
curl_e = ez[i+1][j] - ez[j];
ihy[j] = ihy[j] + fj1[j]*curl_e ;
hy[j] = fi3*hy[j]
+ fi2*.5*(curl_e + ihy[j]);
}
}

}
/* ---- End of the main FDTD loop ---- */

for ( j=1; j < JE; j++ ) {
printf( "%2d ",j);
for ( i=1; i <= IE; i++ ) {
printf( "%4.1f",ez[j]);
}
printf( " \n");
}

/* Write the E field out to a file "Ez" */
fp = fopen( "Ez","w");
for ( j=0; j < JE; j++ ) {
for ( i=0; i < IE; i++ ) {
fprintf( fp,"%6.3f ",ez[j]);
}
fprintf( fp," \n");
}

fclose(fp);

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

}

}

Dear,
Pls define the array fields before the main() function.
I think your program will work.
all the best!

Added after 21 seconds:

Dear,
Pls define the array fields before the main() function.
I think your program will work.
all the best!

the code has no problem, just at beginning of declaration instead of
FILE*fp*fopen();

write only
File*fp

compile the program, there will be no error

run the program and insert your input

thanks

Hi, I need to find this algorithm (FDTD 2D TM MODES with the PML Conditions) in matlab code!
Is there somebody that can help me?!

Thanks!!

pippopamp said:

Hi, I need to find this algorithm (FDTD 2D TM MODES with the PML Conditions) in matlab code!
Is there somebody that can help me?!

Thanks!!

Click to expand...

nobody can help me?