//....................................................
				//  PML for top Hx, j-direction
				//.....................................................
				 jj = nyPML_2 - 2;
				 for(j = Jmax - nyPML_2; j < Jmax-1; ++j) {

					psi_Hxy_2[i][jj][k] = bh_y_2[jj] * psi_Hxy_2[i][jj][k]
						+ ch_y_2[jj] * (Ez[i][j][k] - Ez[i][j+1][k]) / dy;
					Hx[i][j][k] = Hx[i][j][k] + DB * psi_Hxy_2[i][jj][k];
					jj = jj - 1;
		 		  }
				}
   			}

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

			  for(j = 0; j < Jmax-1; ++j) {
				//....................................................
				//  PML for bottom Hx, k-direction
				//................................................
				 for(k = 1; k < nzPML_1; ++k) {

					psi_Hxz_1[i][j][k-1] = bh_z_1[k-1] * psi_Hxz_1[i][j][k-1]
						+ ch_z_1[k-1] * (Ey[i][j][k] - Ey[i][j][k-1]) / dz;
					Hx[i][j][k] = Hx[i][j][k] + DB * psi_Hxz_1[i][j][k-1];
		 		 }
				//....................................................
				//  PML for top Hx, k-direction
				//...............................................
				 kk = nzPML_2 - 2;
				 for(k = Kmax - nzPML_2; k < Kmax-1; ++k) {

					psi_Hxz_2[i][j][kk] = bh_z_2[kk] * psi_Hxz_2[i][j][kk]
						+ ch_z_2[kk] * (Ey[i][j][k] - Ey[i][j][k-1]) / dz;
					Hx[i][j][k] = Hx[i][j][k] + DB * psi_Hxz_2[i][j][kk];
					kk = kk - 1;
		 		 }
			 }
   		   }

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

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

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

					Hy[i][j][k] = DA * Hy[i][j][k] + DB *
					  ((Ez[i+1][j][k] - Ez[i][j][k]) * den_hx[i] +
					  (Ex[i][j][k-1] - Ex[i][j][k]) * den_hz[k] );
			   }
			  }

			  for(j = 0; j < Jmax-1; ++j) {
				//.......................................................
				//  PML for bottom Hy, i-direction
				//.......................................................
				 for(i = 0; i < nxPML_1-1; ++i){

					  psi_Hyx_1[i][j][k] = bh_x_1[i] * psi_Hyx_1[i][j][k]
							+ ch_x_1[i] * (Ez[i+1][j][k] - Ez[i][j][k]) / dx;
					  Hy[i][j][k] = Hy[i][j][k] + DB * psi_Hyx_1[i][j][k];
				 }
				//.........................................................
				//  PML for top Hy, i-direction
				//.........................................................
				 ii = nxPML_2 - 2;
				 for(i = Imax - nxPML_2; i < Imax-1; ++i) {

			  		psi_Hyx_2[ii][j][k] = bh_x_2[ii] * psi_Hyx_2[ii][j][k]
						+ ch_x_2[ii] * (Ez[i+1][j][k] - Ez[i][j][k]) / dx;
					Hy[i][j][k] = Hy[i][j][k] + DB * psi_Hyx_2[ii][j][k];
					ii = ii - 1;
				 }
			  }
		   }

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

			  for(j = 0; j < Jmax-1; ++j) {
				//.......................................................
				//  PML for bottom Hy, k-direction
				//......................................................
				 for(k = 1; k < nzPML_1; ++k) {

					  psi_Hyz_1[i][j][k-1] = bh_z_1[k-1] * psi_Hyz_1[i][j][k-1]
							+ ch_z_1[k-1] * (Ex[i][j][k-1] - Ex[i][j][k]) / dz;
					  Hy[i][j][k] = Hy[i][j][k] + DB * psi_Hyz_1[i][j][k-1];
				 }
				//.......................................................
				//  PML for top Hy, k-direction
				//.........................................................
				 kk = nzPML_2 - 2;
				 for(k = Kmax - nzPML_2; k < Kmax-1; ++k) {

					psi_Hyz_2[i][j][kk] = bh_z_2[kk] * psi_Hyz_2[i][j][kk]
							+ ch_z_2[kk] * (Ex[i][j][k-1] - Ex[i][j][k]) / dz;
					Hy[i][j][k] = Hy[i][j][k] + DB * psi_Hyz_2[i][j][kk];
					kk = kk - 1;
				 }
			 }
		   }

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

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

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

					Hz[i][j][k] = DA * Hz[i][j][k] + DB
						  * ((Ey[i][j][k] - Ey[i+1][j][k]) * den_hx[i] +
						(Ex[i][j+1][k] - Ex[i][j][k]) * den_hy[j]);
				}
			  }

			  for(j = 0; j < Jmax-1; ++j) {
				//..........................................................
				//  PML for bottom Hz, x-direction
				//..........................................................
				 for(i = 0; i < nxPML_1-1; ++i) {

					  psi_Hzx_1[i][j][k] = bh_x_1[i] * psi_Hzx_1[i][j][k]
							+ ch_x_1[i] * (Ey[i][j][k] - Ey[i+1][j][k]) / dx;
					  Hz[i][j][k] = Hz[i][j][k] + DB * psi_Hzx_1[i][j][k];
				 }
				//..........................................................
				//  PML for top Hz, x-direction
				//..........................................................
				 ii = nxPML_2 - 2;
				 for(i = Imax - nxPML_2; i < Imax-1; ++i) {

					  psi_Hzx_2[ii][j][k] = bh_x_2[ii] * psi_Hzx_2[ii][j][k]
							+ ch_x_2[ii] * (Ey[i][j][k] - Ey[i+1][j][k])/dx;
					  Hz[i][j][k] = Hz[i][j][k] + DB * psi_Hzx_2[ii][j][k];
					  ii = ii - 1;
				  }
			  }

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

					psi_Hzy_1[i][j][k] = bh_y_1[j] * psi_Hzy_1[i][j][k]
						+ ch_y_1[j] * (Ex[i][j+1][k] - Ex[i][j][k]) / dy;
					Hz[i][j][k] = Hz[i][j][k] + DB*  psi_Hzy_1[i][j][k];

				 }
				//.........................................................
				//  PML for top Hz, y-direction
				//..........................................................
				 jj = nyPML_2 - 2;
				 for(j = Jmax - nyPML_2; j < Jmax-1; ++j) {

					psi_Hzy_2[i][jj][k] = bh_y_2[jj] * psi_Hzy_2[i][jj][k]
						+ ch_y_2[jj] * (Ex[i][j+1][k] - Ex[i][j][k]) / dy;
					Hz[i][j][k] = Hz[i][j][k] + DB * psi_Hzy_2[i][jj][k];
					jj = jj - 1;
			   }
			  }
		   }

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

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

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

					id = ID1[i][j][k];
					if(id == 1) { // PEC
				
						Ex[i][j][k] = 0;

					} else {

						 Ex[i][j][k] = CA[id] * Ex[i][j][k] + CB[id] *
						  ((Hz[i][j][k] - Hz[i][j-1][k]) * den_ey[j]  +
						  (Hy[i][j][k] - Hy[i][j][k+1]) * den_ez[k] );
					}
			   }
			}

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

					id = ID1[i][j][k];
					psi_Exy_1[i][j][k] = be_y_1[j] * psi_Exy_1[i][j][k]
						+ ce_y_1[j] * (Hz[i][j][k] - Hz[i][j-1][k])/dy;
					Ex[i][j][k] = Ex[i][j][k] + CB[id] * psi_Exy_1[i][j][k];
				 }
				//.............................................................
				//  PML for top Ex, j-direction
				//.............................................................
				 jj = nyPML_2 - 1;
				 for(j = Jmax - nyPML_2; j < Jmax-1; ++j) {

					id = ID1[i][j][k];
					psi_Exy_2[i][jj][k] = be_y_2[jj] * psi_Exy_2[i][jj][k]
						+ ce_y_2[jj] * (Hz[i][j][k] - Hz[i][j-1][k]) / dy;
					Ex[i][j][k] = Ex[i][j][k] + CB[id] * psi_Exy_2[i][jj][k];
					jj = jj - 1;
				 }
			  }
		   }

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

			  for(j = 1; j < Jmax-1; ++j) {
				//.............................................................
				//  PML for bottom Ex, k-direction
				//.............................................................
				 for(k = 0; k < nzPML_1; ++k) {

					id = ID1[i][j][k];
					psi_Exz_1[i][j][k] = be_z_1[k] * psi_Exz_1[i][j][k]
						+ ce_z_1[k] * (Hy[i][j][k] - Hy[i][j][k+1]) / dz;
					Ex[i][j][k] = Ex[i][j][k] + CB[id] * psi_Exz_1[i][j][k];
				 }
				//..............................................................
				//  PML for top Ex, k-direction
				//..............................................................
				 kk = nzPML_2 - 1;
				 for(k = Kmax - nzPML_2 - 1; k < Kmax-1; ++k) {

					id = ID1[i][j][k];
					psi_Exz_2[i][j][kk] = be_z_2[kk] * psi_Exz_2[i][j][kk]
				 		+ ce_z_2[kk] * (Hy[i][j][k] - Hy[i][j][k+1]) / dz;
					Ex[i][j][k] = Ex[i][j][k] + CB[id] * psi_Exz_2[i][j][kk];
					kk = kk - 1;
				 }
			  }
			}

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

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

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

						id = ID2[i][j][k];
						if(id == 1) { // PEC
					
							Ey[i][j][k] = 0;

						} else {

							 Ey[i][j][k] = CA[id] * Ey[i][j][k] + CB[id] *
							 ((Hz[i-1][j][k] - Hz[i][j][k]) * den_ex[i] +
							 (Hx[i][j][k+1] - Hx[i][j][k]) * den_ez[k] );
						}
				   }
			  }

			  for(j = 0; j < Jmax-1; ++j) {
				//...........................................................
				//  PML for bottom Ey, i-direction
				//...........................................................
				 for(i = 1; i < nxPML_1; ++i) {

					  id = ID2[i][j][k];
					  psi_Eyx_1[i][j][k] = be_x_1[i] * psi_Eyx_1[i][j][k]
							+ ce_x_1[i] * (Hz[i-1][j][k] - Hz[i][j][k]) / dx;
					  Ey[i][j][k] = Ey[i][j][k] + CB[id] * psi_Eyx_1[i][j][k];
				 }
				//............................................................
				//  PML for top Ey, i-direction
				//............................................................
				 ii = nxPML_2 - 1;
				 for(i = Imax - nxPML_2; i < Imax-1; ++i) {

					id = ID2[i][j][k];
			  		psi_Eyx_2[ii][j][k] = be_x_2[ii] * psi_Eyx_2[ii][j][k]
						+ ce_x_2[ii] * (Hz[i-1][j][k] - Hz[i][j][k]) / dx;
			  		Ey[i][j][k] = Ey[i][j][k] + CB[id] * psi_Eyx_2[ii][j][k];
					ii = ii - 1;
				 }
			  }
		   }

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

			  for(j = 0; j < Jmax-1; ++j) {
				//...........................................................
				//  PML for bottom Ey, k-direction
				//...........................................................
				 for(k = 0; k < nzPML_1; ++k) {


					  id = ID2[i][j][k];
					  psi_Eyz_1[i][j][k] = be_z_1[k] * psi_Eyz_1[i][j][k]
							+ ce_z_1[k] * (Hx[i][j][k+1] - Hx[i][j][k]) / dz;
					  Ey[i][j][k] = Ey[i][j][k] + CB[id] * psi_Eyz_1[i][j][k];
				 }
				//...........................................................
				//  PML for top Ey, k-direction
				//............................................................
				 kk = nzPML_2 - 1;
				 for(k = Kmax - nzPML_2 - 1; k < Kmax-1; ++k) {

					id = ID2[i][j][k];
					psi_Eyz_2[i][j][kk] = be_z_2[kk] * psi_Eyz_2[i][j][kk]
							+ ce_z_2[kk] * (Hx[i][j][k+1] - Hx[i][j][k]) / dz;
					Ey[i][j][k] = Ey[i][j][k] + CB[id] * psi_Eyz_2[i][j][kk];
					kk = kk - 1;
				 }
			 }
		   }

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

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

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

					id = ID3[i][j][k];
					if(id == 1) { // PEC
				
						Ez[i][j][k] = 0;

					} else {

						Ez[i][j][k] = CA[id] * Ez[i][j][k] + CB[id]
							  * ((Hy[i][j][k] - Hy[i-1][j][k]) * den_ex[i] +
							(Hx[i][j-1][k] - Hx[i][j][k]) * den_ey[j]);
					}
				  }
			  }

			  for(j = 1; j < Jmax-1; ++j) {
				//............................................................
				//  PML for bottom Ez, x-direction
				//.............................................................
				 for(i = 1; i < nxPML_1; ++i) {


					  id = ID3[i][j][k];
					  psi_Ezx_1[i][j][k] = be_x_1[i] * psi_Ezx_1[i][j][k]
							+ ce_x_1[i] * (Hy[i][j][k] - Hy[i-1][j][k]) / dx;
					  Ez[i][j][k] = Ez[i][j][k] + CB[id] * psi_Ezx_1[i][j][k];
				 }
				//............................................................
				//  PML for top Ez, x-direction
				//............................................................
				 ii = nxPML_2 - 1;
				 for(i = Imax - nxPML_2; i < Imax-1; ++i) {

					  id = ID3[i][j][k];
					  psi_Ezx_2[ii][j][k] = be_x_2[ii] * psi_Ezx_2[ii][j][k]
							+ ce_x_2[ii] * (Hy[i][j][k] - Hy[i-1][j][k]) / dx;
					  Ez[i][j][k] = Ez[i][j][k] + CB[id] * psi_Ezx_2[ii][j][k];
					  ii = ii - 1;
				 }
			  }

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

					id = ID3[i][j][k];
					psi_Ezy_1[i][j][k] = be_y_1[j] * psi_Ezy_1[i][j][k]
						+ ce_y_1[j] * (Hx[i][j-1][k] - Hx[i][j][k]) / dy;
					Ez[i][j][k] = Ez[i][j][k] + CB[id] * psi_Ezy_1[i][j][k];
				 }
				//............................................................
				//  PML for top Ez, y-direction
				//............................................................
				 jj = nyPML_2 - 1;
				 for(j = Jmax - nyPML_2; j < Jmax-1; ++j) {

					id = ID3[i][j][k];
					psi_Ezy_2[i][jj][k] = be_y_2[jj] * psi_Ezy_2[i][jj][k]
						+ ce_y_2[jj] * (Hx[i][j-1][k] - Hx[i][j][k]) / dy;
					Ez[i][j][k] = Ez[i][j][k] + CB[id] * psi_Ezy_2[i][jj][k];
					jj = jj - 1;
			   }
			  }
		   }


		//-----------------------------------------------------------
		//   Apply a point source (Soft)
		//-----------------------------------------------------------
	   i = 25;
	   j = 63;
	   k = 12;
	   source = amp * -2.0 * ((n * dt - tO) / tw)
	   				 * exp(-pow(((n * dt - tO) / tw), 2));//Differentiated Gaussian pulse
	 
	   Ez[i][j][k] = Ez[i][j][k] - CB[ID3[i][j][k]] * source;


		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
		//  WRITE TO OUTPUT FILES
		//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

	   if((n % save_modulus) == 0) {

			writeField(n);
   	   }

	}
   	//  END TIME STEP
    printf("Done time-stepping...\n");

}

	 //Builds an object
	void buildObject() {

		//buildSphere();
		buildDipole();
	}

	//Builds a sphere (Sample code - NOt used in this program)
	void buildSphere() {

		double dist;//distance
		double rad = 8;//(double)Imax / 5.0; // sphere radius
      	double sc = (double)Imax / 2.0;//sphere centre
		double rad2 = 0.3;//(double)Imax / 5.0 - 3.0; // sphere radius

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

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

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

					//compute distance form centre to the point i, j, k
            		dist = sqrt((i + 0.5 - sc) * (i + 0.5 - sc) +
								(j + 0.5 - sc) * (j + 0.5 - sc) +
								(k + 0.5 - sc) * (k + 0.5 - sc));

					//if point is within the sphere
					if (dist <= rad) {
					   //set the material at that point
					   yeeCube (i, j, k, 6);

					}
  				}
            }
       }

	}

	//Builds a dipole
	void buildDipole() {

		int centre = (jstart + jend) / 2;

      	for(i = istart; i <= iend; ++i) {

        	for(j = jstart; j <= jend; ++j) {

          		for(k = kstart; k <= kend; ++k) {
					
					if(j != centre) {

						yeeCube (i, j, k, 1);//PEC material
					}
				}
			}
		}

	}

	//creates a dielctric cube (yee cell) made up of the selected material
    void yeeCube (int I, int J,int K, short mType) {

		   //set face 1 (for EX)
           ID1[I][J][K] = mType;
           ID1[I][J][K + 1] = mType;
           ID1[I][J + 1][K + 1] = mType;
           ID1[I][J + 1][K] = mType;

		   //set face 2 (for EY)
           ID2[I][J][K] = mType;
           ID2[I + 1][J][K] = mType;
           ID2[I + 1][J][K + 1] = mType;
           ID2[I][J][K + 1] = mType;

		   //set face 3 (for EZ)
           ID3[I][J][K] = mType;
           ID3[I + 1][J][K] = mType;
           ID3[I + 1][J + 1][K] = mType;
           ID3[I][J + 1][K] = mType;
      }

	//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	// Saving Output Data to files
	//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
	void writeField (int iteration) {

		FILE *ptr;
		char step[10];
		char fileBaseName[] = "E_Field_";
		sprintf(step, "%d", iteration);
		strcat(fileBaseName, step);
		strcat(fileBaseName, ".txt");

		ptr = fopen(fileBaseName,"wt");

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

			for(j = 0 ; j < Jmax-1 ; j++){
				// |E|
				fprintf(ptr, "%f\t", sqrt(pow(Ex[i][j][ksource], 2) +
					pow(Ey[i][j][ksource], 2) + pow( Ez[i][j][ksource], 2)));
			
			//	fprintf(ptr, "%f\t", Ex[i][j][ksource]);//Ex
			//	fprintf(ptr, "%f\t", Ey[i][j][ksource]);//Ey
			//	fprintf(ptr, "%f\t", Ez[i][j][ksource]);//Ez
			//	fprintf(ptr, "%f\t", Hx[i][j][ksource]);//Hx
			//	fprintf(ptr, "%f\t", Hy[i][j][ksource]);//Hy
			//	fprintf(ptr, "%f\t", Hz[i][j][ksource]);//Hz

			// |H|
			//	fprintf(ptr, "%f\t", sqrt(pow(Hx[i][j][ksource], 2) +
			//		pow(Hy[i][j][ksource], 2) + pow( Hz[i][j][ksource], 2)));

			}
			fprintf(ptr, "\n");
		}

		fclose(ptr);

	}

//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
// END OF PROGRAM CPMLFDTD3D
//~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~