📄 radarfdtd.cpp
字号:
{
fprintf(stderr, "Not enough memory available!\n");
exit(1);
}
for(i=0; i<maxReceiver; i++)
{
receiver[i] = (int *) malloc(4 * sizeof(int));
if (receiver[i] == NULL)
{
fprintf(stderr, "Not enough memory available!\n");
exit(1);
}
}
tmpCounter = 0;
for(i=tmpReceiver[h][0]; i<=tmpReceiver[h][3]; i+=tmpReceiver[h][7])
for(j=tmpReceiver[h][1]; j<=tmpReceiver[h][4]; j+=tmpReceiver[h][7])
for(k=tmpReceiver[h][2]; k<=tmpReceiver[h][5]; k+=tmpReceiver[h][7])
{
receiver[tmpCounter][0] = i;
receiver[tmpCounter][1] = j;
receiver[tmpCounter][2] = k;
receiver[tmpCounter][3] = tmpReceiver[h][6]; /* mode */
tmpCounter++;
}
}
else
{
warningStop = 1;
fprintf(stderr, "WARNING: error in receiver box #%d coordinates \n", h);
}
}
for(h=0; h<numberOfReceiverBoxes; h++)
free(tmpReceiver[h]);
free(tmpReceiver);
/* NOW CHECK IF THE PROGRAM CAN CONTINUE BEYOND THIS POINT */
if (warningStop == 0)
{
fprintf(stdout, "Settings seem to be alright. Starting the computation.\n\n");
}
else
{
fprintf(stderr, "Since there is something wrong, please correct your configfile (%s/%s) and restart the program.\n", directory, filename);
}
/*********** Now the computations ***********/
/* Courant-Criterion */
tmpdt = 1 / (LIGHT_SPEED * sqrt(1/(dx*dx) + 1/(dy*dy) + 1/(dz*dz)));
if (dt < 0) dt = -tmpdt * (PRECISION) dt;
/* material constants */
/* set values for vacuum */
materialConstants[0][DTMUDX] = dt/(MU_0*dx);
materialConstants[0][DTEPSDX] = dt/(EPSILON_0*dx);
materialConstants[0][DTMUDY] = dt/(MU_0*dy);
materialConstants[0][DTEPSDY] = dt/(EPSILON_0*dy);
materialConstants[0][DTMUDZ] = dt/(MU_0*dz);
materialConstants[0][DTEPSDZ] = dt/(EPSILON_0*dz);
materialConstants[0][DTSIGEPS] = 1.0;
materialConstants[0][EPSILON] = 1.0;
materialConstants[0][SIGMA] = 0.0;
/* for ideal conductors (this will _not_ be used, but just in case) */
materialConstants[1][DTMUDX] = dt/(MU_0*dx);
materialConstants[1][DTEPSDX] = 0.0;
materialConstants[1][DTMUDY] = dt/(MU_0*dy);
materialConstants[1][DTEPSDY] = 0.0;
materialConstants[1][DTMUDZ] = dt/(MU_0*dz);
materialConstants[1][DTEPSDZ] = 0.0;
materialConstants[1][DTSIGEPS] = 0.0;
materialConstants[1][EPSILON] = 0.0;
materialConstants[1][SIGMA] = 0.0;
/*and for the various other materials */
for(i=0; i<maxMaterials; i++)
{
tmpEpsilon = materialConstants[i+2][0];
tmpSigma = materialConstants[i+2][1];
materialConstants[i+2][DTMUDX] = dt/(MU_0 * dx);
materialConstants[i+2][DTEPSDX] = dt / (dx * EPSILON_0 * tmpEpsilon * (1 + (tmpSigma * dt) / (2.0 * EPSILON_0 * tmpEpsilon)));
materialConstants[i+2][DTMUDY] = dt/(MU_0 * dy);
materialConstants[i+2][DTEPSDY] = dt / (dy * EPSILON_0 * tmpEpsilon * (1 + (tmpSigma * dt) / (2.0 * EPSILON_0 * tmpEpsilon)));
materialConstants[i+2][DTMUDZ] = dt/(MU_0 * dz);
materialConstants[i+2][DTEPSDZ] = dt / (dz * EPSILON_0 * tmpEpsilon * (1 + (tmpSigma * dt) / (2.0 * EPSILON_0 * tmpEpsilon)));
materialConstants[i+2][DTSIGEPS] = (1 - tmpSigma * dt / (2.0 * EPSILON_0 * tmpEpsilon)) / (1 + tmpSigma * dt / (2.0 * EPSILON_0 * tmpEpsilon));
materialConstants[i+2][EPSILON] = tmpEpsilon;
materialConstants[i+2][SIGMA] = tmpSigma;
}
/* Calculate the maximum space - step size */
maxStep = dx;
if (maxStep < dy) maxStep = dy;
if (maxStep < dz) maxStep = dz;
/* Calculate the number of iteration steps needed */
maxIteration = (int) (1.0 + simulationLength / dt);
/* Calculate theoretical reflection from the absorber */
tmpReflection = 0.0;
for(i=0; i<pmlWidth; i++)
tmpReflection += ReturnStretching((i+1) * maxStep, pmlWidth * maxStep, maxStretching, stretchSteepness) * ReturnConductivity((i+1) * maxStep, pmlWidth * maxStep, maxSigma);
theoreticalReflection = exp(-2.0 / (maxEpsilon * EPSILON_0 * LIGHT_SPEED) * tmpReflection * maxStep);
minWavelength = (1.0 + maxStretching) * 3.0 * maxStep;
/***************** Output section ******************/
fprintf(stdout, "\n******************************************************************\n");
fprintf(stdout, VERSION_INFO);
fprintf(stdout, "\n******************************************************************\n");
fprintf(stdout, "\nConfiguration file: %s\n\n", filename);
fprintf(stdout, "Simulation steps: %d (%.4gns), dt = %.4gns, output every %d steps\n\n", maxIteration, 1e9 * (float) maxIteration * dt, 1e9*dt,timeToPlotReceiver);
fprintf(stdout, "Simulation space: %.4f x %.4f x %.4f m砛n", dx*nx, dy*ny, dz*nz);
fprintf(stdout, "Number of cells:\t%d\t%d\t%d\n", nx,ny,nz);
fprintf(stdout, "Cellsizes:\t\t%.4gm\t%.4gm\t%.4gm\n\n", dx,dy,dz);
fprintf(stdout, "Number of absorbing boundary layers: %d\nMaximum conductivity at edge: %.4g, maximum Stretching: %.4g\n\n", pmlWidth, (float) maxSigma, (float) maxStretching);
fprintf(stdout, "\nThe derived theoretical reflection (after Fang and Wu) is: %e,\nand the minimum absorber wavelength is", theoreticalReflection);
PrintCorrectOrder(minWavelength, stdout);
fprintf(stdout, "m (freq: ");
PrintCorrectOrder(LIGHT_SPEED / minWavelength, stdout);
fprintf(stdout, "Hz)\n\n");
fprintf(stdout, "Number of materials: %d\n", (maxMaterials+2));
for(i=0; i<maxMaterials+2; i++)
if (i != 1)
fprintf(stdout, "%d:\tEpsilon: %.6g\tSigma: %.6g\n", i, materialConstants[i][7], materialConstants[i][8]);
else
fprintf(stdout, "1:\tperfect conductor\n");
fprintf(stdout, "\nNumber of tranmsitters: %d\n", maxTransmitter);
fprintf(stdout, "\nNumber of receivers: %d\n\n", maxReceiver);
fprintf(stdout, "Remarks:\nCourant-criterion ");
if (dt < tmpdt)
fprintf(stdout, "OK.\n");
else
fprintf(stdout, "NOT OK! (dt / good) < %6.2g\n", (tmpdt/dt));
fprintf(stdout, "Cellsize: ");
if (maxStep < LIGHT_SPEED/(10.0 * transmitterFrequency * sqrt(maxEpsilon)))
fprintf(stdout, "OK\n");
else
{
fprintf(stdout, "NOT OK! Recommended value is at least:");
PrintCorrectOrder(LIGHT_SPEED/(10 * transmitterFrequency * sqrt(maxEpsilon)), stdout);
fprintf(stdout, "m\n\n");
}
/* close file */
fclose(configFile);
}
PRECISION CradarFDTD::ReturnConductivity(PRECISION position, PRECISION pmlThickness, PRECISION maxConductivity)
{
/* I use the profile discribed in: Fang and Wu, Generlized PML for the Absorpion of Porpagating and Evanescent Waves, IEEE Transactions on Microwave Theory and Techniques, Volume 44, No.12, December 1996 */
PRECISION tmp = sin(M_PI * position / (2.0 * pmlThickness));
return maxConductivity * tmp * tmp;
/* PRECISION tmp;
switch(((int) (0.5+position/pmlThickness)))
{
case 0: tmp= 0.001064;break;
case 1: tmp= 0.001210;break;
case 2: tmp= 0.001296;break;
case 3: tmp= 0.001385;break;
case 4: tmp= 0.001487;break;
case 5: tmp= 0.001606;break;
case 6: tmp= 0.001746;break;
case 7: tmp= 0.001916;break;
case 8: tmp= 0.002126;break;
case 9: tmp= 0.002395;break;
case 10: tmp= 0.002752;break;
case 11: tmp= 0.003257;break;
case 12: tmp= 0.004040;break;
case 13: tmp= 0.005446;break;
case 14: tmp= 0.008747;break;
case 15: tmp= 0.042007;break;
default: tmp= 0.042007;
}
return tmp;*/
}
PRECISION CradarFDTD::ReturnFieldStrength(int x, int y, int z, int mode)
{
PRECISION result = 0.0;
/* basically this is only a big 'if' */
if (x < 0) /* use BottomYZ */
switch (mode) {
case 0: result = sqrt(pmlHXYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlHXYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHXZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlHXZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHYXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlHYXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHYZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlHYZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHZXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlHZXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHZYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlHZYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth]);break;
case 1: result = pmlHXYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHXZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth];break;
case 2: result = pmlHYXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHYZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth];break;
case 3: result = pmlHZXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlHZYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth];break;
case 4: result = pmlEXYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEXZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth];break;
case 5: result = pmlEYXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEYZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth];break;
case 6: result = pmlEZXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEZYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth];break;
default: result = sqrt(pmlEXYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlEXYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEXZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlEXZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEYXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlEYXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEYZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlEYZBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEZXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlEZXBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] + pmlEZYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth] * pmlEZYBottomYZ[x+pmlWidth][y+pmlWidth][z+pmlWidth]);
}
else if (x<nx)
if (z<0)
switch (mode) { /* BottomXY */
case 0: result = sqrt(pmlHXYBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlHXYBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHXZBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlHXZBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHYXBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlHYXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHYZBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlHYZBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHZXBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlHZXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHZYBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlHZYBottomXY[x][y+pmlWidth][z+pmlWidth]);break;
case 1: result = pmlHXYBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHXZBottomXY[x][y+pmlWidth][z+pmlWidth];break;
case 2: result = pmlHYXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHYZBottomXY[x][y+pmlWidth][z+pmlWidth];break;
case 3: result = pmlHZXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlHZYBottomXY[x][y+pmlWidth][z+pmlWidth];break;
case 4: result = pmlEXYBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEXZBottomXY[x][y+pmlWidth][z+pmlWidth];break;
case 5: result = pmlEYXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEYZBottomXY[x][y+pmlWidth][z+pmlWidth];break;
case 6: result = pmlEZXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEZYBottomXY[x][y+pmlWidth][z+pmlWidth];break;
default: result = sqrt(pmlEXYBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlEXYBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEXZBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlEXZBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEYXBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlEYXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEYZBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlEYZBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEZXBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlEZXBottomXY[x][y+pmlWidth][z+pmlWidth] + pmlEZYBottomXY[x][y+pmlWidth][z+pmlWidth] * pmlEZYBottomXY[x][y+pmlWidth][z+pmlWidth]);
}
else
if (z < nz)
if (y < 0)
switch (mode) { /* BottomXZ */
case 0: result = sqrt(pmlHXYBottomXZ[x][y+pmlWidth][z] * pmlHXYBottomXZ[x][y+pmlWidth][z] + pmlHXZBottomXZ[x][y+pmlWidth][z] * pmlHXZBottomXZ[x][y+pmlWidth][z] + pmlHYXBottomXZ[x][y+pmlWidth][z] * pmlHYXBottomXZ[x][y+pmlWidth][z] + pmlHYZBottomXZ[x][y+pmlWidth][z] * pmlHYZBottomXZ[x][y+pmlWidth][z] + pmlHZXBottomXZ[x][y+pmlWidth][z] * pmlHZXBottomXZ[x][y+pmlWidth][z] + pmlHZYBottomXZ[x][y+pmlWidth][z] * pmlHZYBottomXZ[x][y+pmlWidth][z]);break;
case 1: result = pmlHXYBottomXZ[x][y+pmlWidth][z] + pmlHXZBottomXZ[x][y+pmlWidth][z];break;
case 2: result = pmlHYXBottomXZ[x][y+pmlWidth][z] + pmlHYZBottomXZ[x][y+pmlWidth][z];break;
case 3: result = pmlHZXBottomXZ[x][y+pmlWidth][z] + pmlHZYBottomXZ[x][y+pmlWidth][z];break;
case 4: result = pmlEXYBottomXZ[x][y+pmlWidth][z] + pmlEXZBottomXZ[x][y+pmlWidth][z];break;
case 5: result = pmlEYXBottomXZ[x][y+pmlWidth][z] + pmlEYZBottomXZ[x][y+pmlWidth][z];break;
case 6: result = pmlEZXBottomXZ[x][y+pmlWidth][z] + pmlEZYBottomXZ[x][y+pmlWidth][z];break;
default: result = sqrt(pmlEXYBottomXZ[x][y+pmlWidth][z] * pmlEXYBottomXZ[x][y+pmlWidth][z] + pmlEXZBottomXZ[x][y+pmlWidth][z] * pmlEXZBottomXZ[x][y+pmlWidth][z] + pmlEYXBottomXZ[x][y+pmlWidth][z] * pmlEYXBottomXZ[x][y+pmlWidth][z] + pmlEYZBottomXZ[x][y+pmlWidth][z] * pmlEYZBottomXZ[x][y+pmlWidth][z] + pmlEZXBottomXZ[x][y+pmlWidth][z] * pmlEZXBottomXZ[x][y+pmlWidth][z] + pmlEZYBottomXZ[x][y+pmlWidth][z] * pmlEZYBottomXZ[x][y+pmlWidth][z]);
}
else if(y < ny)
switch (mode){ /* Sim Space */
case 0: result = sqrt(hx[x][y][k] * hx[x][y][z] + hy[x][y][k] * hy[x][y][z] + hz[x][y][k] * hz[x][y][z]); break;
case 1: result = hx[x][y][k]; break;
case 2: result = hy[x][y][k]; break;
case 3: result = hz[x][y][k]; break;
case 4: result = ex[x][y][k]; break;
case 5: result = ey[x][y][k]; break;
case 6: result = ez[x][y][k]; break;
default: result = sqrt(ex[x][y][z] * ex[x][y][z] + ey[x][y][z] * ey[x][y][z] + ez[x][y][z] * ez[x][y][z]);
}
else
switch (mode) { /* TopXZ */
case 0: result = sqrt(pmlHXYTopXZ[x][y-ny][z] * pmlHXYTopXZ[x][y-ny][z] + pmlHXZTopXZ[x][y-ny][z] * pmlHXZTopXZ[x][y-ny][z] + pmlHYXTopXZ[x][y-ny][z] * pmlHYXTopXZ[x][y-ny][z] + pmlHYZTopXZ[x][y-ny][z] * pmlHYZTopXZ[x][y-ny][z] + pmlHZXTopXZ[x][y-ny][z] * pmlHZXTopXZ[x][y-ny][z] + pmlHZYTopXZ[x][y-ny][z] * pmlHZYTopXZ[x][y-ny][z]);break;
case 1: result = pmlHXYTopXZ[x][y-ny][z] + pmlHXZTopXZ[x][y-ny][z];break;
case 2: result = pmlHYXTopXZ[x][y-ny][z] + pmlHYZTopXZ[x][y-ny][z];break;
case 3: result = pmlHZXTopXZ[x][y-ny][z] + pmlHZYTopXZ[x][y-ny][z];break;
case 4: result = pmlEXYTopXZ[x][y-ny][z] + pmlEXZTopXZ[x][y-ny][z];break;
case 5: result = pmlEYXTopXZ[x][y-ny][z] + pmlEYZTopXZ[x][y-ny][z];break;
case 6: result = pmlEZXTopXZ[x][y-ny][z] + pmlEZYTopXZ[x][y-ny][z];break;
default: result = sqrt(pmlEXYTopXZ[x][y-ny][z] * pmlEXYTopXZ[x][y-ny][z] + pmlEXZTopXZ[x][y-ny][z] * pmlEXZTopXZ[x][y-ny][z] + pmlEYXTopXZ[x][y-ny][z] * pmlEYXTopXZ[x][y-ny][z] + pmlEYZTopXZ[x][y-ny][z] * pmlEYZTopXZ[x][y-ny][z] + pmlEZXTopXZ[x][y-ny][z] * pmlEZXTopXZ[x][y-ny][z] + pmlEZYTopXZ[x][y-ny][z] * pmlEZYTopXZ[x][y-ny][z]);
}
else
switch (mode) { /* TopXY */
case 0: result = sqrt(pmlHXYTopXY[x][y+pmlWidth][z-nz] * pmlHXYTopXY[x][y+pmlWidth][z-nz] + pmlHXZTopXY[x][y+pmlWidth][z-nz] * pmlHXZTopXY[x][y+pmlWidth][z-nz] + pmlHYXTopXY[x][y+pmlWidth][z-nz] * pmlHYXTopXY[x][y+pmlWidth][z-nz] + pmlHYZTopXY[x][y+pmlWidth][z-nz] * pmlHYZTopXY[x][y+pmlWidth][z-nz] + pmlHZXTopXY[x][y+pmlWidth][z-nz] * pmlHZXTopXY[x][y+pmlWidth][z-nz] + pmlHZYTopXY[x][y+pmlWidth][z-nz] * pmlHZYTopXY[x][y+pmlWidth][z-nz]);break;
case 1: result = pmlHXYTopXY[x][y+pmlWidth][z-nz] + pmlHXZTopXY[x][y+pmlWidth][z-nz];break;
case 2: result = pmlHYXTopXY[x][y+pmlWidth][z-nz] + pmlHYZTopXY[x][y+pmlWidth][z-nz];break;
case 3: result = pmlHZXTopXY[x][y+pmlWidth][z-nz] + pmlHZYTopXY[x][y+pmlWidth][z-nz];break;
case 4: result = pmlEXYTopXY[x][y+pmlWidth][z-nz] + pmlEXZTopXY[x][y+pmlWidth][z-nz];break;
case 5: result = pmlEYXTopXY[x][y+pmlWidth][z-nz] + pmlEYZTopXY[x][y+pmlWidth][z-nz];break;
case 6: result = pmlEZXTopXY[x][y+pmlWidth][z-nz] + pmlEZYTopXY[x][y+pmlWidth][z-nz];break;
default: result = sqrt(pmlEXYTopXY[x][y+pmlWidth][z-nz] * pmlEXYTopXY[x][y+pmlWidth][z-nz] + pmlEXZTopXY[x][y+pmlWidth][z-nz] * pmlEXZTopXY[x][y+pmlWidth][z-nz] + pmlEYXTopXY[x][y+pmlWidth][z-nz] * pmlEYXTopXY[x][y+pmlWidth][z-nz] + pmlEYZTopXY[x][y+pmlWidth][z-nz] * pmlEYZTopXY[x][y+pmlWidth][z-nz] + pmlEZXTopXY[x][y+pmlWidth][z-nz] * pmlEZXTopXY[x][y+pmlWidth][z-nz] + pmlEZYTopXY[x][y+pmlWidth][z-nz] * pmlEZYTopXY[x][y+pmlWidth][z-nz]);
}
else
switch (mode) { /* TopYZ */
case 0: result = sqrt(pmlHXYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlHXYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHXZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlHXZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHYXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlHYXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHYZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlHYZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHZXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlHZXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHZYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlHZYTopYZ[x-nx][y+pmlWidth][z+pmlWidth]);break;
case 1: result = pmlHXYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHXZTopYZ[x-nx][y+pmlWidth][z+pmlWidth];break;
case 2: result = pmlHYXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHYZTopYZ[x-nx][y+pmlWidth][z+pmlWidth];break;
case 3: result = pmlHZXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlHZYTopYZ[x-nx][y+pmlWidth][z+pmlWidth];break;
case 4: result = pmlEXYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEXZTopYZ[x-nx][y+pmlWidth][z+pmlWidth];break;
case 5: result = pmlEYXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEYZTopYZ[x-nx][y+pmlWidth][z+pmlWidth];break;
case 6: result = pmlEZXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEZYTopYZ[x-nx][y+pmlWidth][z+pmlWidth];break;
default: result = sqrt(pmlEXYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlEXYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEXZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlEXZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEYXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlEYXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEYZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlEYZTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEZXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlEZXTopYZ[x-nx][y+pmlWidth][z+pmlWidth] + pmlEZYTopYZ[x-nx][y+pmlWidth][z+pmlWidth] * pmlEZYTopYZ[x-nx][y+pmlWidth][z+pmlWidth]);
}
return result;
}
PRECISION CradarFDTD::ReturnStretching (PRECISION position, PRECISION pmlThickness, PRECISION maxStretching, int stretchSteepness)
{
/* I use the profile discribed in: Fang and Wu, Generlized PML for the Absorpion of Porpagating and Evanescent Waves, IEEE Transactions on Microwave Theory and Techniques, Volume 44, No.12, December 1996 */
/* formula for this profile: s(pos) = 1 + s_max * (pos / thickness)**steepness */
PRECISION tmp = 1.0 + maxStretching * pow(position / pmlThickness, stretchSteepness);
return tmp;
}
void CradarFDTD::SetConductivities(void)
{
/*****************************************************************************************/
/************************************* BottomYZ ******************************************/
/*****************************************************************************************/
for(i=0; i<pmlWidth; i++)
for(j=0; j<ny+2*pmlWidth; j++)
for(k=0; k<nz+2*pmlWidth; k++)
{
pmlSigmaXBottomYZ[i][j][k] = ReturnConductivity((pmlWidth-i) * dx, pmlWidth * dx, maxSigma);
pmlSigmaXStarBottomYZ[i][j][k] = MU_0 / (EPSILON_0 * materialConstants[pmlMaterialBottomYZ[i][j][k]][EPSILON]) * pmlSigmaXBottomYZ[i][j][k];
pmlSXBottomYZ[i][j][k] = ReturnStretching((pmlWidth-i) * 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++)
{
pmlSigmaYBottomYZ[i][j][k] = ReturnConductivity((pmlWidth-j) * 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((pmlWidth-j) * dy, pmlWidth * dy, maxStretching, stretchSteepness);
}
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -