📄 radarfdtd.cpp
字号:
for(i=0; i<pmlWidth; i++)
for(j=ny+pmlWidth; j<ny+2*pmlWidth; j++)
for(k=0; k<nz+2*pmlWidth; k++)
{
pmlSigmaYBottomYZ[i][j][k] = ReturnConductivity((j-pmlWidth-ny+1) * dy, pmlWidth * dy, maxSigma);
pmlSigmaYStarBottomYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomYZ[i][j][k]][EPSILON]) * pmlSigmaYBottomYZ[i][j][k];
pmlSYBottomYZ[i][j][k] = ReturnStretching((j-pmlWidth-ny+1) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
}
/* computing sigmaZ within the edges */
for(i=0; i<pmlWidth; i++)
for(j=0; j<ny+2*pmlWidth; j++)
for(k=0; k<pmlWidth; k++)
{
pmlSigmaZBottomYZ[i][j][k] = ReturnConductivity((pmlWidth-k) * dz, pmlWidth * dz, maxSigma);
pmlSigmaZStarBottomYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomYZ[i][j][k]][EPSILON]) * pmlSigmaZBottomYZ[i][j][k];
pmlSZBottomYZ[i][j][k] = ReturnStretching((pmlWidth-k) * dz, pmlWidth * dz, maxStretching, stretchSteepness);
}
for(i=0; i<pmlWidth; i++)
for(j=0; j<ny+2*pmlWidth; j++)
for(k=nz+pmlWidth; k<nz+2*pmlWidth; k++)
{
pmlSigmaZBottomYZ[i][j][k] = ReturnConductivity((k-pmlWidth-nz+1) * dz, pmlWidth * dz, maxSigma);
pmlSigmaZStarBottomYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomYZ[i][j][k]][EPSILON]) * pmlSigmaZBottomYZ[i][j][k];
pmlSZBottomYZ[i][j][k] = ReturnStretching((k-pmlWidth-nz+1) * dz, pmlWidth * dz, maxStretching, stretchSteepness);
}
/*****************************************************************************************/
/*************************************** TopYZ *******************************************/
/*****************************************************************************************/
for(i=0; i<pmlWidth; i++)
for(j=0; j<ny+2*pmlWidth; j++)
for(k=0; k<nz+2*pmlWidth; k++)
{
pmlSigmaXTopYZ[i][j][k] = ReturnConductivity((i+1) * dx, pmlWidth * dx, maxSigma);
pmlSigmaXStarTopYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopYZ[i][j][k]][EPSILON]) * pmlSigmaXTopYZ[i][j][k];
pmlSXTopYZ[i][j][k] = ReturnStretching((i+1) * dx, pmlWidth * dx, maxStretching, stretchSteepness);
}
/* computing sigmaY within the edges */
for(i=0; i<pmlWidth; i++)
for(j=0; j<pmlWidth; j++)
for(k=0; k<nz+2*pmlWidth; k++)
{
pmlSigmaYTopYZ[i][j][k] = ReturnConductivity((pmlWidth-j) * dy, pmlWidth * dy, maxSigma);
pmlSigmaYStarTopYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopYZ[i][j][k]][EPSILON]) * pmlSigmaYTopYZ[i][j][k];
pmlSYTopYZ[i][j][k] = ReturnStretching((pmlWidth-j) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
}
for(i=0; i<pmlWidth; i++)
for(j=ny+pmlWidth; j<ny+2*pmlWidth; j++)
for(k=0; k<nz+2*pmlWidth; k++)
{
pmlSigmaYTopYZ[i][j][k] = ReturnConductivity((j-pmlWidth-ny+1) * dy, pmlWidth * dy, maxSigma);
pmlSigmaYStarTopYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopYZ[i][j][k]][EPSILON]) * pmlSigmaYTopYZ[i][j][k];
pmlSYTopYZ[i][j][k] = ReturnStretching((j-pmlWidth-ny+1) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
}
/* computing sigmaZ within the edges */
for(i=0; i<pmlWidth; i++)
for(j=0; j<ny+2*pmlWidth; j++)
for(k=0; k<pmlWidth; k++)
{
pmlSigmaZTopYZ[i][j][k] = ReturnConductivity((pmlWidth-k) * dz, pmlWidth * dz, maxSigma);
pmlSigmaZStarTopYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopYZ[i][j][k]][EPSILON]) * pmlSigmaZTopYZ[i][j][k];
pmlSZTopYZ[i][j][k] = ReturnStretching((pmlWidth-k) * dz, pmlWidth * dz, maxStretching, stretchSteepness);
}
for(i=0; i<pmlWidth; i++)
for(j=0; j<ny+2*pmlWidth; j++)
for(k=nz+pmlWidth; k<nz+2*pmlWidth; k++)
{
pmlSigmaZTopYZ[i][j][k] = ReturnConductivity((k-pmlWidth-nz+1) * dz, pmlWidth * dz, maxSigma);
pmlSigmaZStarTopYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopYZ[i][j][k]][EPSILON]) * pmlSigmaZTopYZ[i][j][k];
pmlSZTopYZ[i][j][k] = ReturnStretching((k-pmlWidth-nz+1) * dz, pmlWidth * dz, maxStretching, stretchSteepness);
}
/*****************************************************************************************/
/************************************* BottomXY ******************************************/
/*****************************************************************************************/
for(i=0; i<nx; i++)
for(j=0; j<ny+2*pmlWidth; j++)
for(k=0; k<pmlWidth; k++)
{
pmlSigmaZBottomXY[i][j][k] = ReturnConductivity((pmlWidth-k) * dz, pmlWidth * dz, maxSigma);
pmlSigmaZStarBottomXY[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomXY[i][j][k]][EPSILON]) * pmlSigmaZBottomXY[i][j][k];
pmlSZBottomXY[i][j][k] = ReturnStretching((pmlWidth-k) * dz, pmlWidth * dz, maxStretching, stretchSteepness);
}
/* computing SigmaY within the edges */
for(i=0; i<nx; i++)
for(j=0; j<pmlWidth; j++)
for(k=0; k<pmlWidth; k++)
{
pmlSigmaYBottomXY[i][j][k] = ReturnConductivity((pmlWidth-j) * dy, pmlWidth * dy, maxSigma);
pmlSigmaYStarBottomXY[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomXY[i][j][k]][EPSILON]) * pmlSigmaYBottomXY[i][j][k];
pmlSYBottomXY[i][j][k] = ReturnStretching((pmlWidth-j) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
}
for(i=0; i<nx; i++)
for(j=ny+pmlWidth; j<ny+2*pmlWidth; j++)
for(k=0; k<pmlWidth; k++)
{
pmlSigmaYBottomXY[i][j][k] = ReturnConductivity((j-pmlWidth-ny+1) * dy, pmlWidth * dy, maxSigma);
pmlSigmaYStarBottomXY[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomXY[i][j][k]][EPSILON]) * pmlSigmaYBottomXY[i][j][k];
pmlSYBottomXY[i][j][k] = ReturnStretching((j-pmlWidth-ny+1) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
}
/*****************************************************************************************/
/************************************** TopXY ********************************************/
/*****************************************************************************************/
for(i=0; i<nx; i++)
for(j=0; j<ny+2*pmlWidth; j++)
for(k=0; k<pmlWidth; k++)
{
pmlSigmaZTopXY[i][j][k] = ReturnConductivity((k+1) * dz, pmlWidth * dz, maxSigma);
pmlSigmaZStarTopXY[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopXY[i][j][k]][EPSILON]) * pmlSigmaZTopXY[i][j][k];
pmlSZTopXY[i][j][k] = ReturnStretching((k+1) * dz, pmlWidth * dz, maxStretching, stretchSteepness);
}
/* computing SigmaY within the edges */
for(i=0; i<nx; i++)
for(j=0; j<pmlWidth; j++)
for(k=0; k<pmlWidth; k++)
{
pmlSigmaYTopXY[i][j][k] = ReturnConductivity((pmlWidth-j) * dy, pmlWidth * dy, maxSigma);
pmlSigmaYStarTopXY[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopXY[i][j][k]][EPSILON]) * pmlSigmaYTopXY[i][j][k];
pmlSYTopXY[i][j][k] = ReturnStretching((pmlWidth-j) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
}
for(i=0; i<nx; i++)
for(j=ny+pmlWidth; j<ny+2*pmlWidth; j++)
for(k=0; k<pmlWidth; k++)
{
pmlSigmaYTopXY[i][j][k] = ReturnConductivity((j-pmlWidth-ny+1) * dy, pmlWidth * dy, maxSigma);
pmlSigmaYStarTopXY[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopXY[i][j][k]][EPSILON]) * pmlSigmaYTopXY[i][j][k];
pmlSYTopXY[i][j][k] = ReturnStretching((j-pmlWidth-ny+1) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
}
/*****************************************************************************************/
/************************************ BottomXZ *******************************************/
/*****************************************************************************************/
for(i=0; i<nx; i++)
for(j=0; j<pmlWidth; j++)
for(k=0; k<nz; k++)
{
pmlSigmaYBottomXZ[i][j][k] = ReturnConductivity((pmlWidth-j) * dy, pmlWidth * dy, maxSigma);
pmlSigmaYStarBottomXZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomXZ[i][j][k]][EPSILON]) * pmlSigmaYBottomXZ[i][j][k];
pmlSYBottomXZ[i][j][k] = ReturnStretching((pmlWidth-j) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
}
/*****************************************************************************************/
/************************************** TopXZ ********************************************/
/*****************************************************************************************/
for(i=0; i<nx; i++)
for(j=0; j<pmlWidth; j++)
for(k=0; k<nz; k++)
{
pmlSigmaYTopXZ[i][j][k] = ReturnConductivity((j+1) * dy, pmlWidth * dy, maxSigma);
pmlSigmaYStarTopXZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialTopXZ[i][j][k]][EPSILON]) * pmlSigmaYTopXZ[i][j][k];
pmlSYTopXZ[i][j][k] = ReturnStretching((j+1) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
}
/*
#ifdef DEBUG
for(i=0; i<nx+2*pmlWidth; i++)
{
sprintf(filename, "pmlSigma%05d.dat", (int)i);
/ *sprintf(filename, "current/pic%05d.dat", (int) iteration / timeToPlot);* /
foutput = fopen(filename, "w");
for(j=0; j<ny+2*pmlWidth; j++)
for(k=0; k<nz+2*pmlWidth; k++)
{
if (i<pmlWidth)
{
fprintf(foutput, "%d\t%d\t%.6e\t%.6e\t%.6e\n",j,k,pmlSigmaXBottomYZ[i][j][k], pmlSigmaYBottomYZ[i][j][k], pmlSigmaZBottomYZ[i][j][k]);
}
else if(i<pmlWidth+nx)
{
if (k<pmlWidth)
{
fprintf(foutput, "%d\t%d\t%.6e\t%.6e\t%.6e\n",j,k,pmlSigmaXBottomXY[i-pmlWidth][j][k], pmlSigmaYBottomXY[i-pmlWidth][j][k], pmlSigmaZBottomXY[i-pmlWidth][j][k]);
}
else if (k<pmlWidth+nz)
{
if (j<pmlWidth)
{
fprintf(foutput, "%d\t%d\t%.6e\t%.6e\t%.6e\n",j,k,pmlSigmaXBottomXZ[i-pmlWidth][j][k-pmlWidth], pmlSigmaYBottomXZ[i-pmlWidth][j][k-pmlWidth], pmlSigmaZBottomXZ[i-pmlWidth][j][k-pmlWidth]);
}
else if (j<pmlWidth+ny)
{
fprintf(foutput, "%d\t%d\tNaN\tNan\tNaN\n", j,k);
}
else
{
fprintf(foutput, "%d\t%d\t%.6e\t%.6e\t%.6e\n",j,k,pmlSigmaXTopXZ[i-pmlWidth][j-ny-pmlWidth][k-pmlWidth], pmlSigmaYTopXZ[i-pmlWidth][j-ny-pmlWidth][k-pmlWidth], pmlSigmaZTopXZ[i-pmlWidth][j-ny-pmlWidth][k-pmlWidth]);
}
}
else
{
fprintf(foutput, "%d\t%d\t%.6e\t%.6e\t%.6e\n",j,k,pmlSigmaXTopXY[i-pmlWidth][j][k-nz-pmlWidth], pmlSigmaYTopXY[i-pmlWidth][j][k-nz-pmlWidth], pmlSigmaZTopXY[i-pmlWidth][j][k-nz-pmlWidth]);
}
}
else
{
fprintf(foutput, "%d\t%d\t%.6e\t%.6e\t%.6e\n",j,k,pmlSigmaXTopYZ[i-nx-pmlWidth][j][k], pmlSigmaYTopYZ[i-nx-pmlWidth][j][k], pmlSigmaZTopYZ[i-nx-pmlWidth][j][k]);
}
}
fclose(foutput);
}
#endif*/
}
void CradarFDTD::SetPMLMaterial(void)
{
int i,j,k;
int relPosX, relPosY, relPosZ;
for(i=0; i<nx+2*pmlWidth; i++)
for(j=0; j<ny+2*pmlWidth; j++)
for(k=0; k<nz+2*pmlWidth; k++)
{
if (i<pmlWidth)
relPosX = 0;
else if (i<nx+pmlWidth)
relPosX = i-pmlWidth;
else
relPosX = nx-1;
if (j<pmlWidth)
relPosY = 0;
else if (j<ny+pmlWidth)
relPosY = j-pmlWidth;
else
relPosY = ny-1;
if (k<pmlWidth)
relPosZ = 0;
else if (k<nz+pmlWidth)
relPosZ = k-pmlWidth;
else
relPosZ = nz-1;
if (i<pmlWidth)
pmlMaterialBottomYZ[i][j][k] = material[relPosX][relPosY][relPosZ];
else if (i<nx+pmlWidth)
if (k<pmlWidth)
pmlMaterialBottomXY[i-pmlWidth][j][k] = material[relPosX][relPosY][relPosZ];
else if(k<nz+pmlWidth)
if (j<pmlWidth)
pmlMaterialBottomXZ[i-pmlWidth][j][k-pmlWidth] = material[relPosX][relPosY][relPosZ];
else if (j<ny+pmlWidth)
{}
else
pmlMaterialTopXZ[i-pmlWidth][j-ny-pmlWidth][k-pmlWidth] = material[relPosX][relPosY][relPosZ];
else
pmlMaterialTopXY[i-pmlWidth][j][k-nz-pmlWidth] = material[relPosX][relPosY][relPosZ];
else
pmlMaterialTopYZ[i-nx-pmlWidth][j][k] = material[relPosX][relPosY][relPosZ];
}
/*
#ifdef DEBUG
sprintf(filename, "./material.dat");
foutput = fopen(filename, "w");
for(i=0; i<nx+2*pmlWidth; i++)
for(j=0; j<ny+2*pmlWidth; j++)
for(k=0; k<nz+2*pmlWidth; k++)
{
if (i<pmlWidth)
tmpMaterial = pmlMaterialBottomYZ[i][j][k];
else if (i<nx+pmlWidth)
if (k<pmlWidth)
tmpMaterial = pmlMaterialBottomXY[i-pmlWidth][j][k];
else if(k<nz+pmlWidth)
if (j<pmlWidth)
tmpMaterial = pmlMaterialBottomXZ[i-pmlWidth][j][k-pmlWidth];
else if (j<ny+pmlWidth)
tmpMaterial = -1;
else
tmpMaterial = pmlMaterialTopXZ[i-pmlWidth][j-ny-pmlWidth][k-pmlWidth];
else
tmpMaterial = pmlMaterialTopXY[i-pmlWidth][j][k-nz-pmlWidth];
else
tmpMaterial = pmlMaterialBottomYZ[i-nx-pmlWidth][j][k];
fprintf(foutput, "%d %d %d %d\n", i,j,k,tmpMaterial);
}
fclose(foutput);
#endif*/
}
void CradarFDTD::ShowWhirl(FILE *output, long int counter)
{
switch(counter % 4)
{
case 0: fprintf(output, "\b|");break;
case 1: fprintf(output, "\b/");break;
case 2: fprintf(output, "\b-");break;
default: fprintf(output, "\b\\");
}
}
/* Update the ex values: */
void CradarFDTD::UpdateInnerEx(void)
{
int i,j,k;
PRECISION lowerY, lowerZ;
for(i=0; i<(nx); i++)
for(j=0; j<(ny); j++)
for(k=0; k<(nz); k++)
{
if (j > 0)
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -