diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp
index dfe7687..bab8c41 100644
--- a/FDTD/operator.cpp
+++ b/FDTD/operator.cpp
@@ -521,23 +521,23 @@ void Operator::Calc_ECOperatorPos(int n, unsigned int* pos)
 	unsigned int i = MainOp->SetPos(pos[0],pos[1],pos[2]);
 	if (EC_C[n][i]>0)
 	{
-		GetVV(n,pos[0],pos[1],pos[2]) = (1-dT*EC_G[n][i]/2/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
-		GetVI(n,pos[0],pos[1],pos[2]) = (dT/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]);
+		SetVV(n,pos[0],pos[1],pos[2], (1-dT*EC_G[n][i]/2/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]) );
+		SetVI(n,pos[0],pos[1],pos[2], (dT/EC_C[n][i])/(1+dT*EC_G[n][i]/2/EC_C[n][i]) );
 	}
 	else
 	{
-		GetVV(n,pos[0],pos[1],pos[2]) = 0;
-		GetVI(n,pos[0],pos[1],pos[2]) = 0;
+		SetVV(n,pos[0],pos[1],pos[2], 0 );
+		SetVI(n,pos[0],pos[1],pos[2], 0 );
 	}
 	if (EC_L[n][i]>0)
 	{
-		GetII(n,pos[0],pos[1],pos[2]) = (1-dT*EC_R[n][i]/2/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]);
-		GetIV(n,pos[0],pos[1],pos[2]) = (dT/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]);
+		SetII(n,pos[0],pos[1],pos[2], (1-dT*EC_R[n][i]/2/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]) );
+		SetIV(n,pos[0],pos[1],pos[2], (dT/EC_L[n][i])/(1+dT*EC_R[n][i]/2/EC_L[n][i]) );
 	}
 	else
 	{
-		GetII(n,pos[0],pos[1],pos[2]) = 0;
-		GetIV(n,pos[0],pos[1],pos[2]) = 0;
+		SetII(n,pos[0],pos[1],pos[2], 0 );
+		SetIV(n,pos[0],pos[1],pos[2], 0 );
 	}
 }
 
@@ -632,19 +632,19 @@ void Operator::ApplyElectricBC(bool* dirs)
 			for (pos[nPP]=0;pos[nPP]<numLines[nPP];++pos[nPP])
 			{
 				pos[n]=0;
-				GetVV(nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n];
-				GetVI(nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n];
-				GetVV(nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n];
-				GetVI(nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n];
+				SetVV(nP, pos[0],pos[1],pos[2], GetVV(nP, pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n] );
+				SetVI(nP, pos[0],pos[1],pos[2], GetVI(nP, pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n] );
+				SetVV(nPP,pos[0],pos[1],pos[2], GetVV(nPP,pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n] );
+				SetVI(nPP,pos[0],pos[1],pos[2], GetVI(nPP,pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n] );
 
 				pos[n]=numLines[n]-1;
-				GetVV(n,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1]; // these are outside the FDTD-domain as defined by the main disc
-				GetVI(n,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1]; // these are outside the FDTD-domain as defined by the main disc
+				SetVV(n,  pos[0],pos[1],pos[2], GetVV(n,  pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n+1] ); // these are outside the FDTD-domain as defined by the main disc
+				SetVI(n,  pos[0],pos[1],pos[2], GetVI(n,  pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n+1] ); // these are outside the FDTD-domain as defined by the main disc
 
-				GetVV(nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1];
-				GetVI(nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1];
-				GetVV(nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1];
-				GetVI(nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1];
+				SetVV(nP, pos[0],pos[1],pos[2], GetVV(nP, pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n+1] );
+				SetVI(nP, pos[0],pos[1],pos[2], GetVI(nP, pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n+1] );
+				SetVV(nPP,pos[0],pos[1],pos[2], GetVV(nPP,pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n+1] );
+				SetVI(nPP,pos[0],pos[1],pos[2], GetVI(nPP,pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n+1] );
 			}
 		}
 	}
@@ -663,27 +663,27 @@ void Operator::ApplyMagneticBC(bool* dirs)
 			for (pos[nPP]=0;pos[nPP]<numLines[nPP];++pos[nPP])
 			{
 				pos[n]=0;
-				GetII(n,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n];
-				GetIV(n,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n];
-				GetII(nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n];
-				GetIV(nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n];
-				GetII(nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n];
-				GetIV(nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n];
+				SetII(n,  pos[0],pos[1],pos[2], GetII(n,  pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n] );
+				SetIV(n,  pos[0],pos[1],pos[2], GetIV(n,  pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n] );
+				SetII(nP, pos[0],pos[1],pos[2], GetII(nP, pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n] );
+				SetIV(nP, pos[0],pos[1],pos[2], GetIV(nP, pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n] );
+				SetII(nPP,pos[0],pos[1],pos[2], GetII(nPP,pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n] );
+				SetIV(nPP,pos[0],pos[1],pos[2], GetIV(nPP,pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n] );
 
 				pos[n]=numLines[n]-2;
-				GetII(nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1];
-				GetIV(nP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1];
-				GetII(nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1];
-				GetIV(nPP,pos[0],pos[1],pos[2]) *= (FDTD_FLOAT)!dirs[2*n+1];
+				SetII(nP, pos[0],pos[1],pos[2], GetII(nP, pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n+1] );
+				SetIV(nP, pos[0],pos[1],pos[2], GetIV(nP, pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n+1] );
+				SetII(nPP,pos[0],pos[1],pos[2], GetII(nPP,pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n+1] );
+				SetIV(nPP,pos[0],pos[1],pos[2], GetIV(nPP,pos[0],pos[1],pos[2]) * (FDTD_FLOAT)!dirs[2*n+1] );
 
-				//the last current lines are outside the FDTD domain and cannot be iterated by the FDTD engine
+				// the last current lines are outside the FDTD domain and cannot be iterated by the FDTD engine
 				pos[n]=numLines[n]-1;
-				GetII(n,pos[0],pos[1],pos[2])   = 0;
-				GetIV(n,pos[0],pos[1],pos[2])   = 0;
-				GetII(nP,pos[0],pos[1],pos[2])  = 0;
-				GetIV(nP,pos[0],pos[1],pos[2])  = 0;
-				GetII(nPP,pos[0],pos[1],pos[2]) = 0;
-				GetIV(nPP,pos[0],pos[1],pos[2]) = 0;
+				SetII(n,  pos[0],pos[1],pos[2], 0 );
+				SetIV(n,  pos[0],pos[1],pos[2], 0 );
+				SetII(nP, pos[0],pos[1],pos[2], 0 );
+				SetIV(nP, pos[0],pos[1],pos[2], 0 );
+				SetII(nPP,pos[0],pos[1],pos[2], 0 );
+				SetIV(nPP,pos[0],pos[1],pos[2], 0 );
 			}
 		}
 	}
@@ -1161,8 +1161,8 @@ bool Operator::CalcFieldExcitation()
 								}
 								if (elec->GetExcitType()==1) //hard excite
 								{
-									GetVV(n,pos[0],pos[1],pos[2]) = 0;
-									GetVI(n,pos[0],pos[1],pos[2]) = 0;
+									SetVV(n,pos[0],pos[1],pos[2], 0 );
+									SetVI(n,pos[0],pos[1],pos[2], 0 );
 								}
 							}
 						}
@@ -1203,8 +1203,8 @@ bool Operator::CalcFieldExcitation()
 								}
 								if (elec->GetExcitType()==3) //hard excite
 								{
-									GetII(n,pos[0],pos[1],pos[2]) = 0;
-									GetIV(n,pos[0],pos[1],pos[2]) = 0;
+									SetII(n,pos[0],pos[1],pos[2], 0 );
+									SetIV(n,pos[0],pos[1],pos[2], 0 );
 								}
 							}
 						}
@@ -1266,8 +1266,8 @@ bool Operator::CalcFieldExcitation()
 								}
 								if (elec->GetExcitType()==1) //hard excite
 								{
-									GetVV(n,pos[0],pos[1],pos[2]) = 0;
-									GetVI(n,pos[0],pos[1],pos[2]) = 0;
+									SetVV(n,pos[0],pos[1],pos[2], 0 );
+									SetVI(n,pos[0],pos[1],pos[2], 0 );
 								}
 							}
 						}
@@ -1321,8 +1321,8 @@ void Operator::CalcPEC_Range(unsigned int startX, unsigned int stopX, unsigned i
 					{
 						if (prop->GetType()==CSProperties::METAL) //set to PEC
 						{
-							GetVV(n,pos[0],pos[1],pos[2]) = 0;
-							GetVI(n,pos[0],pos[1],pos[2]) = 0;
+							SetVV(n,pos[0],pos[1],pos[2], 0 );
+							SetVI(n,pos[0],pos[1],pos[2], 0 );
 							++counter[n];
 				//							cerr << "CartOperator::CalcPEC: PEC found at " << pos[0] << " ; "  << pos[1] << " ; " << pos[2] << endl;
 						}
@@ -1359,8 +1359,8 @@ void Operator::CalcPEC_Curves()
 					for (size_t t=0;t<path.dir.size();++t)
 					{
 //						cerr << path.dir.at(t) << " " << path.posPath[0].at(t) << " " << path.posPath[1].at(t) << " " << path.posPath[2].at(t) << endl;
-						GetVV(path.dir.at(t),path.posPath[0].at(t),path.posPath[1].at(t),path.posPath[2].at(t)) = 0;
-						GetVI(path.dir.at(t),path.posPath[0].at(t),path.posPath[1].at(t),path.posPath[2].at(t)) = 0;
+						SetVV(path.dir.at(t),path.posPath[0].at(t),path.posPath[1].at(t),path.posPath[2].at(t), 0 );
+						SetVI(path.dir.at(t),path.posPath[0].at(t),path.posPath[1].at(t),path.posPath[2].at(t), 0 );
 						++m_Nr_PEC[path.dir.at(t)];
 					}
 //					cerr << "found path size: " << path.dir.size() << endl;
diff --git a/FDTD/operator.h b/FDTD/operator.h
index f25efb3..d7ab48c 100644
--- a/FDTD/operator.h
+++ b/FDTD/operator.h
@@ -53,11 +53,17 @@ public:
 
 	virtual bool SetupExcitation(TiXmlElement* Excite, unsigned int maxTS) {return Exc->setupExcitation(Excite,maxTS);};
 
-	inline virtual FDTD_FLOAT& GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vv[n][x][y][z]; }
-	inline virtual FDTD_FLOAT& GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vi[n][x][y][z]; }
+	// the next four functions need to be reimplemented in a derived class
+	inline virtual FDTD_FLOAT GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vv[n][x][y][z]; }
+	inline virtual FDTD_FLOAT GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return vi[n][x][y][z]; }
+	inline virtual FDTD_FLOAT GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return ii[n][x][y][z]; }
+	inline virtual FDTD_FLOAT GetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return iv[n][x][y][z]; }
 
-	inline virtual FDTD_FLOAT& GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return ii[n][x][y][z]; }
-	inline virtual FDTD_FLOAT& GetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return iv[n][x][y][z]; }
+	// the next four functions need to be reimplemented in a derived class
+	inline virtual void SetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { vv[n][x][y][z] = value; }
+	inline virtual void SetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { vi[n][x][y][z] = value; }
+	inline virtual void SetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { ii[n][x][y][z] = value; }
+	inline virtual void SetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { iv[n][x][y][z] = value; }
 
 	virtual void SetBoundaryCondition(int* BCs) {for (int n=0;n<6;++n) m_BC[n]=BCs[n];}
 	virtual void ApplyElectricBC(bool* dirs); //applied by default to all boundaries
diff --git a/FDTD/operator_cylinder.cpp b/FDTD/operator_cylinder.cpp
index ccf9b79..385b29a 100644
--- a/FDTD/operator_cylinder.cpp
+++ b/FDTD/operator_cylinder.cpp
@@ -207,8 +207,8 @@ void Operator_Cylinder::ApplyElectricBC(bool* dirs)
 		{
 			for (pos[2]=0;pos[2]<numLines[2];++pos[2])
 			{
-				GetVV(1,pos[0],pos[1],pos[2]) = 0;
-				GetVI(1,pos[0],pos[1],pos[2]) = 0;
+				SetVV(1,pos[0],pos[1],pos[2], 0 );
+				SetVI(1,pos[0],pos[1],pos[2], 0 );
 			}
 		}
 	}
diff --git a/FDTD/operator_ext_cylinder.cpp b/FDTD/operator_ext_cylinder.cpp
index 4c3fbf9..a5e9458 100644
--- a/FDTD/operator_ext_cylinder.cpp
+++ b/FDTD/operator_ext_cylinder.cpp
@@ -63,7 +63,7 @@ bool Operator_Ext_Cylinder::BuildExtension()
 			C+=inEC[0]*0.5;
 			G+=inEC[1]*0.5;
 		}
-		m_Op->GetVV(2,0,0,pos[2]) = 1;
+		m_Op->SetVV(2,0,0,pos[2], 1);
 		vv_R0[pos[2]] = (1-dT*G/2/C)/(1+dT*G/2/C);
 		vi_R0[pos[2]] = (dT/C)/(1+dT*G/2/C);
 	}
diff --git a/FDTD/operator_ext_upml.cpp b/FDTD/operator_ext_upml.cpp
index 0bfc575..f85f693 100644
--- a/FDTD/operator_ext_upml.cpp
+++ b/FDTD/operator_ext_upml.cpp
@@ -378,8 +378,8 @@ bool Operator_Ext_UPML::BuildExtension()
 							//modify the original operator to perform eq. (7.85) by the main engine (EC-FDTD: equation is multiplied by delta_n)
 							//the engine extension will replace the original voltages with the "voltage flux" (volt*eps0) prior to the voltage updates
 							//after the updates are done the extension will calculate the new voltages (see below) and place them back into the main field domain
-							m_Op->GetVV(n,pos[0],pos[1],pos[2]) = (2*__EPS0__ - kappa_v[nP]*dT) / (2*__EPS0__ + kappa_v[nP]*dT);
-							m_Op->GetVI(n,pos[0],pos[1],pos[2]) = (2*__EPS0__*dT) / (2*__EPS0__ + kappa_v[nP]*dT) * m_Op->GetEdgeLength(n,pos) / m_Op->GetEdgeArea(n,pos);
+							m_Op->SetVV(n,pos[0],pos[1],pos[2], (2*__EPS0__ - kappa_v[nP]*dT) / (2*__EPS0__ + kappa_v[nP]*dT) );
+							m_Op->SetVI(n,pos[0],pos[1],pos[2], (2*__EPS0__*dT) / (2*__EPS0__ + kappa_v[nP]*dT) * m_Op->GetEdgeLength(n,pos) / m_Op->GetEdgeArea(n,pos) );
 
 
 							//operators needed by eq. (7.88) to calculate new voltages from old voltages and old and new "voltage fluxes"
@@ -392,7 +392,7 @@ bool Operator_Ext_UPML::BuildExtension()
 					{
 						//disable upml
 						GetVV(n,loc_pos) = m_Op->GetVV(n,pos[0],pos[1],pos[2]);
-						m_Op->GetVV(n,pos[0],pos[1],pos[2]) = 0;
+						m_Op->SetVV(n,pos[0],pos[1],pos[2], 0 );
 						GetVVFO(n,loc_pos) = 0;
 						GetVVFN(n,loc_pos) = 1;
 					}
@@ -405,8 +405,8 @@ bool Operator_Ext_UPML::BuildExtension()
 							//modify the original operator to perform eq. (7.89) by the main engine (EC-FDTD: equation is multiplied by delta_n)
 							//the engine extension will replace the original currents with the "current flux" (curr*mu0) prior to the current updates
 							//after the updates are done the extension will calculate the new currents (see below) and place them back into the main field domain
-							m_Op->GetII(n,pos[0],pos[1],pos[2]) = (2*__EPS0__ - kappa_i[nP]*dT) / (2*__EPS0__ + kappa_i[nP]*dT);
-							m_Op->GetIV(n,pos[0],pos[1],pos[2]) = (2*__EPS0__*dT) / (2*__EPS0__ + kappa_i[nP]*dT) * m_Op->GetEdgeLength(n,pos,true) / m_Op->GetEdgeArea(n,pos,true);
+							m_Op->SetII(n,pos[0],pos[1],pos[2], (2*__EPS0__ - kappa_i[nP]*dT) / (2*__EPS0__ + kappa_i[nP]*dT) );
+							m_Op->SetIV(n,pos[0],pos[1],pos[2], (2*__EPS0__*dT) / (2*__EPS0__ + kappa_i[nP]*dT) * m_Op->GetEdgeLength(n,pos,true) / m_Op->GetEdgeArea(n,pos,true) );
 
 							//operators needed by eq. (7.90) to calculate new currents from old currents and old and new "current fluxes"
 							GetII(n,loc_pos)   = (2*__EPS0__ - kappa_i[nPP]*dT) / (2*__EPS0__ + kappa_i[nPP]*dT);
@@ -418,7 +418,7 @@ bool Operator_Ext_UPML::BuildExtension()
 					{
 						//disable upml
 						GetII(n,loc_pos) = m_Op->GetII(n,pos[0],pos[1],pos[2]);
-						m_Op->GetII(n,pos[0],pos[1],pos[2]) = 0;
+						m_Op->SetII(n,pos[0],pos[1],pos[2], 0 );
 						GetIIFO(n,loc_pos) = 0;
 						GetIIFN(n,loc_pos) = 1;
 					}
diff --git a/FDTD/operator_sse.h b/FDTD/operator_sse.h
index 7d9687a..d86d089 100644
--- a/FDTD/operator_sse.h
+++ b/FDTD/operator_sse.h
@@ -32,11 +32,15 @@ public:
 
 	virtual void DumpOperator2File(string filename);
 
-	inline virtual FDTD_FLOAT& GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_vv[n][x][y][z%numVectors].f[z/numVectors]; }
-	inline virtual FDTD_FLOAT& GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_vi[n][x][y][z%numVectors].f[z/numVectors]; }
+	inline virtual FDTD_FLOAT GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_vv[n][x][y][z%numVectors].f[z/numVectors]; }
+	inline virtual FDTD_FLOAT GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_vi[n][x][y][z%numVectors].f[z/numVectors]; }
+	inline virtual FDTD_FLOAT GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_ii[n][x][y][z%numVectors].f[z/numVectors]; }
+	inline virtual FDTD_FLOAT GetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_iv[n][x][y][z%numVectors].f[z/numVectors]; }
 
-	inline virtual FDTD_FLOAT& GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_ii[n][x][y][z%numVectors].f[z/numVectors]; }
-	inline virtual FDTD_FLOAT& GetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { return f4_iv[n][x][y][z%numVectors].f[z/numVectors]; }
+	inline virtual void SetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { f4_vv[n][x][y][z%numVectors].f[z/numVectors] = value; }
+	inline virtual void SetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { f4_vi[n][x][y][z%numVectors].f[z/numVectors] = value; }
+	inline virtual void SetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { f4_ii[n][x][y][z%numVectors].f[z/numVectors] = value; }
+	inline virtual void SetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { f4_iv[n][x][y][z%numVectors].f[z/numVectors] = value; }
 
 protected:
 	//! use New() for creating a new Operator
diff --git a/FDTD/operator_sse_compressed.h b/FDTD/operator_sse_compressed.h
index 2928a33..6aeb29e 100644
--- a/FDTD/operator_sse_compressed.h
+++ b/FDTD/operator_sse_compressed.h
@@ -49,11 +49,15 @@ public:
 
 	virtual int CalcECOperator();
 
-	inline virtual FDTD_FLOAT& GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { if (m_Use_Compression) return f4_vv_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetVV(n,x,y,z);}
-	inline virtual FDTD_FLOAT& GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { if (m_Use_Compression) return f4_vi_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetVI(n,x,y,z);}
+	inline virtual FDTD_FLOAT GetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_vv_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetVV(n,x,y,z);}
+	inline virtual FDTD_FLOAT GetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_vi_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetVI(n,x,y,z);}
+	inline virtual FDTD_FLOAT GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_ii_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetII(n,x,y,z);}
+	inline virtual FDTD_FLOAT GetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) const { if (m_Use_Compression) return f4_iv_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetIV(n,x,y,z);}
 
-	inline virtual FDTD_FLOAT& GetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { if (m_Use_Compression) return f4_ii_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetII(n,x,y,z);}
-	inline virtual FDTD_FLOAT& GetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z ) { if (m_Use_Compression) return f4_iv_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors]; else return Operator_sse::GetIV(n,x,y,z);}
+	inline virtual void SetVV( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { if (m_Use_Compression) f4_vv_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors] = value; else Operator_sse::SetVV(n,x,y,z,value);}
+	inline virtual void SetVI( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { if (m_Use_Compression) f4_vi_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors] = value; else Operator_sse::SetVI(n,x,y,z,value);}
+	inline virtual void SetII( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { if (m_Use_Compression) f4_ii_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors] = value; else Operator_sse::SetII(n,x,y,z,value);}
+	inline virtual void SetIV( unsigned int n, unsigned int x, unsigned int y, unsigned int z, FDTD_FLOAT value ) { if (m_Use_Compression) f4_iv_Compressed[n][m_Op_index[x][y][z%numVectors]].f[z/numVectors] = value; else Operator_sse::SetIV(n,x,y,z,value);}
 
 	virtual void ShowStat() const;
 
diff --git a/FDTD/process_efield.cpp b/FDTD/process_efield.cpp
index 464c5e3..b1c189c 100644
--- a/FDTD/process_efield.cpp
+++ b/FDTD/process_efield.cpp
@@ -111,7 +111,7 @@ int ProcessEField::Process()
 				field *= m_weight;
 				for (size_t n=0;n<m_FD_Samples.size();++n)
 				{
-					FD_Values[pol].at(n) += (double)field * std::exp( -2.0 * II * M_PI * m_FD_Samples.at(n) * T );
+					FD_Values[pol].at(n) += (double)field * std::exp( -2.0 * _I * M_PI * m_FD_Samples.at(n) * T );
 				}
 				++m_FD_SampleCount;
 			}
diff --git a/FDTD/process_hfield.cpp b/FDTD/process_hfield.cpp
index c0b5d55..df8acf9 100644
--- a/FDTD/process_hfield.cpp
+++ b/FDTD/process_hfield.cpp
@@ -95,7 +95,7 @@ int ProcessHField::Process()
 				field *= m_weight;
 				for (size_t n=0;n<m_FD_Samples.size();++n)
 				{
-					FD_Values[pol].at(n) += (double)field * std::exp( -2.0 * II * M_PI * m_FD_Samples.at(n) * T );
+					FD_Values[pol].at(n) += (double)field * std::exp( -2.0 * _I * M_PI * m_FD_Samples.at(n) * T );
 				}
 				++m_FD_SampleCount;
 			}
diff --git a/FDTD/processing.h b/FDTD/processing.h
index 2760be2..b6fe330 100644
--- a/FDTD/processing.h
+++ b/FDTD/processing.h
@@ -20,7 +20,7 @@
 
 #include <complex>
 typedef std::complex<double> double_complex;
-#define II double_complex(0.0,1.0)
+#define _I double_complex(0.0,1.0)
 
 #include <iostream>
 #include <fstream>
diff --git a/FDTD/processintegral.cpp b/FDTD/processintegral.cpp
index a690b65..cc2321a 100644
--- a/FDTD/processintegral.cpp
+++ b/FDTD/processintegral.cpp
@@ -77,7 +77,7 @@ int ProcessIntegral::Process()
 			double T = time;
 			for (size_t n=0;n<m_FD_Samples.size();++n)
 			{
-				FD_Values.at(n) += (double)integral * std::exp( -2.0 * II * M_PI * m_FD_Samples.at(n) * T );
+				FD_Values.at(n) += (double)integral * std::exp( -2.0 * _I * M_PI * m_FD_Samples.at(n) * T );
 			}
 			++m_FD_SampleCount;
 			if (m_Flush)
diff --git a/TESTSUITE/enginetests/cavity.m b/TESTSUITE/enginetests/cavity.m
index 61f56e2..3acbc52 100644
--- a/TESTSUITE/enginetests/cavity.m
+++ b/TESTSUITE/enginetests/cavity.m
@@ -6,9 +6,9 @@ function pass = cavity
 CLEANUP = 1;        % if enabled and result is PASS, remove simulation folder
 STOP_IF_FAILED = 1; % if enabled and result is FAILED, stop with error
 
-% engines = {'' '--engine=sse' '--engine=sse-compressed' '--engine=multithreaded' '--engine=sse-compressed-linear'};
-engines = {'--engine=sse-compressed' '--engine=sse-compressed-linear'};
-engines = {'' '--engine=sse-compressed'};
+engines = {'' '--engine=sse' '--engine=sse-compressed' '--engine=multithreaded'};
+% engines = {'--engine=sse-compressed' '--engine=sse-compressed-linear'};
+% engines = {'' '--engine=sse-compressed'};
 
 isOctave = exist('OCTAVE_VERSION','builtin') ~= 0;
 if isOctave
@@ -66,7 +66,7 @@ f_stop = 10e9;
 % setup FDTD parameter
 FDTD = InitFDTD( 1000, 0 );
 FDTD = SetGaussExcite(FDTD,(f_stop-f_start)/2,(f_stop-f_start)/2);
-BC = {'PEC' 'PEC' 'PEC' 'PEC' 'PEC' 'PEC'}; % PEC boundaries
+BC = {'MUR' 'PML_8' 'PEC' 'PEC' 'PEC' 'PEC'}; % PEC boundaries
 FDTD = SetBoundaryCond(FDTD,BC);
 
 % setup CSXCAD geometry