diff --git a/App/App.pro b/App/App.pro
index 3332d775fbd5840a96fd1ffa4a0998c89431b5bc..3c057c0333f6326a565059ea09854cede44e13ee 100644
--- a/App/App.pro
+++ b/App/App.pro
@@ -83,6 +83,7 @@ SOURCES += \
     src/TestPerformance.cpp \
     src/TestRootTree.cpp \
     src/TestRoughness.cpp \
+    src/TestSpecularMagnetic.cpp \
     src/TestSpecularMatrix.cpp \
     src/TestToySimulation.cpp \
     src/TreeEventStructure.cpp \
@@ -143,6 +144,7 @@ HEADERS += \
     inc/TestPerformance.h \
     inc/TestRootTree.h \
     inc/TestRoughness.h \
+    inc/TestSpecularMagnetic.h \
     inc/TestSpecularMatrix.h \
     inc/TestToySimulation.h \
     inc/TreeEventStructure.h \
diff --git a/App/inc/StandardSamples.h b/App/inc/StandardSamples.h
index 79cd1d050829b6be5226b6dac9297a8de108ea92..786208870c1c1bafeeffe892f0d9c6aa63314fe9 100644
--- a/App/inc/StandardSamples.h
+++ b/App/inc/StandardSamples.h
@@ -1,5 +1,5 @@
 // ************************************************************************** //
-//                                                                         
+//
 //  BornAgain: simulate and fit scattering at grazing incidence
 //
 //! @file      App/inc/StandardSamples.h
@@ -30,6 +30,7 @@ ISample *MultilayerOffspecTestcase1a();
 ISample *MultilayerOffspecTestcase1b();
 ISample *MultilayerOffspecTestcase2a();
 ISample *MultilayerOffspecTestcase2b();
+ISample *MultilayerSpecularMagneticTestCase();
 //ISample *IsGISAXS1_CylinderAndPrism();
 //ISample *IsGISAXS2_CylindersMixture();
 //ISample *IsGISAXS3_CylinderDWBA();
diff --git a/App/inc/TestSpecularMagnetic.h b/App/inc/TestSpecularMagnetic.h
new file mode 100644
index 0000000000000000000000000000000000000000..89e881a6a5d4b90ded08c50ba2f2222b78304275
--- /dev/null
+++ b/App/inc/TestSpecularMagnetic.h
@@ -0,0 +1,46 @@
+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      App/inc/TestSpecularMagnetic.h
+//! @brief     Defines class TestSpecularMagnetic.
+//!
+//! @homepage  http://apps.jcns.fz-juelich.de/BornAgain
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2013
+//! @authors   Scientific Computing Group at MLZ Garching
+//! @authors   C. Durniak, G. Pospelov, W. Van Herck, J. Wuttke
+//
+// ************************************************************************** //
+
+#ifndef TESTSPECULARMAGNETIC_H_
+#define TESTSPECULARMAGNETIC_H_
+
+#include "IFunctionalTest.h"
+#include "SpecularMagnetic.h"
+#include "MultiLayer.h"
+#include "OutputData.h"
+
+//! Returns amplitudes for different wave components and polarizations using
+//! the magnetic matrix formalism for a typical multilayer sample and produce
+//! validation plots
+
+class TestSpecularMagnetic : public IFunctionalTest
+{
+public:
+    TestSpecularMagnetic();
+
+    void execute();
+
+private:
+    //! Returns amplitudes vs. alpha_i for several standard samples
+    void test_standard_sample();
+    //! draw results of the test
+    void draw_standard_sample();
+
+    MultiLayer *mp_sample; //!< pointer to multilayer sample
+    OutputData<SpecularMagnetic::MultiLayerCoeff_t  >
+        *mp_coeffs; //!< output data structure
+};
+
+#endif /* TESTSPECULARMAGNETIC_H_ */
diff --git a/App/inc/TestSpecularMatrix.h b/App/inc/TestSpecularMatrix.h
index fe62f84b39cbd752a2a03fae6a3ea13920df7eab..35c8c45edf93710918fd1dbf1bdaf77235b7b51e 100644
--- a/App/inc/TestSpecularMatrix.h
+++ b/App/inc/TestSpecularMatrix.h
@@ -26,12 +26,12 @@
 
 class TestSpecularMatrix : public IFunctionalTest
 {
- public:
+public:
     TestSpecularMatrix();
 
     void execute();
 
- private:
+private:
     //! Returns amplitudes vs. alpha_i for several standard samples
     void test_standard_samples();
     //! draw results of the test
@@ -43,9 +43,8 @@ class TestSpecularMatrix : public IFunctionalTest
     void draw_roughness_set();
 
     MultiLayer *mp_sample; //!< pointer to multilayer sample
-    OutputData<SpecularMatrix::MultiLayerCoeff_t  > *mp_coeffs; //!< output data structure
+    OutputData<SpecularMatrix::MultiLayerCoeff_t  >
+        *mp_coeffs; //!< output data structure
 };
 
-
-
 #endif /* TESTSPECULARMATRIX_H_ */
diff --git a/App/src/FunctionalTestFactory.cpp b/App/src/FunctionalTestFactory.cpp
index b554fffdc4f0d7003eb81b4eee4f78599d730ffc..ea31cb15546de722ead45221ac81c8fb39c4f9d3 100644
--- a/App/src/FunctionalTestFactory.cpp
+++ b/App/src/FunctionalTestFactory.cpp
@@ -48,6 +48,7 @@
 #include "TestPerformance.h"
 #include "TestRootTree.h"
 #include "TestRoughness.h"
+#include "TestSpecularMagnetic.h"
 #include "TestSpecularMatrix.h"
 #include "TestToySimulation.h"
 
@@ -303,5 +304,9 @@ void RegisterFunctionalTests(FunctionalTestFactory *p_test_factory)
         "specularmatrix",
         IFactoryCreateFunction<TestSpecularMatrix, IFunctionalTest>,
         "functional test: specular reflectivity with matrix formalism");
+    p_test_factory->registerItem(
+        "specularmagnetic",
+        IFactoryCreateFunction<TestSpecularMagnetic, IFunctionalTest>,
+        "functional test: specular reflectivity with magnetic matrix formalism");
 }
 
diff --git a/App/src/SampleFactory.cpp b/App/src/SampleFactory.cpp
index 5b7c1a9f83ffa9b0bc5a7a33e93e2afa5d20216b..0bcdff08ebabf16cdb995966320003e8de451e15 100644
--- a/App/src/SampleFactory.cpp
+++ b/App/src/SampleFactory.cpp
@@ -1,5 +1,5 @@
 // ************************************************************************** //
-//                                                                         
+//
 //  BornAgain: simulate and fit scattering at grazing incidence
 //
 //! @file      App/src/SampleFactory.cpp
@@ -32,33 +32,23 @@ SampleFactory::SampleFactory()
 
     // samples used for offspecular reflectivity validation
 
-    //10x2 layers for comparison of diffuse (off-specular) reflectivity with and without cross-correlation in layer's roughnesses
-    registerItem("MultilayerOffspecTestcase1a", StandardSamples::MultilayerOffspecTestcase1a);
-    registerItem("MultilayerOffspecTestcase1b", StandardSamples::MultilayerOffspecTestcase1b);
-
-    // thin layer of air (two different thicknesses) to check diffuse (off-specular) reflectivity
-    registerItem("MultilayerOffspecTestcase2a", StandardSamples::MultilayerOffspecTestcase2a);
-    registerItem("MultilayerOffspecTestcase2b", StandardSamples::MultilayerOffspecTestcase2b);
-
-    // IsGISAXS1 example: cylinder and prism
-    //registerItem("IsGISAXS1_CylinderAndPrism", StandardSamples::IsGISAXS1_CylinderAndPrism);
-
-    // IsGISAXS2 example: mixture of cylindrical particles with two size distribution
-    //registerItem("IsGISAXS2_CylindersMixture", StandardSamples::IsGISAXS2_CylindersMixture);
-
-    // IsGISAXS3 example: cylinder on top of substrate, cylinder in the air, cylinder with size distribution
-//    registerItem("IsGISAXS3_CylinderDWBA", StandardSamples::IsGISAXS3_CylinderDWBA);
-//    registerItem("IsGISAXS3_CylinderBA", StandardSamples::IsGISAXS3_CylinderBA);
-//    registerItem("IsGISAXS3_CylinderBASize", StandardSamples::IsGISAXS3_CylinderBASize);
-
-    // IsGISAXS4 example: cylinders on top of substrate with paracrystal structure factors
-    //registerItem("IsGISAXS4_1DDL", StandardSamples::IsGISAXS4_1DDL);
-    //registerItem("IsGISAXS4_2DDL", StandardSamples::IsGISAXS4_2DDL);
-
-    // IsGISAXS6 example: cylinders with lattice interference function
-//    registerItem("IsGISAXS6_lattice", StandardSamples::IsGISAXS6_lattice);
-//    registerItem("IsGISAXS6_centered", StandardSamples::IsGISAXS6_centered);
-//    registerItem("IsGISAXS6_rotated", StandardSamples::IsGISAXS6_rotated);
+    // 10x2 layers for comparison of diffuse (off-specular) reflectivity with
+    // and without cross-correlation in layer's roughnesses
+    registerItem("MultilayerOffspecTestcase1a",
+            StandardSamples::MultilayerOffspecTestcase1a);
+    registerItem("MultilayerOffspecTestcase1b",
+            StandardSamples::MultilayerOffspecTestcase1b);
+
+    // thin layer of air (two different thicknesses) to check diffuse
+    // (off-specular) reflectivity
+    registerItem("MultilayerOffspecTestcase2a",
+            StandardSamples::MultilayerOffspecTestcase2a);
+    registerItem("MultilayerOffspecTestcase2b",
+            StandardSamples::MultilayerOffspecTestcase2b);
+
+    // 10x2 layers with same index of refraction but opposite magnetization
+    registerItem("MultilayerSpecularMagneticTestCase",
+            StandardSamples::MultilayerSpecularMagneticTestCase);
 
     // IsGISAXS7 example: particle mixture from morphology file
     registerItem("IsGISAXS7_mor", StandardSamples::IsGISAXS7_morphology);
diff --git a/App/src/StandardSamples.cpp b/App/src/StandardSamples.cpp
index 8eeea50a6ca95e5aa44b9d6cbf8190ff2d1d240c..fe0f671d49fb39a6489373a9fd470604760d47af 100644
--- a/App/src/StandardSamples.cpp
+++ b/App/src/StandardSamples.cpp
@@ -1008,4 +1008,52 @@ ISample *StandardSamples::FormFactor_FullSphere()
      return p_multi_layer;
 }
 
+//! Multilayer specular magnetic testcase
+
+ISample *StandardSamples::MultilayerSpecularMagneticTestCase()
+{
+    const IMaterial *mAmbience =
+        MaterialManager::getHomogeneousMaterial
+        ("ambience", 0.0, 0.0 );
+    kvector_t magnetic_field(0.0, 1.0, 0.0);
+    const IMaterial *mPartA =
+        MaterialManager::getHomogeneousMagneticMaterial
+        ("PartA", 5e-6, 0.0, magnetic_field);
+    const IMaterial *mPartB =
+        MaterialManager::getHomogeneousMagneticMaterial
+        ("PartB", 8e-6, 0.0, -magnetic_field );
+    const IMaterial *mSubstrate =
+        MaterialManager::getHomogeneousMaterial
+        ("substrate", 15e-6, 0.0 );
+
+    Layer lAmbience;
+    lAmbience.setMaterial(mAmbience, 0);
+
+    Layer lPartA;
+    lPartA.setMaterial(mPartA, 5.0*Units::nanometer);
+
+    Layer lPartB;
+    lPartB.setMaterial(mPartB, 5.0*Units::nanometer);
+
+    Layer lSubstrate;
+    lSubstrate.setMaterial(mSubstrate, 0);
+
+    MultiLayer *mySample = new MultiLayer;
+
+    // adding layers
+    mySample->addLayer(lAmbience);
+
+    const unsigned nrepetitions = 10;
+    for(unsigned i=0; i<nrepetitions; ++i) {
+        mySample->addLayer(lPartA);
+        mySample->addLayer(lPartB);
+    }
+
+    mySample->addLayer(lSubstrate);
+
+    return mySample;
+}
+
+
+
 
diff --git a/App/src/TestSpecularMagnetic.cpp b/App/src/TestSpecularMagnetic.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..d854f516f96018bbb275a8e00a3d4629c0b9a2d6
--- /dev/null
+++ b/App/src/TestSpecularMagnetic.cpp
@@ -0,0 +1,146 @@
+// ************************************************************************** //
+//
+//  BornAgain: simulate and fit scattering at grazing incidence
+//
+//! @file      App/src/TestSpecularMagnetic.cpp
+//! @brief     Implements class TestSpecularMagnetic.
+//!
+//! @homepage  http://apps.jcns.fz-juelich.de/BornAgain
+//! @license   GNU General Public License v3 or higher (see COPYING)
+//! @copyright Forschungszentrum Jülich GmbH 2013
+//! @authors   Scientific Computing Group at MLZ Garching
+//! @authors   C. Durniak, G. Pospelov, W. Van Herck, J. Wuttke
+//
+// ************************************************************************** //
+
+#include "TestSpecularMagnetic.h"
+
+#include "SampleFactory.h"
+#include "Units.h"
+
+#include <Eigen/Core>
+#include "TCanvas.h"
+#include "TGraph.h"
+#include "DrawHelper.h"
+#include "TH1F.h"
+#include "TLegend.h"
+
+TestSpecularMagnetic::TestSpecularMagnetic()
+: mp_sample(0)
+, mp_coeffs(0)
+{
+    std::cout << "TestSpecularMagnetic::TestSpecularMagnetic() -> Info."
+            << std::endl;
+}
+
+void TestSpecularMagnetic::execute()
+{
+    std::cout << "TestSpecularMagnetic::execute() -> Info." << std::endl;
+
+    // calculate amplitudes for a standard multi-layer sample
+    test_standard_sample();
+}
+
+void TestSpecularMagnetic::test_standard_sample()
+{
+    std::vector<std::string > snames;
+    snames.push_back("MultilayerSpecularMagneticTestCase");
+
+    // loop over standard samples defined in SampleFactory and StandardSamples
+    for(size_t i_sample=0; i_sample<snames.size(); i_sample++){
+        mp_sample = dynamic_cast<MultiLayer *>(SampleFactory::createSample(
+                snames[i_sample]));
+
+        mp_coeffs = new OutputData<SpecularMagnetic::MultiLayerCoeff_t >;
+        mp_coeffs->addAxis(std::string("alpha_i"), 2000, 0.0*Units::degree,
+                2.0*Units::degree);
+        OutputData<SpecularMagnetic::MultiLayerCoeff_t >::iterator it =
+                mp_coeffs->begin();
+        while (it != mp_coeffs->end()) {
+            double alpha_i = mp_coeffs->getValueOfAxis("alpha_i",
+                    it.getIndex());
+            kvector_t kvec;
+            kvec.setLambdaAlphaPhi(1.54*Units::angstrom, -alpha_i, 0.0);
+
+            SpecularMagnetic::MultiLayerCoeff_t coeffs;
+            SpecularMagnetic magneticCalculator;
+            magneticCalculator.execute(*mp_sample, kvec, coeffs);
+
+            *it = coeffs;
+            ++it;
+
+        } // alpha_i
+
+        draw_standard_sample();
+
+        delete mp_sample;
+        delete mp_coeffs;
+    } // i_sample
+}
+
+void TestSpecularMagnetic::draw_standard_sample()
+{
+    static int ncall = 0;
+
+    std::vector<TGraph *> gr_coeff; // different polarization combinations
+    gr_coeff.resize(4);
+    for(size_t i=0; i<4; ++i) {
+        gr_coeff[i] = new TGraph();
+    }
+
+    OutputData<SpecularMagnetic::MultiLayerCoeff_t >::const_iterator it =
+            mp_coeffs->begin();
+    int i_point = 0;
+    while (it != mp_coeffs->end()) {
+        double alpha_i = mp_coeffs->getValueOfAxis("alpha_i", it.getIndex());
+        const SpecularMagnetic::MultiLayerCoeff_t coeffs = *it++;
+
+        // Filling graphics for R coefficients as a function of alpha_i
+        Eigen::Vector2cd rp = coeffs[0].R1plus() + coeffs[0].R2plus();
+        Eigen::Vector2cd rm = coeffs[0].R1min() + coeffs[0].R2min();
+        gr_coeff[0]->SetPoint(i_point, Units::rad2deg(alpha_i),
+                std::abs( (complex_t)rp(0) ) );
+        gr_coeff[1]->SetPoint(i_point, Units::rad2deg(alpha_i),
+                std::abs( (complex_t)rp(1) ) );
+        gr_coeff[2]->SetPoint(i_point, Units::rad2deg(alpha_i),
+                std::abs( (complex_t)rm(0) ) );
+        gr_coeff[3]->SetPoint(i_point++, Units::rad2deg(alpha_i),
+                std::abs( (complex_t)rm(1) ) );
+    }
+
+    // create name of canvas different for each new call of this method
+    std::ostringstream os;
+    os << (ncall++) << std::endl;
+    std::string cname = std::string("c1_test_specular_magnetic")+os.str();
+    TCanvas *c1 = new TCanvas(cname.c_str(),
+            "Fresnel Coefficients in Magnetic Multilayer",1024,768);
+    DrawHelper::SetMagnifier(c1);
+
+    // create subdivision of canvas
+    int ndiv(2);
+    c1->Divide(ndiv,ndiv);
+
+    for(size_t i=0; i<4; ++i) {
+        c1->cd((int)i+1);
+        gPad->SetLogy();
+
+        // calculating histogram limits common for all graphs on given pad
+        double xmin(0), ymin(0), xmax(0), ymax(0);
+        gr_coeff[i]->ComputeRange(xmin, ymin, xmax, ymax);
+        TH1F h1ref("h1ref","h1ref",100, xmin, xmax);
+        h1ref.SetMinimum(1e-6);
+        h1ref.SetMaximum(10);
+        h1ref.SetStats(0);
+        h1ref.SetTitle("");
+        h1ref.GetXaxis()->SetTitle("angle, deg");
+        h1ref.GetYaxis()->SetTitle("|R|");
+        h1ref.DrawCopy();
+
+        TGraph *gr = gr_coeff[i];
+        gr->SetLineColor( kBlue );
+        gr->SetMarkerColor( kBlue );
+        gr->SetMarkerStyle(21);
+        gr->SetMarkerSize(0.2);
+        gr->Draw("pl same");
+    }
+}
diff --git a/Core/Algorithms/src/SpecularMagnetic.cpp b/Core/Algorithms/src/SpecularMagnetic.cpp
index 0ae3e7e9ca3465bea208270d60b090e68c9b8d57..0b83c6170b3d823d287c4ce228ce9df6cbe2f692 100644
--- a/Core/Algorithms/src/SpecularMagnetic.cpp
+++ b/Core/Algorithms/src/SpecularMagnetic.cpp
@@ -42,7 +42,7 @@ void SpecularMagnetic::calculateEigenvalues(const MultiLayer& sample,
                 coeff[i].m_scatt_matrix(1,1) )/2.0;
         coeff[i].lambda(0) = std::sqrt(coeff[i].m_a - coeff[i].m_b_mag);
         coeff[i].lambda(1) = std::sqrt(coeff[i].m_a + coeff[i].m_b_mag);
-        coeff[i].kz = mag_k*coeff[i].lambda;
+        coeff[i].kz = mag_k * coeff[i].lambda;
     }
 }
 
@@ -60,7 +60,7 @@ void SpecularMagnetic::calculateTransferAndBoundary(const MultiLayer& sample,
     coeff[N-1].initializeBottomLayerPhiPsi();
 
     coeff[0].calculateTRMatrices();
-    coeff[0].l = Eigen::Matrix4cd::Identity();
+    coeff[0].l.setIdentity();
     for (int i=(int)N-2; i>0; --i) {
         coeff[i].calculateTRMatrices();
         double t = sample.getLayer(i)->getThickness();
@@ -76,35 +76,34 @@ void SpecularMagnetic::calculateTransferAndBoundary(const MultiLayer& sample,
     // for spin-z polarization in top layer
     if (N>1) {
         // First layer boundary is also top layer boundary:
-        coeff[0].l.setIdentity();
         coeff[0].phi_psi_plus = coeff[1].phi_psi_plus;
         coeff[0].phi_psi_min = coeff[1].phi_psi_min;
         // Normalize all boundary values such that top layer has unit wave
         // amplitude for both spin up and down (and does not contain a
         // transmitted wave amplitude for the opposite polarization)
-        Eigen::Vector2cd T0plusV = coeff[0].T1plus() + coeff[0].T2plus();
-        Eigen::Vector2cd T0minV = coeff[0].T1plus() + coeff[0].T2plus();
-        complex_t cpp, cpm, cmp, cmm;
-        cpp = T0minV(1);
-        cpm = -T0plusV(1);
-        cmp = T0minV(0);
-        cmm = -T0minV(0);
-        Eigen::Vector4cd phipsitemp = cpp * coeff[0].phi_psi_plus +
-                cpm * coeff[0].phi_psi_min;
-        coeff[0].phi_psi_min = cmp * coeff[0].phi_psi_plus +
-                cmm * coeff[0].phi_psi_min;
+        Eigen::Vector2cd T0basisA = coeff[0].T1plus() + coeff[0].T2plus();
+        Eigen::Vector2cd T0basisB = coeff[0].T1min() + coeff[0].T2min();
+        complex_t cpA, cpB, cmA, cmB;
+        cpA = T0basisB(1);
+        cpB = -T0basisA(1);
+        cmA = T0basisB(0);
+        cmB = -T0basisA(0);
+        Eigen::Vector4cd phipsitemp = cpA * coeff[0].phi_psi_plus +
+                cpB * coeff[0].phi_psi_min;
+        coeff[0].phi_psi_min = cmA * coeff[0].phi_psi_plus +
+                cmB * coeff[0].phi_psi_min;
         coeff[0].phi_psi_plus = phipsitemp;
-        complex_t T0plus = coeff[0].phi_psi_plus(2);
-        complex_t T0min = coeff[0].phi_psi_min(3);
-        std::cout << T0plus << std::endl;
-        std::cout << T0min << std::endl;
+        Eigen::Vector2cd T0plusV = coeff[0].T1plus() + coeff[0].T2plus();
+        Eigen::Vector2cd T0minV = coeff[0].T1min() + coeff[0].T2min();
+        complex_t T0plus = T0plusV(0);
+        complex_t T0min = T0minV(1);
         coeff[0].phi_psi_min = coeff[0].phi_psi_min / T0min;
         coeff[0].phi_psi_plus = coeff[0].phi_psi_plus / T0plus;
         for (size_t i=1; i<N; ++i) {
-            phipsitemp = ( cpp * coeff[i].phi_psi_plus +
-                    cpm * coeff[i].phi_psi_min ) / T0plus;
-            coeff[i].phi_psi_min = ( cmp * coeff[i].phi_psi_plus +
-                    cmm * coeff[i].phi_psi_min ) / T0min;
+            phipsitemp = ( cpA * coeff[i].phi_psi_plus +
+                    cpB * coeff[i].phi_psi_min ) / T0plus;
+            coeff[i].phi_psi_min = ( cmA * coeff[i].phi_psi_plus +
+                    cmB * coeff[i].phi_psi_min ) / T0min;
             coeff[i].phi_psi_plus = phipsitemp;
         }
     }
@@ -260,17 +259,17 @@ void SpecularMagnetic::LayerMatrixCoeff::calculateTRWithoutMagnetization()
 
     // R1m:
     R1m.setZero();
-    R1m(0,0) = 0.5;
-    R1m(0,2) = -std::sqrt(m_a)/2.0;
-    R1m(2,0) = -1.0/(2.0*std::sqrt(m_a));
-    R1m(2,2) = 0.5;
+    R1m(1,1) = 0.5;
+    R1m(1,3) = std::sqrt(m_a)/2.0;
+    R1m(3,1) = 1.0/(2.0*std::sqrt(m_a));
+    R1m(3,3) = 0.5;
 
     // T2m:
     T2m.setZero();
-    T2m(1,1) = 0.5;
-    T2m(1,3) = std::sqrt(m_a)/2.0;
-    T2m(3,1) = 1.0/(2.0*std::sqrt(m_a));
-    T2m(3,3) = 0.5;
+    T2m(0,0) = 0.5;
+    T2m(0,2) = -std::sqrt(m_a)/2.0;
+    T2m(2,0) = -1.0/(2.0*std::sqrt(m_a));
+    T2m(2,2) = 0.5;
 
     // R2m:
     R2m.setZero();
diff --git a/Core/Samples/inc/HomogeneousMagneticMaterial.h b/Core/Samples/inc/HomogeneousMagneticMaterial.h
index 12a080449cc9ad6382b758c424abf515c5175f6b..5dc32985d2caf722c86b19b39583f9c9dd6b9f83 100644
--- a/Core/Samples/inc/HomogeneousMagneticMaterial.h
+++ b/Core/Samples/inc/HomogeneousMagneticMaterial.h
@@ -54,6 +54,13 @@ public:
     //! magnetic field and a given wavevector
     virtual Eigen::Matrix2cd getScatteringMatrix(const kvector_t& k) const;
 protected:
+    virtual void print(std::ostream& ostr) const
+    {
+        ostr  << "HomMagMat:" << getName() << "<" << this << ">{ " <<
+                 "R=" << m_refractive_index <<
+                 ", B=" << m_magnetic_field << "}" ;
+    }
+
     kvector_t m_magnetic_field; //!< magnetic field in Tesla
 private:
     //! Function to initialize some private members
diff --git a/Core/Samples/inc/MaterialManager.h b/Core/Samples/inc/MaterialManager.h
index 7ec7854806936277ec7e92fa10dc2d9d2d691266..eba95e275228ffe3eb30615f637007bbcd127c1b 100644
--- a/Core/Samples/inc/MaterialManager.h
+++ b/Core/Samples/inc/MaterialManager.h
@@ -19,7 +19,7 @@
 #include "WinDllMacros.h"
 #include "Exceptions.h"
 #include "ISingleton.h"
-#include "HomogeneousMaterial.h"
+#include "IMaterial.h"
 #include <iostream>
 #include <string>
 #include <map>
@@ -55,8 +55,26 @@ class BA_CORE_API_ MaterialManager: public ISingleton<MaterialManager>
     { return instance().this_getHomogeneousMaterial(
             name, refractive_index_delta, refractive_index_beta); }
 
+    //! Adds magnetic material to database.
+    static const IMaterial *getHomogeneousMagneticMaterial(
+        const std::string& name, const complex_t& refractive_index,
+        const kvector_t &magnetic_field)
+    { return instance().this_getHomogeneousMagneticMaterial(name,
+            refractive_index, magnetic_field); }
+
+    //! Adds magnetic material to database.
+    static const IMaterial *getHomogeneousMagneticMaterial(
+        const std::string& name,
+        double refractive_index_delta,
+        double refractive_index_beta,
+        const kvector_t &magnetic_field)
+    { return instance().this_getHomogeneousMagneticMaterial(
+            name, refractive_index_delta, refractive_index_beta,
+            magnetic_field); }
+
     //! returns number of materials
-    static int getNumberOfMaterials() { return instance().this_getNumberOfMaterials(); }
+    static int getNumberOfMaterials() { return instance().
+            this_getNumberOfMaterials(); }
 
     //! Sends class description to stream.
     friend std::ostream& operator<<(
@@ -97,6 +115,13 @@ class BA_CORE_API_ MaterialManager: public ISingleton<MaterialManager>
     const IMaterial *this_getHomogeneousMaterial(
         const std::string& name,
         double refractive_index_delta, double refractive_index_beta);
+    const IMaterial *this_getHomogeneousMagneticMaterial(
+        const std::string& name, const complex_t& refractive_index,
+        const kvector_t &magnetic_field);
+    const IMaterial *this_getHomogeneousMagneticMaterial(
+        const std::string& name,
+        double refractive_index_delta, double refractive_index_beta,
+        const kvector_t &magnetic_field);
     int this_getNumberOfMaterials() const { return (int)m_materials.size(); }
 
     void check_refractive_index(const complex_t &index);
diff --git a/Core/Samples/src/MaterialManager.cpp b/Core/Samples/src/MaterialManager.cpp
index 6ab618451483b57c4587d1f4b5b4c8ac1876e6f8..e2681ad5e7fa62cb29030abf0ac1000aaa59ae7e 100644
--- a/Core/Samples/src/MaterialManager.cpp
+++ b/Core/Samples/src/MaterialManager.cpp
@@ -14,6 +14,8 @@
 // ************************************************************************** //
 
 #include "MaterialManager.h"
+
+#include "HomogeneousMagneticMaterial.h"
 #include "Exceptions.h"
 #include "MessageService.h"
 #include <boost/thread.hpp>
@@ -86,10 +88,71 @@ const IMaterial *MaterialManager::this_getHomogeneousMaterial(
     double refractive_index_delta,
     double refractive_index_beta)
 {
-    return getHomogeneousMaterial(
+    return this_getHomogeneousMaterial(
         name, complex_t(1.0-refractive_index_delta, refractive_index_beta));
 }
 
+//! Creates magnetic material, and add into database using name of material
+//! as identifier.
+
+const IMaterial* MaterialManager::this_getHomogeneousMagneticMaterial(
+        const std::string& name, const complex_t& refractive_index,
+        const kvector_t& magnetic_field)
+{
+    check_refractive_index(refractive_index);
+
+    static boost::mutex single_mutex;
+    boost::unique_lock<boost::mutex> single_lock( single_mutex );
+    const IMaterial *mat = getMaterial(name);
+    if( mat ) {
+        // check if user is trying to create material
+        // with same name but different parameters
+        const HomogeneousMaterial *old =
+            dynamic_cast<const HomogeneousMaterial *>(mat);
+        if(old->getRefractiveIndex() != refractive_index) {
+            HomogeneousMaterial *non_const_mat =
+                const_cast<HomogeneousMaterial *>(old);
+            non_const_mat->setRefractiveIndex(refractive_index);
+            msglog(MSG::WARNING) <<
+                "MaterialManager::addHomogeneousMagneticMaterial()" <<
+                "-> Redefining refractive index for material '" << name << "'";
+        }
+        const HomogeneousMagneticMaterial *mold =
+            dynamic_cast<const HomogeneousMagneticMaterial *>(mat);
+        if(mold && mold->getMagneticField() != magnetic_field) {
+            HomogeneousMagneticMaterial *non_const_mat =
+                const_cast<HomogeneousMagneticMaterial *>(mold);
+            non_const_mat->setMagneticField(magnetic_field);
+            msglog(MSG::WARNING) <<
+                "MaterialManager::addHomogeneousMagneticMaterial()" <<
+                "-> Redefining magnetic field for material '" << name << "'";
+        }
+        if(!mold) {
+            throw ExistingClassRegistrationException(
+                    "Non-magnetic material with this name"
+                    " was already registered: " + name);
+        }
+        return mat;
+    } else {
+        IMaterial *hmat = new HomogeneousMagneticMaterial(name,
+                refractive_index, magnetic_field);
+        m_materials[name] = hmat;
+        return hmat;
+    }
+}
+
+//! Creates magnetic material, and add into database using name of material
+//! as identifier.
+
+const IMaterial* MaterialManager::this_getHomogeneousMagneticMaterial(
+        const std::string& name, double refractive_index_delta,
+        double refractive_index_beta, const kvector_t& magnetic_field)
+{
+    return this_getHomogeneousMagneticMaterial(
+        name, complex_t(1.0-refractive_index_delta, refractive_index_beta),
+        magnetic_field);
+}
+
 //! Dump this to stream.
 
 void MaterialManager::print(std::ostream& ostr) const