📄 fdtd.cpp
字号:
if (k > HEIGHT)
Ey(i,j,k) = Ey(i,j,k) +
coefE0*((Hx(i,j,k)-Hx(i,j,k-1)) - (Hz(i,j,k)-Hz(i-1,j,k)));
else if (k == HEIGHT)
Ey(i,j,k) = Ey(i,j,k) +
coefE01*((Hx(i,j,k)-Hx(i,j,k-1)) - (Hz(i,j,k)-Hz(i-1,j,k)));
else
Ey(i,j,k) = Ey(i,j,k) +
coefE1*((Hx(i,j,k)-Hx(i,j,k-1)) - (Hz(i,j,k)-Hz(i-1,j,k)));
#ifdef PATCH
/* Zero Ey on metal. */
k=HEIGHT;
/* First the feed strip. */
for (i=LINE_X_START; i<=LINE_X_END; i++)
for (j=0; j<LIMY/2-10; j++)
Ey(i,j,k) = 0.0;
/* Next the patch. */
for (i=PATCH_X_START; i<=PATCH_X_END; i++)
for (j=PATCH_Y_START; j<PATCH_Y_END; j++)
Ey(i,j,k) = 0.0;
#else
/* Zero Ey on metal microstrip. */
k=HEIGHT;
for (i=LINE_X_START; i<=LINE_X_END; i++)
for (j=0; j<LIMY; j++)
Ey(i,j,k) = 0.0;
#endif
#ifdef ABC1
/* ABC on left wall */
/* Ey above substrate */
i = 0;
for (j=0; j<LIMY-1; j++) {
for (k=HEIGHT+1; k<LIMZ; k++) {
Ey(i,j,k) = c10*(Ey(i+2,j,k)+Eyleft(i,j,k,1))
+ c20*(Eyleft(i,j,k,0) + Eyleft(i+2,j,k,0) - Ey(i+1,j,k)
- Eyleft(i+1,j,k,1)) + c30*Eyleft(i+1,j,k,0) - Eyleft(i+2,j,k,1);
for(ii=0; ii<3; ii++) {
Eyleft(ii,j,k,1) = Eyleft(ii,j,k,0);
Eyleft(ii,j,k,0) = Ey(ii,j,k);
}
}
}
/* Ey at substrate */
i=0;
for (j=0; j<LIMY-1; j++) {
k=HEIGHT;
Ey(i,j,k) = c101*(Ey(i+2,j,k)+Eyleft(i,j,k,1))
+ c201*(Eyleft(i,j,k,0) + Eyleft(i+2,j,k,0) - Ey(i+1,j,k)
- Eyleft(i+1,j,k,1)) + c301*Eyleft(i+1,j,k,0) - Eyleft(i+2,j,k,1);
for(ii=0; ii<3; ii++) {
Eyleft(ii,j,k,1) = Eyleft(ii,j,k,0);
Eyleft(ii,j,k,0) = Ey(ii,j,k);
}
}
/* Ey below substrate */
i=0;
for (j=0; j<LIMY-1; j++) {
for (k=1; k<HEIGHT; k++) {
Ey(i,j,k) = c11*(Ey(i+2,j,k)+Eyleft(i,j,k,1))
+ c21*(Eyleft(i,j,k,0) + Eyleft(i+2,j,k,0) - Ey(i+1,j,k)
- Eyleft(i+1,j,k,1)) + c31*Eyleft(i+1,j,k,0) - Eyleft(i+2,j,k,1);
for(ii=0; ii<3; ii++) {
Eyleft(ii,j,k,1) = Eyleft(ii,j,k,0);
Eyleft(ii,j,k,0) = Ey(ii,j,k);
}
}
}
#endif
#ifdef ABC3
/* ABC on right wall */
/* Ey above substrate */
i = LIMX-1;
for (j=0; j<LIMY-1; j++) {
for (k=HEIGHT+1; k<LIMZ; k++) {
Ey(i,j,k) = c10*(Ey(i-2,j,k)+Eyright(i,j,k,1))
+ c20*(Eyright(i,j,k,0) + Eyright(i-2,j,k,0) - Ey(i-1,j,k)
- Eyright(i-1,j,k,1)) + c30*Eyright(i-1,j,k,0) - Eyright(i-2,j,k,1);
for(ii=LIMX-3; ii<LIMX; ii++) {
Eyright(ii,j,k,1) = Eyright(ii,j,k,0);
Eyright(ii,j,k,0) = Ey(ii,j,k);
}
}
}
/* Ey at substrate */
i=LIMX-1;
for (j=0; j<LIMY-1; j++) {
k=HEIGHT;
Ey(i,j,k) = c101*(Ey(i-2,j,k)+Eyright(i,j,k,1))
+ c201*(Eyright(i,j,k,0) + Eyright(i-2,j,k,0) - Ey(i-1,j,k)
- Eyright(i-1,j,k,1)) + c301*Eyright(i-1,j,k,0) - Eyright(i-2,j,k,1);
for(ii=LIMX-3; ii<LIMX; ii++) {
Eyright(ii,j,k,1) = Eyright(ii,j,k,0);
Eyright(ii,j,k,0) = Ey(ii,j,k);
}
}
/* Ey below substrate */
i=LIMX-1;
for (j=0; j<LIMY-1; j++) {
for (k=1; k<HEIGHT; k++) {
Ey(i,j,k) = c11*(Ey(i-2,j,k)+Eyright(i,j,k,1))
+ c21*(Eyright(i,j,k,0) + Eyright(i-2,j,k,0) - Ey(i-1,j,k)
- Eyright(i-1,j,k,1)) + c31*Eyright(i-1,j,k,0) - Eyright(i-2,j,k,1);
for(ii=LIMX-3; ii<LIMX; ii++) {
Eyright(ii,j,k,1) = Eyright(ii,j,k,0);
Eyright(ii,j,k,0) = Ey(ii,j,k);
}
}
}
#endif
#ifdef ABC5
/* ABC at top */
for (i=0; i<LIMX; i++)
for (j=0; j<LIMY-1; j++) {
k=LIMZ-1;
Ey(i,j,k) = c10*(Ey(i,j,k-2)+Eytop(i,j,k,1))
+ c20*(Eytop(i,j,k,0) + Eytop(i,j,k-2,0) - Ey(i,j,k-1)
- Eytop(i,j,k-1,1)) + c30*Eytop(i,j,k-1,0) - Eytop(i,j,k-2,1);
for(kk=LIMZ-3; kk<LIMZ; kk++) {
Eytop(i,j,kk,1) = Eytop(i,j,kk,0);
Eytop(i,j,kk,0) = Ey(i,j,kk);
}
}
#endif
/************ Ez update. *************/
for (i=1; i<LIMX-1; i++)
for (j=1; j<LIMY-1; j++)
for (k=0; k<LIMZ-1; k++)
if (k > HEIGHT)
Ez(i,j,k) = Ez(i,j,k)
+ coefE0*((Hy(i,j,k)-Hy(i-1,j,k)) - (Hx(i,j,k)-Hx(i,j-1,k)));
else
Ez(i,j,k) = Ez(i,j,k)
+ coefE1*((Hy(i,j,k)-Hy(i-1,j,k)) - (Hx(i,j,k)-Hx(i,j-1,k)));
#ifdef ABC4
if (ntime<SWITCH_SRC) {
#endif
/* Ez update on "near" wall. */
j = 0;
for (i=1; i<LIMX; i++)
for (k=0; k<LIMZ; k++)
if (k > HEIGHT)
Ez(i,j,k) = Ez(i,j,k)
+ coefE0*((Hy(i,j,k)-Hy(i-1,j,k)) - 2.0*Hx(i,j,k));
else
Ez(i,j,k) = Ez(i,j,k)
+ coefE1*((Hy(i,j,k)-Hy(i-1,j,k)) - 2.0*Hx(i,j,k));
/* source at y=0 wall of computational domain */
holdez = gaussian(ntime, cdtds, PPW)/3.0;
j=0;
for (i=LINE_X_START; i<=LINE_X_END; i++)
for (k=0; k<HEIGHT; k++)
Ez(i,j,k) = holdez;
/* Ez(i,j,k) = Ez(i,j,k) + holdez; */
#ifdef ABC4
}
#endif
#ifdef ABC1
/* ABC on left wall */
/* Ez above substrate */
i = 0;
for (j=0; j<LIMY-1; j++) {
for (k=HEIGHT; k<LIMZ; k++) {
Ez(i,j,k) = c10*(Ez(i+2,j,k)+Ezleft(i,j,k,1))
+ c20*(Ezleft(i,j,k,0) + Ezleft(i+2,j,k,0) - Ez(i+1,j,k)
- Ezleft(i+1,j,k,1)) + c30*Ezleft(i+1,j,k,0) - Ezleft(i+2,j,k,1);
for(ii=0; ii<3; ii++) {
Ezleft(ii,j,k,1) = Ezleft(ii,j,k,0);
Ezleft(ii,j,k,0) = Ez(ii,j,k);
}
}
}
/* Ez below substrate */
i=0;
for (j=0; j<LIMY-1; j++) {
for (k=0; k<HEIGHT; k++) {
Ez(i,j,k) = c11*(Ez(i+2,j,k)+Ezleft(i,j,k,1))
+ c21*(Ezleft(i,j,k,0) + Ezleft(i+2,j,k,0) - Ez(i+1,j,k)
- Ezleft(i+1,j,k,1)) + c31*Ezleft(i+1,j,k,0) - Ezleft(i+2,j,k,1);
for(ii=0; ii<3; ii++) {
Ezleft(ii,j,k,1) = Ezleft(ii,j,k,0);
Ezleft(ii,j,k,0) = Ez(ii,j,k);
}
}
}
#endif
#ifdef ABC2
/* ABC at far wall */
/* Ez above substrate */
for (i=0; i<LIMX-1; i++) {
j=LIMY-1;
for (k=HEIGHT; k<LIMZ; k++) {
Ez(i,j,k) = c10*(Ez(i,j-2,k)+Ezfar(i,j,k,1))
+ c20*(Ezfar(i,j,k,0) + Ezfar(i,j-2,k,0) - Ez(i,j-1,k)
- Ezfar(i,j-1,k,1)) + c30*Ezfar(i,j-1,k,0) - Ezfar(i,j-2,k,1);
for(jj=LIMY-3; jj<LIMY; jj++) {
Ezfar(i,jj,k,1) = Ezfar(i,jj,k,0);
Ezfar(i,jj,k,0) = Ez(i,jj,k);
}
}
}
/* Ez below substrate */
for (i=0; i<LIMX-1; i++) {
j=LIMY-1;
for (k=0; k<HEIGHT; k++) {
Ez(i,j,k) = c11*(Ez(i,j-2,k)+Ezfar(i,j,k,1))
+ c21*(Ezfar(i,j,k,0) + Ezfar(i,j-2,k,0) - Ez(i,j-1,k)
- Ezfar(i,j-1,k,1)) + c31*Ezfar(i,j-1,k,0) - Ezfar(i,j-2,k,1);
for(jj=LIMY-3; jj<LIMY; jj++) {
Ezfar(i,jj,k,1) = Ezfar(i,jj,k,0);
Ezfar(i,jj,k,0) = Ez(i,jj,k);
}
}
}
#endif
#ifdef ABC4
/* ABC at near wall -- only apply after source introduced */
if (ntime >= SWITCH_SRC) {
/* Ez above substrate */
for (i=0; i<LIMX-1; i++) {
j=0;
for (k=HEIGHT; k<LIMZ; k++) {
Ez(i,j,k) = Eznear(i,j+1,k,0) + c40*(Ez(i,j+1,k) - Ez(i,j,k));
/* Ez(i,j,k) = c10*(Ez(i,j+2,k)+Eznear(i,j,k,1))
+ c20*(Eznear(i,j,k,0) + Eznear(i,j+2,k,0) - Ez(i,j+1,k)
- Eznear(i,j+1,k,1)) + c30*Eznear(i,j+1,k,0) - Eznear(i,j+2,k,1);
*/
for(jj=2; jj>=0; jj--) {
Eznear(i,jj,k,1) = Eznear(i,jj,k,0);
Eznear(i,jj,k,0) = Ez(i,jj,k);
}
}
}
/* Ez below substrate */
for (i=0; i<LIMX-1; i++) {
j=0;
for (k=0; k<HEIGHT; k++) {
Ez(i,j,k) = Eznear(i,j+1,k,0) + c41*(Ez(i,j+1,k) - Ez(i,j,k));
/*
Ez(i,j,k) = c11*(Ez(i,j+2,k)+Eznear(i,j,k,1))
+ c21*(Eznear(i,j,k,0) + Eznear(i,j+2,k,0) - Ez(i,j+1,k)
- Eznear(i,j+1,k,1)) + c31*Eznear(i,j+1,k,0) - Eznear(i,j+2,k,1);
*/
for(jj=2; jj>=0; jj--) {
Eznear(i,jj,k,1) = Eznear(i,jj,k,0);
Eznear(i,jj,k,0) = Ez(i,jj,k);
}
}
}
}
#endif
#ifdef ABC3
/* ABC on right wall */
/* Ez above substrate */
i = LIMX-1;
for (j=0; j<LIMY-1; j++) {
for (k=HEIGHT; k<LIMZ; k++) {
Ez(i,j,k) = c10*(Ez(i-2,j,k)+Ezright(i,j,k,1))
+ c20*(Ezright(i,j,k,0) + Ezright(i-2,j,k,0) - Ez(i-1,j,k)
- Ezright(i-1,j,k,1)) + c30*Ezright(i-1,j,k,0) - Ezright(i-2,j,k,1);
for(ii=LIMX-3; ii<LIMX; ii++) {
Ezright(ii,j,k,1) = Ezright(ii,j,k,0);
Ezright(ii,j,k,0) = Ez(ii,j,k);
}
}
}
/* Ez below substrate */
i=LIMX-1;
for (j=0; j<LIMY-1; j++) {
for (k=0; k<HEIGHT; k++) {
Ez(i,j,k) = c11*(Ez(i-2,j,k)+Ezright(i,j,k,1))
+ c21*(Ezright(i,j,k,0) + Ezright(i-2,j,k,0) - Ez(i-1,j,k)
- Ezright(i-1,j,k,1)) + c31*Ezright(i-1,j,k,0) - Ezright(i-2,j,k,1);
for(ii=LIMX-3; ii<LIMX; ii++) {
Ezright(ii,j,k,1) = Ezright(ii,j,k,0);
Ezright(ii,j,k,0) = Ez(ii,j,k);
}
}
}
#endif
/************ Hx update. ************/
for (i=0; i<LIMX; i++)
for (j=0; j<LIMY-1; j++)
for (k=0; k<LIMZ-1; k++)
Hx(i,j,k) = Hx(i,j,k)
+ coefH*((Ey(i,j,k+1)-Ey(i,j,k)) - (Ez(i,j+1,k)-Ez(i,j,k)));
/************ Hy update. ************/
for (i=0; i<LIMX-1; i++)
for (j=0; j<LIMY; j++)
for (k=0; k<LIMZ-1; k++)
Hy(i,j,k) = Hy(i,j,k)
+ coefH*((Ez(i+1,j,k)-Ez(i,j,k)) - (Ex(i,j,k+1)-Ex(i,j,k)));
/************ Hz update. ************/
for (i=0; i<LIMX-1; i++)
for (j=0; j<LIMY-1; j++)
for (k=0; k<LIMZ; k++)
Hz(i,j,k) = Hz(i,j,k)
+ coefH*((Ex(i,j+1,k)-Ex(i,j,k)) - (Ey(i+1,j,k)-Ey(i,j,k)));
fprintf(obs,"%f\n",Ez((LINE_X_START+LINE_X_END)/2,PATCH_Y_START-10,2));
}
fclose(obs);
return 0;
}
/* ------------------------- end of main ------------------------- */
/* ########################## gaussian ########################### */
/* Gaussian pulse.
*/
double gaussian(int ntime, double cdtds, int ippw) {
double arg, time;
time = (double)(ntime);
arg = pow(M_PI*(cdtds*time/(double)(ippw)-1.),2);
/* If the argument is greater than 70.0, the result is so small we
may as well skip calculating the exponential and return zero. */
if (arg > 70.0) {
return 0.0;
} else {
return exp(-arg);
}
}
/* ------------------------ end of gaussian ------------------------ */
/* ########################### init ############################# */
void init() {
int i, j, k;
/* Initialize all fields. */
/* Ex. */
for (k=0;k<LIMZ;k++)
for (j=0;j<LIMY;j++)
for (i=0;i<LIMX;i++)
Ex(i,j,k) = 0.;
/* Ey. */
for (k=0;k<LIMZ;k++)
for (j=0;j<LIMY;j++)
for (i=0;i<LIMX;i++)
Ey(i,j,k) = 0.;
/* Ez. */
for (k=0;k<LIMZ;k++)
for (j=0;j<LIMY;j++)
for (i=0;i<LIMX;i++)
Ez(i,j,k) = 0.;
/* Hx. */
for (k=0;k<LIMZ;k++)
for (j=0;j<LIMY;j++)
for (i=0;i<LIMX;i++)
Hx(i,j,k) = 0.;
/* Hy. */
for (i=0;i<LIMX;i++)
for (j=0;j<LIMY; j++)
for (k=0; k<LIMZ; k++)
Hy(i,j,k) = 0.;
/* Hz. */
for (i=0; i<LIMX; i++)
for (j=0; j<LIMY; j++)
for (k=0; k<LIMZ; k++)
Hz(i,j,k) = 0.;
/* ABC arrays */
for (i=0; i<LIMX; i++)
for(j=LIMY-3; j<LIMY; j++)
for(k=0; k<LIMZ; k++) {
Exfar(i,j,k,0) = 0.0;
Exfar(i,j,k,1) = 0.0;
Ezfar(i,j,k,0) = 0.0;
Ezfar(i,j,k,1) = 0.0;
Exnear(i,j,k,0) = 0.0;
Exnear(i,j,k,1) = 0.0;
Eznear(i,j,k,0) = 0.0;
Eznear(i,j,k,1) = 0.0;
}
for (i=0; i<3; i++)
for(j=0; j<LIMY; j++)
for(k=0; k<LIMZ; k++) {
Eyleft(i,j,k,0) = 0.0;
Eyleft(i,j,k,1) = 0.0;
Ezleft(i,j,k,0) = 0.0;
Ezleft(i,j,k,1) = 0.0;
}
for (i=0; i<LIMX; i++)
for(j=0; j<LIMY; j++)
for(k=0; k<3; k++) {
Extop(i,j,k,0) = 0.0;
Extop(i,j,k,1) = 0.0;
Eytop(i,j,k,0) = 0.0;
Eytop(i,j,k,1) = 0.0;
}
for (i=LIMX-3; i<LIMX; i++)
for(j=0; j<LIMY; j++)
for(k=0; k<LIMZ; k++) {
Eyright(i,j,k,0) = 0.0;
Eyright(i,j,k,1) = 0.0;
Ezright(i,j,k,0) = 0.0;
Ezright(i,j,k,1) = 0.0;
}
}
/* ------------------------- end of init ------------------------ */⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -