/* ******************************************************************** wave2d.c : 2-D Wave Equation FDTD TE code with Split Field PML ******************************************************************** A 2-D Wave Equation FDTD TE implementation of Berenger's Split Field PML ABC version 1.0, 5/22/2008, Doug Neubauer This program implements a FDTD of the Maxwell's equation using the "wave equation" formulation, instead of the standard "Yee" curl equations. In addition, it also implements a "wave equation" version of Berenger's Split Field PML ABC. This program builds upon a previous work: Aoyagi, Lee, Mittra; "A Hybrid Yee Algorithm/Scalar-Wave Equation Approach"; IEEE Transactions on Microwave Theory and Techniques, vol. 41, no. 9, September 1993 which used the wave equation, but not the PML ABC. (As it was written before the Berenger PML ABC was created.) This program is also adapted from the earlier standard 2-D Curl Equation code of Dr. Susan C. Hagness (See Taflove2000) The simulation results were compared with the simulation results of the standard "Yee" curl equations version. The results were essentially identical (to within round off errors (around 12 or 13 decimal places)) This program probably has no practical value as implementing the wave equation version of the PML is much more memory intensive and compute intensive than the standard Yee algorithm. One of the penalties of the wave equation is that it is necessary to keep around an old value of Ex and Ey for every point on the grid. This program is provided as an existence proof. Implementing the PML as a wave equation can be done! (But you probably wouldn't want to!!) Note: no attempt was made to Optimize the algorithm, as it was only intended as a proof of concept. Therefore, for simplicity, the pml equations were also used in the main grid as well. This results in a slower running program, but an easier one to debug. The point of the wave equation is to "decouple" the H and E fields. This can actually be done, but with the penalty of more complicated field equations and additional memory requirements. This program's algorithm uses Ex and Ey only (Hz is used only for display purposes, ie to compare with earlier simulations. such as the Yee Pml program) Note: The PML split fields are different than in the standard PML algorithm. In the standard 2-D algorithm the split fields are: Hzx and Hzy. In the wave equation version the split fields are: Exxa, Exxb, Exya, Exyb, Eyxa, Eyxb, Eyya, Eyyb. Much more complicated! and more memory needed as well! ------------- This program implements the finite-difference time-domain solution of Maxwell's wave equations over a two-dimensional Cartesian space lattice comprised of uniform square grid cells. To illustrate the algorithm, a 6-cm-diameter metal cylindrical scatterer in free space is modeled. The source excitation is an Infinite Magnetic Line-source Gaussian pulse with a carrier frequency of 5 GHz. The grid resolution (dx = 3 mm) was chosen to provide 20 samples per wavelength at the center frequency of the pulse (which in turn provides approximately 10 samples per wavelength at the high end of the excitation spectrum, around 10 GHz). The computational domain is truncated using the perfectly matched layer (PML) absorbing boundary conditions. The formulation used in this code is based on the original split-field Berenger PML. The PML regions are labeled as shown in the following diagram: xSize-1,ySize-1 / ----------------------------------------------. | | BACK PML | | ---------------------------------------------- |L | | R| |E | | I| |F | | G| |T | | H| | | MAIN GRID | T| |P | | | |M | | P| |L | | M| | | | L| ---------------------------------------------- |. | FRONT PML | | /---------------------------------------------- 0,0 Below: Detailed view of the Yee grid... (where N=xSize, M=ySize, see below) Note: an extra column of ey on right edge and an extra row of ex on back edge Note: ey at x=0 and x=N are PEC and ex at y=0 and y=M are PEC. (0,M) (N-1,M) ___ex___ ___ex___ ___ex___ ___ex___ .. ___ex___ ___ex___ ___ex___ ......... ......... ......... ......... ......... ......... ......... | .| .| .| . | .| .| .| | 0,M-1 .| .| .| . | .| .| N-1,M-1.|(N,M-1) ey hz .ey hz .ey hz .ey hz . ey hz .ey hz .ey hz .ey | .| .| .| . | .| .| .| |___ex___.|___ex___.|___ex___.|___ex___. |___ex___.|___ex___.|___ex___.| ......... ......... ......... ........... ......... ......... ......... | .| .| .| . | .| .| .| | 0,M-2 .| .| .| . | .| .| N-1,M-2.| ey hz .ey hz .ey hz .ey hz . ey hz .ey hz .ey hz .ey | .| .| .| . | .| .| .| |___ex___.|___ex___.|___ex___.|___ex___...|___ex___.|___ex___.|___ex___.| . . . . . . . . ......... ......... ......... ........... ......... ......... ......... | .| .| .| . | .| .| .| | 0,2 .| .| .| . | .| .| N-1,2 .| ey hz .ey hz .ey hz .ey hz . ey hz .ey hz .ey hz .ey | .| .| .| . | .| .| .| |___ex___.|___ex___.|___ex___.|___ex___. |___ex___.|___ex___.|___ex___.| ......... ......... ......... ........... ......... ......... ......... | .| .| .| . | .| .| .| | 0,1 .| .| .| . | .| .| N-1,1 .| ey hz .ey hz .ey hz .ey hz . ey hz .ey hz .ey hz .ey | .| .| .| . | .| .| .| |___ex___.|___ex___.|___ex___.|___ex___. |___ex___.|___ex___.|___ex___.| ......... ......... ......... ........... ......... ......... ......... | .| .| .| . | .| .| .| | 0,0 .| 1,0 .| 2,0 .| 3.0 . | N-3,0 .| N-2,0 .| N-1,0 .|(N,0) ey hz .ey hz .ey hz .ey hz . ey hz .ey hz .ey hz .ey | .| .| .| . | .| .| .| |___ex___.|___ex___.|___ex___.|___ex___...|___ex___.|___ex___.|___ex___.| PML Reflection Simulation Results: ------------------------------------------ #layer r0 gradingOrder dB ------ ------ ------------ ----- 8 1.0e-7 2 -74.4 8 1.0e-7 3 -81.2 <<== use this one 8 1.0e-7 4 -81.15 10 1.0e-7 3 -90.3 10 1.0e-8 3 6 1.0e-7 3 -71.8 4 1.0e-7 3 -51.8 Results are the same as the standard E/H curl equations split-field. Also a Comparison of the cylinder simulation was identical to the standard E/H curl equations split-field simulation. In the spirit of the ToyFdtd programs, a point was made to try to heavily comment the source code. ********************************************************************** */ #include #include #include #define ABCSIZECONSTANT (8) // thickness of PML region #define MEDIACONSTANT (2) // number of different medias, ie 2: vacuum, metallicCylinder #define NUMBEROFITERATIONCONSTANT (400) // was 300 #define BIGLINESIZE (8192) // for creating an output filename #define SOURCESIZE (4) // hack for source // These global variables are used by ExecuteFdtd() int xSize ; // number of total grid cells in x-direction (left+main+right) int ySize ; // number of total grid cells in y-direction (front+main+back) int maximumIteration; // how long to run the simulation int xSource, ySource; // location of z-directed hard source double **ex; // the fields double **ey; // "" double **exxa; // split-fields for pml double **exxb; // "" double **exya; // "" double **exyb; // "" double **exxaOld; // "" double **exxbOld; // "" double **exyaOld; // "" double **exybOld; // "" double **eyxa; // "" double **eyxb; // "" double **eyya; // "" double **eyyb; // "" double **eyxaOld; // "" double **eyxbOld; // "" double **eyyaOld; // "" double **eyybOld; // "" double **hz; // "" used only for plotting double **caex; // fdtd/pml coefficents double **cbex; // "" double **caey; // "" double **cbey; // "" double **dahzx; // pml coefficient double **dbhzx; // "" double **dahzy; // "" double **dbhzy; // "" // some variables to Implement the source // The source is a Hz Magnetic infinite line source located at x=Xsource+0.5, y=ySource+0.5 // The source is actually implemented without using Hz. Instead 4 equations for the 4 Ex,Ey field points // around the Hz source point is used. double cbSource,dbSource; // coefficients for the source equation double currentSourceValue, oldSourceValue; // hold the current and old values of the source double eOldSource[SOURCESIZE]; // hold old values of Ex and Ey in the region around the source double sourceValue[NUMBEROFITERATIONCONSTANT]; // holds the pre-calculated values of the source, for the run of the simulation // standard C memory allocation for 2-D array double **AllocateMemory (int imax, int jmax, double initialValue) { int i,j; double **pointer; pointer = (double **)malloc(imax * sizeof(double *)); for (i = 0; i < imax; i++) { pointer[i] = (double *)malloc(jmax * sizeof(double)); for (j = 0; j < jmax; j++) { pointer[i][j] = initialValue; } /* jForLoop */ } /* iForLoop */ return(pointer) ; } // standard C memory allocation for 1-D array double *AllocateMemory1D (int size, double initialValue) { int j; double *pointer; pointer = (double *)malloc(size * sizeof(double)); for (j = 0; j < size; j++) { pointer[j] = initialValue; } /* jForLoop */ return(pointer) ; } void InitializeFdtd () { double mediaPermittivity[MEDIACONSTANT] = {1.0, 1.0}; // eps, index=0 is for vacuum, index=1 is for the metallic cylinder double mediaConductivity[MEDIACONSTANT] = {0.0, 1.0e+7}; // sig, double mediaPermeability[MEDIACONSTANT] = {1.0, 1.0}; // mur double mediaResistivity[MEDIACONSTANT] = {0.0, 0.0}; // sim double mediaCa[MEDIACONSTANT]; double mediaCb[MEDIACONSTANT]; double mediaDa[MEDIACONSTANT]; double mediaDb[MEDIACONSTANT]; double pi,speedOfLight,magneticPermeability0,electricalPermittivity0,frequency,wavelength,angularFrequency; double dx,dt,reflectionCoefficient0,gradingOrder,temporary; int i,j, media, xSizeMain,ySizeMain; int abcSize ; double cylinderDiameter, cylinderRadius, temporaryi,temporaryj,distance2 ; int cylinderXCenter,cylinderYCenter; double x,x1,x2; double electricalConductivityMaximum, boundaryWidth, gradientConductivity, gradientResistivity, boundaryFactor; double gradientCa1[ABCSIZECONSTANT]; double gradientCb1[ABCSIZECONSTANT]; double gradientDa1[ABCSIZECONSTANT]; double gradientDb1[ABCSIZECONSTANT]; double rtau, tau, delay; // char ch; //*********************************************************************** // Fundamental constants //*********************************************************************** pi = (acos(-1.0)); speedOfLight = 2.99792458e8; //speed of light in free space (meters/second) magneticPermeability0 = 4.0 * pi * 1.0e-7; //permeability of free space electricalPermittivity0 = 1.0 / (speedOfLight * speedOfLight * magneticPermeability0); //permittivity of free space frequency = 5.0e+9; //center frequency of source excitation (Hz) wavelength = speedOfLight / frequency ; //center wavelength of source excitation angularFrequency = 2.0 * pi * frequency; //center frequency in radians //*********************************************************************** // Grid parameters //*********************************************************************** xSizeMain = 100; // number of main grid cells in x-direction // ySizeMain = 100; // number of main grid cells in y-direction ySizeMain = 50; // number of main grid cells in y-direction (cylinder simulation) abcSize = ABCSIZECONSTANT; // thickness of PML region xSize = xSizeMain + 2 * abcSize; // number of total grid cells in x-direction ySize = ySizeMain + 2 * abcSize; // number of total grid cells in y-direction // xSource = (xSizeMain / 2) + abcSize; //location of z-directed hard source // ySource = (ySizeMain / 2) + abcSize; //location of z-directed hard source xSource = 15 + abcSize; //location of z-directed hard source (cylinder simulation) ySource = ySize / 2; //location of z-directed hard source dx = 3.0e-3; //space increment of square lattice (meters) dt = dx / (2.0 * speedOfLight); //time step, seconds, courant limit, Taflove1995 page 177 maximumIteration = NUMBEROFITERATIONCONSTANT; //total number of time steps reflectionCoefficient0 = 1.0e-7; // for PML, Nikolova part4 p.25 gradingOrder = 3; // for PML, (m) was 2; optimal values: 2 <= m <= 6, Nikolova part4 p.29 //*********************************************************************** // Material parameters //*********************************************************************** media = MEDIACONSTANT; // number of different medias, ie 2: vacuum, metallicCylinder //*********************************************************************** // Wave excitation //*********************************************************************** #if 0 for (i = 0; i < maximumIteration; i++) { temporary = (double )i; sourceValue[i] = sin( angularFrequency * temporary * dt); // simple sine wave } /* iForLoop */ #endif #if 1 rtau = 160.0e-12; tau = rtau / dt; delay = 3 * tau; for (i = 0; i < maximumIteration; i++) { sourceValue[i] = 0.0; } /* iForLoop */ for (i = 0; i < (int )(7.0 * tau); i++) { // Gaussian pulse with a carrier frequency of 5 GHz temporary = (double )i - delay; sourceValue[i] = sin( angularFrequency * (temporary) * dt) * exp(-( (temporary * temporary)/(tau * tau) ) ); } /* iForLoop */ #endif //*********************************************************************** // Field arrays //*********************************************************************** ex = AllocateMemory(xSize, ySize + 1, 0.0 ); // 1 extra in y direction for pec ey = AllocateMemory(xSize + 1,ySize, 0.0 ); // 1 extra in x direction for pec exxa = AllocateMemory(xSize, ySize, 0.0 ); // split-fields for pml (for simplicity of the equations, use them also in the main grid) exxb = AllocateMemory(xSize, ySize, 0.0 ); exya = AllocateMemory(xSize, ySize, 0.0 ); exyb = AllocateMemory(xSize, ySize, 0.0 ); exxaOld = AllocateMemory(xSize, ySize, 0.0 ); exxbOld = AllocateMemory(xSize, ySize, 0.0 ); exyaOld = AllocateMemory(xSize, ySize, 0.0 ); exybOld = AllocateMemory(xSize, ySize, 0.0 ); eyxa = AllocateMemory(xSize, ySize, 0.0 ); eyxb = AllocateMemory(xSize, ySize, 0.0 ); eyya = AllocateMemory(xSize, ySize, 0.0 ); eyyb = AllocateMemory(xSize, ySize, 0.0 ); eyxaOld = AllocateMemory(xSize, ySize, 0.0 ); eyxbOld = AllocateMemory(xSize, ySize, 0.0 ); eyyaOld = AllocateMemory(xSize, ySize, 0.0 ); eyybOld = AllocateMemory(xSize, ySize, 0.0 ); hz = AllocateMemory(xSize, ySize, 0.0 ); // for plotting and display for (i = 0; i < SOURCESIZE; i++) { eOldSource[i] = 0.0; // initialize the old source values to zero } /* iForLoop */ //*********************************************************************** // Media coefficients //*********************************************************************** for (i = 0; i < media; i++) { temporary = dt * mediaConductivity[i] / (2.0 * electricalPermittivity0 * mediaPermittivity[i] ); // Taflove1995 p.67 mediaCa[i] = (1.0 - temporary) / (1.0 + temporary); // ditto mediaCb[i] = dt / (electricalPermittivity0 * mediaPermittivity[i] * dx * (1.0 + temporary)); // ditto temporary = dt * mediaResistivity[i] / (2.0 * magneticPermeability0 * mediaPermeability[i]); // ditto mediaDa[i] = (1.0 - temporary) / (1.0 + temporary); // ditto mediaDb[i] = dt / (magneticPermeability0 * mediaPermeability[i] * dx * (1.0 + temporary)); // ditto } /* iForLoop */ // for the source (standard fdtd coefficients in vacuum) cbSource = dt / (electricalPermittivity0 * dx); dbSource = dt / (magneticPermeability0 * dx); //*********************************************************************** // Grid Coefficients //*********************************************************************** // Initialize entire grid to free space caex = AllocateMemory(xSize, ySize, mediaCa[0]); // note: don't need to allocate for pec region, as it is not evaluated cbex = AllocateMemory(xSize, ySize, mediaCb[0]); // also: Initialize the entire grid to vacuum. caey = AllocateMemory(xSize, ySize, mediaCa[0] ); cbey = AllocateMemory(xSize, ySize, mediaCb[0] ); dahzx = AllocateMemory(xSize, ySize, mediaDa[0] ); dbhzx = AllocateMemory(xSize, ySize, mediaDb[0] ); dahzy = AllocateMemory(xSize, ySize, mediaDa[0] ); dbhzy = AllocateMemory(xSize, ySize, mediaDb[0] ); // Add metal cylinder #if 1 cylinderDiameter = 20; // diameter of cylinder: 6 cm cylinderRadius = cylinderDiameter / 2.0; // radius of cylinder: 3 cm cylinderXCenter = (4 * xSizeMain) / 5 + abcSize; // i-coordinate of cylinder's center cylinderYCenter = ySize / 2; // j-coordinate of cylinder's center for (i = 0; i < xSize; i++) { for (j = 0; j < ySize; j++) { temporaryi = (double )(i - cylinderXCenter) ; temporaryj = (double )(j - cylinderYCenter) ; distance2 = (temporaryi + 0.5) * (temporaryi + 0.5) + (temporaryj) * (temporaryj); if (distance2 <= (cylinderRadius * cylinderRadius)) { caex[i][j] = mediaCa[1]; cbex[i][j] = mediaCb[1]; } /* if */ // This looks tricky! Why can't caey/cbey use the same 'if' statement as caex/cbex above ?? distance2 = (temporaryj + 0.5) * (temporaryj + 0.5) + (temporaryi) * (temporaryi); if (distance2 <= (cylinderRadius * cylinderRadius)) { caey[i][j] = mediaCa[1]; cbey[i][j] = mediaCb[1]; } /* if */ } /* jForLoop */ } /* iForLoop */ #endif //*********************************************************************** // Fill the PML regions --- (Caution...Here there be Tygers!) //*********************************************************************** // The most important part of the PML fdtd simulation is getting the // PML Coefficients correct. Which requires getting the correct PML gradient and // positioning the coefficients correctly on the x-y grid. // ALERT: It is possible to make a mistake here, and yet the simulation may appear to // be working properly. However a detailed analysis of reflections off the PML // will show they may be (much) larger than those for a correctly designed PML. boundaryWidth = (double )abcSize * dx; // width of PML region (in mm) // SigmaMaximum, using polynomial grading (Nikolova part 4, p.30), rmax=reflectionMax in percent electricalConductivityMaximum = -log(reflectionCoefficient0) * (gradingOrder + 1.0) * electricalPermittivity0 * speedOfLight / (2.0 * boundaryWidth); // boundaryFactor comes from the polynomial grading equation: sigma_x = sigmaxMaximum * (x/d)^m, where d=width of PML, m=gradingOrder, sigmaxMaximum = electricalConductivityMaximum (Nikolova part4, p.28) // IMPORTANT: The conductivity (sigma) must use the "average" value at each mesh point as follows: // sigma_x = sigma_Maximum/dx * Integral_from_x0_to_x1 of (x/d)^m dx, where x0=currentx-0.5, x1=currentx+0.5 (Nikolova part 4, p.32) // integrating gives: sigma_x = (sigmaMaximum / (dx * d^m * m+1)) * ( x1^(m+1) - x0^(m+1) ) (Nikolova part 4, p.32) // the first part is "boundaryFactor", so, sigma_x = boundaryFactor * ( x1^(m+1) - x0^(m+1) ) (Nikolova part 4, p.32) // note: it's not exactly clear what the term eps[0] is for. It's probably to cover the case in which eps[0] is not equal to one (ie the main grid area next to the pml boundary is not vacuum) boundaryFactor = mediaPermittivity[0] * electricalConductivityMaximum / ( dx * (pow(boundaryWidth,gradingOrder)) * (gradingOrder + 1)); // build the gradient // caution: if the gradient is built improperly, the PML will not function correctly for (i = 0, x = 0.0; i < abcSize; i++, x++) { // 0=border between pml and vacuum // even: for ex and ey x1 = (x + 0.5) * dx; // upper bounds for point i x2 = (x - 0.5) * dx; // lower bounds for point i if (i == 0) { gradientConductivity = boundaryFactor * (pow(x1,(gradingOrder+1)) ); // polynomial grading (special case: on the edge, 1/2 = pml, 1/2 = vacuum) } /* if */ else { gradientConductivity = boundaryFactor * (pow(x1,(gradingOrder+1)) - pow(x2,(gradingOrder+1)) ); // polynomial grading } /* else */ gradientCa1[i] = exp (-gradientConductivity * dt / (electricalPermittivity0 * mediaPermittivity[0]) ); // exponential time step, Taflove1995 p.77,78 gradientCb1[i] = (1.0 - gradientCa1[i]) / (gradientConductivity * dx); // ditto, but note sign change from Taflove1995 // odd: for hzx and hzy x1 = (x + 1.0) * dx; // upper bounds for point i x2 = (x + 0.0) * dx; // lower bounds for point i gradientConductivity = boundaryFactor * (pow(x1,(gradingOrder+1)) - pow(x2,(gradingOrder+1)) ); // polynomial grading gradientResistivity = gradientConductivity * (magneticPermeability0 / (electricalPermittivity0 * mediaPermittivity[0]) ); // Taflove1995 p.182 (for no reflection: sigmaM = sigmaE * mu0/eps0) gradientDa1[i] = exp(-gradientResistivity * dt / magneticPermeability0); // exponential time step, Taflove1995 p.77,78 gradientDb1[i] = (1.0 - gradientDa1[i]) / (gradientResistivity * dx); // ditto, but note sign change from Taflove1995 } /* iForLoop */ // ex --- front/back for (j = 0; j < abcSize; j++) { // do coefficients for ex for (i = 0; i < xSize; i++) { // do coefficients for ex for front and back regions caex[i][abcSize - j] = gradientCa1[j]; // front cbex[i][abcSize - j] = gradientCb1[j]; caex[i][ySize - abcSize + j] = gradientCa1[j]; // back cbex[i][ySize - abcSize + j] = gradientCb1[j]; dahzy[i][abcSize - j - 1] = gradientDa1[j]; // front, (note that the index is offset by 1 from caex) dbhzy[i][abcSize - j - 1] = gradientDb1[j]; // ( ditto ) dahzy[i][ySize - abcSize + j] = gradientDa1[j]; // back dbhzy[i][ySize - abcSize + j] = gradientDb1[j]; // } /* iForLoop */ } /* jForLoop */ // ey --- left/right for (i = 0; i < abcSize; i++) { // do coefficients for ey for (j = 0; j < ySize; j++) { // do coefficients for ey for left and right regions caey[abcSize - i][j] = gradientCa1[i]; // left cbey[abcSize - i][j] = gradientCb1[i]; caey[xSize - abcSize + i][j] = gradientCa1[i]; // right cbey[xSize - abcSize + i][j] = gradientCb1[i]; dahzx[abcSize - i - 1][j] = gradientDa1[i]; // left, (note that the index is offset by 1 from caey) dbhzx[abcSize - i - 1][j] = gradientDb1[i]; // (ditto) dahzx[xSize - abcSize + i][j] = gradientDa1[i]; // right dbhzx[xSize - abcSize + i][j] = gradientDb1[i]; // } /* iForLoop */ } /* jForLoop */ // all done with Initialization! //*********************************************************************** // Print variables (diagnostic) //*********************************************************************** printf("# pi:%16.10g\n",pi); printf("# speedOfLight:%16.10g\n",speedOfLight); printf("# magneticPermeability0:%16.10g\n",magneticPermeability0); printf("# electricalPermittivity0:%16.10g\n",electricalPermittivity0); printf("# frequency:%16.10g\n",frequency); printf("# wavelength:%16.10g\n",wavelength); printf("# angularFrequency:%16.10g\n",angularFrequency); printf("# dx:%16.10g\n",dx); printf("# dt:%16.10g\n",dt); printf("# reflectionCoefficient0:%16.10g\n",reflectionCoefficient0); printf("# gradingOrder:%16.10g\n",gradingOrder); printf("# xSizeMain:%d\n",xSizeMain); printf("# ySizeMain:%d\n",ySizeMain); printf("# abcSize:%d\n",abcSize); printf("# xSize:%d\n",xSize); printf("# ySize:%d\n",ySize); printf("# xSource:%d\n",xSource); printf("# ySource:%d\n",ySource); printf("# maximumIteration:%d\n",maximumIteration); #if 0 for (i = 0; i < maximumIteration; i++) { printf("%d: source: %16.10g\n",i,sourceValue[i]); } /* iForLoop */ // print main grid geometry printf("total grid:\n"); for (i = 0; i < xSize; i++) { printf("%3d: ",i); for (j = 0; j < ySize; j++) { ch = '.'; if (caex[i][j] != mediaCa[0]) { ch = 'x'; } /* if */ if (caey[i][j] != mediaCa[0]) { if (ch == '.') { ch = 'y'; } /* if */ else { ch = 'B'; } /* else */ } /* if */ if ((i == xSource) && (j == ySource)) { ch = 'S'; } /* if */ printf("%c",ch); } /* jForLoop */ printf("\n"); } /* iForLoop */ printf("\n"); exit(0); #endif } void EvaluateFdtd (double minimumValue, double maximumValue) { int n,i,j; double temporary; FILE *filePointer ; // for plotting int plottingInterval,iValue; // for plotting double scaleValue; // for plotting char filename[BIGLINESIZE] ; // for plotting int centerx = 50+ABCSIZECONSTANT ; // for printing int centery = 25+ABCSIZECONSTANT ; // "" double eSourceNew[SOURCESIZE]; // hack for source (temporarily holds the new source values) //*********************************************************************** // BEGIN TIME-STEPPING LOOP //*********************************************************************** plottingInterval = 0; for (n = 0; n < maximumIteration; n++) { // iteration loop fprintf(stderr,"n:%d\n",n); //*********************************************************************** // Update electric fields (EX and EY) (including pml) , step 1 //*********************************************************************** // Note: the new value in temporarily saved in the "old value", and updated // in step 3. for (i = 0; i < xSize; i++) { for (j = 1; j < ySize; j++) { // j=0 = pec, so don't evaluate exxaOld[i][j] = (caex[i][j] + dahzx[i][j]) * exxa[i][j] - (caex[i][j] * dahzx[i][j]) * exxaOld[i][j] + (cbex[i][j] * dbhzx[i][j]) * (ey[i][j] - ey[i+1][j]); exxbOld[i][j] = (caex[i][j] + dahzx[i][j-1]) * exxb[i][j] - (caex[i][j] * dahzx[i][j-1]) * exxbOld[i][j] - (cbex[i][j] * dbhzx[i][j-1]) * (ey[i][j-1] - ey[i+1][j-1]); exyaOld[i][j] = (caex[i][j] + dahzy[i][j]) * exya[i][j] - (caex[i][j] * dahzy[i][j]) * exyaOld[i][j] + (cbex[i][j] * dbhzy[i][j]) * (ex[i][j+1] - ex[i][j]); exybOld[i][j] = (caex[i][j] + dahzy[i][j-1]) * exyb[i][j] - (caex[i][j] * dahzy[i][j-1]) * exybOld[i][j] - (cbex[i][j] * dbhzy[i][j-1]) * (ex[i][j] - ex[i][j-1]); } /* jForLoop */ } /* iForLoop */ for (i = 1; i < xSize; i++) { // i=0 = pec, so don't evaluate for (j = 0; j < ySize; j++) { eyxaOld[i][j] = (caey[i][j] + dahzx[i-1][j]) * eyxa[i][j] - (caey[i][j] * dahzx[i-1][j]) * eyxaOld[i][j] + (cbey[i][j] * dbhzx[i-1][j]) * (ey[i-1][j] - ey[i][j]); eyxbOld[i][j] = (caey[i][j] + dahzx[i][j]) * eyxb[i][j] - (caey[i][j] * dahzx[i][j]) * eyxbOld[i][j] - (cbey[i][j] * dbhzx[i][j]) * (ey[i][j] - ey[i+1][j]); eyyaOld[i][j] = (caey[i][j] + dahzy[i-1][j]) * eyya[i][j] - (caey[i][j] * dahzy[i-1][j]) * eyyaOld[i][j] + (cbey[i][j] * dbhzy[i-1][j]) * (ex[i-1][j+1] - ex[i-1][j]); eyybOld[i][j] = (caey[i][j] + dahzy[i][j]) * eyyb[i][j] - (caey[i][j] * dahzy[i][j]) * eyybOld[i][j] - (cbey[i][j] * dbhzy[i][j]) * (ex[i][j+1] - ex[i][j]); } /* jForLoop */ } /* iForLoop */ //*********************************************************************** // Update source electric fields (EX and EY) , step 2 //*********************************************************************** // This is much more complicated than the standard Yee algorithm using Hz // It is necessary to calculate four equations for the four Ex,Ey points // surrounding the Hz Magnetic infinite line source. i = xSource ; j = ySource; temporary = cbSource * (currentSourceValue - oldSourceValue); eSourceNew[0] = (2.0 - cbSource * dbSource) * ex[i][j] - eOldSource[0] + temporary + (cbSource * dbSource) * ( ex[i][j-1] - ey[i][j-1] + ey[i+1][j-1] ); eSourceNew[1] = (2.0 - cbSource * dbSource) * ex[i][j+1] - eOldSource[1] - temporary + (cbSource * dbSource) * ( ex[i][j+2] + ey[i][j+1] - ey[i+1][j+1] ); eSourceNew[2] = (2.0 - cbSource * dbSource) * ey[i][j] - eOldSource[2] - temporary + (cbSource * dbSource) * ( ey[i-1][j] - ex[i-1][j] + ex[i-1][j+1] ); eSourceNew[3] = (2.0 - cbSource * dbSource) * ey[i+1][j] - eOldSource[3] + temporary + (cbSource * dbSource) * ( ey[i+2][j] + ex[i+1][j] - ex[i+1][j+1] ); oldSourceValue = currentSourceValue; currentSourceValue = sourceValue[n]; // update ex/ey Old Source Values (must do at this point, as ex and ey will be overwritten in step 3) eOldSource[0] = ex[i][j] ; eOldSource[1] = ex[i][j+1]; eOldSource[2] = ey[i][j] ; eOldSource[3] = ey[i+1][j]; #if 0 // printout, test of area around source, compare with yee pml simulation printf("ex[i][j] : %16.10g\n", eSourceNew[0]); printf("ex[i][j+1]: %16.10g\n", eSourceNew[1]); printf("ey[i][j] : %16.10g\n", eSourceNew[2]); printf("ey[i+1][j]: %16.10g\n", eSourceNew[3]); printf("temporary: %16.10g\n", temporary); printf("oldSourceValue: %16.10g\n", oldSourceValue); printf("currentSourceValue: %16.10g\n", currentSourceValue); printf("hz[i][j] : %16.10g\n", hz[i][j]); printf("hz[i][j+1]: %16.10g\n", hz[i][j+1]); printf("hz[i+1][j]: %16.10g\n", hz[i+1][j]); printf("hz[i-1][j]: %16.10g\n", hz[i-1][j]); #endif //*********************************************************************** // Update electric fields (EX and EY) , step 3 //*********************************************************************** for (i = 0; i < xSize; i++) { for (j = 1; j < ySize; j++) { // j=0 = pec, so don't evaluate temporary = exxaOld[i][j]; exxaOld[i][j] = exxa[i][j]; exxa[i][j] = temporary; temporary = exxbOld[i][j]; exxbOld[i][j] = exxb[i][j]; exxb[i][j] = temporary; temporary = exyaOld[i][j]; exyaOld[i][j] = exya[i][j]; exya[i][j] = temporary; temporary = exybOld[i][j]; exybOld[i][j] = exyb[i][j]; exyb[i][j] = temporary; ex[i][j] = exxa[i][j] + exxb[i][j] + exya[i][j] + exyb[i][j]; // combine the split fields to form the actual field. } /* jForLoop */ } /* iForLoop */ for (i = 1; i < xSize; i++) { // i=0 = pec, so don't evaluate for (j = 0; j < ySize; j++) { temporary = eyxaOld[i][j]; eyxaOld[i][j] = eyxa[i][j]; eyxa[i][j] = temporary; temporary = eyxbOld[i][j]; eyxbOld[i][j] = eyxb[i][j]; eyxb[i][j] = temporary; temporary = eyyaOld[i][j]; eyyaOld[i][j] = eyya[i][j]; eyya[i][j] = temporary; temporary = eyybOld[i][j]; eyybOld[i][j] = eyyb[i][j]; eyyb[i][j] = temporary; ey[i][j] = eyxa[i][j] + eyxb[i][j] + eyya[i][j] + eyyb[i][j]; // combine the split fields } /* jForLoop */ } /* iForLoop */ // fixup area around source (as it just got overwritten in step 3 above) i = xSource ; j = ySource; ex[i][j] = eSourceNew[0]; ex[i][j+1] = eSourceNew[1]; ey[i][j] = eSourceNew[2]; ey[i+1][j] = eSourceNew[3]; //*********************************************************************** // Plot fields //*********************************************************************** #if 1 // generate Hz for plotting and printing for (j = ABCSIZECONSTANT; j < (ySize-ABCSIZECONSTANT); j++) { for (i = ABCSIZECONSTANT; i < (xSize-ABCSIZECONSTANT); i++) { hz[i][j] += dbSource * ( ex[i][j+1] - ex[i][j] + ey[i][j] - ey[i+1][j] ); // hack for plotting hz } /* xForLoop */ } /* yForLoop */ hz[xSource][ySource] = sourceValue[n]; // hack to Visually match previous simulation programs, such as Yee PML. if (plottingInterval == 0) { plottingInterval = 2; // plot a field (based somewhat on the ToyFdtd BlockOfBricks (bob) output format) // at each plottingInterval, print out a file of field data (field values normalized to a range between 0 and 255). sprintf(filename, "c_%06d.bob", n); filePointer = fopen(filename, "wb") ; if (filePointer == 0) { fprintf(stderr, "Difficulty opening c_%06d.bob", n); exit(1); } /* if */ scaleValue = 256.0 / (maximumValue - minimumValue); for (j = 0; j < ySize; j++) { for (i = 0; i < xSize; i++) { temporary = hz[i][j]; temporary = (temporary - minimumValue) * scaleValue; iValue = (int )(temporary) ; if (iValue < 0) { iValue = 0; } /* if */ if (iValue > 255) { iValue = 255; } /* if */ putc( iValue, filePointer); } /* xForLoop */ } /* yForLoop */ fclose(filePointer); } /* if */ plottingInterval--; #endif //*********************************************************************** // Print fields //*********************************************************************** printf("{ %d\n",n); for (i = 0; i < 8; i++) { // printout some fields printf("eX_%d_%d: %16.10g\n", i,50, ex[ centerx-48 + i][centery]); printf("eY_%d_%d: %16.10g\n", i,50, ey[ centerx-48 + i][centery]); printf("hZ_%d_%d: %16.10g\n", i,50, hz[ centerx-48 + i][centery]); } /* iForLoop */ // printf("eX_%d_%d: %16.10g\n", 50,25, ex[ 50][25]); // printf("eY_%d_%d: %16.10g\n", 50,25, ey[ 50][25]); printf("hZ_%d_%d: %16.10g\n", 50,25, hz[ 50][25]); printf("}\n"); } /* iteration forLoop */ } int main () { InitializeFdtd(); EvaluateFdtd(-0.1,0.1); return(0) ; }