S
sixfu
Guest
Ik ben het leren FDTD via Sullivan's boek --- Elektromagnetische simulatieprofiel Met de FDTD methode, en ik kan niet het juiste resultaat van de code opgenomen in het boek voor fd2d_3.2.c.Ook ik gedownload heb deze code van zijn website, die niet werkt ook.Kan een lichaam mij helpen bij het uit!
Thanks a lot!
Hieronder staat de broncode van het boek:
************************************************** *****
/ * Fd2d_3.2.c.2D TM-programma met de 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, gespreid, polsslag;
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 grootte * /
dt = ddx/6e8; / * Tijd stappen * /
epsz = 8.8e-12;
pi = 3,14159;
/ * Initialize de 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");
)
/ * Bereken de 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 ( "Aantal PML cellen ->");
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;
verspreiding = 12,0;
T = 0;
nsteps = 1;
while (nsteps> 0) (
printf ( "nsteps ->");
scanf ( "% d", & nsteps);
printf ( "% d \ n", nsteps);
voor (n = 1; n <= nsteps; n ) (
T = T 1;
/ * ---- Start van de Algemene FDTD lus ---- * /
/ * Bereken het gebied Dz * /
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]);
)
)
/ * Sinusvormige Bron * /
/ * Puls = sin (2 * pi * 1500 * 1e6 * dt * T); * /
puls = exp (-. 5 * pow ((T-t0) / uitbreiden, 2.));
dz [ic] [jc] = impuls;
/ * Bereken het gebied Ez * /
/ * Laat de Ez randen op 0, als onderdeel van de PML * /
for (j = 1; j <JE-1; j ) (
for (i = 1; i <IE-1; i ) (
ez [j] = g [j] * dz [j];
)
)
printf ( "% 3f% 6.2f \ n", T, ez [ic] [jc]);
/ * Bereken de Hx gebied * /
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]);
)
)
/ * Bereken de Hy gebied * /
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]);
)
)
)
/ * ---- Einde van de belangrijkste FDTD lus ---- * /
for (j = 1; j <JE; j ) (
printf ( "% 2d", j);
for (i = 1; i <= IE; i ) (
printf ( "% 4.1f", EZ [j]);
)
printf ( "\ n");
)
/ * Schrijf de E gebied uit tot een bestand "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);
)
)
Thanks a lot!
Hieronder staat de broncode van het boek:
************************************************** *****
/ * Fd2d_3.2.c.2D TM-programma met de 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, gespreid, polsslag;
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 grootte * /
dt = ddx/6e8; / * Tijd stappen * /
epsz = 8.8e-12;
pi = 3,14159;
/ * Initialize de 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");
)
/ * Bereken de 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 ( "Aantal PML cellen ->");
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;
verspreiding = 12,0;
T = 0;
nsteps = 1;
while (nsteps> 0) (
printf ( "nsteps ->");
scanf ( "% d", & nsteps);
printf ( "% d \ n", nsteps);
voor (n = 1; n <= nsteps; n ) (
T = T 1;
/ * ---- Start van de Algemene FDTD lus ---- * /
/ * Bereken het gebied Dz * /
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]);
)
)
/ * Sinusvormige Bron * /
/ * Puls = sin (2 * pi * 1500 * 1e6 * dt * T); * /
puls = exp (-. 5 * pow ((T-t0) / uitbreiden, 2.));
dz [ic] [jc] = impuls;
/ * Bereken het gebied Ez * /
/ * Laat de Ez randen op 0, als onderdeel van de PML * /
for (j = 1; j <JE-1; j ) (
for (i = 1; i <IE-1; i ) (
ez [j] = g [j] * dz [j];
)
)
printf ( "% 3f% 6.2f \ n", T, ez [ic] [jc]);
/ * Bereken de Hx gebied * /
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]);
)
)
/ * Bereken de Hy gebied * /
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]);
)
)
)
/ * ---- Einde van de belangrijkste FDTD lus ---- * /
for (j = 1; j <JE; j ) (
printf ( "% 2d", j);
for (i = 1; i <= IE; i ) (
printf ( "% 4.1f", EZ [j]);
)
printf ( "\ n");
)
/ * Schrijf de E gebied uit tot een bestand "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);
)
)