📄 fdtd_3d_pml_ey.cpp
字号:
#include "fdtd_3d_pml_ey.h"
///////////////////////////////////////////////////////////////////////////////////////
//Calculate Ey field
///////////////////////////////////////////////////////////////////////////////////////
void calc_Ey_1(double ***Ey, double ***Hx, double ***Hz, double ***Fy_1, double ***Gy_1,
long ***Ind, double *K_a, double *K_b, double *K_Gy_a_1, double *K_Gy_b_1,
double *K_Ey_a_1, double *K_Ey_b_1, double *K_Ey_c_1, double *K_Ey_d_1,
double inv_dx, double inv_dz, long nPML_x_1, long nPML_y_1, long nPML_z_1)
{
long i, j, k;
double Fy_r, Gy_r;
//1
for (i = 1; i < nPML_x_1; i++)
{
for (j = 0; j < nPML_y_1; j++)
{
for (k = 1; k < nPML_z_1; k++)
{
Fy_r = Fy_1[i][j][k];
Fy_1[i][j][k] = K_a[Ind[i][j][k]]*Fy_1[i][j][k] +
K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz -
(Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );
Gy_r = Gy_1[i][j][k];
Gy_1[i][j][k] = K_Gy_a_1[k]*Gy_1[i][j][k] +
K_Gy_b_1[k]*(Fy_1[i][j][k] - Fy_r);
Ey[i][j][k] = K_Ey_a_1[i]*Ey[i][j][k] +
K_Ey_b_1[i]*( K_Ey_c_1[j]*Gy_1[i][j][k] + K_Ey_d_1[j]*Gy_r);
}
}
}
}
void calc_Ey_2(double ***Ey, double ***Hx, double ***Hz, double ***Fy_2, double ***Gy_2,
long ***Ind, double *K_a, double *K_b, double *K_Gy_a_1, double *K_Gy_b_1,
double *K_Ey_c_1, double *K_Ey_d_1, double inv_dx, double inv_dz,
long nPML_x_1, long nx_MIN_nPML_x_2, long nPML_y_1, long nPML_z_1,
double inv_TwoOREp0)
{
long i, j, k, ia;
double Fy_r, Gy_r;
//2
ia = 0;
for (i = nPML_x_1; i < nx_MIN_nPML_x_2; i++)
{
for (j = 0; j < nPML_y_1; j++)
{
for (k = 1; k < nPML_z_1; k++)
{
Fy_r = Fy_2[ia][j][k];
Fy_2[ia][j][k] = K_a[Ind[i][j][k]]*Fy_2[ia][j][k] +
K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz -
(Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );
Gy_r = Gy_2[ia][j][k];
Gy_2[ia][j][k] = K_Gy_a_1[k]*Gy_2[ia][j][k] +
K_Gy_b_1[k]*(Fy_2[ia][j][k] - Fy_r);
Ey[i][j][k] += inv_TwoOREp0*( K_Ey_c_1[j]*Gy_2[ia][j][k] +
K_Ey_d_1[j]*Gy_r);
}
}
ia++;
}
}
void calc_Ey_3(double ***Ey, double ***Hx, double ***Hz, double ***Fy_3, double ***Gy_3,
long ***Ind, double *K_a, double *K_b, double *K_Gy_a_1, double *K_Gy_b_1,
double *K_Ey_a_2, double *K_Ey_b_2, double *K_Ey_c_1, double *K_Ey_d_1,
double inv_dx, double inv_dz, long nx_MIN_nPML_x_2, long nx_MIN_1,
long nPML_y_1, long nPML_z_1, long nPML_x_2_MIN_1)
{
long i, j, k, ii, ia;
double Fy_r, Gy_r;
//3
ia = 0;
ii = nPML_x_2_MIN_1;
for (i = nx_MIN_nPML_x_2; i < nx_MIN_1; i++)
{
for (j = 0; j < nPML_y_1; j++)
{
for (k = 1; k < nPML_z_1; k++)
{
Fy_r = Fy_3[ia][j][k];
Fy_3[ia][j][k] = K_a[Ind[i][j][k]]*Fy_3[ia][j][k] +
K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz -
(Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );
Gy_r = Gy_3[ia][j][k];
Gy_3[ia][j][k] = K_Gy_a_1[k]*Gy_3[ia][j][k] +
K_Gy_b_1[k]*(Fy_3[ia][j][k] - Fy_r);
Ey[i][j][k] = K_Ey_a_2[ii]*Ey[i][j][k] +
K_Ey_b_2[ii]*( K_Ey_c_1[j]*Gy_3[ia][j][k] + K_Ey_d_1[j]*Gy_r);
}
}
ia++;
ii--;
}
}
void calc_Ey_4(double ***Ey, double ***Hx, double ***Hz, double ***Fy_4, double ***Gy_4,
long ***Ind, double *K_a, double *K_b, double *K_Gy_a_1, double *K_Gy_b_1,
double *K_Ey_a_1, double *K_Ey_b_1, double inv_dx, double inv_dz,
long nPML_x_1, long nPML_y_1, long ny_MIN_1_MIN_nPML_y_2, long nPML_z_1,
double TwoOREp0)
{
long i, j, k, ja;
double Fy_r, Gy_r;
//4
for (i = 1; i < nPML_x_1; i++)
{
ja = 0;
for (j = nPML_y_1; j < ny_MIN_1_MIN_nPML_y_2; j++)
{
for (k = 1; k < nPML_z_1; k++)
{
Fy_r = Fy_4[i][ja][k];
Fy_4[i][ja][k] = K_a[Ind[i][j][k]]*Fy_4[i][ja][k] +
K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz -
(Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );
Gy_r = Gy_4[i][ja][k];
Gy_4[i][ja][k] = K_Gy_a_1[k]*Gy_4[i][ja][k] +
K_Gy_b_1[k]*(Fy_4[i][ja][k] - Fy_r);
Ey[i][j][k] = K_Ey_a_1[i]*Ey[i][j][k] +
K_Ey_b_1[i]*TwoOREp0*( Gy_4[i][ja][k] - Gy_r);
}
ja++;
}
}
}
void calc_Ey_5(double ***Ey, double ***Hx, double ***Hz, double ***Fy_5, long ***Ind,
double *K_a, double *K_b, double *K_Gy_a_1, double *K_Gy_b_1,
double inv_dx, double inv_dz, long nPML_x_1, long nx_MIN_nPML_x_2,
long nPML_y_1, long ny_MIN_1_MIN_nPML_y_2, long nPML_z_1)
{
long i, j, k, ia, ja;
double Fy_r;
//5
ia = 0;
for (i = nPML_x_1; i < nx_MIN_nPML_x_2; i++)
{
ja = 0;
for (j = nPML_y_1; j < ny_MIN_1_MIN_nPML_y_2; j++)
{
for (k = 1; k < nPML_z_1; k++)
{
Fy_r = Fy_5[ia][ja][k];
Fy_5[ia][ja][k] = K_a[Ind[i][j][k]]*Fy_5[ia][ja][k] +
K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz -
(Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );
Ey[i][j][k] = K_Gy_a_1[k]*Ey[i][j][k] +
K_Gy_b_1[k]*(Fy_5[ia][ja][k] - Fy_r);
}
ja++;
}
ia++;
}
}
void calc_Ey_6(double ***Ey, double ***Hx, double ***Hz, double ***Fy_6, double ***Gy_6,
long ***Ind, double *K_a, double *K_b, double *K_Gy_a_1, double *K_Gy_b_1,
double *K_Ey_a_2, double *K_Ey_b_2, double inv_dx, double inv_dz,
long nx_MIN_nPML_x_2, long nx_MIN_1, long nPML_y_1,
long ny_MIN_1_MIN_nPML_y_2, long nPML_z_1, long nPML_x_2_MIN_1,
double TwoOREp0)
{
long i, j, k, ii, ia, ja;
double Fy_r, Gy_r;
//6
ia = 0;
ii = nPML_x_2_MIN_1;
for (i = nx_MIN_nPML_x_2; i < nx_MIN_1; i++)
{
ja = 0;
for (j = nPML_y_1; j < ny_MIN_1_MIN_nPML_y_2; j++)
{
for (k = 1; k < nPML_z_1; k++)
{
Fy_r = Fy_6[ia][ja][k];
Fy_6[ia][ja][k] = K_a[Ind[i][j][k]]*Fy_6[ia][ja][k] +
K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz -
(Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );
Gy_r = Gy_6[ia][ja][k];
Gy_6[ia][ja][k] = K_Gy_a_1[k]*Gy_6[ia][ja][k] +
K_Gy_b_1[k]*(Fy_6[ia][ja][k] - Fy_r);
Ey[i][j][k] = K_Ey_a_2[ii]*Ey[i][j][k] +
K_Ey_b_2[ii]*TwoOREp0*( Gy_6[ia][ja][k] - Gy_r);
}
ja++;
}
ia++;
ii--;
}
}
void calc_Ey_7(double ***Ey, double ***Hx, double ***Hz, double ***Fy_7, double ***Gy_7,
long ***Ind, double *K_a, double *K_b, double *K_Gy_a_1, double *K_Gy_b_1,
double *K_Ey_a_1, double *K_Ey_b_1, double *K_Ey_c_2, double *K_Ey_d_2,
double inv_dx, double inv_dz, long nPML_x_1, long ny_MIN_1_MIN_nPML_y_2,
long ny_MIN_1, long nPML_z_1, long nPML_y_2_MIN_1)
{
long i, j, k, jj, ja;
double Fy_r, Gy_r;
//7
for (i = 1; i < nPML_x_1; i++)
{
ja = 0;
jj = nPML_y_2_MIN_1;
for (j = ny_MIN_1_MIN_nPML_y_2; j < ny_MIN_1; j++)
{
for (k = 1; k < nPML_z_1; k++)
{
Fy_r = Fy_7[i][ja][k];
Fy_7[i][ja][k] = K_a[Ind[i][j][k]]*Fy_7[i][ja][k] +
K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz -
(Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );
Gy_r = Gy_7[i][ja][k];
Gy_7[i][ja][k] = K_Gy_a_1[k]*Gy_7[i][ja][k] +
K_Gy_b_1[k]*(Fy_7[i][ja][k] - Fy_r);
Ey[i][j][k] = K_Ey_a_1[i]*Ey[i][j][k] +
K_Ey_b_1[i]*( K_Ey_c_2[jj]*Gy_7[i][ja][k] + K_Ey_d_2[jj]*Gy_r);
}
ja++;
jj--;
}
}
}
void calc_Ey_8(double ***Ey, double ***Hx, double ***Hz, double ***Fy_8, double ***Gy_8,
long ***Ind, double *K_a, double *K_b, double *K_Gy_a_1, double *K_Gy_b_1,
double *K_Ey_c_2, double *K_Ey_d_2, double inv_dx, double inv_dz,
long nPML_x_1, long nx_MIN_nPML_x_2, long ny_MIN_1_MIN_nPML_y_2,
long ny_MIN_1, long nPML_z_1, long nPML_y_2_MIN_1, double inv_TwoOREp0)
{
long i, j, k, jj, ia, ja;
double Fy_r, Gy_r;
//8
ia = 0;
for (i = nPML_x_1; i < nx_MIN_nPML_x_2; i++)
{
ja = 0;
jj = nPML_y_2_MIN_1;
for (j = ny_MIN_1_MIN_nPML_y_2; j < ny_MIN_1; j++)
{
for (k = 1; k < nPML_z_1; k++)
{
Fy_r = Fy_8[ia][ja][k];
Fy_8[ia][ja][k] = K_a[Ind[i][j][k]]*Fy_8[ia][ja][k] +
K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz -
(Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );
Gy_r = Gy_8[ia][ja][k];
Gy_8[ia][ja][k] = K_Gy_a_1[k]*Gy_8[ia][ja][k] +
K_Gy_b_1[k]*(Fy_8[ia][ja][k] - Fy_r);
Ey[i][j][k] += inv_TwoOREp0*( K_Ey_c_2[jj]*Gy_8[ia][ja][k] +
K_Ey_d_2[jj]*Gy_r);
}
ja++;
jj--;
}
ia++;
}
}
void calc_Ey_9(double ***Ey, double ***Hx, double ***Hz, double ***Fy_9, double ***Gy_9,
long ***Ind, double *K_a, double *K_b, double *K_Gy_a_1, double *K_Gy_b_1,
double *K_Ey_a_2, double *K_Ey_b_2, double *K_Ey_c_2, double *K_Ey_d_2,
double inv_dx, double inv_dz, long nx_MIN_nPML_x_2, long nx_MIN_1,
long ny_MIN_1_MIN_nPML_y_2, long ny_MIN_1, long nPML_z_1, long nPML_x_2_MIN_1,
long nPML_y_2_MIN_1)
{
long i, j, k, ii, jj, ia, ja;
double Fy_r, Gy_r;
//9
ia = 0;
ii = nPML_x_2_MIN_1;
for (i = nx_MIN_nPML_x_2; i < nx_MIN_1; i++)
{
ja = 0;
jj = nPML_y_2_MIN_1;
for (j = ny_MIN_1_MIN_nPML_y_2; j < ny_MIN_1; j++)
{
for (k = 1; k < nPML_z_1; k++)
{
Fy_r = Fy_9[ia][ja][k];
Fy_9[ia][ja][k] = K_a[Ind[i][j][k]]*Fy_9[ia][ja][k] +
K_b[Ind[i][j][k]]*( (Hx[i][j][k] - Hx[i][j][k-1])*inv_dz -
(Hz[i][j][k] - Hz[i-1][j][k])*inv_dx );
Gy_r = Gy_9[ia][ja][k];
Gy_9[ia][ja][k] = K_Gy_a_1[k]*Gy_9[ia][ja][k] +
K_Gy_b_1[k]*(Fy_9[ia][ja][k] - Fy_r);
Ey[i][j][k] = K_Ey_a_2[ii]*Ey[i][j][k] +
K_Ey_b_2[ii]*( K_Ey_c_2[jj]*Gy_9[ia][ja][k] + K_Ey_d_2[jj]*Gy_r);
}
ja++;
jj--;
}
ia++;
ii--;
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -