?? fdtd_3d_pml_hy.cpp
字號:
#include "fdtd_3d_pml_hy.h"
///////////////////////////////////////////////////////////////////////////////////////
//Calculate Hy field
///////////////////////////////////////////////////////////////////////////////////////
void calc_Hy_1(double ***Hy, double ***Ex, double ***Ez, double ***By_1, long ***Ind,
double *mu_r, double *K_By_a_1, double *K_By_b_1, double *K_Hy_a_1,
double *K_Hy_b_1, double *K_Hy_c_1, double *K_Hy_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 By_r;
//1
for (i = 0; i < nPML_x_1; i++)
{
for (j = 0; j < nPML_y_1; j++)
{
for (k = 0; k < nPML_z_1; k++)
{
By_r = By_1[i][j][k];
By_1[i][j][k] = K_By_a_1[k]*By_1[i][j][k] +
K_By_b_1[k]*( (Ez[i+1][j][k] - Ez[i][j][k])*inv_dx -
(Ex[i][j][k+1] - Ex[i][j][k])*inv_dz );
Hy[i][j][k] = K_Hy_a_1[i]*Hy[i][j][k] +
K_Hy_b_1[i]*( K_Hy_c_1[j]*By_1[i][j][k] +
K_Hy_d_1[j]*By_r )/mu_r[Ind[i][j][k]];
}
}
}
}
void calc_Hy_2(double ***Hy, double ***Ex, double ***Ez, double ***By_2, long ***Ind,
double *mu_r, double *K_By_a_1, double *K_By_b_1, double *K_Hy_c_1,
double *K_Hy_d_1, double inv_dx, double inv_dz, long nPML_x_1,
long nx_MIN_1_MIN_nPML_x_2, long nPML_y_1, long nPML_z_1, double inv_TwoOREp0)
{
long i, j, k, ia;
double By_r;
//2
ia = 0;
for (i = nPML_x_1; i < nx_MIN_1_MIN_nPML_x_2; i++)
{
for (j = 0; j < nPML_y_1; j++)
{
for (k = 0; k < nPML_z_1; k++)
{
By_r = By_2[ia][j][k];
By_2[ia][j][k] = K_By_a_1[k]*By_2[ia][j][k] +
K_By_b_1[k]*( (Ez[i+1][j][k] - Ez[i][j][k])*inv_dx -
(Ex[i][j][k+1] - Ex[i][j][k])*inv_dz );
Hy[i][j][k] += inv_TwoOREp0*( K_Hy_c_1[j]*By_2[ia][j][k] +
K_Hy_d_1[j]*By_r )/mu_r[Ind[i][j][k]];
}
}
ia++;
}
}
void calc_Hy_3(double ***Hy, double ***Ex, double ***Ez, double ***By_3, long ***Ind,
double *mu_r, double *K_By_a_1, double *K_By_b_1, double *K_Hy_a_2,
double *K_Hy_b_2, double *K_Hy_c_1, double *K_Hy_d_1, double inv_dx,
double inv_dz, long nx_MIN_1_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, ia, ii;
double By_r;
//3
ia = 0;
ii = nPML_x_2_MIN_1;
for (i = nx_MIN_1_MIN_nPML_x_2; i < nx_MIN_1; i++)
{
for (j = 0; j < nPML_y_1; j++)
{
for (k = 0; k < nPML_z_1; k++)
{
By_r = By_3[ia][j][k];
By_3[ia][j][k] = K_By_a_1[k]*By_3[ia][j][k] +
K_By_b_1[k]*( (Ez[i+1][j][k] - Ez[i][j][k])*inv_dx -
(Ex[i][j][k+1] - Ex[i][j][k])*inv_dz );
Hy[i][j][k] = K_Hy_a_2[ii]*Hy[i][j][k] +
K_Hy_b_2[ii]*( K_Hy_c_1[j]*By_3[ia][j][k] +
K_Hy_d_1[j]*By_r )/mu_r[Ind[i][j][k]];
}
}
ia++;
ii--;
}
}
void calc_Hy_4(double ***Hy, double ***Ex, double ***Ez, double ***By_4, long ***Ind,
double *mu_r, double *K_By_a_1, double *K_By_b_1, double *K_Hy_a_1,
double *K_Hy_b_1, double inv_dx, double inv_dz, long nPML_x_1, long nPML_y_1,
long ny_MIN_nPML_y_2, long nPML_z_1, double TwoOREp0DIVMu0)
{
long i, j, k, ja;
double By_r;
//4
for (i = 0; i < nPML_x_1; i++)
{
ja = 0;
for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
{
for (k = 0; k < nPML_z_1; k++)
{
By_r = By_4[i][ja][k];
By_4[i][ja][k] = K_By_a_1[k]*By_4[i][ja][k] +
K_By_b_1[k]*( (Ez[i+1][j][k] - Ez[i][j][k])*inv_dx -
(Ex[i][j][k+1] - Ex[i][j][k])*inv_dz );
Hy[i][j][k] = K_Hy_a_1[i]*Hy[i][j][k] +
K_Hy_b_1[i]*TwoOREp0DIVMu0*( By_4[i][ja][k] - By_r )/mu_r[Ind[i][j][k]];
}
ja++;
}
}
}
void calc_Hy_5(double ***Hy, double ***Ex, double ***Ez, long ***Ind, double *mu_r,
double *K_By_a_1, double *K_By_b_1, double inv_dxORmu_0, double inv_dzORmu_0,
long nPML_x_1, long nx_MIN_1_MIN_nPML_x_2, long nPML_y_1, long ny_MIN_nPML_y_2,
long nPML_z_1)
{
long i, j, k;
//5
for (i = nPML_x_1; i < nx_MIN_1_MIN_nPML_x_2; i++)
{
for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
{
for (k = 0; k < nPML_z_1; k++)
{
Hy[i][j][k] = K_By_a_1[k]*Hy[i][j][k] +
K_By_b_1[k]/mu_r[Ind[i][j][k]]*
( (Ez[i+1][j][k] - Ez[i][j][k])*inv_dxORmu_0 -
(Ex[i][j][k+1] - Ex[i][j][k])*inv_dzORmu_0 );
}
}
}
}
void calc_Hy_6(double ***Hy, double ***Ex, double ***Ez, double ***By_6, long ***Ind,
double *mu_r, double *K_By_a_1, double *K_By_b_1, double *K_Hy_a_2,
double *K_Hy_b_2, double inv_dx, double inv_dz, long nx_MIN_1_MIN_nPML_x_2,
long nx_MIN_1, long nPML_y_1, long ny_MIN_nPML_y_2, long nPML_z_1,
long nPML_x_2_MIN_1, double TwoOREp0DIVMu0)
{
long i, j, k, ia, ja, ii;
double By_r;
//6
ia = 0;
ii = nPML_x_2_MIN_1;
for (i = nx_MIN_1_MIN_nPML_x_2; i < nx_MIN_1; i++)
{
ja = 0;
for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
{
for (k = 0; k < nPML_z_1; k++)
{
By_r = By_6[ia][ja][k];
By_6[ia][ja][k] = K_By_a_1[k]*By_6[ia][ja][k] +
K_By_b_1[k]*( (Ez[i+1][j][k] - Ez[i][j][k])*inv_dx -
(Ex[i][j][k+1] - Ex[i][j][k])*inv_dz );
Hy[i][j][k] = K_Hy_a_2[ii]*Hy[i][j][k] +
K_Hy_b_2[ii]*TwoOREp0DIVMu0*( By_6[ia][ja][k] - By_r )/mu_r[Ind[i][j][k]];
}
ja++;
}
ia++;
ii--;
}
}
void calc_Hy_7(double ***Hy, double ***Ex, double ***Ez, double ***By_7, long ***Ind,
double *mu_r, double *K_By_a_1, double *K_By_b_1, double *K_Hy_a_1,
double *K_Hy_b_1, double *K_Hy_c_2, double *K_Hy_d_2, double inv_dx,
double inv_dz, long nPML_x_1, long ny_MIN_nPML_y_2, long ny, long nPML_z_1,
long nPML_y_2_MIN_1)
{
long i, j, k, ja, jj;
double By_r;
//7
for (i = 0; i < nPML_x_1; i++)
{
ja = 0;
jj = nPML_y_2_MIN_1;
for (j = ny_MIN_nPML_y_2; j < ny; j++)
{
for (k = 0; k < nPML_z_1; k++)
{
By_r = By_7[i][ja][k];
By_7[i][ja][k] = K_By_a_1[k]*By_7[i][ja][k] +
K_By_b_1[k]*( (Ez[i+1][j][k] - Ez[i][j][k])*inv_dx -
(Ex[i][j][k+1] - Ex[i][j][k])*inv_dz );
Hy[i][j][k] = K_Hy_a_1[i]*Hy[i][j][k] +
K_Hy_b_1[i]*( K_Hy_c_2[jj]*By_7[i][ja][k] +
K_Hy_d_2[jj]*By_r )/mu_r[Ind[i][j][k]];
}
ja++;
jj--;
}
}
}
void calc_Hy_8(double ***Hy, double ***Ex, double ***Ez, double ***By_8, long ***Ind,
double *mu_r, double *K_By_a_1, double *K_By_b_1, double *K_Hy_c_2,
double *K_Hy_d_2, double inv_dx, double inv_dz, long nPML_x_1,
long nx_MIN_1_MIN_nPML_x_2, long ny_MIN_nPML_y_2, long ny, long nPML_z_1,
long nPML_y_2_MIN_1, double inv_TwoOREp0)
{
long i, j, k, ia, ja, jj;
double By_r;
//8
ia = 0;
for (i = nPML_x_1; i < nx_MIN_1_MIN_nPML_x_2; i++)
{
ja = 0;
jj = nPML_y_2_MIN_1;
for (j = ny_MIN_nPML_y_2; j < ny; j++)
{
for (k = 0; k < nPML_z_1; k++)
{
By_r = By_8[ia][ja][k];
By_8[ia][ja][k] = K_By_a_1[k]*By_8[ia][ja][k] +
K_By_b_1[k]*( (Ez[i+1][j][k] - Ez[i][j][k])*inv_dx -
(Ex[i][j][k+1] - Ex[i][j][k])*inv_dz );
Hy[i][j][k] += inv_TwoOREp0*( K_Hy_c_2[jj]*By_8[ia][ja][k] +
K_Hy_d_2[jj]*By_r )/mu_r[Ind[i][j][k]];
}
ja++;
jj--;
}
ia++;
}
}
void calc_Hy_9(double ***Hy, double ***Ex, double ***Ez, double ***By_9, long ***Ind,
double *mu_r, double *K_By_a_1, double *K_By_b_1, double *K_Hy_a_2,
double *K_Hy_b_2, double *K_Hy_c_2, double *K_Hy_d_2, double inv_dx,
double inv_dz, long nx_MIN_1_MIN_nPML_x_2, long nx_MIN_1, long ny_MIN_nPML_y_2,
long ny, long nPML_z_1, long nPML_x_2_MIN_1, long nPML_y_2_MIN_1)
{
long i, j, k, ia, ja, ii, jj;
double By_r;
//9
ia = 0;
ii = nPML_x_2_MIN_1;
for (i = nx_MIN_1_MIN_nPML_x_2; i < nx_MIN_1; i++)
{
ja = 0;
jj = nPML_y_2_MIN_1;
for (j = ny_MIN_nPML_y_2; j < ny; j++)
{
for (k = 0; k < nPML_z_1; k++)
{
By_r = By_9[ia][ja][k];
By_9[ia][ja][k] = K_By_a_1[k]*By_9[ia][ja][k] +
K_By_b_1[k]*( (Ez[i+1][j][k] - Ez[i][j][k])*inv_dx -
(Ex[i][j][k+1] - Ex[i][j][k])*inv_dz );
Hy[i][j][k] = K_Hy_a_2[ii]*Hy[i][j][k] +
K_Hy_b_2[ii]*( K_Hy_c_2[jj]*By_9[ia][ja][k] +
K_Hy_d_2[jj]*By_r )/mu_r[Ind[i][j][k]];
}
ja++;
jj--;
}
ia++;
ii--;
}
}
void calc_Hy_10(double ***Hy, double ***Ex, double ***Ez, double ***By_10, long ***Ind,
double *mu_r, double *K_Hy_a_1, double *K_Hy_b_1, double *K_Hy_c_1,
double *K_Hy_d_1, double dtDIVdx, double dtDIVdz, long nPML_x_1, long nPML_y_1,
long nPML_z_1, long nz_MIN_1_MIN_nPML_z_2)
{
long i, j, k, ka;
double By_r;
//10
for (i = 0; i < nPML_x_1; i++)
{
for (j = 0; j < nPML_y_1; j++)
{
ka = 0;
for (k = nPML_z_1; k < nz_MIN_1_MIN_nPML_z_2; k++)
{
By_r = By_10[i][j][ka];
By_10[i][j][ka] += ( (Ez[i+1][j][k] - Ez[i][j][k])*dtDIVdx -
(Ex[i][j][k+1] - Ex[i][j][k])*dtDIVdz );
Hy[i][j][k] = K_Hy_a_1[i]*Hy[i][j][k] +
K_Hy_b_1[i]*( K_Hy_c_1[j]*By_10[i][j][ka] +
K_Hy_d_1[j]*By_r )/mu_r[Ind[i][j][k]];
ka++;
}
}
}
}
void calc_Hy_11(double ***Hy, double ***Ex, double ***Ez, double ***By_11, long ***Ind,
double *mu_r, double *K_Hy_c_1, double *K_Hy_d_1, double dtDIVdx,
double dtDIVdz, long nPML_x_1, long nx_MIN_1_MIN_nPML_x_2, long nPML_y_1,
long nPML_z_1, long nz_MIN_1_MIN_nPML_z_2, double inv_TwoOREp0)
{
long i, j, k, ia, ka;
double By_r;
//11
ia = 0;
for (i = nPML_x_1; i < nx_MIN_1_MIN_nPML_x_2; i++)
{
for (j = 0; j < nPML_y_1; j++)
{
ka = 0;
for (k = nPML_z_1; k < nz_MIN_1_MIN_nPML_z_2; k++)
{
By_r = By_11[ia][j][ka];
By_11[ia][j][ka] += ( (Ez[i+1][j][k] - Ez[i][j][k])*dtDIVdx -
(Ex[i][j][k+1] - Ex[i][j][k])*dtDIVdz );
Hy[i][j][k] += inv_TwoOREp0*( K_Hy_c_1[j]*By_11[ia][j][ka] +
K_Hy_d_1[j]*By_r )/mu_r[Ind[i][j][k]];
ka++;
}
}
ia++;
}
}
void calc_Hy_12(double ***Hy, double ***Ex, double ***Ez, double ***By_12, long ***Ind,
double *mu_r, double *K_Hy_a_2, double *K_Hy_b_2, double *K_Hy_c_1,
double *K_Hy_d_1, double dtDIVdx, double dtDIVdz, long nx_MIN_1_MIN_nPML_x_2,
long nx_MIN_1, long nPML_y_1, long nPML_z_1, long nz_MIN_1_MIN_nPML_z_2,
long nPML_x_2_MIN_1)
{
long i, j, k, ia, ka, ii;
double By_r;
//12
ia = 0;
ii = nPML_x_2_MIN_1;
for (i = nx_MIN_1_MIN_nPML_x_2; i < nx_MIN_1; i++)
{
for (j = 0; j < nPML_y_1; j++)
{
ka = 0;
for (k = nPML_z_1; k < nz_MIN_1_MIN_nPML_z_2; k++)
{
By_r = By_12[ia][j][ka];
By_12[ia][j][ka] += ( (Ez[i+1][j][k] - Ez[i][j][k])*dtDIVdx -
(Ex[i][j][k+1] - Ex[i][j][k])*dtDIVdz );
Hy[i][j][k] = K_Hy_a_2[ii]*Hy[i][j][k] +
K_Hy_b_2[ii]*( K_Hy_c_1[j]*By_12[ia][j][ka] +
K_Hy_d_1[j]*By_r )/mu_r[Ind[i][j][k]];
ka++;
}
}
ia++;
ii--;
}
}
void calc_Hy_13(double ***Hy, double ***Ex, double ***Ez, long ***Ind, double *mu_r,
double *K_Hy_a_1, double *K_Hy_b_1, double TwoOREp0DIVMu0ORdtDIVdx,
double TwoOREp0DIVMu0ORdtDIVdz, long nPML_x_1, long nPML_y_1,
long ny_MIN_nPML_y_2, long nPML_z_1, long nz_MIN_1_MIN_nPML_z_2)
{
long i, j, k;
//13
for (i = 0; i < nPML_x_1; i++)
{
for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
{
for (k = nPML_z_1; k < nz_MIN_1_MIN_nPML_z_2; k++)
{
Hy[i][j][k] = K_Hy_a_1[i]*Hy[i][j][k] + K_Hy_b_1[i]/mu_r[Ind[i][j][k]]*
( (Ez[i+1][j][k] - Ez[i][j][k])*TwoOREp0DIVMu0ORdtDIVdx -
(Ex[i][j][k+1] - Ex[i][j][k])*TwoOREp0DIVMu0ORdtDIVdz );
}
}
}
}
void calc_Hy_14(double ***Hy, double ***Ex, double ***Ez, long ***Ind, double *mu_r,
double dtDIVMu0DIVdx, double dtDIVMu0DIVdz, long nx_a,
long nx_b, long ny_a, long ny_b, long nz_a, long nz_b)
{
long i, j, k;
//14 Outside of the PML region
for (i = nx_a; i < nx_b; i++)
{
for (j = ny_a; j < ny_b; j++)
{
for (k = nz_a; k < nz_b; k++)
{
Hy[i][j][k] += ( (Ez[i+1][j][k] - Ez[i][j][k])*dtDIVMu0DIVdx -
(Ex[i][j][k+1] - Ex[i][j][k])*dtDIVMu0DIVdz )/mu_r[Ind[i][j][k]];
}
}
}
}
void calc_Hy_15(double ***Hy, double ***Ex, double ***Ez, long ***Ind, double *mu_r,
double *K_Hy_a_2, double *K_Hy_b_2, double TwoOREp0DIVMu0ORdtDIVdx,
double TwoOREp0DIVMu0ORdtDIVdz, long nx_MIN_1_MIN_nPML_x_2, long nx_MIN_1,
long nPML_y_1, long ny_MIN_nPML_y_2, long nPML_z_1, long nz_MIN_1_MIN_nPML_z_2,
long nPML_x_2_MIN_1)
{
long i, j, k, ii;
//15
ii = nPML_x_2_MIN_1;
for (i = nx_MIN_1_MIN_nPML_x_2; i < nx_MIN_1; i++)
{
for (j = nPML_y_1; j < ny_MIN_nPML_y_2; j++)
{
for (k = nPML_z_1; k < nz_MIN_1_MIN_nPML_z_2; k++)
?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -