/******************************************************************* Analysis.cpp -- Hertzian dipole with Gaussian pulse source (Courtesy of Maxima!) ********************************************************************/ /* Copyright (C) 2003-2006 Doug Neubauer * * This program is free software; you can redistribute it and/or * modify it under the terms of the GNU General Public License * as published by the Free Software Foundation. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU General Public License * along with this program; if not, write to the Free Software * Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. */ #include #include #include #include const double PIEVALUE = (acos(-1.0)); const double speedOfLight = (2.99792458e8); /* m/second */ const int BIGLINESIZE = (8192); FILE *fpdebug = 0; class World { public: World () { angularFrequency = 0.0 ; electricalPermittivity0 = 0.0 ; dx = 0.0 ; dTime = 0.0 ; waveNumber = 0.0 ; } void World::InitializeWorld ( int numberOfTimeSteps1 ); void World::EvaluateWorld ( ); double World::AnalyticalSolution (double sourceValue, int xIndex,int yIndex,int zIndex, char *s, int nearField); private: int numberOfTimeSteps; double impedance0; double dTime; double angularFrequency ; double electricalPermittivity0; double dx; double waveNumber; /* aka K = angularFrequency/speedOfLight */ }; #if 1 const int MAXIMUMSIZEX = (400); const int MAXIMUMSIZEY = (60); const int MAXIMUMSIZEZ = (60); const double region_dx = (0.000999308); // meters const double region_frequencyMaximum = (1e+010); // hertz // test area: x:0.399723 y:0.0599585 z:0.0599585 meters // wavelength: 0.0299792 // WAVELENGTHPERDX: 30 // speedOfLight: 2.99792e+008 #endif #if 0 // dx = wavelength/20 const int MAXIMUMSIZEX = (240); const int MAXIMUMSIZEY = (40); const int MAXIMUMSIZEZ = (40); const double region_dx = (0.00149896); // meters const double region_frequencyMaximum = (1e+010); // hertz // test area: x:0.359751 y:0.0599585 z:0.0599585 meters #endif void World::InitializeWorld ( int numberOfTimeSteps1 ) { double wavelength ; double magneticPermeability0; double frequency; numberOfTimeSteps = numberOfTimeSteps1; frequency = 10.0e9; /* Hertz */ if (frequency > region_frequencyMaximum) { fprintf(stderr,"error: maximum frequency exceeded: %g > %g\n",frequency,region_frequencyMaximum); exit(1); } /* if */ angularFrequency = 2.0 * PIEVALUE * frequency; /* radians */ wavelength = speedOfLight / frequency; /* meters */ dx = region_dx ; // meters dTime = dx / (speedOfLight * 2.0); /* seconds, courant limit, Taflove page 177 (?? but dx/dt =2C, velocity > C !!!) */ magneticPermeability0 = 4.0e-7 * PIEVALUE ; impedance0 = magneticPermeability0 * speedOfLight ; electricalPermittivity0 = 1.0 / (impedance0 * speedOfLight) ; /* values for analytical solution */ waveNumber = 2.0 * PIEVALUE / wavelength; // k , waveNumber fprintf(fpdebug,"numberOfTimeSteps: %d\n",numberOfTimeSteps); fprintf(fpdebug,"dimensions: x:%d y:%d z:%d (voxels)\n",MAXIMUMSIZEX,MAXIMUMSIZEY,MAXIMUMSIZEZ); fprintf(fpdebug,"frequency: %g (Hz)\n",frequency); fprintf(fpdebug,"magneticPermeability0: %g (Henrys/meter)\n",magneticPermeability0); fprintf(fpdebug,"impedance0: %g (ohms)\n",impedance0); fprintf(fpdebug,"electricalPermittivity0: %g (farads/meter)\n",electricalPermittivity0); fprintf(fpdebug,"wavelength: %g (meters)\n",wavelength); fprintf(fpdebug,"doubleSize:%d intSize:%d char:%d\n",sizeof(double),sizeof(int),sizeof(char)); fprintf(fpdebug,"dx: %g (meters)\n",dx); fprintf(fpdebug,"dTime: %g (seconds)\n",dTime); fprintf(fpdebug,"angularFrequency: %g (radians)\n",angularFrequency); fprintf(fpdebug,"waveNumber: %g\n",waveNumber); fprintf(fpdebug,"speedOfLight: %g (meters/second)\n",speedOfLight); fprintf(fpdebug,"PIEVALUE: %g\n",PIEVALUE); fprintf(fpdebug,"\n"); fflush(fpdebug); } #if 0 // dipole is at {19, 19, 19} const double xSource = (19.0); const double ySource = (19.0); const double zSource = (19.0); const double amplitude = (-1.0e-9 * 1.734); #endif #if 1 // dipole is at {29, 29, 29} const double xSource = (29.0); const double ySource = (29.0); const double zSource = (29.0); // const double amplitude = (-1.0e-9 * 0.487); //sine wave const double amplitude = (1.0e-9 * 1.22); // gp #endif double World::AnalyticalSolution (double currentTime, int xIndex,int yIndex,int zIndex, char *s, int nearField) { double x,y,z,r,theta,phi,c,k,omega ; double retarded, ez; double d = 1.0e-9; // meters double q0 = amplitude ; double t0 = 84.0 * dTime; // was 56 double td = 24.0 * dTime; // was 16 double pi,t,epsilon0; double temporary1,temporary2; // hertzian dipole antenna x = (double )(xIndex) - xSource; y = (double )(yIndex) - ySource; z = (double )(zIndex) - zSource; r = sqrt(x * x + y * y + z * z) ; if (r == 0.0) { return(0.0); } /* if */ if (y == 0.0) { /* fix Microsoft bug */ phi = 0.0; } /* if */ else { phi = atan2(y,x) ; // phi GriffthsEd notation } /* else */ theta = acos(z / r); // theta r *= dx; // get the units right!!!!! meters c = speedOfLight; pi = PIEVALUE; t = currentTime; epsilon0 = electricalPermittivity0; // temporary1 = (t - t0) / td; // Gaussian pulse // temporary2 = 0.1 * exp(-(temporary1 * temporary1)); // printf("gp: %g\n",temporary2); // printf("dgp: %g\n",temporary2 * (-2.0 * temporary1) ); retarded = (currentTime - r / speedOfLight) ; // retarded time calculation // retarded = angularFrequency * (currentTime - r / speedOfLight) ; // retarded time calculation // printf("retarded:%g currentTime:%g r:%g speedOfLight:%g\n",retarded,currentTime,r,speedOfLight); if (retarded < 0.0) { printf("aez_%d_%d_%d: %g\n" ,xIndex,yIndex,zIndex ,0.0); return(0.0); // wave front hasn't reached this point yet. } /* if */ k = waveNumber; omega = angularFrequency; #if 0 // Source: Jz(t)=sine wave, q(t)=cosine wave which should match FDTD Jz(t)=sine wave ez = (d / (4 * pi * pow(c , 2) * epsilon0 * r) * cos( omega * (t - r / c) ) * pow(omega , 2) * q0 + sin( theta ) / (4 * pi * epsilon0 * r) * (d / (2 * c * r) * (1 - d / (2 * r) * cos( theta )) * omega * q0 * sin( theta ) * sin( omega * (t - r / c * (1 + d / (2 * r) * cos( theta )) ) ) + d / (2 * c * r) * (1 + d / (2 * r) * cos( theta )) * omega * q0 * sin( theta ) * sin( omega * (t - r / c * (1 - d / (2 * r) * cos( theta )) ) ) - d / (2 * pow(r , 2) ) * cos( omega * (t - r / c * (1 + d / (2 * r) * cos( theta )) ) ) * q0 * sin( theta ) - d / (2 * pow(r , 2) ) * cos( omega * (t - r / c * (1 - d / (2 * r) * cos( theta )) ) ) * q0 * sin( theta )) - cos( theta ) / (4 * pi * epsilon0) * (omega / r * (1 - d / (2 * r) * cos( theta )) * (d / (2 * c * r) * cos( theta ) - (1 + d / (2 * r) * cos( theta )) / c ) * q0 * sin( omega * (t - r / c * (1 + d / (2 * r) * cos( theta )) ) ) + q0 / pow(r , 2) * (1 - d / (2 * r) * cos( theta )) * cos( omega * (t - r / c * (1 + d / (2 * r) * cos( theta )) ) ) - d / (2 * pow(r , 3) ) * cos( theta ) * cos( omega * (t - r / c * (1 + d / (2 * r) * cos( theta )) ) ) * q0 - omega / r * (1 + d / (2 * r) * cos( theta )) * (-(1 - d / (2 * r) * cos( theta )) / c - d / (2 * c * r) * cos( theta )) * q0 * sin( omega * (t - r / c * (1 - d / (2 * r) * cos( theta )) ) ) - q0 / pow(r , 2) * (1 + d / (2 * r) * cos( theta )) * cos( omega * (t - r / c * (1 - d / (2 * r) * cos( theta )) ) ) - d / (2 * pow(r , 3) ) * cos( theta ) * cos( omega * (t - r / c * (1 - d / (2 * r) * cos( theta )) ) ) * q0) ) ; /* 510 */ #endif #if 0 // sine wave from Maxima ez = (d / (4 * pi * pow(c , 2) * epsilon0 * r) * pow(omega , 2) * q0 * sin( omega * (t - r / c) ) + sin( theta ) / (4 * pi * epsilon0 * r) * (-d / (2 * pow(r , 2) ) * q0 * sin( theta ) * sin( omega * (t - r / c * (1 + d / (2 * r) * cos( theta )) ) ) - d / (2 * c * r) * (1 - d / (2 * r) * cos( theta )) * cos( omega * (t - r / c * (1 + d / (2 * r) * cos( theta )) ) ) * omega * q0 * sin( theta ) - d / (2 * pow(r , 2) ) * q0 * sin( theta ) * sin( omega * (t - r / c * (1 - d / (2 * r) * cos( theta )) ) ) - d / (2 * c * r) * (1 + d / (2 * r) * cos( theta )) * cos( omega * (t - r / c * (1 - d / (2 * r) * cos( theta )) ) ) * omega * q0 * sin( theta )) - cos( theta ) / (4 * pi * epsilon0) * (omega / r * (1 + d / (2 * r) * cos( theta )) * (-(1 - d / (2 * r) * cos( theta )) / c - d / (2 * c * r) * cos( theta )) * cos( omega * (t - r / c * (1 - d / (2 * r) * cos( theta )) ) ) * q0 + q0 / pow(r , 2) * (1 - d / (2 * r) * cos( theta )) * sin( omega * (t - r / c * (1 + d / (2 * r) * cos( theta )) ) ) - d / (2 * pow(r , 3) ) * cos( theta ) * q0 * sin( omega * (t - r / c * (1 + d / (2 * r) * cos( theta )) ) ) - omega / r * (1 - d / (2 * r) * cos( theta )) * cos( omega * (t - r / c * (1 + d / (2 * r) * cos( theta )) ) ) * (d / (2 * c * r) * cos( theta ) - (1 + d / (2 * r) * cos( theta )) / c ) * q0 - q0 / pow(r , 2) * (1 + d / (2 * r) * cos( theta )) * sin( omega * (t - r / c * (1 - d / (2 * r) * cos( theta )) ) ) - d / (2 * pow(r , 3) ) * cos( theta ) * q0 * sin( omega * (t - r / c * (1 - d / (2 * r) * cos( theta )) ) )) ) ; /* 511 */ #endif #if 1 // gaussian pulse ez = (q0 * sin( theta ) * (d * (1 - d * cos( theta ) / (2 * r) ) * (-r * (d * cos( theta ) / (2 * r) + 1) / c - t0 + t) * exp((-(pow((-r * (d * cos( theta ) / (2 * r) + 1) / c - t0 + t) , (2) ) / pow(td , (2) ) ) ) ) * sin( theta ) / (c * r * pow(td , (2) ) ) - d * exp((-(pow((-r * (d * cos( theta ) / (2 * r) + 1) / c - t0 + t) , (2) ) / pow(td , (2) ) ) ) ) * sin( theta ) / (2 * pow(r , (2) ) ) + d * (d * cos( theta ) / (2 * r) + 1) * (-r * (1 - d * cos( theta ) / (2 * r) ) / c - t0 + t) * exp((-(pow((-r * (1 - d * cos( theta ) / (2 * r) ) / c - t0 + t) , (2) ) / pow(td , (2) ) ) ) ) * sin( theta ) / (c * r * pow(td , (2) ) ) - d * exp((-(pow((-r * (1 - d * cos( theta ) / (2 * r) ) / c - t0 + t) , (2) ) / pow(td , (2) ) ) ) ) * sin( theta ) / (2 * pow(r , (2) ) ) ) / (4 * pi * epsilon0 * r) - q0 * cos( theta ) * (2 * (1 - d * cos( theta ) / (2 * r) ) * (d * cos( theta ) / (2 * c * r) - (d * cos( theta ) / (2 * r) + 1) / c ) * (-r * (d * cos( theta ) / (2 * r) + 1) / c - t0 + t) * exp((-(pow((-r * (d * cos( theta ) / (2 * r) + 1) / c - t0 + t) , (2) ) / pow(td , (2) ) ) ) ) / (r * pow(td , (2) ) ) + (1 - d * cos( theta ) / (2 * r) ) * exp((-(pow((-r * (d * cos( theta ) / (2 * r) + 1) / c - t0 + t) , (2) ) / pow(td , (2) ) ) ) ) / pow(r , (2) ) - d * cos( theta ) * exp((-(pow((-r * (d * cos( theta ) / (2 * r) + 1) / c - t0 + t) , (2) ) / pow(td , (2) ) ) ) ) / (2 * pow(r , (3) ) ) - 2 * (d * cos( theta ) / (2 * r) + 1) * (-(1 - d * cos( theta ) / (2 * r) ) / c - d * cos( theta ) / (2 * c * r) ) * (-r * (1 - d * cos( theta ) / (2 * r) ) / c - t0 + t) * exp((-(pow((-r * (1 - d * cos( theta ) / (2 * r) ) / c - t0 + t) , (2) ) / pow(td , (2) ) ) ) ) / (r * pow(td , (2) ) ) - (d * cos( theta ) / (2 * r) + 1) * exp((-(pow((-r * (1 - d * cos( theta ) / (2 * r) ) / c - t0 + t) , (2) ) / pow(td , (2) ) ) ) ) / pow(r , (2) ) - d * cos( theta ) * exp((-(pow((-r * (1 - d * cos( theta ) / (2 * r) ) / c - t0 + t) , (2) ) / pow(td , (2) ) ) ) ) / (2 * pow(r , (3) ) ) ) / (4 * pi * epsilon0) + d * q0 * exp((-(pow((-t0 + t - r / c ) , (2) ) / pow(td , (2) ) ) ) ) / (2 * pi * pow(c , (2) ) * epsilon0 * r * pow(td , (2) ) ) - d * q0 * pow((-t0 + t - r / c ) , (2) ) * exp((-(pow((-t0 + t - r / c ) , (2) ) / pow(td , (2) ) ) ) ) / (pi * pow(c , (2) ) * epsilon0 * r * pow(td , (4) ) ) ) ; /* 0 */ #endif // printf("; (x:%g y:%g z:%g voxels) r:%g(meters) phi:%g theta:%g\n",x,y,z,r,phi,theta); printf("aez_%d_%d_%d: %g\n" ,xIndex,yIndex,zIndex , ez); return(ez) ; } char *printIndicators = "*OXABCDEFGHJKLMNPSTUVWYZ123456789abcdefghjkmnprstuvwyz"; const int NUMBEROFPRINTPOSITIONS = (18); // 22 typedef struct _PrintPosition{ int x,y,z ; double maximumex; double maximumey; double maximumez; } PrintPosition ; PrintPosition printPositions[NUMBEROFPRINTPOSITIONS] = { {30,29,29, -1.0e100, -1.0e100, -1.0e100}, {50,29,29, -1.0e100, -1.0e100, -1.0e100}, {70,29,29, -1.0e100, -1.0e100, -1.0e100}, {90,29,29, -1.0e100, -1.0e100, -1.0e100}, {110,29,29, -1.0e100, -1.0e100, -1.0e100}, {130,29,29, -1.0e100, -1.0e100, -1.0e100}, {150,29,29, -1.0e100, -1.0e100, -1.0e100}, {170,29,29, -1.0e100, -1.0e100, -1.0e100}, {190,29,29, -1.0e100, -1.0e100, -1.0e100}, {210,29,29, -1.0e100, -1.0e100, -1.0e100}, {230,29,29, -1.0e100, -1.0e100, -1.0e100}, {250,29,29, -1.0e100, -1.0e100, -1.0e100}, {270,29,29, -1.0e100, -1.0e100, -1.0e100}, {290,29,29, -1.0e100, -1.0e100, -1.0e100}, {310,29,29, -1.0e100, -1.0e100, -1.0e100}, {330,29,29, -1.0e100, -1.0e100, -1.0e100}, {350,29,29, -1.0e100, -1.0e100, -1.0e100}, {370,29,29, -1.0e100, -1.0e100, -1.0e100}, #if 0 {30,19,19, -1.0e100, -1.0e100, -1.0e100}, {50,19,19, -1.0e100, -1.0e100, -1.0e100}, {70,19,19, -1.0e100, -1.0e100, -1.0e100}, {90,19,19, -1.0e100, -1.0e100, -1.0e100}, {110,19,19, -1.0e100, -1.0e100, -1.0e100}, {130,19,19, -1.0e100, -1.0e100, -1.0e100}, {150,19,19, -1.0e100, -1.0e100, -1.0e100}, {170,19,19, -1.0e100, -1.0e100, -1.0e100}, {190,19,19, -1.0e100, -1.0e100, -1.0e100}, {210,19,19, -1.0e100, -1.0e100, -1.0e100}, {230,19,19, -1.0e100, -1.0e100, -1.0e100}, {250,19,19, -1.0e100, -1.0e100, -1.0e100}, {270,19,19, -1.0e100, -1.0e100, -1.0e100}, {290,19,19, -1.0e100, -1.0e100, -1.0e100}, {310,19,19, -1.0e100, -1.0e100, -1.0e100}, {330,19,19, -1.0e100, -1.0e100, -1.0e100}, {350,19,19, -1.0e100, -1.0e100, -1.0e100}, {370,19,19, -1.0e100, -1.0e100, -1.0e100}, {19,19,19, -1.0e100, -1.0e100, -1.0e100}, {21,19,19, -1.0e100, -1.0e100, -1.0e100}, {25,19,19, -1.0e100, -1.0e100, -1.0e100}, {35,19,19, -1.0e100, -1.0e100, -1.0e100}, {55,19,19, -1.0e100, -1.0e100, -1.0e100}, {19,21,19, -1.0e100, -1.0e100, -1.0e100}, {19,25,19, -1.0e100, -1.0e100, -1.0e100}, {19,35,19, -1.0e100, -1.0e100, -1.0e100}, {21,21,21, -1.0e100, -1.0e100, -1.0e100}, {25,25,25, -1.0e100, -1.0e100, -1.0e100}, {35,35,35, -1.0e100, -1.0e100, -1.0e100}, {25,15,30, -1.0e100, -1.0e100, -1.0e100}, {15,13,19, -1.0e100, -1.0e100, -1.0e100}, {29,27,13, -1.0e100, -1.0e100, -1.0e100}, {29,13,13, -1.0e100, -1.0e100, -1.0e100}, {29,13,25, -1.0e100, -1.0e100, -1.0e100}, {29,25,13, -1.0e100, -1.0e100, -1.0e100}, {29,25,25, -1.0e100, -1.0e100, -1.0e100}, { 9,13,13, -1.0e100, -1.0e100, -1.0e100}, { 9,13,25, -1.0e100, -1.0e100, -1.0e100}, { 9,25,13, -1.0e100, -1.0e100, -1.0e100}, { 9,25,25, -1.0e100, -1.0e100, -1.0e100}, #endif }; void PrintHeader () { int i ; for (i = 0; i < NUMBEROFPRINTPOSITIONS; i++) { printf(";data: \"aez_%d_%d_%d\" \"aez(%d,%d,%d)\" %c %g %g\n", printPositions[i].x,printPositions[i].y,printPositions[i].z, printPositions[i].x,printPositions[i].y,printPositions[i].z, printIndicators[i], -printPositions[i].maximumez,printPositions[i].maximumez); } /* iForLoop */ } void World::EvaluateWorld ( ) { double t,ez; int i,j; int xIndex,yIndex,zIndex; for (t = 0.0 * dTime,i = 0; i < numberOfTimeSteps; t += dTime,i++) { // fprintf(stderr,"%d\n",i); fprintf(fpdebug,"\n---------------------------\nTime:%g:\n",t); printf("\n{ %d\n", i ); for (j = 0; j < NUMBEROFPRINTPOSITIONS; j++) { xIndex = printPositions[j].x; yIndex = printPositions[j].y; zIndex = printPositions[j].z; ez = AnalyticalSolution(t,xIndex,yIndex,zIndex, "??", 0); ez = fabs(ez); // printf("; ez:%g\n",ez); if (ez > printPositions[j].maximumez) { printPositions[j].maximumez = ez; } /* if */ } /* iForLoop */ printf("}\n"); } /* forLoop */ PrintHeader(); } void main () { World testFixture ; int numberOfTimeSteps = 4000; fpdebug = stdout ; testFixture.InitializeWorld( numberOfTimeSteps ) ; testFixture.EvaluateWorld( ); fclose(fpdebug); }