diff --git a/Common/operator_base.cpp b/Common/operator_base.cpp
index bea1581..9ba853b 100644
--- a/Common/operator_base.cpp
+++ b/Common/operator_base.cpp
@@ -86,12 +86,13 @@ void Operator_Base::SetMaterialStoreFlags(int type, bool val)
 }
 
 
-void Operator_Base::SetBackgroundMaterial(double epsR, double mueR, double kappa, double sigma)
+void Operator_Base::SetBackgroundMaterial(double epsR, double mueR, double kappa, double sigma, double density)
 {
 	SetBackgroundEpsR(epsR);
 	SetBackgroundMueR(mueR);
 	SetBackgroundKappa(kappa);
 	SetBackgroundSigma(sigma);
+	SetBackgroundDensity(density);
 }
 
 void Operator_Base::SetBackgroundEpsR(double val)
@@ -133,3 +134,14 @@ void Operator_Base::SetBackgroundSigma(double val)
 	}
 	m_BG_sigma=val;
 }
+
+
+void Operator_Base::SetBackgroundDensity(double val)
+{
+	if (val<0)
+	{
+		cerr << __func__ << ": Warning, a mass density <0 it not supported, skipping" << endl;
+		return;
+	}
+	m_BG_density=val;
+}
diff --git a/Common/operator_base.h b/Common/operator_base.h
index 9a9fdfa..9180656 100644
--- a/Common/operator_base.h
+++ b/Common/operator_base.h
@@ -109,7 +109,7 @@ public:
 	virtual double GetDiscMaterial(int type, int ny, const unsigned int pos[3]) const = 0;
 
 	//! Set the background material (default is vacuum)
-	virtual void SetBackgroundMaterial(double epsR=0, double mueR=0, double kappa=0, double sigma=0);
+	virtual void SetBackgroundMaterial(double epsR=0, double mueR=0, double kappa=0, double sigma=0, double density=0);
 
 	//! Get background rel. electric permittivity
 	double GetBackgroundEpsR() const {return m_BG_epsR;}
@@ -131,6 +131,11 @@ public:
 	//! Set background magnetic conductivity (artificial)
 	void SetBackgroundSigma(double val);
 
+	//! Get background mass density
+	double GetBackgroundDensity() const {return m_BG_density;}
+	//! Set background mass density
+	void SetBackgroundDensity(double val);
+
 protected:
 	Operator_Base();
 
@@ -155,6 +160,7 @@ protected:
 	double m_BG_mueR;
 	double m_BG_kappa;
 	double m_BG_sigma;
+	double m_BG_density;
 
 	CoordinateSystem m_MeshType;
 	unsigned int numLines[3];
diff --git a/FDTD/operator.cpp b/FDTD/operator.cpp
index 2ab5dcd..ecf1a52 100644
--- a/FDTD/operator.cpp
+++ b/FDTD/operator.cpp
@@ -812,6 +812,7 @@ bool Operator::SetGeometryCSX(ContinuousStructure* geo)
 	SetBackgroundMueR(bg_mat->GetMue());
 	SetBackgroundKappa(bg_mat->GetKappa());
 	SetBackgroundSigma(bg_mat->GetSigma());
+	SetBackgroundDensity(0);
 
 	CSRectGrid* grid=CSX->GetGrid();
 	return SetupCSXGrid(CSRectGrid::Clone(grid));
@@ -1225,6 +1226,8 @@ double Operator::GetMaterial(int ny, const double* coords, int MatType, bool mar
 			return mat->GetMueWeighted(ny,coords);
 		case 3:
 			return mat->GetSigmaWeighted(ny,coords);
+		case 4:
+			return mat->GetDensityWeighted(coords);
 		default:
 			cerr << "Operator::GetMaterial: Error: unknown material type" << endl;
 			return 0;
@@ -1241,6 +1244,8 @@ double Operator::GetMaterial(int ny, const double* coords, int MatType, bool mar
 		return GetBackgroundMueR();
 	case 3:
 		return GetBackgroundSigma();
+	case 4:
+		return GetBackgroundDensity();
 	default:
 		cerr << "Operator::GetMaterial: Error: unknown material type" << endl;
 		return 0;