TOP OF MAIN.c /////////////////////////////////////////////////////////////////////////////////////// /// AUTHOR: Bill Klemens /// (Some of Code is Modified From Code Distributed by Professor J.B. Schneider) /// DATE: January 28th, 2007 /// BRIEF: 1)Simulation of pulse normally incident on water. /// 2)FDTD Implementation of 1D Grid w/ additive Ricker wavelet w/120ppw. /// 3)Debye Material and PLRC Implementation are used for H20. /// /// PURPOSE: This is the ninth and final homework assignment for EE535, /// Fall semester of 2005. /// /// /// REVISION HISTORY: /// Date Completed /// /////////////////////////////////////////////////////////////////////////////////////// // SYSTEM INCLUDES: // #include #include #include #include "Extra.h" #include "cpml_tez.h" Snip Snip ***************** Extra.C //////////////////////////////////////////////////////////////////////////////// /// \author Bill Klemens /// \date September 21st, 2005 /// \brief FDTD Implementation of TMz grid w/ Ricker Excitation /// /// PURPOSE: A longer description and purpose for this file. /// /// /// REVISION HISTORY: /// date Comment on what was changed. /// /////////////////////////////////////////////////////////////////////////////// #include "Extra.h" #include #include #include #define M_PI 3.141592653589793238462643383 static double cdtds_Extra; static double ppw_Extra; // Functions /////////////////////////////////////////////////////////////////////// /// Initialize Coefficients /// \pre /// \post /// \param Time;Courant 'Sc'; Np, points per wavelength /// \return /// \throw (omit if the function throws no exceptions, or repeat for every exception) /////////////////////////////////////////////////////////////////////// void init_Extra(double the_cdtds, double the_ppw) { cdtds_Extra = the_cdtds; ppw_Extra = the_ppw; } /////////////////////////////////////////////////////////////////////// /// Round: Returns Number Rounded up or Down. /// \pre /// \post /// \param /// \return /// \throw (omit if the function throws no exceptions, or repeat for every exception) /////////////////////////////////////////////////////////////////////// int Round ( double the_temp ) { double temp, convert; int final, hold; temp = the_temp; hold = (int)temp; convert = (double)hold + 0.5; if( temp >= convert ) { final = hold + 1; } else { final = hold; } return final; } /////////////////////////////////////////////////////////////////////// /// Ricker Wavelet. /// \pre /// \post /// \param Time;Courant 'Sc'; Np, points per wavelength /// \return /// \throw (omit if the function throws no exceptions, or repeat for every exception) /////////////////////////////////////////////////////////////////////// double ricker_ADD(double time) { double delay = 1.0; // delay = multiple of inverse of most // energetic frequency. double pi = 3.141592;//Math::PI; double arg; double arg_max = 70.0; // Arguments beyond this are assumed to // yield zero. This allows us to avoid // calling the exponential when result // will be effectively zero anyway. arg = pow(pi *( (cdtds_Extra * time)/ppw_Extra - delay),2); if (arg > arg_max) { return 0.0; } else { return (1.0 - 2.0*arg) * exp(-arg); } } double quasi(double time) { double func; double temp; double pi = 3.14159265359; temp = (2 * pi * time) / 30.0; //temp = (2 * pi ) / 30.0; func = (1-exp(-time/20)) * sin(temp); return func; } /////////////////////////////////////////////////////////////////////// /// WRTRAW: Used to write raw data. /// \pre /// \post /// \param Time;Courant 'Sc'; Np, points per wavelength /// \return /// \throw (omit if the function throws no exceptions, or repeat for every exception) /////////////////////////////////////////////////////////////////////// void wrtraw(double *field, int size_x, int size_y, int frame, char *basename) { int n , m; float dim1 = (float)size_x; float dim2 = (float)size_y; float tmp; char filename[100]; FILE *out; /* construct output file name and open it */ sprintf(filename,"%s.%d",basename,frame); out = fopen(filename,"wb"); /* write dimensions to output file */ fwrite(&dim1,sizeof(float),1,out); fwrite(&dim2,sizeof(float),1,out); /* write remaining data */ for (n = size_y-1; n >= 0; n--) { for (m = 0; m < size_x; m++) { tmp = (float)field[ m * size_y + n ]; // store data in a float fwrite( &tmp, sizeof(float), 1, out ); // write the float } } fclose(out); // close file return; } /////////////////////////////////////////////////////////////////////// /// WRTRAW: Used to write raw data. /// \pre /// \post /// \param Time;Courant 'Sc'; Np, points per wavelength /// \return /// \throw (omit if the function throws no exceptions, or repeat for every exception) /////////////////////////////////////////////////////////////////////// void wrtraw2(double *field, int size_x, int size_y) { int n , m; float dim1 = (float)size_x; float dim2 = (float)size_y; float tmp; char ez120[20] = "EZ_120.0"; FILE *out; /* construct output file name and open it */ out = fopen(ez120,"wb"); /* write dimensions to output file */ fwrite(&dim1,sizeof(float),1,out); fwrite(&dim2,sizeof(float),1,out); /* write remaining data */ for (n = size_y-1; n >= 0; n--) { for (m = 0; m < size_x; m++) { tmp = (float)field[ m * size_y + n ]; // store data in a float fwrite( &tmp, sizeof(float), 1, out ); // write the float } } fclose(out); // close file return; } ********************TOP OF CPML-TEz.h #include #include #include "cpml_tez.h" #include "Extra.h" // contains ALLOC_2D /* macros for accessing array elements */ #define Hz(I,J) hz[(I)*(size_y-1) + (J)] #define Ex(I,J) ex[(I)*size_y + (J)] #define Ey(I,J) ey[(I)*(size_y-1) + (J)] /* macros for accessing PML array */ #define Psi_hzy0(I,J) psi_hzy0[(I)*npml + (J)] #define Psi_hzy1(I,J) psi_hzy1[(I)*npml + (J-(size_y-npml-1))] #define Psi_hzx0(I,J) psi_hzx0[(I)*(size_y-1) + (J)] #define Psi_hzx1(I,J) psi_hzx1[(I-(size_x-npml-1))*(size_y-1) + (J)] #define Psi_eyx0(I,J) psi_eyx0[(I)*(size_y-1) + (J)] #define Psi_eyx1(I,J) psi_eyx1[(I-(size_x-npml))*(size_y-1) + (J)] #define Psi_exy0(I,J) psi_exy0[(I)*npml + (J)] #define Psi_exy1(I,J) psi_exy1[(I)*npml + (J-(size_y-npml))] /************** global variables ***************/ static double *kappa_ex, *kappa_hx, *kappa_ey, *kappa_hy, *be_x, *bh_x, *ce_x, *ch_x, *be_y, *bh_y, *ce_y, *ch_y, *psi_exy0, *psi_exy1, *psi_eyx0, *psi_eyx1, SNIP SNIP