cpmlfdtd3d.c

FDTD

	/*******************************************************************
	//     3-D FDTD code with CPML absorbing boundary conditions - Dipole
	//******************************************************************
	//
	//     This C code implements the finite-difference time-domain
	//     solution of Maxwell's curl equations over a three-dimensional
	//     Cartesian space lattice.  The grid is terminated by CPML ABCs.
	//     The thickness of the PML in each Cartesian direction can be 
	//     varied independently.
	//
	//******************************************************************/
	//Header files (Libraries to be included)
	#include<stdio.h>
	#include<math.h>
	#include<string.h>
	#include<stdlib.h>

	//  Fundamental Constants (MKS units)
    double pi = 3.14159265358979;
    double C = 2.99792458E8;
    double muO;
    double epsO;

	//  Specify Material Relative Permittivity and Conductivity
    double epsR = 1.0;//free space

	//  Specify Number of Time Steps and Grid Size Parameters
    int nMax = 500;  // total number of time steps
    
	// grid size corresponding to the number of Ez field components
    int Imax = 51;
    int Jmax = 126;
    int Kmax = 26;

	//  Specify Grid Cell Size in Each Direction and Calculate the
  	//  Resulting Courant-Stable Time Step
    double dx = 1.0E-3;
    double dy = 1.0E-3;
    double dz = 1.0E-3; // cell size in each direction
    // time step increment
    double dt;

	//  Specify the Impulsive Source (Differentiated Gaussian) parameters
    double tw = 53.0E-12;//pulse width
    double tO;//delay
    double source;//Differentiated Gaussian source
	double amp = 1000;// Amplitude

    //Specify the Time Step at which the data has to be saved for Visualization
    int save_modulus = 10;

	//  Specify the dipole Boundaries(A cuboidal rode- NOT as a cylinder)
    int istart, iend, jstart;
    int  jend, kstart, kend;

	//Output recording point
    int ksource = 12;

	//  Specify the CPML Thickness in Each Direction (Value of Zero
	//  Corresponds to No PML, and the Grid is Terminated with a PEC)
    // PML thickness in each direction
    int nxPML_1, nxPML_2, nyPML_1;
    int nyPML_2, nzPML_1, nzPML_2;

	//  Specify the CPML Order and Other Parameters:
    int m = 3, ma = 1;

    double sig_x_max;
    double sig_y_max;
    double sig_z_max;
    double alpha_x_max;
    double alpha_y_max;
    double alpha_z_max;
    double kappa_x_max;
    double kappa_y_max;
    double kappa_z_max;
	
	//Loop indices
    int i,j,ii,jj,k,kk,n;


	// H & E Field components
    float ***Hx;
    float ***Hy;
    float ***Hz;
    float ***Ex;
    float ***Ey;
    float ***Ez;

	short ***ID1;//medium definition array for Ex
	short ***ID2;//medium definition array for Ey
	short ***ID3;//medium definition array for Ez

	//  CPML components (Taflove 3rd Edition, Chapter 7)
    float ***psi_Ezx_1;
    float ***psi_Ezx_2;
    float ***psi_Hyx_1;
    float ***psi_Hyx_2;
    float ***psi_Ezy_1;
    float ***psi_Ezy_2;
    float ***psi_Hxy_1;
    float ***psi_Hxy_2;
    float ***psi_Hxz_1;
    float ***psi_Hxz_2;
    float ***psi_Hyz_1;
    float ***psi_Hyz_2;
    float ***psi_Exz_1;
    float ***psi_Exz_2;
    float ***psi_Eyz_1;
    float ***psi_Eyz_2;
    float ***psi_Hzx_1;
    float ***psi_Eyx_1;
    float ***psi_Hzx_2;
    float ***psi_Eyx_2;
    float ***psi_Hzy_1;
    float ***psi_Exy_1;
    float ***psi_Hzy_2;
    float ***psi_Exy_2;

    float *be_x_1, *ce_x_1, *alphae_x_PML_1, *sige_x_PML_1, *kappae_x_PML_1;
    float *bh_x_1, *ch_x_1, *alphah_x_PML_1, *sigh_x_PML_1, *kappah_x_PML_1;
    float *be_x_2, *ce_x_2, *alphae_x_PML_2, *sige_x_PML_2, *kappae_x_PML_2;
    float *bh_x_2, *ch_x_2, *alphah_x_PML_2, *sigh_x_PML_2, *kappah_x_PML_2;
    float *be_y_1, *ce_y_1, *alphae_y_PML_1, *sige_y_PML_1, *kappae_y_PML_1;
    float *bh_y_1, *ch_y_1, *alphah_y_PML_1, *sigh_y_PML_1, *kappah_y_PML_1;
    float *be_y_2, *ce_y_2, *alphae_y_PML_2, *sige_y_PML_2, *kappae_y_PML_2;
    float *bh_y_2, *ch_y_2, *alphah_y_PML_2, *sigh_y_PML_2, *kappah_y_PML_2;
    float *be_z_1, *ce_z_1, *alphae_z_PML_1, *sige_z_PML_1, *kappae_z_PML_1;
    float *bh_z_1, *ch_z_1, *alphah_z_PML_1, *sigh_z_PML_1, *kappah_z_PML_1;
    float *be_z_2, *ce_z_2, *alphae_z_PML_2, *sige_z_PML_2, *kappae_z_PML_2;
    float *bh_z_2, *ch_z_2, *alphah_z_PML_2, *sigh_z_PML_2, *kappah_z_PML_2;

	// denominators for the update equations
    float *den_ex;
    float *den_hx;
    float *den_ey;
    float *den_hy;
    float *den_ez;
    float *den_hz;

	//Max number of materials allowed
	int numMaterials = 50;

    //permittivity, permeability and conductivity of diffrent materials
	double *epsilon;
	double *mu;
	double *sigma;

	//E field update coefficients
	float *CA;
	float *CB;

	//H field update coefficients
    float DA;
    float DB;
	
	//Function prototype definitions
	void initialize();//Memeory initialization
	void setUp();//Coefficients, parameters etc will get computed
	void initializeCPML();//CPML coefficient computation
	void compute();//E & H Field update equation
	void buildObject();//Creates the object geometry
	void yeeCube (int, int,int, short);//Sets material properties to a cell
	void writeField(int);//Writes output
	void buildSphere();//Builds a spherical object
	void buildDipole();//Builds a dipole

	//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	void main() {

		initialize();
		setUp();
		buildObject();
		initializeCPML();
		compute();

	}

	//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	void initialize() {

		muO = 4.0 * pi * 1.0E-7;
		epsO = 1.0 / (C * C * muO);
		 
		//PML Layers (10 layers)
		nxPML_1 = 11;
		nxPML_2 = 11;
		nyPML_1 = 11;
		nyPML_2 = 11;
		nzPML_1 = 11;
		nzPML_2 = 11;

		//Dynamic memory allocation
	    epsilon = (double *)malloc((numMaterials) * sizeof(double));

		for(i = 0; i < numMaterials; i++) {

			epsilon[i] = epsO;
		}

	    mu = (double *)malloc((numMaterials) * sizeof(double));

		for(i = 0; i < numMaterials; i++) {

			mu[i] = muO;
		}

	    sigma = (double *)malloc((numMaterials) * sizeof(double));

		for(i = 0; i < numMaterials; i++) {

			sigma[i] = 0.0;
		}

	    CA = (float *)malloc((numMaterials) * sizeof(float));

		for(i = 0; i < numMaterials; i++) {

			CA[i] = 0.0;
		}

	    CB = (float *)malloc((numMaterials) * sizeof(float));

		for(i = 0; i < numMaterials; i++) {

			CB[i] = 0.0;
		}

	    Ez = (float ***)malloc(Imax * sizeof(float **));

	    for(i = 0; i < Imax; i++) {

		    Ez[i] = (float **)malloc(Jmax * sizeof(float *));

		    for(j = 0; j < Jmax; j++) {

				Ez[i][j] = (float *)malloc(Kmax * sizeof(float));

				for(k = 0; k < Kmax; k++) {

					Ez[i][j][k] = 0.0;
				}
			}
		}

	    Ey = (float ***)malloc((Imax) * sizeof(float **));

	    for(i = 0; i < Imax; i++) {

		    Ey[i] = (float **)malloc((Jmax-1) * sizeof(float *));

		    for(j = 0; j < Jmax -1; j++) {

				Ey[i][j] = (float *)malloc((Kmax-1) * sizeof(float));

				for(k = 0; k < Kmax-1; k++) {

					Ey[i][j][k] = 0.0;
				}
			}
		}

	    Ex = (float ***)malloc((Imax-1) * sizeof(float **));

	    for(i = 0; i < Imax-1; i++) {

		    Ex[i] = (float **)malloc(Jmax * sizeof(float *));

		    for(j = 0; j < Jmax; j++) {

				Ex[i][j] = (float *)malloc((Kmax-1) * sizeof(float));

				for(k = 0; k < Kmax-1; k++) {

					Ex[i][j][k] = 0.0;
				}
			}
		}

	    Hx = (float ***)malloc(Imax * sizeof(float **));

	    for(i = 0; i < Imax; i++) {

		    Hx[i] = (float **)malloc((Jmax-1) * sizeof(float *));

		    for(j = 0; j < Jmax-1; j++) {

				Hx[i][j] = (float *)malloc(Kmax * sizeof(float));

				for(k = 0; k < Kmax; k++) {

					Hx[i][j][k] = 0.0;
				}
			}
		}

	    Hy = (float ***)malloc((Imax-1) * sizeof(float **));

	    for(i = 0; i < Imax-1; i++) {

		    Hy[i] = (float **)malloc(Jmax * sizeof(float *));

		    for(j = 0; j < Jmax; j++) {

				Hy[i][j] = (float *)malloc(Kmax * sizeof(float));

				for(k = 0; k < Kmax; k++) {

					Hy[i][j][k] = 0.0;
				}
			}
		}

	    Hz = (float ***)malloc((Imax-1) * sizeof(float **));

	    for(i = 0; i < Imax-1; i++) {

		    Hz[i] = (float **)malloc((Jmax-1) * sizeof(float *));

		    for(j = 0; j < Jmax-1; j++) {

				Hz[i][j] = (float *)malloc((Kmax-1) * sizeof(float));

				for(k = 0; k < Kmax-1; k++) {

					Hz[i][j][k] = 0.0;
				}
			}
		}

	    ID1 = (short ***)malloc(Imax * sizeof(short **));

	    for(i = 0; i < Imax; i++) {

		    ID1[i] = (short **)malloc(Jmax * sizeof(short *));

		    for(j = 0; j < Jmax; j++) {

				ID1[i][j] = (short *)malloc(Kmax * sizeof(short));

				for(k = 0; k < Kmax; k++) {

					ID1[i][j][k] = 0;
				}
			}
		}

	    ID2 = (short ***)malloc(Imax * sizeof(short **));

	    for(i = 0; i < Imax; i++) {

		    ID2[i] = (short **)malloc(Jmax * sizeof(short *));

		    for(j = 0; j < Jmax; j++) {

				ID2[i][j] = (short *)malloc(Kmax * sizeof(short));

				for(k = 0; k < Kmax; k++) {

					ID2[i][j][k] = 0;
				}
			}
		}

	    ID3 = (short ***)malloc(Imax * sizeof(short **));

	    for(i = 0; i < Imax; i++) {

		    ID3[i] = (short **)malloc(Jmax * sizeof(short *));

		    for(j = 0; j < Jmax; j++) {

				ID3[i][j] = (short *)malloc(Kmax * sizeof(short));

				for(k = 0; k < Kmax; k++) {

					ID3[i][j][k] = 0;
				}
			}
		}

	    psi_Ezx_1 = (float ***)malloc(nxPML_1 * sizeof(float **));

	    for(i = 0; i < nxPML_1; i++) {

		    psi_Ezx_1[i] = (float **)malloc(Jmax * sizeof(float *));

		    for(j = 0; j < Jmax; j++) {

				psi_Ezx_1[i][j] = (float *)malloc(Kmax * sizeof(float));

				for(k = 0; k < Kmax; k++) {

					psi_Ezx_1[i][j][k] = 0.0;
				}
			}
		}

	    psi_Ezx_2 = (float ***)malloc(nxPML_2 * sizeof(float **));

	    for(i = 0; i < nxPML_2; i++) {

		    psi_Ezx_2[i] = (float **)malloc(Jmax * sizeof(float *));

		    for(j = 0; j < Jmax; j++) {

				psi_Ezx_2[i][j] = (float *)malloc(Kmax * sizeof(float));

				for(k = 0; k < Kmax; k++) {

					psi_Ezx_2[i][j][k] = 0.0;
				}
			}
		}

	    psi_Hyx_1 = (float ***)malloc((nxPML_1-1) * sizeof(float **));

	    for(i = 0; i < nxPML_1-1; i++) {

		    psi_Hyx_1[i] = (float **)malloc(Jmax * sizeof(float *));

		    for(j = 0; j < Jmax; j++) {

				psi_Hyx_1[i][j] = (float *)malloc(Kmax * sizeof(float));

				for(k = 0; k < Kmax; k++) {

					psi_Hyx_1[i][j][k] = 0.0;
				}
			}
		}

	    psi_Hyx_2 = (float ***)malloc((nxPML_2-1) * sizeof(float **));

	    for(i = 0; i < nxPML_1-1; i++) {

		    psi_Hyx_2[i] = (float **)malloc(Jmax * sizeof(float *));

		    for(j = 0; j < Jmax; j++) {

				psi_Hyx_2[i][j] = (float *)malloc(Kmax * sizeof(float));

				for(k = 0; k < Kmax; k++) {

					psi_Hyx_2[i][j][k] = 0.0;
				}
			}
		}

	    psi_Ezy_1 = (float ***)malloc(Imax * sizeof(float **));

	    for(i = 0; i < Imax; i++) {

		    psi_Ezy_1[i] = (float **)malloc(nyPML_1 * sizeof(float *));

		    for(j = 0; j < nyPML_1; j++) {

				psi_Ezy_1[i][j] = (float *)malloc(Kmax * sizeof(float));

				for(k = 0; k < Kmax; k++) {

					psi_Ezy_1[i][j][k] = 0.0;
				}
			}
		}

	    psi_Ezy_2 = (float ***)malloc(Imax * sizeof(float **));

	    for(i = 0; i < Imax; i++) {

		    psi_Ezy_2[i] = (float **)malloc(nyPML_2 * sizeof(float *));

		    for(j = 0; j < nyPML_2; j++) {

				psi_Ezy_2[i][j] = (float *)malloc(Kmax * sizeof(float));

				for(k = 0; k < Kmax; k++) {

					psi_Ezy_2[i][j][k] = 0.0;
				}
			}
		}

	    psi_Hxy_1 = (float ***)malloc(Imax * sizeof(float **));

	    for(i = 0; i < Imax; i++) {

		    psi_Hxy_1[i] = (float **)malloc((nyPML_1-1) * sizeof(float *));

		    for(j = 0; j < nyPML_1-1; j++) {

				psi_Hxy_1[i][j] = (float *)malloc(Kmax * sizeof(float));

				for(k = 0; k < Kmax; k++) {

					psi_Hxy_1[i][j][k] = 0.0;
				}
			}
		}

	    psi_Hxy_2 = (float ***)malloc(Imax * sizeof(float **));

	    for(i = 0; i < Imax; i++) {

		    psi_Hxy_2[i] = (float **)malloc((nyPML_2-1) * sizeof(float *));

		    for(j = 0; j < nyPML_2-1; j++) {

				psi_Hxy_2[i][j] = (float *)malloc(Kmax * sizeof(float));

				for(k = 0; k < Kmax; k++) {

					psi_Hxy_2[i][j][k] = 0.0;
				}
			}
		}

	    psi_Hxz_1 = (float ***)malloc(Imax * sizeof(float **));

	    for(i = 0; i < Imax; i++) {

		    psi_Hxz_1[i] = (float **)malloc((Jmax-1) * sizeof(float *));

		    for(j = 0; j < Jmax; j++) {

				psi_Hxz_1[i][j] = (float *)malloc((nzPML_1-1) * sizeof(float));

				for(k = 0; k < nzPML_1-1; k++) {

					psi_Hxz_1[i][j][k] = 0.0;
				}
			}
		}

	    psi_Hxz_2 = (float ***)malloc(Imax * sizeof(float **));

	    for(i = 0; i < Imax; i++) {

		    psi_Hxz_2[i] = (float **)malloc((Jmax-1) * sizeof(float *));

		    for(j = 0; j < Jmax; j++) {

				psi_Hxz_2[i][j] = (float *)malloc((nzPML_2-1) * sizeof(float));

				for(k = 0; k < nzPML_2-1; k++) {

					psi_Hxz_2[i][j][k] = 0.0;
				}
			}