📄 fdtd3d.txt
字号:
/* Fd3d_4.2.c 3D FDTD with PML added */
# include <math.h>
# include <stdlib.h>
# include <stdio.h>
#define IE 60
#define JE 60
#define KE 60
#define ia 7
#define ja 7
#define ka 7
#define NFREQS 3
main ()
{
float dx[IE][JE][KE],dy[IE][JE][KE],dz[IE][JE][KE];
float ex[IE][JE][KE],ey[IE][JE][KE],ez[IE][JE][KE];
float hx[IE][JE][KE],hy[IE][JE][KE],hz[IE][JE][KE];
float ix[IE][JE][KE],iy[IE][JE][KE],iz[IE][JE][KE];
float gax[IE][JE][KE],gay[IE][JE][KE],gaz[IE][JE][KE];
float gbx[IE][JE][KE],gby[IE][JE][KE],gbz[IE][JE][KE];
int l,m,n,i,j,k,ic,jc,kc,nsteps,n_pml;
float ddx,dt,T,epsz,muz,pi,eaf,npml,Eav;
int ib,jb,kb;
float xn,xxn,xnum,xd,curl_e;
float t0,spread,pulse;
FILE *fp, *fopen();
float idxl[ia][JE][KE],idxh[ia][JE][KE];
float ihxl[ia][JE][KE],ihxh[ia][JE][KE];
float idyl[IE][ja][KE],idyh[IE][ja][KE];
float ihyl[IE][ja][KE],ihyh[IE][ja][KE];
float idzl[IE][JE][ka],idzh[IE][JE][ka];
float ihzl[IE][JE][ka],ihzh[IE][JE][ka];
int ixh, jyh, kzh;
float gi0[IE],fi0[IE];
float gj0[JE],fj0[JE];
float gk0[KE],fk0[KE];
float gi1[IE],gi2[IE],gi3[IE];
float gj1[JE],gj2[JE],gj3[JE];
float gk1[KE],gk2[KE],gk3[KE];
float fi1[IE],fi2[IE],fi3[IE];
float fj1[JE],fj2[JE],fj3[JE];
float fk1[KE],fk2[KE],fk3[KE];
float curl_h,curl_d;
ic = IE/2 ;
jc = JE/2 ;
kc = KE/2 ;
ib = IE - ia - 1;
jb = JE - ja - 1;
kb = KE - ka - 1;
pi = 3.14159;
epsz = 8.8e-12;
muz = 4*pi*1.e-7;
ddx = .01; /* Cell size */
dt = ddx/6e8; /* Time steps */
/* Initialize the arrays */
for ( j=0; j < JE; j++ ) {
for ( k=0; k < KE; k++ ) {
for ( i=0; i < IE; i++ ) {
ex[i][j][k]= 0.0 ;
ey[i][j][k]= 0.0 ;
ez[i][j][k]= 0.0 ;
dx[i][j][k]= 0.0 ;
dy[i][j][k]= 0.0 ;
dz[i][j][k]= 0.0 ;
hx[i][j][k]= 0.0 ;
hy[i][j][k]= 0.0 ;
hz[i][j][k]= 0.0 ;
ix[i][j][k]= 0.0 ;
iy[i][j][k]= 0.0 ;
iz[i][j][k]= 0.0 ;
gax[i][j][k]= 1. ;
gay[i][j][k]= 1. ;
gaz[i][j][k]= 1. ;
} } }
/* Specify the dipole */
for ( k=kc-20; k < kc+21; k++ ) {
gaz[ic][jc][k] = 0.;
}
gaz[ic][jc][kc] = 1.0;
/* Boundary Conditions */
for ( i=0; i < ia; i++ ) {
for ( j=0; j < JE; j++ ) {
for ( k=0; k < KE; k++ ) {
idxl[i][j][k] = 0.0;
idxh[i][j][k] = 0.0;
ihxl[i][j][k] = 0.0;
ihxh[i][j][k] = 0.0;
} } }
for ( i=0; i < IE; i++ ) {
for ( j=0; j < ja; j++ ) {
for ( k=0; k < KE; k++ ) {
idyl[i][j][k] = 0.0;
idyh[i][j][k] = 0.0;
ihyl[i][j][k] = 0.0;
ihyh[i][j][k] = 0.0;
} } }
for ( i=0; i < IE; i++ ) {
for ( j=0; j < JE; j++ ) {
for ( k=0; k < ka; k++ ) {
idzl[i][j][k] = 0.0;
idzh[i][j][k] = 0.0;
ihzl[i][j][k] = 0.0;
ihzh[i][j][k] = 0.0;
} } }
for ( i=0; i < IE; i++ ) {
gi0[i] = 0.;
fi0[i] = 0.;
gi1[i] = 0.;
fi1[i] = 0.;
gi2[i] = 1.;
fi2[i] = 1.;
gi3[i] = 1.;
fi3[i] = 1.;
}
for ( j=0; j < JE; j++ ) {
gj0[j] = 0.;
fj0[j] = 0.;
gj1[j] = 0.;
fj1[j] = 0.;
gj2[j] = 1.;
fj2[j] = 1.;
gj3[j] = 1.;
fj3[j] = 1.;
}
for ( k=0; k < IE; k++ ) {
gk0[k] = 0.;
fk0[k] = 0.;
gk1[k] = 0.;
fk1[k] = 0.;
gk2[k] = 1.;
fk2[k] = 1.;
gk3[k] = 1.;
fk3[k] = 1.;
}
printf( "npml --> ");
scanf("%f", &npml);
printf("%f \n", npml);
n_pml = npml;
for ( i=0; i < n_pml; i++ ) {
xxn = (npml-i)/npml;
xn = .33*pow(xxn,3.);
printf( "%d xn = %8.4f xn = %8.4f \n",
i,xxn,xn);
fi0[i] = 1.;
fi0[IE-i-1] = 1.;
fi1[i] = xn;
fi1[IE-i-1] = xn;
gi2[i] = 1./(1.+xn);
gi2[IE-i-1] = 1./(1.+xn);
gi3[i] = (1.-xn)/(1.+xn);
gi3[IE-i-1] = (1.-xn)/(1.+xn);
xxn = (npml-i-.5)/npml;
xn = .33*pow(xxn,3.);
gi0[i] = 1.;
gi0[IE-i-2] = 1.;
gi1[i] = xn;
gi1[IE-i-2] = xn;
fi2[i] = 1./(1.+xn);
fi2[IE-i-2] = 1./(1.+xn);
fi3[i] = (1.-xn)/(1.+xn);
fi3[IE-i-2] = (1.-xn)/(1.+xn);
}
printf( "f \n");
for ( i=0; i < IE; i++ ) {
printf( "%2d %6.4f %6.4f %6.4f %6.4f\n",
i,fi0[i],fi1[i],gi2[i],gi3[i]);
printf( " %6.4f %6.4f %6.4f %6.4f\n",
gi0[i],gi1[i],fi2[i],fi3[i]);
}
for ( j=0; j < n_pml; j++ ) {
xxn = (npml-j)/npml;
xn = .33*pow(xxn,3.);
fj0[j] = 1.;
fj0[JE-j-1] = 1.;
fj1[j] = xn;
fj1[JE-j-1] = xn;
gj2[j] = 1./(1.+xn);
gj2[JE-j-1] = 1./(1.+xn);
gj3[j] = (1.-xn)/(1.+xn);
gj3[JE-j-1] = (1.-xn)/(1.+xn);
xxn = (npml-j-.5)/npml;
xn = .33*pow(xxn,3.);
gj0[j] = 1.;
gj0[JE-j-2] = 1.;
gj1[j] = xn;
gj1[JE-j-2] = xn;
fj2[j] = 1./(1.+xn);
fj2[JE-j-2] = 1./(1.+xn);
fj3[j] = (1.-xn)/(1.+xn);
fj3[JE-j-2] = (1.-xn)/(1.+xn);
}
printf( "fj & gj \n");
for ( j=0; j < JE; j++ ) {
printf( "%2d %6.4f %6.4f %6.4f %6.4f\n",
j,fj0[j],fj1[j],gj2[j],gj3[j]);
printf( " %6.4f %6.4f %6.4f %6.4f\n",
gj0[j],gj1[j],fj2[j],fj3[j]);
}
for ( k=0; k < n_pml; k++ ) {
xxn = (npml-k)/npml;
xn = .33*pow(xxn,3.);
fk0[k] = 1.;
fk0[KE-k-1] = 1.;
fk1[k] = xn;
fk1[KE-k-1] = xn;
gk2[k] = 1./(1.+xn);
gk2[KE-k-1] = 1./(1.+xn);
gk3[k] = (1.-xn)/(1.+xn);
gk3[KE-k-1] = (1.-xn)/(1.+xn);
xxn = (npml-k-.5)/npml;
xn = .33*pow(xxn,3.);
gk0[k] = 1.;
gk0[KE-k-2] = 1.;
gk1[k] = xn;
gk1[KE-k-2] = xn;
fk2[k] = 1./(1.+xn);
fk2[KE-k-2] = 1./(1.+xn);
fk3[k] = (1.-xn)/(1.+xn);
fk3[KE-k-2] = (1.-xn)/(1.+xn);
}
printf( "fk & gk \n");
for ( k=0; k < KE; k++ ) {
printf( "%2d %6.4f %6.4f %6.4f %6.4f\n",
k,fk0[k],fk1[k],gk2[k],gk3[k]);
printf( " %6.4f %6.4f %6.4f %6.4f\n",
gk0[k],gk1[k],fk2[k],fk3[k]);
}
/* Write the gaz out to a file "gaz" */
fp = fopen( "gaz","w");
for ( k=0; k < KE; k++ ) {
for ( i=0; i < IE; i++ ) {
fprintf( fp,"%7.4f ",gaz[i][jc][k]);
}
fprintf( fp," \n");
}
fclose(fp);
t0 = 40.0;
spread = 10.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 Dx field */
for ( i=1; i < ia; i++ ) {
for ( j=1; j < JE; j++ ) {
for ( k=1; k < KE; k++ ) {
curl_h = ( hz[i][j][k] - hz[i][j-1][k]
- hy[i][j][k] + hy[i][j][k-1]) ;
idxl[i][j][k] = idxl[i][j][k] + gi0[i]*curl_h;
dx[i][j][k] = gj3[j]*gk3[k]*dx[i][j][k]
+ gj2[j]*gk2[k]*.5*(curl_h + gi1[i]*idxl[i][j][k]);
} } }
for ( i=ia; i <= ib; i++ ) {
for ( j=1; j < JE; j++ ) {
for ( k=1; k < KE; k++ ) {
curl_h = ( hz[i][j][k] - hz[i][j-1][k]
- hy[i][j][k] + hy[i][j][k-1]) ;
dx[i][j][k] = gj3[j]*gk3[k]*dx[i][j][k]
+ gj2[j]*gk2[k]*.5*curl_h ;
} } }
for ( i=ib+1; i < IE; i++ ) {
ixh = i - ib - 1;
for ( j=1; j < JE; j++ ) {
for ( k=1; k < KE; k++ ) {
curl_h = ( hz[i][j][k] - hz[i][j-1][k]
- hy[i][j][k] + hy[i][j][k-1]) ;
idxh[ixh][j][k] = idxh[ixh][j][k] + gi0[i]*curl_h;
dx[i][j][k] = gj3[j]*gk3[k]*dx[i][j][k]
+ gj2[j]*gk2[k]*.5*(curl_h + gi1[i]*idxh[ixh][j][k]);
} } }
/* Calculate the Dy field */
for ( i=1; i < IE; i++ ) {
for ( j=1; j < ja; j++ ) {
for ( k=1; k < KE; k++ ) {
curl_h = ( hx[i][j][k] - hx[i][j][k-1]
- hz[i][j][k] + hz[i-1][j][k]) ;
idyl[i][j][k] = idyl[i][j][k] + gj0[j]*curl_h;
dy[i][j][k] = gi3[i]*gk3[k]*dy[i][j][k]
+ gi2[i]*gk2[k]*.5*( curl_h + gj1[j]*idyl[i][j][k]);
} } }
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -