📄 cpmlfdtd3d.c
字号:
//....................................................
// PML for top Hx, j-direction
//.....................................................
jj = nyPML_2 - 2;
for(j = Jmax - nyPML_2; j < Jmax-1; ++j) {
psi_Hxy_2[i][jj][k] = bh_y_2[jj] * psi_Hxy_2[i][jj][k]
+ ch_y_2[jj] * (Ez[i][j][k] - Ez[i][j+1][k]) / dy;
Hx[i][j][k] = Hx[i][j][k] + DB * psi_Hxy_2[i][jj][k];
jj = jj - 1;
}
}
}
for(i = 0; i < Imax-1; ++i) {
for(j = 0; j < Jmax-1; ++j) {
//....................................................
// PML for bottom Hx, k-direction
//................................................
for(k = 1; k < nzPML_1; ++k) {
psi_Hxz_1[i][j][k-1] = bh_z_1[k-1] * psi_Hxz_1[i][j][k-1]
+ ch_z_1[k-1] * (Ey[i][j][k] - Ey[i][j][k-1]) / dz;
Hx[i][j][k] = Hx[i][j][k] + DB * psi_Hxz_1[i][j][k-1];
}
//....................................................
// PML for top Hx, k-direction
//...............................................
kk = nzPML_2 - 2;
for(k = Kmax - nzPML_2; k < Kmax-1; ++k) {
psi_Hxz_2[i][j][kk] = bh_z_2[kk] * psi_Hxz_2[i][j][kk]
+ ch_z_2[kk] * (Ey[i][j][k] - Ey[i][j][k-1]) / dz;
Hx[i][j][k] = Hx[i][j][k] + DB * psi_Hxz_2[i][j][kk];
kk = kk - 1;
}
}
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// UPDATE Hy
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for(k = 1; k < Kmax-1; ++k) {
for(i = 0; i < Imax-1; ++i) {
for(j = 0; j < Jmax-1; ++j) {
Hy[i][j][k] = DA * Hy[i][j][k] + DB *
((Ez[i+1][j][k] - Ez[i][j][k]) * den_hx[i] +
(Ex[i][j][k-1] - Ex[i][j][k]) * den_hz[k] );
}
}
for(j = 0; j < Jmax-1; ++j) {
//.......................................................
// PML for bottom Hy, i-direction
//.......................................................
for(i = 0; i < nxPML_1-1; ++i){
psi_Hyx_1[i][j][k] = bh_x_1[i] * psi_Hyx_1[i][j][k]
+ ch_x_1[i] * (Ez[i+1][j][k] - Ez[i][j][k]) / dx;
Hy[i][j][k] = Hy[i][j][k] + DB * psi_Hyx_1[i][j][k];
}
//.........................................................
// PML for top Hy, i-direction
//.........................................................
ii = nxPML_2 - 2;
for(i = Imax - nxPML_2; i < Imax-1; ++i) {
psi_Hyx_2[ii][j][k] = bh_x_2[ii] * psi_Hyx_2[ii][j][k]
+ ch_x_2[ii] * (Ez[i+1][j][k] - Ez[i][j][k]) / dx;
Hy[i][j][k] = Hy[i][j][k] + DB * psi_Hyx_2[ii][j][k];
ii = ii - 1;
}
}
}
for(i = 0; i < Imax-1; ++i) {
for(j = 0; j < Jmax-1; ++j) {
//.......................................................
// PML for bottom Hy, k-direction
//......................................................
for(k = 1; k < nzPML_1; ++k) {
psi_Hyz_1[i][j][k-1] = bh_z_1[k-1] * psi_Hyz_1[i][j][k-1]
+ ch_z_1[k-1] * (Ex[i][j][k-1] - Ex[i][j][k]) / dz;
Hy[i][j][k] = Hy[i][j][k] + DB * psi_Hyz_1[i][j][k-1];
}
//.......................................................
// PML for top Hy, k-direction
//.........................................................
kk = nzPML_2 - 2;
for(k = Kmax - nzPML_2; k < Kmax-1; ++k) {
psi_Hyz_2[i][j][kk] = bh_z_2[kk] * psi_Hyz_2[i][j][kk]
+ ch_z_2[kk] * (Ex[i][j][k-1] - Ex[i][j][k]) / dz;
Hy[i][j][k] = Hy[i][j][k] + DB * psi_Hyz_2[i][j][kk];
kk = kk - 1;
}
}
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// UPDATE Hz
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for(k = 0; k < Kmax-1; ++k) {
for(i = 0; i < Imax-1; ++i) {
for(j = 0; j< Jmax-1; ++j) {
Hz[i][j][k] = DA * Hz[i][j][k] + DB
* ((Ey[i][j][k] - Ey[i+1][j][k]) * den_hx[i] +
(Ex[i][j+1][k] - Ex[i][j][k]) * den_hy[j]);
}
}
for(j = 0; j < Jmax-1; ++j) {
//..........................................................
// PML for bottom Hz, x-direction
//..........................................................
for(i = 0; i < nxPML_1-1; ++i) {
psi_Hzx_1[i][j][k] = bh_x_1[i] * psi_Hzx_1[i][j][k]
+ ch_x_1[i] * (Ey[i][j][k] - Ey[i+1][j][k]) / dx;
Hz[i][j][k] = Hz[i][j][k] + DB * psi_Hzx_1[i][j][k];
}
//..........................................................
// PML for top Hz, x-direction
//..........................................................
ii = nxPML_2 - 2;
for(i = Imax - nxPML_2; i < Imax-1; ++i) {
psi_Hzx_2[ii][j][k] = bh_x_2[ii] * psi_Hzx_2[ii][j][k]
+ ch_x_2[ii] * (Ey[i][j][k] - Ey[i+1][j][k])/dx;
Hz[i][j][k] = Hz[i][j][k] + DB * psi_Hzx_2[ii][j][k];
ii = ii - 1;
}
}
for(i = 0; i < Imax-1; ++i) {
//........................................................
// PML for bottom Hz, y-direction
//.........................................................
for(j = 0; j < nyPML_1-1; ++j) {
psi_Hzy_1[i][j][k] = bh_y_1[j] * psi_Hzy_1[i][j][k]
+ ch_y_1[j] * (Ex[i][j+1][k] - Ex[i][j][k]) / dy;
Hz[i][j][k] = Hz[i][j][k] + DB* psi_Hzy_1[i][j][k];
}
//.........................................................
// PML for top Hz, y-direction
//..........................................................
jj = nyPML_2 - 2;
for(j = Jmax - nyPML_2; j < Jmax-1; ++j) {
psi_Hzy_2[i][jj][k] = bh_y_2[jj] * psi_Hzy_2[i][jj][k]
+ ch_y_2[jj] * (Ex[i][j+1][k] - Ex[i][j][k]) / dy;
Hz[i][j][k] = Hz[i][j][k] + DB * psi_Hzy_2[i][jj][k];
jj = jj - 1;
}
}
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// UPDATE Ex
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for(k = 0; k < Kmax-1; ++k) {
for(i = 0; i < Imax-1; ++i) {
for(j = 1; j < Jmax-1; ++j) {
id = ID1[i][j][k];
if(id == 1) { // PEC
Ex[i][j][k] = 0;
} else {
Ex[i][j][k] = CA[id] * Ex[i][j][k] + CB[id] *
((Hz[i][j][k] - Hz[i][j-1][k]) * den_ey[j] +
(Hy[i][j][k] - Hy[i][j][k+1]) * den_ez[k] );
}
}
}
for(i = 0; i < Imax-1; ++i) {
//..............................................................
// PML for bottom Ex, j-direction
//..............................................................
for(j = 1; j < nyPML_1; ++j) {
id = ID1[i][j][k];
psi_Exy_1[i][j][k] = be_y_1[j] * psi_Exy_1[i][j][k]
+ ce_y_1[j] * (Hz[i][j][k] - Hz[i][j-1][k])/dy;
Ex[i][j][k] = Ex[i][j][k] + CB[id] * psi_Exy_1[i][j][k];
}
//.............................................................
// PML for top Ex, j-direction
//.............................................................
jj = nyPML_2 - 1;
for(j = Jmax - nyPML_2; j < Jmax-1; ++j) {
id = ID1[i][j][k];
psi_Exy_2[i][jj][k] = be_y_2[jj] * psi_Exy_2[i][jj][k]
+ ce_y_2[jj] * (Hz[i][j][k] - Hz[i][j-1][k]) / dy;
Ex[i][j][k] = Ex[i][j][k] + CB[id] * psi_Exy_2[i][jj][k];
jj = jj - 1;
}
}
}
for(i = 0; i < Imax-1; ++i) {
for(j = 1; j < Jmax-1; ++j) {
//.............................................................
// PML for bottom Ex, k-direction
//.............................................................
for(k = 0; k < nzPML_1; ++k) {
id = ID1[i][j][k];
psi_Exz_1[i][j][k] = be_z_1[k] * psi_Exz_1[i][j][k]
+ ce_z_1[k] * (Hy[i][j][k] - Hy[i][j][k+1]) / dz;
Ex[i][j][k] = Ex[i][j][k] + CB[id] * psi_Exz_1[i][j][k];
}
//..............................................................
// PML for top Ex, k-direction
//..............................................................
kk = nzPML_2 - 1;
for(k = Kmax - nzPML_2 - 1; k < Kmax-1; ++k) {
id = ID1[i][j][k];
psi_Exz_2[i][j][kk] = be_z_2[kk] * psi_Exz_2[i][j][kk]
+ ce_z_2[kk] * (Hy[i][j][k] - Hy[i][j][k+1]) / dz;
Ex[i][j][k] = Ex[i][j][k] + CB[id] * psi_Exz_2[i][j][kk];
kk = kk - 1;
}
}
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// UPDATE Ey
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for(k = 0; k < Kmax-1; ++k) {
for(i = 1; i < Imax-1; ++i) {
for(j = 0; j < Jmax-1; ++j) {
id = ID2[i][j][k];
if(id == 1) { // PEC
Ey[i][j][k] = 0;
} else {
Ey[i][j][k] = CA[id] * Ey[i][j][k] + CB[id] *
((Hz[i-1][j][k] - Hz[i][j][k]) * den_ex[i] +
(Hx[i][j][k+1] - Hx[i][j][k]) * den_ez[k] );
}
}
}
for(j = 0; j < Jmax-1; ++j) {
//...........................................................
// PML for bottom Ey, i-direction
//...........................................................
for(i = 1; i < nxPML_1; ++i) {
id = ID2[i][j][k];
psi_Eyx_1[i][j][k] = be_x_1[i] * psi_Eyx_1[i][j][k]
+ ce_x_1[i] * (Hz[i-1][j][k] - Hz[i][j][k]) / dx;
Ey[i][j][k] = Ey[i][j][k] + CB[id] * psi_Eyx_1[i][j][k];
}
//............................................................
// PML for top Ey, i-direction
//............................................................
ii = nxPML_2 - 1;
for(i = Imax - nxPML_2; i < Imax-1; ++i) {
id = ID2[i][j][k];
psi_Eyx_2[ii][j][k] = be_x_2[ii] * psi_Eyx_2[ii][j][k]
+ ce_x_2[ii] * (Hz[i-1][j][k] - Hz[i][j][k]) / dx;
Ey[i][j][k] = Ey[i][j][k] + CB[id] * psi_Eyx_2[ii][j][k];
ii = ii - 1;
}
}
}
for(i = 1; i < Imax-1; ++i) {
for(j = 0; j < Jmax-1; ++j) {
//...........................................................
// PML for bottom Ey, k-direction
//...........................................................
for(k = 0; k < nzPML_1; ++k) {
id = ID2[i][j][k];
psi_Eyz_1[i][j][k] = be_z_1[k] * psi_Eyz_1[i][j][k]
+ ce_z_1[k] * (Hx[i][j][k+1] - Hx[i][j][k]) / dz;
Ey[i][j][k] = Ey[i][j][k] + CB[id] * psi_Eyz_1[i][j][k];
}
//...........................................................
// PML for top Ey, k-direction
//............................................................
kk = nzPML_2 - 1;
for(k = Kmax - nzPML_2 - 1; k < Kmax-1; ++k) {
id = ID2[i][j][k];
psi_Eyz_2[i][j][kk] = be_z_2[kk] * psi_Eyz_2[i][j][kk]
+ ce_z_2[kk] * (Hx[i][j][k+1] - Hx[i][j][k]) / dz;
Ey[i][j][k] = Ey[i][j][k] + CB[id] * psi_Eyz_2[i][j][kk];
kk = kk - 1;
}
}
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// UPDATE Ez
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
for(k = 1; k < Kmax-1; ++k) {
for(i = 1; i < Imax-1; ++i) {
for(j = 1; j < Jmax-1; ++j) {
id = ID3[i][j][k];
if(id == 1) { // PEC
Ez[i][j][k] = 0;
} else {
Ez[i][j][k] = CA[id] * Ez[i][j][k] + CB[id]
* ((Hy[i][j][k] - Hy[i-1][j][k]) * den_ex[i] +
(Hx[i][j-1][k] - Hx[i][j][k]) * den_ey[j]);
}
}
}
for(j = 1; j < Jmax-1; ++j) {
//............................................................
// PML for bottom Ez, x-direction
//.............................................................
for(i = 1; i < nxPML_1; ++i) {
id = ID3[i][j][k];
psi_Ezx_1[i][j][k] = be_x_1[i] * psi_Ezx_1[i][j][k]
+ ce_x_1[i] * (Hy[i][j][k] - Hy[i-1][j][k]) / dx;
Ez[i][j][k] = Ez[i][j][k] + CB[id] * psi_Ezx_1[i][j][k];
}
//............................................................
// PML for top Ez, x-direction
//............................................................
ii = nxPML_2 - 1;
for(i = Imax - nxPML_2; i < Imax-1; ++i) {
id = ID3[i][j][k];
psi_Ezx_2[ii][j][k] = be_x_2[ii] * psi_Ezx_2[ii][j][k]
+ ce_x_2[ii] * (Hy[i][j][k] - Hy[i-1][j][k]) / dx;
Ez[i][j][k] = Ez[i][j][k] + CB[id] * psi_Ezx_2[ii][j][k];
ii = ii - 1;
}
}
for(i = 1; i < Imax-1; ++i) {
//..........................................................
// PML for bottom Ez, y-direction
//..........................................................
for(j = 1; j < nyPML_1; ++j) {
id = ID3[i][j][k];
psi_Ezy_1[i][j][k] = be_y_1[j] * psi_Ezy_1[i][j][k]
+ ce_y_1[j] * (Hx[i][j-1][k] - Hx[i][j][k]) / dy;
Ez[i][j][k] = Ez[i][j][k] + CB[id] * psi_Ezy_1[i][j][k];
}
//............................................................
// PML for top Ez, y-direction
//............................................................
jj = nyPML_2 - 1;
for(j = Jmax - nyPML_2; j < Jmax-1; ++j) {
id = ID3[i][j][k];
psi_Ezy_2[i][jj][k] = be_y_2[jj] * psi_Ezy_2[i][jj][k]
+ ce_y_2[jj] * (Hx[i][j-1][k] - Hx[i][j][k]) / dy;
Ez[i][j][k] = Ez[i][j][k] + CB[id] * psi_Ezy_2[i][jj][k];
jj = jj - 1;
}
}
}
//-----------------------------------------------------------
// Apply a point source (Soft)
//-----------------------------------------------------------
i = 25;
j = 63;
k = 12;
source = amp * -2.0 * ((n * dt - tO) / tw)
* exp(-pow(((n * dt - tO) / tw), 2));//Differentiated Gaussian pulse
Ez[i][j][k] = Ez[i][j][k] - CB[ID3[i][j][k]] * source;
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// WRITE TO OUTPUT FILES
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if((n % save_modulus) == 0) {
writeField(n);
}
}
// END TIME STEP
printf("Done time-stepping...\n");
}
//Builds an object
void buildObject() {
//buildSphere();
buildDipole();
}
//Builds a sphere (Sample code - NOt used in this program)
void buildSphere() {
double dist;//distance
double rad = 8;//(double)Imax / 5.0; // sphere radius
double sc = (double)Imax / 2.0;//sphere centre
double rad2 = 0.3;//(double)Imax / 5.0 - 3.0; // sphere radius
for(i = 0; i < Imax; ++i) {
for(j = 0; j < Jmax; ++j) {
for(k = 0; k < Kmax; ++k) {
//compute distance form centre to the point i, j, k
dist = sqrt((i + 0.5 - sc) * (i + 0.5 - sc) +
(j + 0.5 - sc) * (j + 0.5 - sc) +
(k + 0.5 - sc) * (k + 0.5 - sc));
//if point is within the sphere
if (dist <= rad) {
//set the material at that point
yeeCube (i, j, k, 6);
}
}
}
}
}
//Builds a dipole
void buildDipole() {
int centre = (jstart + jend) / 2;
for(i = istart; i <= iend; ++i) {
for(j = jstart; j <= jend; ++j) {
for(k = kstart; k <= kend; ++k) {
if(j != centre) {
yeeCube (i, j, k, 1);//PEC material
}
}
}
}
}
//creates a dielctric cube (yee cell) made up of the selected material
void yeeCube (int I, int J,int K, short mType) {
//set face 1 (for EX)
ID1[I][J][K] = mType;
ID1[I][J][K + 1] = mType;
ID1[I][J + 1][K + 1] = mType;
ID1[I][J + 1][K] = mType;
//set face 2 (for EY)
ID2[I][J][K] = mType;
ID2[I + 1][J][K] = mType;
ID2[I + 1][J][K + 1] = mType;
ID2[I][J][K + 1] = mType;
//set face 3 (for EZ)
ID3[I][J][K] = mType;
ID3[I + 1][J][K] = mType;
ID3[I + 1][J + 1][K] = mType;
ID3[I][J + 1][K] = mType;
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// Saving Output Data to files
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
void writeField (int iteration) {
FILE *ptr;
char step[10];
char fileBaseName[] = "E_Field_";
sprintf(step, "%d", iteration);
strcat(fileBaseName, step);
strcat(fileBaseName, ".txt");
ptr = fopen(fileBaseName,"wt");
for(i = 0 ; i < Imax-1 ; i++) {
for(j = 0 ; j < Jmax-1 ; j++){
// |E|
fprintf(ptr, "%f\t", sqrt(pow(Ex[i][j][ksource], 2) +
pow(Ey[i][j][ksource], 2) + pow( Ez[i][j][ksource], 2)));
// fprintf(ptr, "%f\t", Ex[i][j][ksource]);//Ex
// fprintf(ptr, "%f\t", Ey[i][j][ksource]);//Ey
// fprintf(ptr, "%f\t", Ez[i][j][ksource]);//Ez
// fprintf(ptr, "%f\t", Hx[i][j][ksource]);//Hx
// fprintf(ptr, "%f\t", Hy[i][j][ksource]);//Hy
// fprintf(ptr, "%f\t", Hz[i][j][ksource]);//Hz
// |H|
// fprintf(ptr, "%f\t", sqrt(pow(Hx[i][j][ksource], 2) +
// pow(Hy[i][j][ksource], 2) + pow( Hz[i][j][ksource], 2)));
}
fprintf(ptr, "\n");
}
fclose(ptr);
}
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// END OF PROGRAM CPMLFDTD3D
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -