diff --git a/FDTD/engine.cpp b/FDTD/engine.cpp
index 3371e2f..7aa7d8a 100644
--- a/FDTD/engine.cpp
+++ b/FDTD/engine.cpp
@@ -50,8 +50,8 @@ void Engine::Init()
 {
 	Reset();
 	numTS = 0;
-	volt = Create_N_3DArray(numLines);
-	curr = Create_N_3DArray(numLines);
+	volt = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	curr = Create_N_3DArray<FDTD_FLOAT>(numLines);
 
 	file_et.open( "et" );
 
diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp
index f1b7090..c4aa9b4 100644
--- a/FDTD/operator.cpp
+++ b/FDTD/operator.cpp
@@ -319,8 +319,7 @@ void Operator::DumpOperator2File(string filename)
 	}
 
 	cout << "Dumping PEC information to vtk file: " << filename << " ..." ;
-
-	FDTD_FLOAT**** exc = Create_N_3DArray(numLines);
+	FDTD_FLOAT**** exc = Create_N_3DArray<FDTD_FLOAT>(numLines);
 	if (Exc) {
 		for (unsigned int n=0;n<Exc->E_Count;++n)
 			exc[Exc->E_dir[n]][Exc->E_index[0][n]][Exc->E_index[1][n]][Exc->E_index[2][n]] = Exc->E_amp[n];
@@ -351,7 +350,7 @@ void Operator::DumpPEC2File( string filename )
 
 	cout << "Dumping PEC information to vtk file: " << filename << " ..." ;
 
-	FDTD_FLOAT**** pec = Create_N_3DArray( numLines );
+	FDTD_FLOAT**** pec = Create_N_3DArray<FDTD_FLOAT>( numLines );
 	unsigned int pos[3];
 
 	for (pos[0]=0; pos[0]<numLines[0]; pos[0]++) {
@@ -392,10 +391,10 @@ void Operator::DumpMaterial2File(string filename)
 	unsigned int pos[3];
 	double inMat[4];
 
-	epsilon = Create3DArray( numLines);
-	mue = Create3DArray( numLines);
-	kappa = Create3DArray( numLines);
-	sigma = Create3DArray( numLines);
+	epsilon = Create3DArray<FDTD_FLOAT>( numLines);
+	mue = Create3DArray<FDTD_FLOAT>( numLines);
+	kappa = Create3DArray<FDTD_FLOAT>( numLines);
+	sigma = Create3DArray<FDTD_FLOAT>( numLines);
 	for (pos[0]=0;pos[0]<numLines[0];++pos[0])
 	{
 		for (pos[1]=0;pos[1]<numLines[1];++pos[1])
@@ -459,10 +458,10 @@ void Operator::InitOperator()
 	Delete_N_3DArray(vi,numLines);
 	Delete_N_3DArray(iv,numLines);
 	Delete_N_3DArray(ii,numLines);
-	vv = Create_N_3DArray(numLines);
-	vi = Create_N_3DArray(numLines);
-	iv = Create_N_3DArray(numLines);
-	ii = Create_N_3DArray(numLines);
+	vv = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	vi = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	iv = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	ii = Create_N_3DArray<FDTD_FLOAT>(numLines);
 }
 
 void Operator::InitExcitation()
diff --git a/FDTD/operator_sse.cpp b/FDTD/operator_sse.cpp
index 70a7bff..e9866ef 100644
--- a/FDTD/operator_sse.cpp
+++ b/FDTD/operator_sse.cpp
@@ -93,16 +93,16 @@ void Operator_sse::DumpOperator2File(string filename)
 		return;
 	}
 
-	FDTD_FLOAT**** exc = Create_N_3DArray(numLines);
+	FDTD_FLOAT**** exc = Create_N_3DArray<FDTD_FLOAT>(numLines);
 	if (Exc) {
 		for (unsigned int n=0;n<Exc->E_Count;++n)
 			exc[Exc->E_dir[n]][Exc->E_index[0][n]][Exc->E_index[1][n]][Exc->E_index[2][n]] = Exc->E_amp[n];
 	}
 
-	vv = Create_N_3DArray(numLines);
-	vi = Create_N_3DArray(numLines);
-	iv = Create_N_3DArray(numLines);
-	ii = Create_N_3DArray(numLines);
+	vv = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	vi = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	iv = Create_N_3DArray<FDTD_FLOAT>(numLines);
+	ii = Create_N_3DArray<FDTD_FLOAT>(numLines);
 
 	unsigned int pos[3];
 	for (pos[0]=0;pos[0]<numLines[0];++pos[0])
diff --git a/FDTD/processfields_td.cpp b/FDTD/processfields_td.cpp
index c0cb0cb..10cb3dc 100644
--- a/FDTD/processfields_td.cpp
+++ b/FDTD/processfields_td.cpp
@@ -34,7 +34,7 @@ void ProcessFieldsTD::DumpNodeInterpol(string filename)
 	if (m_DumpType==H_FIELD_DUMP)
 	{
 		//create array
-		FDTD_FLOAT**** H_T = Create_N_3DArray(numLines);
+		FDTD_FLOAT**** H_T = Create_N_3DArray<FDTD_FLOAT>(numLines);
 		unsigned int pos[3] = {start[0],start[1],start[2]};
 		unsigned int OpPos[3];
 		double delta;
@@ -95,7 +95,7 @@ void ProcessFieldsTD::DumpNodeInterpol(string filename)
 	if (m_DumpType==E_FIELD_DUMP)
 	{
 		//create array
-		FDTD_FLOAT**** E_T = Create_N_3DArray(numLines);
+		FDTD_FLOAT**** E_T = Create_N_3DArray<FDTD_FLOAT>(numLines);
 		unsigned int pos[3] = {start[0],start[1],start[2]};
 		unsigned int OpPos[3];
 		unsigned int OpPosUp[3];
@@ -154,7 +154,7 @@ void ProcessFieldsTD::DumpCellInterpol(string filename)
 	if (m_DumpType==E_FIELD_DUMP)
 	{
 		//create array
-		FDTD_FLOAT**** E_T = Create_N_3DArray(numDLines);
+		FDTD_FLOAT**** E_T = Create_N_3DArray<FDTD_FLOAT>(numDLines);
 		unsigned int pos[3] = {start[0],start[1],start[2]};
 		unsigned int OpPos[3];
 		double delta;
@@ -215,7 +215,7 @@ void ProcessFieldsTD::DumpCellInterpol(string filename)
 	if (m_DumpType==1)
 	{
 		//create array
-		FDTD_FLOAT**** H_T = Create_N_3DArray(numDLines);
+		FDTD_FLOAT**** H_T = Create_N_3DArray<FDTD_FLOAT>(numDLines);
 		unsigned int pos[3] = {start[0],start[1],start[2]};
 		unsigned int OpPos[3];
 		unsigned int OpPosUp[3];
@@ -277,7 +277,7 @@ void ProcessFieldsTD::DumpNoInterpol(string filename)
 	if (m_DumpType==E_FIELD_DUMP)
 	{
 		//create array
-		FDTD_FLOAT**** E_T = Create_N_3DArray(numLines);
+		FDTD_FLOAT**** E_T = Create_N_3DArray<FDTD_FLOAT>(numLines);
 		for (pos[0]=0;pos[0]<numLines[0];++pos[0])
 		{
 			OpPos[0]=start[0]+pos[0]*subSample[0];
@@ -322,7 +322,7 @@ void ProcessFieldsTD::DumpNoInterpol(string filename)
 	if (m_DumpType==H_FIELD_DUMP)
 	{
 		//create array
-		FDTD_FLOAT**** H_T = Create_N_3DArray(numLines);
+		FDTD_FLOAT**** H_T = Create_N_3DArray<FDTD_FLOAT>(numLines);
 		for (pos[0]=0;pos[0]<numLines[0];++pos[0])
 		{
 			OpPos[0]=start[0]+pos[0]*subSample[0];
diff --git a/tools/array_ops.cpp b/tools/array_ops.cpp
index 2150090..1926811 100644
--- a/tools/array_ops.cpp
+++ b/tools/array_ops.cpp
@@ -57,62 +57,6 @@ void Delete2DArray(FDTD_FLOAT** array, const unsigned int* numLines)
 	delete[] array;
 }
 
-FDTD_FLOAT*** Create3DArray(const unsigned int* numLines)
-{
-	FDTD_FLOAT*** array=NULL;
-	unsigned int pos[3];
-	array = new FDTD_FLOAT**[numLines[0]];
-	for (pos[0]=0;pos[0]<numLines[0];++pos[0])
-	{
-		array[pos[0]] = new FDTD_FLOAT*[numLines[1]];
-		for (pos[1]=0;pos[1]<numLines[1];++pos[1])
-		{
-			array[pos[0]][pos[1]] = new FDTD_FLOAT[numLines[2]];
-			for (pos[2]=0;pos[2]<numLines[2];++pos[2])
-			{
-				array[pos[0]][pos[1]][pos[2]] = 0;
-			}
-		}
-	}
-	return array;
-}
-
-void Delete3DArray(FDTD_FLOAT*** array, const unsigned int* numLines)
-{
-	if (array==NULL) return;
-	unsigned int pos[3];
-	for (pos[0]=0;pos[0]<numLines[0];++pos[0])
-	{
-		for (pos[1]=0;pos[1]<numLines[1];++pos[1])
-		{
-			delete[] array[pos[0]][pos[1]];
-		}
-		delete[] array[pos[0]];
-	}
-	delete[] array;
-}
-
-FDTD_FLOAT**** Create_N_3DArray(const unsigned int* numLines)
-{
-	FDTD_FLOAT**** array=NULL;
-	array = new FDTD_FLOAT***[3];
-	for (int n=0;n<3;++n)
-	{
-		array[n]=Create3DArray(numLines);
-	}
-	return array;
-}
-
-void Delete_N_3DArray(FDTD_FLOAT**** array, const unsigned int* numLines)
-{
-	if (array==NULL) return;
-	for (int n=0;n<3;++n)
-	{
-		Delete3DArray(array[n],numLines);
-	}
-	delete[] array;
-}
-
 void Dump_N_3DArray2File(ostream &file, FDTD_FLOAT**** array, unsigned int* numLines)
 {
 	unsigned int pos[3];
diff --git a/tools/array_ops.h b/tools/array_ops.h
index 53921e7..3644b2b 100644
--- a/tools/array_ops.h
+++ b/tools/array_ops.h
@@ -38,21 +38,84 @@ union f4vector
 FDTD_FLOAT** Create2DArray(const unsigned int* numLines);
 void Delete2DArray(FDTD_FLOAT** array, const unsigned int* numLines);
 
-FDTD_FLOAT*** Create3DArray(const unsigned int* numLines);
-void Delete3DArray(FDTD_FLOAT*** array, const unsigned int* numLines);
+void Dump_N_3DArray2File(ostream &file, FDTD_FLOAT**** array, const unsigned int* numLines);
 
-inline FDTD_FLOAT& Access_N_3DArray(FDTD_FLOAT**** array, unsigned int n, unsigned int* pos)
+
+//
+// templates
+//
+template <typename T>
+inline T& Access_N_3DArray(T**** array, unsigned int n, unsigned int* pos)
 {
 	return array[n][pos[0]][pos[1]][pos[2]];
 }
-inline FDTD_FLOAT& Access_N_3DArray(FDTD_FLOAT**** array, unsigned int n, unsigned int x, unsigned int y, unsigned int z )
+
+template <typename T>
+inline T& Access_N_3DArray(T**** array, unsigned int n, unsigned int x, unsigned int y, unsigned int z )
 {
 	return array[n][x][y][z];
 }
-FDTD_FLOAT**** Create_N_3DArray(const unsigned int* numLines);
-void Delete_N_3DArray(FDTD_FLOAT**** array, const unsigned int* numLines);
 
-void Dump_N_3DArray2File(ostream &file, FDTD_FLOAT**** array, const unsigned int* numLines);
+template <typename T>
+T*** Create3DArray(const unsigned int* numLines)
+{
+	T*** array=NULL;
+	unsigned int pos[3];
+	array = new T**[numLines[0]];
+	for (pos[0]=0;pos[0]<numLines[0];++pos[0])
+	{
+		array[pos[0]] = new T*[numLines[1]];
+		for (pos[1]=0;pos[1]<numLines[1];++pos[1])
+		{
+			array[pos[0]][pos[1]] = new T[numLines[2]];
+			for (pos[2]=0;pos[2]<numLines[2];++pos[2])
+			{
+				array[pos[0]][pos[1]][pos[2]] = 0;
+			}
+		}
+	}
+	return array;
+}
+
+template <typename T>
+T**** Create_N_3DArray(const unsigned int* numLines)
+{
+	T**** array=NULL;
+	array = new T***[3];
+	for (int n=0;n<3;++n)
+	{
+		array[n]=Create3DArray<T>( numLines );
+	}
+	return array;
+}
+
+template <typename T>
+void Delete3DArray(T*** array, const unsigned int* numLines)
+{
+	if (!array) return;
+	unsigned int pos[3];
+	for (pos[0]=0;pos[0]<numLines[0];++pos[0])
+	{
+		for (pos[1]=0;pos[1]<numLines[1];++pos[1])
+		{
+			delete[] array[pos[0]][pos[1]];
+		}
+		delete[] array[pos[0]];
+	}
+	delete[] array;
+}
+
+template <typename T>
+void Delete_N_3DArray(T**** array, const unsigned int* numLines)
+{
+	if (!array) return;
+	for (int n=0;n<3;++n)
+	{
+		Delete3DArray<T>(array[n],numLines);
+	}
+	delete[] array;
+}
+
 
 
 void Delete1DArray_v4sf(f4vector* array);