bh_z_1[i] = 0.0;
		}

	    ch_z_1 = (float *)malloc((nzPML_1-1) * sizeof(float));
		for(i = 0; i < nzPML_1-1; i++) {

			ch_z_1[i] = 0.0;
		}

	    alphah_z_PML_1 = (float *)malloc((nzPML_1-1) * sizeof(float));
		for(i = 0; i < nzPML_1-1; i++) {

			alphah_z_PML_1[i] = 0.0;
		}

	    sigh_z_PML_1 = (float *)malloc((nzPML_1-1) * sizeof(float));
		for(i = 0; i < nzPML_1-1; i++) {

			sigh_z_PML_1[i] = 0.0;
		}

	    kappah_z_PML_1 = (float *)malloc((nzPML_1-1) * sizeof(float));
		for(i = 0; i < nzPML_1-1; i++) {

			kappah_z_PML_1[i] = 0.0;
		}

	    be_z_2 = (float *)malloc((nzPML_2) * sizeof(float));
		for(i = 0; i < nzPML_2; i++) {

			be_z_2[i] = 0.0;
		}

	    ce_z_2 = (float *)malloc((nzPML_2) * sizeof(float));
		for(i = 0; i < nzPML_2; i++) {

			ce_z_2[i] = 0.0;
		}

	    alphae_z_PML_2 = (float *)malloc((nzPML_2) * sizeof(float));
		for(i = 0; i < nzPML_2; i++) {

			alphae_z_PML_2[i] = 0.0;
		}

	    sige_z_PML_2 = (float *)malloc((nzPML_2) * sizeof(float));
		for(i = 0; i < nzPML_2; i++) {

			sige_z_PML_2[i] = 0.0;
		}

	    kappae_z_PML_2 = (float *)malloc((nzPML_2) * sizeof(float));
		for(i = 0; i < nzPML_2; i++) {

			kappae_z_PML_2[i] = 0.0;
		}

	    bh_z_2 = (float *)malloc((nzPML_2-1) * sizeof(float));
		for(i = 0; i < nzPML_2-1; i++) {

			bh_z_2[i] = 0.0;
		}

	    ch_z_2 = (float *)malloc((nzPML_2-1) * sizeof(float));
		for(i = 0; i < nzPML_2-1; i++) {

			ch_z_2[i] = 0.0;
		}

	    alphah_z_PML_2 = (float *)malloc((nzPML_2-1) * sizeof(float));
		for(i = 0; i < nzPML_2-1; i++) {

			alphah_z_PML_2[i] = 0.0;
		}

	    sigh_z_PML_2 = (float *)malloc((nzPML_2-1) * sizeof(float));
		for(i = 0; i < nzPML_2-1; i++) {

			sigh_z_PML_2[i] = 0.0;
		}


	    kappah_z_PML_2 = (float *)malloc((nzPML_2-1) * sizeof(float));
		for(i = 0; i < nzPML_1-1; i++) {

			kappah_z_PML_2[i] = 0.0;
		}
	
	 }

	//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	void setUp() {

		//Time step
		dt = 0.99 / (C * sqrt(1.0 / (dx * dx) + 1.0 / (dy * dy) +
							  1.0 / (dz * dz)));
		//delay
		tO = 4.0 * tw;

		//  Specify the dipole size 
		istart = 24;
		iend = 26;
		jstart = 55;
		jend = 71;
		kstart = 11;
		kend = 13;

		//Material properties
		//Location '0' is for free space and '1' is for PEC
		epsilon[2] = 4.0 * epsO;
		sigma[2] = 0.005;
		epsilon[3] = 8.0 * epsO;
		sigma[3] = 3.96E7;// aluminum
		epsilon[4] = 9.5 * epsO;
		sigma[4] = 5.76E7;//copper
		epsilon[5] = 9.0 * epsO;
		sigma[5] = 2e6;//steel
		epsilon[6] = 2.1 * epsO;
		sigma[6] = 7.8e-4;//teflon
		epsilon[7] = 81 * epsO;
		sigma[7] = 1e-2;//water

		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		//  COMPUTING FIELD UPDATE EQUATION COEFFICIENTS
		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		DA = 1.0;
		DB = dt / muO;

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

			CA[i] = (1.0 - sigma[i] * dt / (2.0 * epsilon[i])) /
				   (1.0 + sigma[i] * dt / (2.0 * epsilon[i]));
			CB[i] = (dt / (epsilon[i])) /
				   (1.0 + sigma[i] * dt / (2.0 * epsilon[i]));

		}

		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		//  PML parameters
		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

		sig_x_max = 0.75 * (0.8 * (m + 1) / (dx * sqrt(muO / (epsO * epsR))));
		sig_y_max = 0.75 * (0.8 * (m + 1) / (dy * sqrt(muO / (epsO * epsR))));
		sig_z_max = 0.75 * (0.8 * (m + 1) / (dz * sqrt(muO / (epsO * epsR))));
		alpha_x_max = 0.24;
		alpha_y_max = alpha_x_max;
		alpha_z_max = alpha_x_max;
		kappa_x_max = 15.0;
		kappa_y_max = kappa_x_max;
		kappa_z_max = kappa_x_max;
		printf("\nTIme step = %e", dt);
		printf("\n Number of steps = %d", nMax);
		printf("\n Total Simulation time = %e Seconds", nMax * dt);

	 }

	  /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		SET CPML PARAMETERS IN EACH DIRECTION
		~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/

	   void initializeCPML() {

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

			sige_x_PML_1[i] = sig_x_max * pow(((nxPML_1 - 1 - i)
							  / (nxPML_1 - 1.0)), m);
			alphae_x_PML_1[i] = alpha_x_max * pow(((float)i
			 							    / (nxPML_1 - 1.0)), ma);
			kappae_x_PML_1[i] = 1.0 + (kappa_x_max - 1.0)*
							    pow((nxPML_1 - 1- i) / (nxPML_1 - 1.0), m);
			be_x_1[i] = exp(-(sige_x_PML_1[i] / kappae_x_PML_1[i] +
									   alphae_x_PML_1[i]) * dt / epsO);

			if ((sige_x_PML_1[i] == 0.0) &&
			   (alphae_x_PML_1[i] == 0.0) && (i == nxPML_1 - 1)) {

			   ce_x_1[i] = 0.0;

			} else {

			   ce_x_1[i] = sige_x_PML_1[i] * (be_x_1[i] - 1.0) /
				  (sige_x_PML_1[i] + kappae_x_PML_1[i] * alphae_x_PML_1[i])
				  / kappae_x_PML_1[i];
			}
		 }

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

		  sigh_x_PML_1[i] = sig_x_max * pow(((nxPML_1 - 1 - i - 0.5)
		  							  / (nxPML_1-1.0)), m);
		  alphah_x_PML_1[i] = alpha_x_max * pow(((i + 1 -0.5)
		  								  / (nxPML_1-1.0)), ma);
		  kappah_x_PML_1[i] = 1.0 + (kappa_x_max - 1.0) *
					pow(((nxPML_1 - 1 - i - 0.5) / (nxPML_1 - 1.0)), m);
		  bh_x_1[i] = exp(-(sigh_x_PML_1[i] / kappah_x_PML_1[i] +
									 alphah_x_PML_1[i]) * dt / epsO);
		  ch_x_1[i] = sigh_x_PML_1[i] * (bh_x_1[i] - 1.0) /
					  (sigh_x_PML_1[i] + kappah_x_PML_1[i] * alphah_x_PML_1[i])
					  / kappah_x_PML_1[i];
		}

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

			sige_x_PML_2[i] = sig_x_max * pow(((nxPML_2 - 1 - i)
							            / (nxPML_2 - 1.0)), m);
			alphae_x_PML_2[i] = alpha_x_max * pow(((float)i
											/ (nxPML_2 - 1.0)), ma);
			kappae_x_PML_2[i] = 1.0 + (kappa_x_max - 1.0)*
							pow((nxPML_2 - 1- i) / (nxPML_2 - 1.0), m);
			be_x_2[i] = exp(-(sige_x_PML_2[i] / kappae_x_PML_2[i] +
									   alphae_x_PML_2[i]) * dt / epsO);

			if ((sige_x_PML_2[i] == 0.0) &&
			   (alphae_x_PML_2[i] == 0.0) && (i == nxPML_2 - 1)) {

			   ce_x_2[i] = 0.0;

			} else {

			   ce_x_2[i] = sige_x_PML_2[i] * (be_x_2[i] - 1.0) /
				  (sige_x_PML_2[i] + kappae_x_PML_2[i] * alphae_x_PML_2[i])
				  / kappae_x_PML_2[i];
			}

	   }

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

		  sigh_x_PML_2[i] = sig_x_max * pow(((nxPML_2 - 1 - i - 0.5)
		  							  / (nxPML_2-1.0)), m);
		  alphah_x_PML_2[i] = alpha_x_max * pow(((i + 1 -0.5)
		  								  / (nxPML_2-1.0)), ma);
		  kappah_x_PML_2[i] = 1.0 + (kappa_x_max - 1.0) *
						pow(((nxPML_2 - 1 - i - 0.5) / (nxPML_2 - 1.0)), m);
		  bh_x_2[i] = exp(-(sigh_x_PML_2[i] / kappah_x_PML_2[i] +
									 alphah_x_PML_2[i]) * dt / epsO);
		  ch_x_2[i] = sigh_x_PML_2[i] * (bh_x_2[i] - 1.0) /
					  (sigh_x_PML_2[i] + kappah_x_PML_2[i] * alphah_x_PML_2[i])
					  / kappah_x_PML_2[i];
		}

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

		   sige_y_PML_1[j] = sig_y_max * pow(((nyPML_1 - 1 - j)
		   							   / (nyPML_1 - 1.0)), m);
		   alphae_y_PML_1[j] = alpha_y_max * pow(((float)j
		   									/ (nyPML_1 - 1.0)), ma);
		   kappae_y_PML_1[j] = 1.0 + (kappa_y_max - 1.0)*
							pow((nyPML_1 - 1- j) / (nyPML_1 - 1.0), m);
		   be_y_1[j] = exp(-(sige_y_PML_1[j] / kappae_y_PML_1[j] +
									  alphae_y_PML_1[j]) * dt / epsO);

		   if ((sige_y_PML_1[j] == 0.0) &&
			  (alphae_y_PML_1[j] == 0.0) && (j == nyPML_1 - 1)) {

			  ce_y_1[j] = 0.0;

		   } else {

			  ce_y_1[j] = sige_y_PML_1[j] * (be_y_1[j] - 1.0) /
				 (sige_y_PML_1[j] + kappae_y_PML_1[j] * alphae_y_PML_1[j])
				 / kappae_y_PML_1[j];
		   }
		}

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

			 sigh_y_PML_1[j] = sig_y_max * pow(((nyPML_1 - 1 - j - 0.5)
			 							 / (nyPML_1-1.0)), m);
			 alphah_y_PML_1[j] = alpha_y_max * pow(((j + 1 -0.5)
			 								 / (nyPML_1-1.0)), ma);
			 kappah_y_PML_1[j] = 1.0 + (kappa_y_max - 1.0) *
							pow(((nyPML_1 - 1 - j - 0.5) / (nyPML_1 - 1.0)), m);
			 bh_y_1[j] = exp(-(sigh_y_PML_1[j] / kappah_y_PML_1[j] +
										alphah_y_PML_1[j]) * dt / epsO);
			 ch_y_1[j] = sigh_y_PML_1[j] * (bh_y_1[j] - 1.0) /
						 (sigh_y_PML_1[j] + kappah_y_PML_1[j] * alphah_y_PML_1[j])
						 / kappah_y_PML_1[j];
		}

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

		   sige_y_PML_2[j] = sig_y_max * pow(((nyPML_2 - 1 - j)
		   							   / (nyPML_2 - 1.0)), m);
		   alphae_y_PML_2[j] = alpha_y_max * pow(((float)j
		   								   / (nyPML_2 - 1.0)), ma);
		   kappae_y_PML_2[j] = 1.0 + (kappa_y_max - 1.0)*
							pow((nyPML_2 - 1- j) / (nyPML_2 - 1.0), m);
		   be_y_2[j] = exp(-(sige_y_PML_2[j] / kappae_y_PML_2[j] +
									  alphae_y_PML_2[j]) * dt / epsO);

		   if ((sige_y_PML_2[j] == 0.0) &&
			  (alphae_y_PML_2[j] == 0.0) && (j == nyPML_2 - 1)) {

			  ce_y_2[j] = 0.0;

		   } else {

			  ce_y_2[j] = sige_y_PML_2[j] * (be_y_2[j] - 1.0) /
				 (sige_y_PML_2[j] + kappae_y_PML_2[j] * alphae_y_PML_2[j])
				 / kappae_y_PML_2[j];
		   }
	  }

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

		 sigh_y_PML_2[j] = sig_y_max * pow(((nyPML_2 - 1 - j - 0.5)
		 							 / (nyPML_2-1.0)), m);
		 alphah_y_PML_2[j] = alpha_y_max * pow(((j + 1 -0.5)
		 								 / (nyPML_2-1.0)), ma);
		 kappah_y_PML_2[j] = 1.0 + (kappa_y_max - 1.0) *
					pow(((nyPML_2 - 1 - j - 0.5) / (nyPML_2 - 1.0)), m);
		 bh_y_2[j] = exp(-(sigh_y_PML_2[j] / kappah_y_PML_2[j] +
									alphah_y_PML_2[j]) * dt / epsO);
		 ch_y_2[j] = sigh_y_PML_2[j] * (bh_y_2[j] - 1.0) /
					 (sigh_y_PML_2[j] + kappah_y_PML_2[j] * alphah_y_PML_2[j])
					 / kappah_y_PML_2[j];
	   }

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

		   sige_z_PML_1[k] = sig_z_max * pow(((nzPML_1 - 1 - k)
		   							   / (nzPML_1 - 1.0)), m);
		   alphae_z_PML_1[k] = alpha_z_max * pow(((float)k
		   								   / (nzPML_1 - 1.0)), ma);
		   kappae_z_PML_1[k] = 1.0 + (kappa_z_max - 1.0)*
							pow((nzPML_1 - 1- k) / (nzPML_1 - 1.0), m);
		   be_z_1[k] = exp(-(sige_z_PML_1[k] / kappae_z_PML_1[k] +
									  alphae_z_PML_1[k]) * dt / epsO);

		   if ((sige_z_PML_1[k] == 0.0) &&
			  (alphae_z_PML_1[k] == 0.0) && (k == nzPML_1 - 1)) {

			  ce_z_1[k] = 0.0;

		   } else {

			  ce_z_1[k] = sige_z_PML_1[k] * (be_z_1[k] - 1.0) /
				 (sige_z_PML_1[k] + kappae_z_PML_1[k] * alphae_z_PML_1[k])
				 / kappae_z_PML_1[k];
		   }
		}

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

		 sigh_z_PML_1[k] = sig_z_max * pow(((nzPML_1 - 1 - k - 0.5)
		 							 / (nzPML_1-1.0)), m);
		 alphah_z_PML_1[k] = alpha_z_max * pow(((k + 1 -0.5)
		 								 / (nzPML_1-1.0)), ma);
		 kappah_z_PML_1[k] = 1.0 + (kappa_z_max - 1.0) *
						pow(((nzPML_1 - 1 - k - 0.5) / (nzPML_1 - 1.0)), m);
		 bh_z_1[k] = exp(-(sigh_z_PML_1[k] / kappah_z_PML_1[k] +
									alphah_z_PML_1[k]) * dt / epsO);
		 ch_z_1[k] = sigh_z_PML_1[k] * (bh_z_1[k] - 1.0) /
					 (sigh_z_PML_1[k] + kappah_z_PML_1[k] * alphah_z_PML_1[k])
					 / kappah_z_PML_1[k];
	   }

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

		   sige_z_PML_2[k] = sig_z_max * pow(((nzPML_2 - 1 - k)
		   							   / (nzPML_2 - 1.0)), m);
		   alphae_z_PML_2[k] = alpha_z_max * pow(((float)k
		   						           / (nzPML_2 - 1.0)), ma);
		   kappae_z_PML_2[k] = 1.0 + (kappa_z_max - 1.0)*
							   pow((nzPML_2 - 1- k) / (nzPML_2 - 1.0), m);
		   be_z_2[k] = exp(-(sige_z_PML_2[k] / kappae_z_PML_2[k] +
									  alphae_z_PML_2[k]) * dt / epsO);

		   if ((sige_z_PML_2[k] == 0.0) &&
			  (alphae_z_PML_2[k] == 0.0) && (k == nzPML_2 - 1)) {

			  ce_z_2[k] = 0.0;

		   } else {

			  ce_z_2[k] = sige_z_PML_2[k] * (be_z_2[k] - 1.0) /
				 (sige_z_PML_2[k] + kappae_z_PML_2[k] * alphae_z_PML_2[k])
				 / kappae_z_PML_2[k];
		   }
		}

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

		 sigh_z_PML_2[k] = sig_z_max * pow(((nzPML_2 - 1 - k - 0.5)
		 							 / (nzPML_2-1.0)), m);
		 alphah_z_PML_2[k] = alpha_z_max * pow(((k + 1 -0.5)
		 							 / (nzPML_2-1.0)), ma);
		 kappah_z_PML_2[k] = 1.0 + (kappa_z_max - 1.0) *
					pow(((nzPML_2 - 1 - k - 0.5) / (nzPML_2 - 1.0)), m);
		 bh_z_2[k] = exp(-(sigh_z_PML_2[k] / kappah_z_PML_2[k] +
									alphah_z_PML_2[k]) * dt / epsO);
		 ch_z_2[k] = sigh_z_PML_2[k] * (bh_z_2[k] - 1.0) /
					 (sigh_z_PML_2[k] + kappah_z_PML_2[k] * alphah_z_PML_2[k])
					 / kappah_z_PML_2[k];
   }

	//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	//  DENOMINATORS FOR FIELD UPDATES
	//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	   ii = nxPML_2 - 2;

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

		  if (i < nxPML_1-1) {

			 den_hx[i] = 1.0 / (kappah_x_PML_1[i] * dx);

		  } else if (i >= Imax - nxPML_2) {

			 den_hx[i] = 1.0 / (kappah_x_PML_2[ii] * dx);
			 ii = ii - 1;

		  } else {

			 den_hx[i] = 1.0 / dx;
  		  }
	   }

	   jj = nyPML_2 - 2;

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

		  if (j < nyPML_1-1) {

			 den_hy[j] = 1.0 / (kappah_y_PML_1[j] * dy);

		  } else if (j >= Jmax - nyPML_2) {

			 den_hy[j] = 1.0 / (kappah_y_PML_2[jj] * dy);
			 jj = jj - 1;

		  } else {

			 den_hy[j] = 1.0 / dy;
  		  }
	   }

	   kk = nzPML_2 - 2;

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

		  if (k < nzPML_1) {

			 den_hz[k] = 1.0 / (kappah_z_PML_1[k-1] * dz);

		  } else if (k >= Kmax - nzPML_2) {

			 den_hz[k] = 1.0 / (kappah_z_PML_2[kk] * dz);
			 kk = kk - 1;

		  } else {

			 den_hz[k] = 1.0 / dz;
		  }
	   }

	   ii = nxPML_2 - 1;

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

		  if (i < nxPML_1) {

			 den_ex[i] = 1.0 / (kappae_x_PML_1[i] * dx);

		  } else if (i >= Imax - nxPML_2) {

			 den_ex[i] = 1.0 / (kappae_x_PML_2[ii] * dx);
			 ii = ii - 1;

		  } else{

			 den_ex[i] = 1.0 / dx;
		  }
	   }

	   jj = nyPML_2 - 1;

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

		  if (j < nyPML_1) {

			 den_ey[j] = 1.0 / (kappae_y_PML_1[j] * dy);

		  } else if (j >= Jmax - nyPML_2) {

			 den_ey[j] = 1.0 / (kappae_y_PML_2[jj] * dy);
			 jj = jj - 1;

		  } else {

			 den_ey[j] = 1.0 / dy;
		  }
	   }

	   kk = nzPML_2 - 1;

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

		  if (k < nzPML_1) {

			 den_ez[k] = 1.0 / (kappae_z_PML_1[k] * dz);

		  }else if (k >= Kmax - 1 - nzPML_2) {

			 den_ez[k] = 1.0 / (kappae_z_PML_2[kk] * dz);
			 kk = kk - 1;

		  } else {

			 den_ez[k] = 1.0 / dz;
		  }
	   }
	  }

	void compute() {

		short id;
		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		//  BEGIN TIME STEP
		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		printf("\nBegin time-stepping...\n");

		for(n = 1; n <= nMax; ++n) {

		   printf("Ez at time step %d at (25, 63, 12) :  %f\n", n, Ez[25][63][12]);

			//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
			//  UPDATE Hx
			//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		   for(k = 1; k < Kmax-1; ++k) {

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

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

				  Hx[i][j][k] = DA * Hx[i][j][k] + DB *
					((Ez[i][j][k] - Ez[i][j+1][k]) * den_hy[j]  +
					  (Ey[i][j][k] - Ey[i][j][k-1]) * den_hz[k] );
	   		   }
			  }

			  for(i = 0; i < Imax-1; ++i) {
				//...............................................
				//  PML for bottom Hx, j-direction
				//...............................................
				 for(j = 0; j < nyPML_1-1; ++j) {

					psi_Hxy_1[i][j][k] = bh_y_1[j] * psi_Hxy_1[i][j][k]
						+ ch_y_1[j] * (Ez[i][j][k] - Ez[i][j+1][k]) / dy;
					Hx[i][j][k] = Hx[i][j][k] + DB * psi_Hxy_1[i][j][k];
		 		  }