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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#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 <stdio.h>
#include <stdlib.h>
#include <math.h>
#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 <stdio.h>
#include <math.h>
#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