diff --git a/C++/basicClasses/libBasicClasses.cc b/C++/basicClasses/libBasicClasses.cc
index e0fb4870f0a6f031caea0e27d020b199348396a4..c7f578126f38d439169b6ec6cb3de13b6000c1a0 100644
--- a/C++/basicClasses/libBasicClasses.cc
+++ b/C++/basicClasses/libBasicClasses.cc
@@ -246,6 +246,7 @@ void declare_ngData_class(py::module &m, std::string name, std::string suffix=""
     .def("GetMipsAsMatrix", &ngData_class<T, DIM>::getMipsAsMatrix)
     .def("GetRotatedValues", &ngData_class<T, DIM>::GetRotatedValues, py::arg("angle"), py::arg("axis")=2, "rotate all values around the axis in the center")
     .def_property("autoUpdate", &ngData_class<T, DIM>::GetAutoUpdate, &ngData_class<T, DIM>::SetAutoUpdate, "Set the additional field to mark the object to be updated")
+    .def_property("derivative", &ngData_class<T, DIM>::GetDerivative, &ngData_class<T, DIM>::SetDerivative, "Set the derivative used for AutoDiff")
     .def("copy", &ngData_class<T, DIM>::copy)
     .def_property_readonly("numEntries", &ngData_class<T, DIM>::GetNumEntries, "number of allocated objects of type T")
     .def_property("vec", &ngData_class<T, DIM>::GetVec, &ngData_class<T, DIM>::SetVec , "vec")
diff --git a/C++/basicClasses/ngData_class.cc b/C++/basicClasses/ngData_class.cc
index ad08037bc7f254e37b183cfd98f1833715442a04..c75d6fbda3081308898a2435e13d1ca96624b273 100644
--- a/C++/basicClasses/ngData_class.cc
+++ b/C++/basicClasses/ngData_class.cc
@@ -479,6 +479,17 @@ void ngData_class<T, DIM_MESH>::SetAutoUpdate(const bool val){
     autoUpdate = val;
 }
 
+template<typename T, unsigned int DIM_MESH>
+shared_ptr<ngData_class<T, DIM_MESH>> ngData_class<T, DIM_MESH>::GetDerivative() const{
+    return derivative;
+}
+
+template<typename T, unsigned int DIM_MESH>
+void ngData_class<T, DIM_MESH>::SetDerivative(const shared_ptr<ngData_class<T, DIM_MESH>> val){
+    derivative = val;
+}
+
+
 
 template<typename T, unsigned int DIM_MESH>
 void ngData_class<T, DIM_MESH>::Evaluate(const BaseMappedIntegrationPoint& ip, FlatVector<Complex> result) const{
@@ -590,7 +601,15 @@ void ngData_class<T, DIM_MESH>::Evaluate (const BaseMappedIntegrationRule & ir,
 
                         // double eps = 1e-5;
                         double val = this->Evaluate(ir[i]);
-                        double dval2 = this->Evaluate(ir[i]);
+
+                        
+                        double dval2;
+                        if (derivative == nullptr){
+                            dval2 = this->Evaluate(ir[i]);
+                        }
+                        else{
+                            dval2 = derivative->Evaluate(ir[i]);
+                        }
                         // cout << "dval = " << dval << " =?= " << dval2 << endl;
 
                         AutoDiff<1,double> res(val);
diff --git a/C++/basicClasses/ngData_class.h b/C++/basicClasses/ngData_class.h
index 580c63a9b79baf501902c4f0f20aed6af0720255..6ef54a3b786c2ddb310e006df1eaf1edddcc0eae 100644
--- a/C++/basicClasses/ngData_class.h
+++ b/C++/basicClasses/ngData_class.h
@@ -48,6 +48,8 @@ class ngData_class : public CoefficientFunction{
         std::map<std::pair<const unsigned int, const unsigned int>, T> data;
         unsigned int num = 0;
 
+        shared_ptr<ngData_class<T, DIM_MESH>> derivative = nullptr;
+
 
         
     public:
@@ -136,6 +138,9 @@ class ngData_class : public CoefficientFunction{
         [[nodiscard]] bool GetAutoUpdate() const;
         void SetAutoUpdate(const bool val);
 
+        [[nodiscard]] shared_ptr<ngData_class<T, DIM_MESH>> GetDerivative() const;
+        void SetDerivative(const shared_ptr<ngData_class<T, DIM_MESH>> val);
+
         template<typename iT, unsigned int iDIM_MESH>
         friend std::ostream& operator << (std::ostream& o, const ngData_class<iT, iDIM_MESH>& m);
         
diff --git a/C++/preisach/libPreisachScalar.cc b/C++/preisach/libPreisachScalar.cc
index ebc7a4228aaa5ab461e5084f3667a652215c3243..55f32a3edf81fcdc6c90fbb4e7f33d313430bd4a 100644
--- a/C++/preisach/libPreisachScalar.cc
+++ b/C++/preisach/libPreisachScalar.cc
@@ -137,8 +137,8 @@ void ExportReclibPreisachScalar(py::module &main_m){
     .def("addB",  overload_cast_<double>()(&preisach::addB), "add a B value")
     .def("addH",  overload_cast_<const std::vector<double>&>()(&preisach::addH), "add a sequence of H values")
     .def("addB",  overload_cast_<const std::vector<double>&>()(&preisach::addB), "add a sequence of B calues")
-    .def("pilotH",  &preisach::pilotH, "check what would change by adding this H value. results in getPilots. \n return value: [B, mu, mu_diff]")
-    .def("pilotB",  &preisach::pilotB, "check what would change by adding this B value. results in getPilots. \n return value: [H, mu, mu_diff]")
+    .def("pilotH",  &preisach::pilotH, "check what would change by adding this H value. results in getPilot. \n return value: [B, mu, mu_diff]")
+    .def("pilotB",  &preisach::pilotB, "check what would change by adding this B value. results in getPilot. \n return value: [H, mu, mu_diff]")
     .def("pilotB_Honly",  &preisach::pilotB_Honly, "check what would change by adding this B value. results in getPilots")
     // .def("pilotB_quantised",  &preisach::pilotB_quantised, "check what would change when adding this (quasi-quantised) B value.")
     
@@ -155,7 +155,7 @@ void ExportReclibPreisachScalar(py::module &main_m){
     .def("getNu",  &preisach::getNu, "get the current reluctivity. If H = 0 -> mu = 0")
     .def("getNuDiff",  &preisach::getNuDiff, "get the current differential reluctivity")
     .def("getPointList",  &preisach::getPointList, "get the current history.")
-    .def("getPilots", &preisach::getPilots, "get the containers which hold the values of the pilot function.\n return values [H, B, mu, mu_diff]")
+    .def("getPilot", &preisach::getPilot, "get the containers which hold the values of the pilot function.\n return values [H, B, mu, mu_diff]")
     .def("setNegativeSaturation", &preisach::setNegativeSaturation, "set all hysterons to zero. -> output value is the negative saturation. -Bmax")
     .def("demagnetise", (void (preisach::*)())&preisach::demagnetise , "optimale demagnetise the preisach model. With perfectDemag=False -> demag. according to the discretisation of the Preisach plane.")
     .def("getInitialMagnetisationCurve", &preisach::getInitialMagnetisationCurve, "get the inital magnetisation curve of the material. The Accuracy depends on the discretisation of the Preisach plane.\n Use .GetData() to get the pure data. ")
@@ -208,6 +208,7 @@ void ExportReclibPreisachScalar(py::module &main_m){
     .def_property("dim", &EverettMatrix::getDIM, &EverettMatrix::setDIM, "for differnent distributions of the SPMs different adaptions of the EV are necessary.")
     .def_property_readonly("size", &EverettMatrix::getSize, "number of elements in the container")
     .def_property("NA", &EverettMatrix::getNA, &EverettMatrix::setNA, "number of discretisation steps of the Preisach plane in one axis ")
+    .def_property("Transpose", &EverettMatrix::getTranspose, &EverettMatrix::setTranspose, "set the matrix to be transposed")
     .def_property_readonly("isNonLin", &EverettMatrix::getIsNonLin, "check if the EV is set to be nonlienar")
     .def_property_readonly("isInverse", &EverettMatrix::getIsInverse, "check if the EV an inverse")
     .def_property_readonly("maxH", &EverettMatrix::getMaxH, "model limits")
diff --git a/C++/preisach/libPreisachVector.cc b/C++/preisach/libPreisachVector.cc
index a9114b3ad332ed482f3208a5842d03fa9a0a6964..146b5ed5ad614191e90316ba2195fc0eb4520569 100644
--- a/C++/preisach/libPreisachVector.cc
+++ b/C++/preisach/libPreisachVector.cc
@@ -183,7 +183,7 @@ void declare_preisachVector(py::module &m) {
     .def_property("differential", &preisachVector<DIM>::getDifferential, &preisachVector<DIM>::setDifferential, "set the parameter differential. If false, no differential permeability is calculated. It is needed to calculate the interpolation of the output and. ")
     .def_property("conventionalMatDiff", &preisachVector<DIM>::getConventionalMatDiff, &preisachVector<DIM>::setConventionalMatDiff, "set the parameter conventionalMatDiff. if false, the differential are calculated based on pre-compoted Matrices")
     .def_property("NumThreads", &preisachVector<DIM>::getNumThreads, &preisachVector<DIM>::setNumThreads, "set or get the numer of used threads for updating the scalar preisach models.")
-    .def("getPilots", ( pVpilots_t<DIM> &(preisachVector<DIM>::*)()) &preisachVector<DIM>::getPilots, \
+    .def("getPilot", &preisachVector<DIM>::getPilot, \
     "get the containers which hold the values of the pilot function.\n return values [vec_tD H, vec_tD B, mat3x3 mu, mat3x3  mu_diff]\t to get By -> getPilots().second.second; mu_xy  = getPilots().third.xy")
     .def_property_readonly("numBytes", &preisachVector<DIM>::getNumBytes, "approximation of the current memory demand")
     .def("__str__", &preisachVector<DIM>::toString)
diff --git a/C++/preisach/ngPreisachScalar.cc b/C++/preisach/ngPreisachScalar.cc
index bade50eb07d32254223c66a57cd98cf23a0d32cc..d8f2b148685973be54a8b8f2cfdac3afdd123df6 100644
--- a/C++/preisach/ngPreisachScalar.cc
+++ b/C++/preisach/ngPreisachScalar.cc
@@ -180,7 +180,7 @@ void ngPreisachScalar<DIM>::UpdatePast(const shared_ptr<CoefficientFunction>& pt
                 if(ptrFieldStrength != nullptr && ptrFieldStrength->GetAutoUpdate())ptrFieldStrength->Set(el, ip, preisachScalarPoints[el][ip]->getH_interpolated()); 
                 if(ptrMu != nullptr && ptrMu->GetAutoUpdate())ptrMu->Set(el, ip, preisachScalarPoints[el][ip]->getMu()); 
                 if(ptrMuDiff != nullptr && ptrMuDiff->GetAutoUpdate())ptrMuDiff->Set(el, ip, preisachScalarPoints[el][ip]->getMuDiff()); 
-                if(ptrMuDiffDiff != nullptr && ptrMuDiffDiff->GetAutoUpdate())ptrMuDiffDiff->Set(el, ip, preisachScalarPoints[el][ip]->calcMatDiffDiff(false)); 
+                if(ptrMatDiffDiff != nullptr && ptrMatDiffDiff->GetAutoUpdate())ptrMatDiffDiff->Set(el, ip, preisachScalarPoints[el][ip]->calcMatDiffDiff(false)); 
                 if(ptrNu != nullptr && ptrNu->GetAutoUpdate())ptrNu->Set(el, ip, preisachScalarPoints[el][ip]->getNu()); 
                 if(ptrNuDiff != nullptr && ptrNuDiff->GetAutoUpdate())ptrNuDiff->Set(el, ip, preisachScalarPoints[el][ip]->getNuDiff()); 
            
@@ -246,7 +246,7 @@ void ngPreisachScalar<DIM>::Update(const shared_ptr<CoefficientFunction>& ptrIn)
     for(int el = 0; el != mips.size(); ++el){
         double Muold = 0;
         tripleD tD_pilot;
-        quadD qD_pilot;
+        shared_ptr<preisach> pilot;
         unsigned int ip;
 
         for(ip = 0; ip != mips[el].size(); ++ip){
@@ -264,7 +264,7 @@ void ngPreisachScalar<DIM>::Update(const shared_ptr<CoefficientFunction>& ptrIn)
                     }      
                 }
 
-                qD_pilot = preisachScalarPoints[el][ip]->getPilots(); // pilotHBMuMuDiff
+                pilot = preisachScalarPoints[el][ip]->getPilot(); // pilotHBMuMuDiff
                 // cout << "update" << el << " "  << ip << " " << qD_pilot.second << endl;
 
 
@@ -276,30 +276,20 @@ void ngPreisachScalar<DIM>::Update(const shared_ptr<CoefficientFunction>& ptrIn)
                     Muold = MU0;
                 }
 
-                if(ptrFieldStrength != nullptr && ptrFieldStrength->GetAutoUpdate())ptrFieldStrength->Set(el, ip, qD_pilot.first);
-                if(ptrFluxDens != nullptr && ptrFluxDens->GetAutoUpdate())ptrFluxDens->Set(el, ip, qD_pilot.second) ;
+                if(ptrFieldStrength != nullptr && ptrFieldStrength->GetAutoUpdate())ptrFieldStrength->Set(el, ip, pilot->getH_interpolated());
+                if(ptrFluxDens != nullptr && ptrFluxDens->GetAutoUpdate())ptrFluxDens->Set(el, ip, pilot->getB_interpolated()) ;
 	            
 
-                if(!this->ptrEverettMatrix->getIsInverse()){
-                    if(ptrMu != nullptr && ptrMu->GetAutoUpdate())ptrMu->Set(el, ip, qD_pilot.third);
-                    if(ptrMuDiff != nullptr && ptrMuDiff->GetAutoUpdate())ptrMuDiff->Set(el, ip, qD_pilot.fourth);
-                    if(ptrMuDiffDiff != nullptr && ptrMuDiffDiff->GetAutoUpdate())ptrMuDiffDiff->Set(el, ip, preisachScalarPoints[el][ip]->calcMatDiffDiff(true)); 
-                    if(ptrNu != nullptr && ptrNu->GetAutoUpdate())ptrNu->Set(el, ip, (qD_pilot.third == 0 ? 0 :1/qD_pilot.third));
-                    if(ptrNuDiff != nullptr && ptrNuDiff->GetAutoUpdate())ptrNuDiff->Set(el, ip, 1/qD_pilot.fourth);
-                }
-                else{
-                    if(ptrMu != nullptr && ptrMu->GetAutoUpdate())ptrMu->Set(el, ip, (qD_pilot.third == 0 ? 0 :1/qD_pilot.third));
-                    if(ptrMuDiff != nullptr && ptrMuDiff->GetAutoUpdate())ptrMuDiff->Set(el, ip, 1/qD_pilot.fourth);
-                    if(ptrMuDiffDiff != nullptr && ptrMuDiffDiff->GetAutoUpdate())throw std::runtime_error("not implemented yet"); 
-                    if(ptrNu != nullptr && ptrNu->GetAutoUpdate())ptrNu->Set(el, ip, qD_pilot.third);
-                    if(ptrNuDiff != nullptr && ptrNuDiff->GetAutoUpdate())ptrNuDiff->Set(el, ip, qD_pilot.fourth);
-
-                }
+                if(ptrMu != nullptr && ptrMu->GetAutoUpdate())ptrMu->Set(el, ip, pilot->getMu());
+                if(ptrMuDiff != nullptr && ptrMuDiff->GetAutoUpdate())ptrMuDiff->Set(el, ip, pilot->getMuDiff());
+                if(ptrMatDiffDiff != nullptr && ptrMatDiffDiff->GetAutoUpdate())ptrMatDiffDiff->Set(el, ip, preisachScalarPoints[el][ip]->calcMatDiffDiff(true)); 
+                if(ptrNu != nullptr && ptrNu->GetAutoUpdate())ptrNu->Set(el, ip, pilot->getNu());
+                if(ptrNuDiff != nullptr && ptrNuDiff->GetAutoUpdate())ptrNuDiff->Set(el, ip, pilot->getNuDiff());
 
 
                 // save values for Error calculations
 
-                Mu = qD_pilot.third;   
+                Mu = pilot->getMat();   
                 
             
                 // calc Error
@@ -793,12 +783,13 @@ const shared_ptr<ngScalar_class<double, DIM>> ngPreisachScalar<DIM>::GetMuDiff(c
 template <unsigned int DIM>
 const shared_ptr<ngScalar_class<double, DIM>> ngPreisachScalar<DIM>::GetMatDiffDiff(const bool autoupdate){
     // diffdiff_mu
-    if(ptrMuDiffDiff == nullptr){
-        ptrMuDiffDiff = std::make_shared<ngScalar_class<double, DIM>>(mesh, intrule, mips);
-        ptrMuDiffDiff->SetAutoUpdate(true);
+    if(ptrMatDiffDiff == nullptr){
+        ptrMatDiffDiff = std::make_shared<ngScalar_class<double, DIM>>(mesh, intrule, mips);
+        ptrMatDiffDiff->SetAutoUpdate(true);        
         Update();
+
     }
-    return ptrMuDiffDiff;
+    return ptrMatDiffDiff;
 }
 
 template <unsigned int DIM>
diff --git a/C++/preisach/ngPreisachScalar.h b/C++/preisach/ngPreisachScalar.h
index 260e08115f4f472aa8204ecca85decf8155f4639..3ed1cb69a9da3dd3f0cf2dde29cca1208023c3af 100644
--- a/C++/preisach/ngPreisachScalar.h
+++ b/C++/preisach/ngPreisachScalar.h
@@ -141,7 +141,7 @@ class ngPreisachScalar:public CoefficientFunction
         shared_ptr<ngScalar_class<double, DIM>> ptrFieldStrength;
         shared_ptr<ngScalar_class<double, DIM>> ptrMu;
         shared_ptr<ngScalar_class<double, DIM>> ptrMuDiff;
-        shared_ptr<ngScalar_class<double, DIM>> ptrMuDiffDiff;
+        shared_ptr<ngScalar_class<double, DIM>> ptrMatDiffDiff;
         shared_ptr<ngScalar_class<double, DIM>> ptrNu;
         shared_ptr<ngScalar_class<double, DIM>> ptrNuDiff;
         shared_ptr<ngScalar_class<double, DIM>> ptrEnergyDens;
diff --git a/C++/preisach/ngPreisachVector.cc b/C++/preisach/ngPreisachVector.cc
index 674abcea8005fd269e2f84d28ebe9f09bc88eb2c..dac9757d7d5b6567747d33de3fd46bcb1a834e7d 100644
--- a/C++/preisach/ngPreisachVector.cc
+++ b/C++/preisach/ngPreisachVector.cc
@@ -243,7 +243,7 @@ void ngPreisachVector<DIM_MESH, DIM_PREISACH>::Update(const shared_ptr<Coefficie
     // update all preisach points
     parFor(mips.size(), [&](int el) {
         unsigned int ip;
-        pVpilots_t<DIM_PREISACH> val_pilots;
+        shared_ptr<preisachVector<DIM_PREISACH>> pilot;
         double tmp[DIM_PREISACH] = {0};
         FlatVector<double> fV_in(DIM_PREISACH, tmp);
         vec_t<double, DIM_PREISACH> tD_in;
@@ -271,27 +271,25 @@ void ngPreisachVector<DIM_MESH, DIM_PREISACH>::Update(const shared_ptr<Coefficie
                     }      
                 }
 
-                val_pilots = preisachVectorPoints[el][ip]->getPilots();
+                pilot = preisachVectorPoints[el][ip]->getPilot();
                 // cout << "update containers" << endl;
                 // save values for Error calculations
 
 
-                if(ptrFieldStrength != nullptr && ptrFieldStrength->GetAutoUpdate())ptrFieldStrength->Set(el, ip, val_pilots.first);
-                if(ptrFluxDens != nullptr && ptrFluxDens->GetAutoUpdate())ptrFluxDens->Set(el, ip, val_pilots.second) ;
+                if(ptrFieldStrength != nullptr && ptrFieldStrength->GetAutoUpdate())ptrFieldStrength->Set(el, ip, pilot->getH());
+                if(ptrFluxDens != nullptr && ptrFluxDens->GetAutoUpdate())ptrFluxDens->Set(el, ip, pilot->getB()) ;
 
-                if(!this->ptrEverettMatrix->getIsInverse()){
-                    if(ptrMu != nullptr && ptrMu->GetAutoUpdate())ptrMu->Set(el, ip, val_pilots.third);
-                    if(ptrMuDiff != nullptr && ptrMuDiff->GetAutoUpdate())ptrMuDiff->Set(el, ip, val_pilots.fourth);
-                    if(ptrNu != nullptr && ptrNu->GetAutoUpdate())ptrNu->Set(el, ip, (val_pilots.third == 0 ? 0 :val_pilots.third.elementInverse()));
-                    if(ptrNuDiff != nullptr && ptrNuDiff->GetAutoUpdate())ptrNuDiff->Set(el, ip, val_pilots.fourth.elementInverse());
-                }
-                else{
-                    if(ptrMu != nullptr && ptrMu->GetAutoUpdate())ptrMu->Set(el, ip, (val_pilots.third == 0 ? 0 :val_pilots.third.elementInverse()));
-                    if(ptrMuDiff != nullptr && ptrMuDiff->GetAutoUpdate())ptrMuDiff->Set(el, ip, val_pilots.fourth.elementInverse());
-                    if(ptrNu != nullptr && ptrNu->GetAutoUpdate())ptrNu->Set(el, ip, val_pilots.third);
-                    if(ptrNuDiff != nullptr && ptrNuDiff->GetAutoUpdate())ptrNuDiff->Set(el, ip, val_pilots.fourth);
+                if(ptrMu != nullptr && ptrMu->GetAutoUpdate())ptrMu->Set(el, ip, pilot->getMu());
+                if(ptrMuDiff != nullptr && ptrMuDiff->GetAutoUpdate())ptrMuDiff->Set(el, ip, pilot->getMuDiff());
+                if(ptrNu != nullptr && ptrNu->GetAutoUpdate())ptrNu->Set(el, ip, pilot->getNu());
+                if(ptrNuDiff != nullptr && ptrNuDiff->GetAutoUpdate())ptrNuDiff->Set(el, ip, pilot->getNuDiff());
+                // else{
+                //     if(ptrMu != nullptr && ptrMu->GetAutoUpdate())ptrMu->Set(el, ip, (val_pilots.third == 0 ? 0 :val_pilots.third.elementInverse()));
+                //     if(ptrMuDiff != nullptr && ptrMuDiff->GetAutoUpdate())ptrMuDiff->Set(el, ip, val_pilots.fourth.elementInverse());
+                //     if(ptrNu != nullptr && ptrNu->GetAutoUpdate())ptrNu->Set(el, ip, val_pilots.third);
+                //     if(ptrNuDiff != nullptr && ptrNuDiff->GetAutoUpdate())ptrNuDiff->Set(el, ip, val_pilots.fourth);
 
-                }
+                // }
 
 
             }
diff --git a/C++/preisach/preisachScalar/everett.cc b/C++/preisach/preisachScalar/everett.cc
index 78d519d489cafd5654508992f44d4a37a8bc864f..23b1a9a8ac76d1f00426f6bb0436ee41a4c45692 100644
--- a/C++/preisach/preisachScalar/everett.cc
+++ b/C++/preisach/preisachScalar/everett.cc
@@ -44,14 +44,26 @@ namespace everett_ns{
 
 
             
-
-            
+        DIM = o.DIM;
         name=o.name;
-        isNonLin = o.getIsNonLin();
-        isInverse = o.getIsInverse();
-        DIM = o.getDIM();
 
+        isNonLin = o.isNonLin;
+        isInverse = o.isInverse;
+        initialMagnetisationCurveSplineOrder = o.initialMagnetisationCurveSplineOrder;
+        smoothMagnetisationCurveSplineOrder = o.initialMagnetisationCurveSplineOrder;
+
+
+        dE_rising = o.dE_rising;
+        dE_falling = o.dE_falling;
+        transpose = o.transpose;
+
+        // std::shared_ptr<KL> currentMagCurve = nullptr;
         initialMagCurve = nullptr;
+        initialMagCurve2D = nullptr;
+        initialMagCurve3D = nullptr;
+        smoothMagCurve = nullptr;
+
+
     }
 
     // double& EverettMatrix::operator[](const std::pair<double, double>& o){
@@ -60,15 +72,17 @@ namespace everett_ns{
 
 
     std::size_t EverettMatrix::getKey(const double alpha, const double beta) const{
-        
-        // std::size_t h = std::hash<int>{}(xi);
-        // hash_combine(h, yi);
-        // return h;
+        if(transpose){
+            std::size_t h = std::hash<double>{}(-beta);
+            hash_combine(h, -alpha);
+            return h;
+        }
+        else{
+            std::size_t h = std::hash<double>{}(alpha);
+            hash_combine(h, beta);
+            return h;
+        }
 
-        std::size_t h = std::hash<double>{}(alpha);
-        hash_combine(h, beta);
-        return h;
-        
     }
 
     std::vector<std::vector<double>> EverettMatrix::getData() const{
@@ -2013,10 +2027,10 @@ namespace everett_ns{
                 alpha = discInputValues[i];
                 if(i == this->NA -1 && j == this->NA - 1){
                     // corner top right
-                    dE_rising[order-1]->addValue(alpha, beta, MU0);
+                    dE_rising[order-1]->addValue(alpha, beta, (order == 1 ? MU0 : 0));
                 }
                 else if(i == this->NA-1){
-                    // top edge
+                    // top edge -> backward
                     alpha_l = discInputValues[i-1];
                     dE_rising[order-1]->addValue(alpha, beta, (cE->getValue(alpha, beta) - cE->getValue(alpha_l, beta))/(alpha-alpha_l));
                 }
@@ -2074,26 +2088,26 @@ namespace everett_ns{
                 alpha = discInputValues[i];
                 if(i == 0 && j == 0){
                     // corner bottom left
-                    dE_falling[order-1]->addValue(alpha, beta, MU0);
+                    dE_falling[order-1]->addValue(alpha, beta, (order == 1 ? MU0 : 0));
                 }
                 else if(j == 0){
-                    // left edge
+                    // left edge -> forward diff
                     beta_u = discInputValues[1];
-                    dE_falling[order-1]->addValue(alpha, beta, -(cE->getValue(alpha, beta_u) - cE->getValue(alpha, beta))/(beta_u-beta));
+                    dE_falling[order-1]->addValue(alpha, beta, (cE->getValue(alpha, beta_u) - cE->getValue(alpha, beta))/(beta_u-beta));
                 }
                 // else if(i == j){
                 //     // diagonal - dirty fix
                 //     dE_falling[order-1]->addValue(alpha, beta, dE_falling[order-1]->getValue(alpha, discInputValues[j-1]));
                 // }
-                else if (discretisation == "linear"){
+                else if (discretisation == "linear" && j < NA - 1){
                     beta_l = discInputValues[j-1];
                     beta_u = discInputValues[j+1];
-                    dE_falling[order-1]->addValue(alpha, beta, -(cE->getValue(alpha, beta_u) - cE->getValue(alpha, beta_l))/(beta_u-beta_l));
+                    dE_falling[order-1]->addValue(alpha, beta, (cE->getValue(alpha, beta_u) - cE->getValue(alpha, beta_l))/(beta_u-beta_l));
                 }
                 else{
-                    // inner
+                    // inner - > backward
                     beta_l = discInputValues[j-1];
-                    dE_falling[order-1]->addValue(alpha, beta, -(cE->getValue(alpha, beta) - cE->getValue(alpha, beta_l))/(beta-beta_l));
+                    dE_falling[order-1]->addValue(alpha, beta, (cE->getValue(alpha, beta) - cE->getValue(alpha, beta_l))/(beta-beta_l));
                 }
             }
         }
@@ -2621,6 +2635,13 @@ namespace everett_ns{
         isNonLin = in;
     }
 
+    bool EverettMatrix::getTranspose() const{
+        return transpose;
+    }
+    void EverettMatrix::setTranspose(const bool in){
+        transpose = in;
+    }
+
     bool EverettMatrix::getIsInverse() const{
         return isInverse;
     }
diff --git a/C++/preisach/preisachScalar/everett.h b/C++/preisach/preisachScalar/everett.h
index 6f2b82fa9c03f519ad593c5acdf7a4e789bc5244..817db741ff93449ec66ca59f5e4eb9f6b13dd1ab 100644
--- a/C++/preisach/preisachScalar/everett.h
+++ b/C++/preisach/preisachScalar/everett.h
@@ -172,6 +172,9 @@ namespace everett_ns{
             bool getIsInverse() const;
             void setIsInverse(const bool);
 
+            bool getTranspose() const;
+            void setTranspose(const bool);
+
             double getSumOverAllElements() const;
 
             void setEvalVersion(const std::string& o);
@@ -244,8 +247,7 @@ namespace everett_ns{
 
             std::vector<std::shared_ptr<Everett_Empty>> dE_rising;
             std::vector<std::shared_ptr<Everett_Empty>> dE_falling;
-
-
+            bool transpose = false;
 
             // std::shared_ptr<KL> currentMagCurve = nullptr;
             std::shared_ptr<KL> initialMagCurve = nullptr;
diff --git a/C++/preisach/preisachScalar/preisach.cc b/C++/preisach/preisachScalar/preisach.cc
index 82c4b3cfd5372387aec2f7ec384041b83d5d8a94..0247587cdc4c4864cf3f407d5ba38283cb771d63 100644
--- a/C++/preisach/preisachScalar/preisach.cc
+++ b/C++/preisach/preisachScalar/preisach.cc
@@ -10,9 +10,9 @@ namespace preisachScalar_ns{
 
 
     // constructor
-    preisach::preisach(const shared_ptr<EverettMatrix>& E, const bool bDemagnetise, const bool interpolate, const bool differential, const bool conventionalMatDiff_in):
+    preisach::preisach(const shared_ptr<EverettMatrix>& E, const bool bDemagnetise, const bool interpolate, const bool differential, const bool conventionalMatDiff_in, const bool isPilot):
         E{E}, interpolate{interpolate}, differential{differential}, maxInputValIn{E->getMaxInputVal()}, maxOutputValOut{E->getMaxOutputVal()},
-         maxValueForOutside{E->getMaxInputVal()}, matDiffForOutside{CHECK_E_FW_INV(MU0, 1/MU0)}, conventionalMatDiff{conventionalMatDiff_in}{
+         maxValueForOutside{E->getMaxInputVal()}, matDiffForOutside{CHECK_E_FW_INV(MU0, 1/MU0)}, conventionalMatDiff{conventionalMatDiff_in}, isPilot{isPilot}{
 
         if(E == nullptr)
             throw std::runtime_error("invalid EverettMatrix pointer: nullptr");
@@ -36,13 +36,20 @@ namespace preisachScalar_ns{
             setNegativeSaturation();
         }
 
+
+        if(! isPilot){
+            pilot = make_shared<preisach>(E, bDemagnetise, interpolate, differential, conventionalMatDiff, true);
+        }
+
+
     }
 
 
 
     // copy const 
-    preisach::preisach(const preisach &o){
-
+    preisach::preisach(const preisach &o):
+         E{o.E}, interpolate{o.interpolate}, differential{o.differential}, maxInputValIn{o.getMaxInputVal()}, maxOutputValOut{o.getMaxOutputVal()},
+         maxValueForOutside{o.E->getMaxInputVal()}, matDiffForOutside{CHECK_E_FW_INV(MU0, 1/MU0)}, conventionalMatDiff{o.conventionalMatDiff}, isPilot{o.isPilot}{
         // cout << "enter abs preisachScalar" << endl;
 
         outputVal = o.outputVal;
@@ -62,44 +69,31 @@ namespace preisachScalar_ns{
         prevInOutMatMatDiffRising = o.prevInOutMatMatDiffRising;
         prevInOutMatMatDiffFalling = o.prevInOutMatMatDiffFalling;
 
-        E = o.E;
-        interpolate = o.interpolate;
-        differential = o.differential;
-
-        maxInputValIn = o.maxInputValIn;
-        maxOutputValOut = o.maxOutputValOut;
 
         oKL = nullptr;
         bAutoUpdateKL = o.bAutoUpdateKL;
 
-        usePreviousForMuDiffOnly = o.getUsePreviousForMuDiffOnly();
-        usePerfectDemag = o.getUsePerfectDemag();
+        usePreviousForMuDiffOnly = o.usePreviousForMuDiffOnly;
+        usePerfectDemag = o.usePerfectDemag;
 
-        pilotHBMatMatDiff = o.pilotHBMatMatDiff;
-
-        maxValueForOutside = o.maxValueForOutside;
-        matDiffForOutside = o.matDiffForOutside;
-
-        conventionalMatDiff = o.conventionalMatDiff;
+        // pilotHBMatMatDiff = o.pilotHBMatMatDiff;
 
+        if(! o.isPilot){
+            pilot = make_shared<preisach>(E, true, interpolate, differential, conventionalMatDiff, true);
+        }
     }
 
     preisach::preisach(const shared_ptr<preisach> &pD):preisach{*pD}{
         // cout << "enter rel preisachScalar" << endl;
-        
-
     }
     // move
     preisach::preisach(const preisach &&pD){
 
         UNUSED(pD);
         cout << "abs move const " << __PRETTY_FUNCTION__ << endl;
-
-
     }
     preisach::preisach(const shared_ptr<preisach> &&pD):preisach{*pD}{
         cout << "ptr move const " << __PRETTY_FUNCTION__ << endl;
-
     }
     // copy assignment
     preisach& preisach::operator =(const preisach& o){
@@ -135,7 +129,7 @@ namespace preisachScalar_ns{
         usePreviousForMuDiffOnly = o.getUsePreviousForMuDiffOnly();
         usePerfectDemag = o.getUsePerfectDemag();
 
-        pilotHBMatMatDiff = o.pilotHBMatMatDiff;
+        // pilotHBMatMatDiff = o.pilotHBMatMatDiff;
 
 
         maxValueForOutside = o.maxValueForOutside;
@@ -143,10 +137,19 @@ namespace preisachScalar_ns{
 
         conventionalMatDiff = o.conventionalMatDiff;
         
-
+        
+        if(! o.isPilot){
+            pilot = make_shared<preisach>(E, true, interpolate, differential, conventionalMatDiff, true);
+        }
+        
         return *this;
     }
 
+    shared_ptr<preisach> preisach::getPilot() const{
+        
+        return pilot;
+    }
+
     preisach& preisach::operator =(const shared_ptr<preisach>& pD){
         *this = *pD;
         return *this;
@@ -332,9 +335,9 @@ namespace preisachScalar_ns{
         // // calculate the last diffmu of initcurve for values outside of the model 
         // matDiffForOutside = E->getInitialMagnetisationCurve()->Differentiate(getMaxInputVal());
 
-        auto tmp(*this);
-        tmp.demagnetise();
-        matDiffForOutside = tmp.pilotForward(tmp.getE()->getMaxInputVal()*0.99).fourth;
+        *pilot = *this;
+        pilot->demagnetise();
+        matDiffForOutside = pilotForward(getE()->getMaxInputVal()*0.99)->mat_diff;;
         maxValueForOutside = val;       
     }
     double preisach::getMaxValueForOutside()const{
@@ -446,38 +449,25 @@ namespace preisachScalar_ns{
     *  --------------------------------------------------------------------------- */
 
     double preisach::pilotForward_Bonly(double H) const{ 
-        return pilotForward(H).second;
+        return pilotForward(H)->getB_interpolated();
     }
 
     //{B, mu, mu_diff}
-    quadD preisach::pilotForward(double H) const{ 
+    shared_ptr<preisach> preisach::pilotForward(double H) const{ 
     //{B, mu, mu_diff }
-        
-        preisach tmpPoint(*this);
 
+        (*getPilot()) = *this;
+        cout << "bef_pilot" << endl;
+        pilot->addForward(H);
         // cout << tmpPoint << endl;
-        // cout << *this << endl;
-        tmpPoint.addForward(H);
-        // cout << tmpPoint << endl;
-        // cout << tmpPoint.getH() <<  " " << H << "pilot" << endl;
-        if(interpolate){
-            pilotHBMatMatDiff.first = tmpPoint.getH_interpolated();
-            pilotHBMatMatDiff.second = tmpPoint.getB_interpolated();
-            pilotHBMatMatDiff.third = tmpPoint.getMat();
-            pilotHBMatMatDiff.fourth = tmpPoint.getMatDiff();
-        }
-        else{
-            pilotHBMatMatDiff.first = tmpPoint.getH();
-            pilotHBMatMatDiff.second = tmpPoint.getB();
-            pilotHBMatMatDiff.third = tmpPoint.getMat();
-            pilotHBMatMatDiff.fourth = tmpPoint.getMatDiff();
-        }
+        // cout << pilot->getH() <<  " " << H << "pilot" << endl;
+
             
-        return pilotHBMatMatDiff;
+        return pilot;
 
     }
 
-    quadD preisach::pilotInverse(double B/*, const unsigned int iNmax, const double dErr*/) const{ 
+    shared_ptr<preisach> preisach::pilotInverse(double B/*, const unsigned int iNmax, const double dErr*/) const{ 
     //{B, mu, mu_diff }
         
         // check if simply the same
@@ -488,34 +478,23 @@ namespace preisachScalar_ns{
         //         return  quadD(inputVal,outputVal, mat, mat_diff);
         // }
 
-        preisach tmpPoint(*this);
+        
+        (*getPilot()) = *this;
 
         // cout << tmpPoint << endl;
         // cout << *this << endl;
-        tmpPoint.addInverse(B);
+        pilot->addInverse(B);
         // cout << tmpPoint << endl;
         // cout << tmpPoint.getB() <<  " " << B << "pilot" << endl;
 
 
-        // Save the values in pilotVector
-        if(interpolate){
-            pilotHBMatMatDiff.first = tmpPoint.getH_interpolated();
-            pilotHBMatMatDiff.second = tmpPoint.getB_interpolated();
-            pilotHBMatMatDiff.third = tmpPoint.getMat();
-            pilotHBMatMatDiff.fourth = tmpPoint.getMatDiff();
-        }
-        else{
-            pilotHBMatMatDiff.first = tmpPoint.getH();
-            pilotHBMatMatDiff.second = tmpPoint.getB();
-            pilotHBMatMatDiff.third = tmpPoint.getMat();
-            pilotHBMatMatDiff.fourth = tmpPoint.getMatDiff();
-        }
+
             
-        return pilotHBMatMatDiff;
+        return pilot;
 
     }
     double preisach::pilotInverse_Honly(double b_pilot) const{
-        return pilotInverse(b_pilot).first;
+        return pilotInverse(b_pilot)->getH_interpolated();
     }
     double preisach::pilotInverse_quantised(double b_pilot) const{
         double H;
@@ -557,7 +536,7 @@ namespace preisachScalar_ns{
         do{
             In_current = (dTop + dBtm)/2;
 
-            Output_new = CHECK_E_FW_INV(pilotForward(In_current).second, pilotForward(In_current).first);
+            Output_new = pilotForward(In_current)->getOutputVal_interpolated();
 
 
             dErrOld = dErr;
@@ -1161,6 +1140,7 @@ namespace preisachScalar_ns{
                         B_u = E->getDE_rising(order)->getValue(points.back().first, points[points.size()-2].second);
                     }
                     else{
+                        
                         B_u = E->getDE_rising(order)->getValue(points.back().first, -maxInputValIn);
                     }
                 }
@@ -1171,11 +1151,11 @@ namespace preisachScalar_ns{
                     }
                     // B_u = E->getDE_falling(order)->getValue(points.back().first, points.back().second);
                     if (points.size() > 2){
-                        B_u = E->getDE_falling(order)->getValue(points[points.size()-2].first, points.back().second);
+                        B_u = -E->getDE_falling(order)->getValue(points[points.size()-2].first, points.back().second);
                     }
                     else{
-                        B_u = E->getDE_falling(order)->getValue(maxInputValIn, points.back().second);
-
+                        // start with falling
+                        B_u = -E->getDE_falling(order)->getValue(maxInputValIn, points.back().second);
                     }
                     
                     
@@ -1325,17 +1305,26 @@ namespace preisachScalar_ns{
 
     double preisach::calcMatDiffDiff(const bool bPilot){
 
-        preisach p(*this);
-        auto InS = this->E->quantiseInputVal((bPilot ? CHECK_E_FW_INV(pilotHBMatMatDiff.first, pilotHBMatMatDiff.second): inputVal));
-        InS = this->E->quantiseInputVal((rising ? InS.fourth : InS.second));
-        p.addForward((rising ? InS.fourth : InS.second));
+        // quadD pilotHBMatMatDiff_save = pilotHBMatMatDiff;
+        // double in_WP = (bPilot ? CHECK_E_FW_INV(pilotHBMatMatDiff.first, pilotHBMatMatDiff.second): inputVal_interpolated);
+        // quadD in_WP_s = this->E->quantiseInputVal(in_WP);
+        // double matDiff_WP = (bPilot ? pilotHBMatMatDiff.fourth : mat_diff);
+        // double in_pilot = (rising ? in_WP_s.fourth : in_WP_s.second);
 
-        if(conventionalMatDiff){            
-            return (p.getMatDiff() - prevDiffValInOutMatMatDiff.fourth)/(InS.fourth - prevDiffValInOutMatMatDiff.first);
-        }
-        else{
-            return calcMatDiff(InS, 2);
-        }
+        double matdiffdiff;
+        // if(conventionalMatDiff){            
+        //     UNUSED(pilotForward(in_pilot));
+        //     matdiffdiff = (pilotHBMatMatDiff.fourth - matDiff_WP)/(in_pilot - in_WP);
+            
+        //     pilotHBMatMatDiff = pilotHBMatMatDiff_save;
+        
+        // }
+        // else{
+        //     matdiffdiff = calcMatDiff(in_WP_s, 2);
+        //     cout << bPilot << " " << in_WP_s.first << " " << matdiffdiff << endl;
+        // }
+
+        return matdiffdiff;
     }
 
     // double preisach::calcMuDiff_secondOrder(const double dH){
@@ -1678,7 +1667,7 @@ namespace preisachScalar_ns{
         // all steps of the nonlin discreteInputValues
         points.clear();
         prevDiffValInOutMatMatDiff = quadD(E->getMaxInputVal(), E->getMaxOutputVal(), MU0, MU0); // not sure if MU0 is correct, same as some lines bellow
-        pilotHBMatMatDiff = quadD{};
+        pilot->demagnetise();
         
         rising = true;
 
@@ -1948,7 +1937,8 @@ namespace preisachScalar_ns{
         differential = b;
         if(b == false){
             mat_diff = 0;
-            pilotHBMatMatDiff.fourth = 0;
+            if(pilot != nullptr)
+                pilot->mat_diff = 0;
         }
     }
 
@@ -1994,9 +1984,6 @@ namespace preisachScalar_ns{
         return inputVal_interpolated;
     }
 
-    const quadD& preisach::getPilots() const{
-        return pilotHBMatMatDiff;
-    }
     // quadD& preisach::getPilots(){
     //     return pilotHBMatMatDiff;
     // }
@@ -2116,10 +2103,10 @@ namespace preisachScalar_ns{
         o << "conventionalMatDiff:\t" << pD.conventionalMatDiff << endl;
         o << "maxValueForOutside:\t" << pD.getMaxValueForOutside() << endl;
         o << "matDiffForOutside:\t" << pD.matDiffForOutside << endl;
-        o << "pilots:\n\tH: " << pD.getPilots().first << endl;
-        o << "\tB: " << pD.getPilots().second << endl;
-        o << "\tmat: " << pD.getPilots().third << endl;
-        o << "\tmat_diff: " << pD.getPilots().fourth << endl;
+        o << "pilots:\n\tH: " << pD.getPilot()->getH_interpolated() << endl;
+        o << "\tB: " << pD.getPilot()->getB_interpolated() << endl;
+        o << "\tmat: " << pD.getPilot()->getMat() << endl;
+        o << "\tmat_diff: " << pD.getPilot()->getMatDiff() << endl;
         o << "===========================" << endl;
 
         return o;
diff --git a/C++/preisach/preisachScalar/preisach.h b/C++/preisach/preisachScalar/preisach.h
index efe0a196a8ef5b138f285d70457cffd63227f703..5853f0f63a0a49f9ee60907a00cdb456815a9f21 100644
--- a/C++/preisach/preisachScalar/preisach.h
+++ b/C++/preisach/preisachScalar/preisach.h
@@ -63,7 +63,7 @@ namespace preisachScalar_ns{
             preisach(const preisach &&);
             preisach(const shared_ptr<preisach> &&);
             
-            preisach(const shared_ptr<EverettMatrix>& E, const bool bDemagnetise=true, const  bool interpolate=true, const bool differential=true, const bool conventionalMatDiff=true);        
+            preisach(const shared_ptr<EverettMatrix>& E, const bool bDemagnetise=true, const  bool interpolate=true, const bool differential=true, const bool conventionalMatDiff=true, const bool NoPilot=false);        
             preisach& operator= (const preisach&);
             preisach& operator= (const shared_ptr<preisach>&);
             
@@ -113,8 +113,6 @@ namespace preisachScalar_ns{
 
             [[nodiscard]] double getMat() const;
             [[nodiscard]] double getMatDiff() const;
-            [[nodiscard]] const quadD& getPilots() const;
-            // [[nodiscard]] quadD& getPilots(); 
 
 
             [[nodiscard]] double getMaxB() const{return CHECK_E_FW_INV(getMaxOutputVal(), getMaxInputVal());}
@@ -159,13 +157,13 @@ namespace preisachScalar_ns{
             void setConventionalMatDiff(const bool b=true);
 
 
-            quadD pilotH(double h_pilot) const {return CHECK_E_FW_INV(pilotForward(h_pilot), pilotInverse(h_pilot));} //{H, B, mu, mu_diff}
+            shared_ptr<preisach> pilotH(double h_pilot) const {return CHECK_E_FW_INV(pilotForward(h_pilot), pilotInverse(h_pilot));} //{H, B, mu, mu_diff}
             double pilotH_Bonly(double h_pilot) const {return CHECK_E_FW_INV(pilotForward_Bonly(h_pilot), pilotInverse_Honly(h_pilot));}  //B
-            quadD pilotB(double b_pilot) const  {return CHECK_E_FW_INV(pilotInverse(b_pilot), pilotForward(b_pilot));}//{H, B, mu, mu_diff}
+            shared_ptr<preisach> pilotB(double b_pilot) const  {return CHECK_E_FW_INV(pilotInverse(b_pilot), pilotForward(b_pilot));}//{H, B, mu, mu_diff}
             double pilotB_Honly(double b_pilot) const {return CHECK_E_FW_INV(pilotInverse_Honly(b_pilot), pilotForward_Bonly(b_pilot));} //B
 
-            [[nodiscard]] quadD pilotForward(double h_pilot) const ; //{H, B, mu, mu_diff}
-            [[nodiscard]] quadD pilotInverse(double b_pilot) const ;
+            [[nodiscard]] shared_ptr<preisach> pilotForward(double h_pilot) const ; //{H, B, mu, mu_diff}
+            [[nodiscard]] shared_ptr<preisach> pilotInverse(double b_pilot) const ;
 
             
             // KL
@@ -211,6 +209,8 @@ namespace preisachScalar_ns{
             [[nodiscard]] std::shared_ptr<Everett_Empty> getInverseEverettMatrix(const std::vector<double> vecB) const;
             [[nodiscard]] std::shared_ptr<Everett_Empty> getInverseEverettMatrix() const;
             
+
+            shared_ptr<preisach> getPilot() const;
             [[nodiscard]] double calcMatDiffDiff(const bool bPilot);
 
 
@@ -243,17 +243,18 @@ namespace preisachScalar_ns{
             double mat = MU0;
             double mat_diff = MU0;
             bool conventionalMatDiff;
-
+            const bool isPilot = false;
 
 
             // history
             PointList_t points;      // M, m, B 
             
-
-            quadD pilotHBMatMatDiff;  
+            // quadD pilotHBMatMatDiff;  
             quadD prevInOutMatMatDiffRising;
             quadD prevInOutMatMatDiffFalling;
-            
+
+            shared_ptr<preisach> pilot = nullptr;
+
 
             [[nodiscard]] double interpolateOutputVal(const double) const; // returns B and dMu
             [[nodiscard]] double interpolateInputVal(const double) const; // returns B and dMu
@@ -304,8 +305,8 @@ namespace preisachScalar_ns{
 
         
             [[nodiscard]] double pilotInverse_Honly(double b_pilot) const ;
-            [[nodiscard]] double  pilotInverse_quantised(double b_pilot) const;
-            [[nodiscard]] double  pilotInverse_cont(double b_pilot, double eps=0) const;
+            [[nodiscard]] double pilotInverse_quantised(double b_pilot) const;
+            [[nodiscard]] double pilotInverse_cont(double b_pilot, double eps=0) const;
             [[nodiscard]] double pilotInverse_BS(double b_pilot) const ;
 
             [[nodiscard]] double getMaxInputVal() const;
diff --git a/C++/preisach/preisachScalar/unit_test.cc b/C++/preisach/preisachScalar/unit_test.cc
index 7c2a13ba257ac25a7e797ca3e83de8d6c77c00b7..887d43a511721f00bd0adaad525c1dac171fddbc 100644
--- a/C++/preisach/preisachScalar/unit_test.cc
+++ b/C++/preisach/preisachScalar/unit_test.cc
@@ -298,9 +298,9 @@ bool preisachScalar_ns::unit_test_specificEverett(std::shared_ptr<EverettMatrix>
         double nu_old = pA.getNu();
         double muDiff_old = pA.getMuDiff();
         double nuDiff_old = pA.getNuDiff();
-        assert_equal(pA.pilotForward(ptrE->getMaxInputVal()).second, maxB);
+        assert_equal(pA.pilotForward(ptrE->getMaxInputVal())->getOutputVal_interpolated(), maxB);
         // assert_equal(B_old, pA.getB());
-        assert_almostequal(pA.pilotInverse(maxB).first, ptrE->getMaxInputVal(), 1e-6);
+        assert_almostequal(pA.pilotInverse(maxB)->getInputVal_interpolated(), ptrE->getMaxInputVal(), 1e-6);
         assert_almostequal(OutputVal_old, pA.getOutputVal(), 1e-6);
         assert_almostequal(OutputVal_old_interpolatd, pA.getOutputVal_interpolated(), 1e-6);
         assert_almostequal(mu_old, pA.getMu(), 1e-6);
@@ -315,13 +315,13 @@ bool preisachScalar_ns::unit_test_specificEverett(std::shared_ptr<EverettMatrix>
     B_old = pA.getOutputVal();
     
     
-    assert_equal(pA.pilotForward(ptrE->getMaxInputVal()).second, ptrE->getMaxB());
-    assert_notequal(pA.pilotForward(ptrE->getMaxInputVal()).third, pA.getMat());
-    assert_notequal(pA.pilotForward(ptrE->getMaxInputVal()*0.5).fourth, pA.getMatDiff());
+    assert_equal(pA.pilotForward(ptrE->getMaxInputVal())->getOutputVal_interpolated(), ptrE->getMaxB());
+    assert_notequal(pA.pilotForward(ptrE->getMaxInputVal())->getMat(), pA.getMat());
+    assert_notequal(pA.pilotForward(ptrE->getMaxInputVal()*0.5)->getMatDiff(), pA.getMatDiff());
     assert_equal(B_old, pA.getOutputVal());
-    assert_almostequal(pA.pilotInverse(ptrE->getMaxOutputVal()).first, ptrE->getMaxH(), 1);
-    assert_notequal(pA.pilotInverse(ptrE->getMaxOutputVal()).third, pA.getMat());
-    assert_notequal(pA.pilotInverse(ptrE->getMaxOutputVal()*0.5).fourth, pA.getMatDiff());
+    assert_almostequal(pA.pilotInverse(ptrE->getMaxOutputVal())->getInputVal_interpolated(), ptrE->getMaxH(), 1);
+    assert_notequal(pA.pilotInverse(ptrE->getMaxOutputVal())->getMat(), pA.getMat());
+    assert_notequal(pA.pilotInverse(ptrE->getMaxOutputVal()*0.5)->getMatDiff(), pA.getMatDiff());
     assert_equal(B_old, pA.getOutputVal());
 
     
@@ -404,7 +404,7 @@ bool preisachScalar_ns::unit_test_specificEverett(std::shared_ptr<EverettMatrix>
     cout << "check if pilotH functions has same result as addH function" << endl;
     pA.demagnetise();
     for(i = 0; i!= 100; ++i){
-        double B = pA.pilotH(i*1.0/99.0*maxH).second;
+        double B = pA.pilotH(i*1.0/99.0*maxH)->getOutputVal_interpolated();
         pA.addH(i*1.0/99.0*maxH);
         assert_almostequal(B, pA.getB_interpolated(), 1e-9);
     }
@@ -412,7 +412,7 @@ bool preisachScalar_ns::unit_test_specificEverett(std::shared_ptr<EverettMatrix>
     cout << "check if pilotB functions has same result as addB function" << endl;
     pA.demagnetise();
     for(i = 0; i!= 100; ++i){
-        double H = pA.pilotB(i*1.0/99.0*maxB).first;
+        double H = pA.pilotB(i*1.0/99.0*maxB)->getInputVal_interpolated();
         pA.addB(i*1.0/99.0*maxB);
 
         assert_almostequal(H, pA.getH_interpolated(), 3e-3);
@@ -432,15 +432,15 @@ bool preisachScalar_ns::unit_test_specificEverett(std::shared_ptr<EverettMatrix>
         double H_2 = KL_data[KL_data.size()/2 + 6].first;
         double H_3 = KL_data[KL_data.size()/2 + 8].first;
 
-        assert_equal(p.pilotH(H_1).second, KL->EvaluateH(H_1));
+        assert_equal(p.pilotH(H_1)->getOutputVal_interpolated(), KL->EvaluateH(H_1));
         p.addH(H_2);
-        assert_notequal(p.pilotH(H_1).second, KL->EvaluateH(H_1));
+        assert_notequal(p.pilotH(H_1)->getOutputVal_interpolated(), KL->EvaluateH(H_1));
         p.updateKL();
-        assert_equal(p.pilotH(H_1).second, KL->EvaluateH(H_1));
+        assert_equal(p.pilotH(H_1)->getOutputVal_interpolated(), KL->EvaluateH(H_1));
         p.demagnetise();
         p.setAutoUpdateKL(true);
         p.addH(H_3);
-        assert_equal(p.pilotH(H_1).second, KL->EvaluateH(H_1));
+        assert_equal(p.pilotH(H_1)->getOutputVal_interpolated(), KL->EvaluateH(H_1));
 
     }
 
diff --git a/C++/preisach/preisachVector/preisachVector.cc b/C++/preisach/preisachVector/preisachVector.cc
index 10a545f5cc473034fadd14939a9ec0c71005ec65..63e50de9142258d32a5b9f26c45651809635814a 100644
--- a/C++/preisach/preisachVector/preisachVector.cc
+++ b/C++/preisach/preisachVector/preisachVector.cc
@@ -96,7 +96,7 @@ namespace preisachVector_ns{
     // Constructor
     template<unsigned int DIM>
     preisachVector<DIM>::preisachVector(const shared_ptr<EverettMatrix>& E, const shared_ptr<distributions>& dist,
-                    const bool interpolate, const bool differential, const std::vector<double> rotationalParameters):E{E}, dist{dist}, rotationalParameters{ rotationalParameters}{
+                    const bool interpolate, const bool differential, const std::vector<double> rotationalParameters, bool isPilot):E{E}, dist{dist}, rotationalParameters{ rotationalParameters}, isPilot{isPilot}{
         
         if(dist == nullptr || E == nullptr){
             throw std::invalid_argument("dist or E are invalid shared pointers");
@@ -130,6 +130,10 @@ namespace preisachVector_ns{
         calcMat();
         calcMatDiff();
 
+        if(!isPilot){
+            pilot = make_shared<preisachVector>(E, dist, interpolate, differential, rotationalParameters, true);
+        }
+
     }
 
     // copy const
@@ -204,7 +208,7 @@ namespace preisachVector_ns{
         mat_diff = o.getMatDiff();
         mat = o.getMat();
 
-        pilotHBMuMuDiff = o.getPilots();
+        
         prevDiffValInOutMatMatDiff = o.getPrevDiffValInOutMatMatDiff();
 
         
@@ -214,6 +218,12 @@ namespace preisachVector_ns{
             this->push_back(make_shared<preisachScalar>(o[i]));        
         }
 
+        isPilot = o.isPilot;
+
+        if(!isPilot){
+            pilot = make_shared<preisachVector>(E, dist, (*this)[0]->getInterpolate(), (*this)[0]->getDifferential(), rotationalParameters, true);
+        }
+
         return *this;
     }
 
@@ -285,189 +295,61 @@ namespace preisachVector_ns{
     }
 
     template<unsigned int DIM>
-    const vec_t<double, DIM> preisachVector<DIM>::calcOutputValue( const bool bPilot){
-        if(!bPilot){
-            calcOutputValue(OutputValue, false);
-
-            return OutputValue;
-        }
-        else{
-            vec_t<double, DIM> Output_pilot;
-            calcOutputValue(Output_pilot, true);
-            return Output_pilot;
-        }
-    }
-    template<unsigned int DIM>
-    void preisachVector<DIM>::calcOutputValue(vec_t<double, DIM>& B, const bool bPilot) const{
+    const vec_t<double, DIM> preisachVector<DIM>::calcOutputValue(){
         
-        B.reset();
+        OutputValue.reset();
         vec_t<double, DIM> tmp;
         unsigned int  i;
 
         // normal use
-        if(!bPilot){
-            // if constexpr(DIM == 2){
-                
-            //     vec_t<double, DIM> tmp_o;
-
-            //     for(i = 0; i != this->size(); ++i){
-            //         tmp_o=tmp;
-
-            //         tmp.first() = (*this)[i]->getXYZW().x * (*this)[i]->getXYZW().weight * (*this)[i]->getOutputVal_interpolated();
-            //         tmp.second() = (*this)[i]->getXYZW().y* (*this)[i]->getXYZW().weight * (*this)[i]->getOutputVal_interpolated();
-
-            //         if(i == 0 || i == this->size() -1){
-            //             B += tmp;
-            //         }
-            //         else{
-            //             B += 0.5*(tmp + tmp_o);
-            //         }
-            //     }
-                
-            // }
-            // else{
-                // cout << "add" << endl;
-            for(i = 0; i != this->size(); ++i){
-
-                tmp.first() = (*this)[i]->getXYZW().x * (*this)[i]->getXYZW().weight * (*this)[i]->getOutputVal_interpolated();
-                if constexpr(DIM > 1){
-                    tmp.second() = (*this)[i]->getXYZW().y* (*this)[i]->getXYZW().weight * (*this)[i]->getOutputVal_interpolated();
-                }
-                if constexpr(DIM > 2){
-                    tmp.third() = (*this)[i]->getXYZW().z * (*this)[i]->getXYZW().weight * (*this)[i]->getOutputVal_interpolated();
-                }
-                B += tmp;
+        
+        for(i = 0; i != this->size(); ++i){
 
-                // cout << i << " " << (*this)[i]->getOutputVal_interpolated() << " " <<  (*this)[i]->getPilots().second << endl;
+            tmp.first() = (*this)[i]->getXYZW().x * (*this)[i]->getXYZW().weight * (*this)[i]->getOutputVal_interpolated();
+            if constexpr(DIM > 1){
+                tmp.second() = (*this)[i]->getXYZW().y* (*this)[i]->getXYZW().weight * (*this)[i]->getOutputVal_interpolated();
             }
-        }
-        // use Perfect input: pilot Function
-        else{
-            // cout << "pilot" << endl;
-            for(i = 0; i != this->size(); ++i){
-
-                if(!E->getIsInverse()){
-                    tmp.first() = (*this)[i]->getXYZW().x * (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().second;
-                    if constexpr(DIM > 1){
-                        tmp.second() = (*this)[i]->getXYZW().y* (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().second;
-                    }
-                    if constexpr(DIM > 2){
-                        tmp.third() = (*this)[i]->getXYZW().z * (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().second;
-                    }
-                }
-                else{
-                    tmp.first() = (*this)[i]->getXYZW().x * (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().first;
-                    if constexpr(DIM > 1){
-                        tmp.second() = (*this)[i]->getXYZW().y* (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().first;
-                    }
-                    if constexpr(DIM > 2){
-                        tmp.third() = (*this)[i]->getXYZW().z * (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().first;
-                    }
-                }
-
-                if(std::isnan(tmp.first()) || std::isnan(tmp.second()) || std::isnan(tmp.third())){
-                        cout << __PRETTY_FUNCTION__ << endl;
-                        cout << __FILE__ << ":" << __LINE__<< " calc B is not good"  <<endl;
-                        cout << B << endl;
-
-                        cout << i << endl;
-                        cout << tmp << endl;
-                        cout << (*this)[i]->getPilots().first << endl;
-                        cout << (*this)[i]->getPilots().second << endl;
-                        throw std::runtime_error("in pilot Function for PV: Bx or By or Bz is nan");
-                }
-
-                B += tmp;
-
-                // cout << i << " " << (*this)[i]->getB() <<  " "<< (*this)[i]->getPilots().second << endl;
+            if constexpr(DIM > 2){
+                tmp.third() = (*this)[i]->getXYZW().z * (*this)[i]->getXYZW().weight * (*this)[i]->getOutputVal_interpolated();
             }
+            OutputValue += tmp;
+
+            // cout << i << " " << (*this)[i]->getOutputVal_interpolated() << " " <<  (*this)[i]->getPilots().second << endl;
         }
+        
+        OutputValue *= normOutputValue;
 
-        B *= normOutputValue;
+        return OutputValue;
     }
 
 
     template<unsigned int DIM>
-    const vec_t<double, DIM> preisachVector<DIM>::calcInputValue( const bool bPilot) {
-        
-        if(!bPilot){
-            calcInputValue(InputValue, bPilot);
-            return InputValue;
-        }
-        else{
-            vec_t<double, DIM> Input_pilot;
-            calcInputValue(Input_pilot, bPilot);
-            return Input_pilot;
-        }
-    }
-    template<unsigned int DIM>
-    void preisachVector<DIM>::calcInputValue( vec_t<double, DIM>& H, const bool bPilot) const{
+    const vec_t<double, DIM> preisachVector<DIM>::calcInputValue(){
         
 
-        H.reset();
+        InputValue.reset();
         vec_t<double, DIM> tmp;
         unsigned int  i;
-
+        
         // normal use
-        if(!bPilot){
-            for(i = 0; i != this->size(); ++i){
-                tmp.first() = (*this)[i]->getXYZW().x * (*this)[i]->getXYZW().weight * (*this)[i]->getInputVal_interpolated();
-                if constexpr(DIM > 1){
-                    tmp.second() = (*this)[i]->getXYZW().y* (*this)[i]->getXYZW().weight * (*this)[i]->getInputVal_interpolated();
-                }
-                if constexpr(DIM > 2){
-                    tmp.third() = (*this)[i]->getXYZW().z * (*this)[i]->getXYZW().weight * (*this)[i]->getInputVal_interpolated();
-                }
-                H += tmp;
+        for(i = 0; i != this->size(); ++i){
+            tmp.first() = (*this)[i]->getXYZW().x * (*this)[i]->getXYZW().weight * (*this)[i]->getInputVal_interpolated();
+            if constexpr(DIM > 1){
+                tmp.second() = (*this)[i]->getXYZW().y* (*this)[i]->getXYZW().weight * (*this)[i]->getInputVal_interpolated();
             }
-        }
-        // use Perfect input: pilot Function
-        else{
-            for(i = 0; i != this->size(); ++i){
-
-
-                if(!E->getIsInverse()){
-                    tmp.first() = (*this)[i]->getXYZW().x * (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().first;
-                    if constexpr(DIM > 1){
-                        tmp.second() = (*this)[i]->getXYZW().y* (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().first;
-                    }
-                    if constexpr(DIM > 2){
-                        tmp.third() = (*this)[i]->getXYZW().z * (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().first;
-                    }
-                }
-                else{
-                    tmp.first() = (*this)[i]->getXYZW().x * (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().second;
-                    if constexpr(DIM > 1){
-                        tmp.second() = (*this)[i]->getXYZW().y* (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().second;
-                    }
-                    if constexpr(DIM > 2){
-                        tmp.third() = (*this)[i]->getXYZW().z * (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().second;
-                    }
-                }
-                
-                H += tmp;
+            if constexpr(DIM > 2){
+                tmp.third() = (*this)[i]->getXYZW().z * (*this)[i]->getXYZW().weight * (*this)[i]->getInputVal_interpolated();
             }
+            InputValue += tmp;
         }
 
-        H *= normInputValue;
-    }
-
-    template<unsigned int DIM>
-    const mat_t<double, DIM, DIM> preisachVector<DIM>::calcMat(const bool bPilot){
-        if(!bPilot){
-            calcMat(mat, bPilot);
-            return mat;
-        }
-        else{
-            mat_t<double, DIM, DIM> mat_pilot;
-            calcMat(mat_pilot, bPilot);
-            return mat_pilot;
-        }
-        
+        InputValue *= normInputValue;
 
+        return InputValue;
     }
+
     template<unsigned int DIM>
-    void preisachVector<DIM>::calcMat(mat_t<double, DIM, DIM>& mat,  const bool bPilot) const {
+    const mat_t<double, DIM, DIM>  preisachVector<DIM>::calcMat() {
         
 
         mat_t<double, DIM, DIM> tmp;
@@ -476,108 +358,54 @@ namespace preisachVector_ns{
         // double phi, theta;
         quadD R;
         // normal use
-        if(!bPilot){
-            if constexpr(DIM == 3){
-                for(i = 0; i != this->size(); ++i){
-                    R = (*this)[i]->getXYZW();
-
-                    tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.x*R.x;
-                    tmp.xy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.x*R.y;
-                    tmp.xz() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.x*R.z;
-                    tmp.yy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.y*R.y;
-                    tmp.yz() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.y*R.z;
-                    tmp.zz() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.z*R.z;
-                    tmp.yx() = tmp.xy();
-                    tmp.zx() = tmp.xz();
-                    tmp.zy() = tmp.yz();
-                    
-                    mat += tmp;
-                }
-            }
-            else if constexpr(DIM == 2){
-                for(i = 0; i != this->size(); ++i){
-                    R = (*this)[i]->getXYZW();
+        if constexpr(DIM == 3){
+            for(i = 0; i != this->size(); ++i){
+                R = (*this)[i]->getXYZW();
 
-                    tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.x*R.x;
-                    tmp.xy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.x*R.y;
-                    tmp.yy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.y*R.y;
-                    tmp.yx() = tmp.xy();
+                tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.x*R.x;
+                tmp.xy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.x*R.y;
+                tmp.xz() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.x*R.z;
+                tmp.yy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.y*R.y;
+                tmp.yz() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.y*R.z;
+                tmp.zz() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.z*R.z;
+                tmp.yx() = tmp.xy();
+                tmp.zx() = tmp.xz();
+                tmp.zy() = tmp.yz();
+                
+                mat += tmp;
+            }
+        }
+        else if constexpr(DIM == 2){
+            for(i = 0; i != this->size(); ++i){
+                R = (*this)[i]->getXYZW();
 
-                    mat += tmp;
-                }
+                tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.x*R.x;
+                tmp.xy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.x*R.y;
+                tmp.yy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.y*R.y;
+                tmp.yx() = tmp.xy();
 
+                mat += tmp;
             }
-            else{
-                for(i = 0; i != this->size(); ++i){
-                    R = (*this)[i]->getXYZW();
-                    tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.x*R.x;
-                                
-                    mat += tmp;
-                }
 
-            }
         }
-        // use Perfect input: pilot Function
         else{
-            if constexpr(DIM == 3){
-                for(i = 0; i != this->size(); ++i){
-
-                    R = (*this)[i]->getXYZW();
-
-                    tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().third * R.x*R.x;
-                    tmp.xy() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().third * R.x*R.y;
-                    tmp.xz() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().third * R.x*R.z;
-                    tmp.yy() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().third * R.y*R.y;
-                    tmp.yz() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().third * R.y*R.z;
-                    tmp.zz() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().third * R.z*R.z;
-                    tmp.yx() = tmp.xy();
-                    tmp.zx() = tmp.xz();
-                    tmp.zy() = tmp.yz();
-                    
-                    mat += tmp;
-                    
-                }
+            for(i = 0; i != this->size(); ++i){
+                R = (*this)[i]->getXYZW();
+                tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getMat() * R.x*R.x;
+                            
+                mat += tmp;
             }
-            if constexpr(DIM == 2){
-                for(i = 0; i != this->size(); ++i){
 
-                    R = (*this)[i]->getXYZW();
-
-                    tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().third * R.x*R.x;
-                    tmp.xy() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().third * R.x*R.y;
-                    tmp.yy() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().third * R.y*R.y;
-                    tmp.yx() = tmp.xy();                
-                    mat += tmp;
-                    
-                }
-            }
-            if constexpr(DIM == 1){
-                for(i = 0; i != this->size(); ++i){
-                    R = (*this)[i]->getXYZW();
-                    tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().third * R.x*R.x;
-                    mat += tmp;
-                    
-                }
-            }
         }
-
+        
         mat *= (normOutputValue);
-    }
 
-    template<unsigned int DIM>
-    const mat_t<double, DIM, DIM> preisachVector<DIM>::calcMatDiff(const bool bPilot){
-        if(!bPilot){
-            calcMatDiff(mat_diff, bPilot);
-            return mat_diff;
-        }
-        else{
-            mat_t<double, DIM, DIM> mat_diff_pilot;
-            calcMatDiff(mat_diff_pilot, bPilot);
-            return mat_diff_pilot;
-        }
+        return mat;
     }
+
+
     template<unsigned int DIM>
-    void preisachVector<DIM>::calcMatDiff(mat_t<double, DIM, DIM>& mat_diff,  const bool bPilot) const{
+    const mat_t<double, DIM, DIM>  preisachVector<DIM>::calcMatDiff(){
         
 
         mat_t<double, DIM, DIM> tmp;
@@ -586,104 +414,54 @@ namespace preisachVector_ns{
         // double phi, theta;
         quadD R;
         // normal use
-        if(!bPilot){
-            if constexpr(DIM == 3){
-                for(i = 0; i != this->size(); ++i){
-                    R = (*this)[i]->getXYZW();
-
-                    tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.x*R.x;
-                    tmp.xy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.x*R.y;
-                    tmp.xz() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.x*R.z;
-                    tmp.yy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.y*R.y;
-                    tmp.yz() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.y*R.z;
-                    tmp.zz() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.z*R.z;
-                    tmp.yx() = tmp.xy();
-                    tmp.zx() = tmp.xz();
-                    tmp.zy() = tmp.yz();
-                    
-                    mat_diff += tmp;
-                }
-            }
-            else if constexpr(DIM == 2){
-                for(i = 0; i != this->size(); ++i){
-                    R = (*this)[i]->getXYZW();
+        if constexpr(DIM == 3){
+            for(i = 0; i != this->size(); ++i){
+                R = (*this)[i]->getXYZW();
 
-                    tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.x*R.x;
-                    tmp.xy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.x*R.y;
-                    tmp.yy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.y*R.y;
-                    tmp.yx() = tmp.xy();
+                tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.x*R.x;
+                tmp.xy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.x*R.y;
+                tmp.xz() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.x*R.z;
+                tmp.yy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.y*R.y;
+                tmp.yz() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.y*R.z;
+                tmp.zz() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.z*R.z;
+                tmp.yx() = tmp.xy();
+                tmp.zx() = tmp.xz();
+                tmp.zy() = tmp.yz();
+                
+                mat_diff += tmp;
+            }
+        }
+        else if constexpr(DIM == 2){
+            for(i = 0; i != this->size(); ++i){
+                R = (*this)[i]->getXYZW();
 
-                    mat_diff += tmp;
-                }
+                tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.x*R.x;
+                tmp.xy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.x*R.y;
+                tmp.yy() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.y*R.y;
+                tmp.yx() = tmp.xy();
 
+                mat_diff += tmp;
             }
-            else{
-                for(i = 0; i != this->size(); ++i){
-                    R = (*this)[i]->getXYZW();
-                    tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.x*R.x;
-                                
-                    mat_diff += tmp;
-                }
 
-            }
         }
-        // use Perfect input: pilot Function
         else{
-            if constexpr(DIM == 3){
-                for(i = 0; i != this->size(); ++i){
-
-                    R = (*this)[i]->getXYZW();
-
-                    tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().fourth * R.x*R.x;
-                    tmp.xy() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().fourth * R.x*R.y;
-                    tmp.xz() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().fourth * R.x*R.z;
-                    tmp.yy() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().fourth * R.y*R.y;
-                    tmp.yz() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().fourth * R.y*R.z;
-                    tmp.zz() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().fourth * R.z*R.z;
-                    tmp.yx() = tmp.xy();
-                    tmp.zx() = tmp.xz();
-                    tmp.zy() = tmp.yz();
-                    
-                    mat_diff += tmp;
-                    
-                }
+            for(i = 0; i != this->size(); ++i){
+                R = (*this)[i]->getXYZW();
+                tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getMatDiff() * R.x*R.x;
+                            
+                mat_diff += tmp;
             }
-            if constexpr(DIM == 2){
-                for(i = 0; i != this->size(); ++i){
 
-                    R = (*this)[i]->getXYZW();
-
-                    tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().fourth * R.x*R.x;
-                    tmp.xy() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().fourth * R.x*R.y;
-                    tmp.yy() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().fourth * R.y*R.y;
-                    tmp.yx() = tmp.xy();                
-                    mat_diff += tmp;
-                    
-                }
-            }
-            if constexpr(DIM == 1){
-                for(i = 0; i != this->size(); ++i){
-                    R = (*this)[i]->getXYZW();
-                    tmp.xx() = (*this)[i]->getXYZW().weight * (*this)[i]->getPilots().fourth * R.x*R.x;
-                    mat_diff += tmp;
-                    
-                }
-            }
         }
-
+        
         mat_diff *= (normOutputValue);
-    }
 
-    template<unsigned int DIM>
-    const mat_t<double, DIM, DIM> preisachVector<DIM>::calcMatDiffDiff(const bool bPilot){
-        mat_t<double, DIM, DIM> tmp;
-        calcMatDiffDiff(tmp, bPilot);
-        return tmp;
+        return mat_diff;
     }
     
     template<unsigned int DIM>
-    void preisachVector<DIM>::calcMatDiffDiff(mat_t<double, DIM, DIM>& mat_diffdiff,  const bool bPilot) const{
-
+    const mat_t<double, DIM, DIM> preisachVector<DIM>::calcMatDiffDiff(){
+        mat_t<double, DIM, DIM> mat_diffdiff;
         mat_t<double, DIM, DIM> tmp;
         double mat_diffdiff_scalar;
         mat_diffdiff.clear();
@@ -691,6 +469,7 @@ namespace preisachVector_ns{
         // double phi, theta;
         quadD R;
         // normal use
+        bool bPilot = true;
 
         if constexpr(DIM == 3){
             for(i = 0; i != this->size(); ++i){
@@ -735,6 +514,8 @@ namespace preisachVector_ns{
         
 
         mat_diffdiff *= (normOutputValue);
+
+        return mat_diffdiff;
     }
 
     template<unsigned int DIM>
@@ -813,7 +594,7 @@ namespace preisachVector_ns{
         if(iRec >= iRecMax){
             return std::make_pair(PhiRecursive, pilotIn);
         }
-        vec_t<double, DIM> pilotOut = CHECK_E_FW_INV(pilotForward(pilotIn).second, pilotForward(pilotIn).first); 
+        vec_t<double, DIM> pilotOut = pilotForward(pilotIn)->getOutputVal();
 
         double PhiRem = 0;
         
@@ -913,28 +694,22 @@ namespace preisachVector_ns{
     const mat_t<double, DIM, DIM>& preisachVector<DIM>::getMat() const{
         return mat;
     }
-    template<unsigned int DIM>
-    vec_t<double, DIM> preisachVector<DIM>::addForward(const vec_t<double, DIM> H){
-        return addForward(H, false);
-    }
+    
     template<unsigned int DIM>
     vec_t<double, DIM> preisachVector<DIM>::addForward(double H_x, double H_y, double H_z){
         // cout << H_x << " " << H_y << " " << H_z << endl;
-        return addForward(vec_t<double, DIM>{H_x, H_y, H_z}, false);
+        return addForward(vec_t<double, DIM>{H_x, H_y, H_z});
     }
 
     template<unsigned int DIM>
     vec_t<double, DIM> preisachVector<DIM>::addForward(const std::initializer_list<double> H){
         vec_t<double, DIM> o(H);
-        return addForward(o, false);
-    }
-    template<unsigned int DIM>
-    vec_t<double, DIM> preisachVector<DIM>::addForward(const std::initializer_list<double> H, bool bPilot){
-        vec_t<double, DIM> o(H);
-        return addForward(o, bPilot);
+        return addForward(o);
     }
+    
     template<unsigned int DIM>
-    vec_t<double, DIM> preisachVector<DIM>::addForward(const vec_t<double, DIM> newVal_in, bool bPilot){
+    vec_t<double, DIM> preisachVector<DIM>::addForward(const vec_t<double, DIM> newVal_in){
+        
 
         unsigned int i;
         double alpha ;
@@ -1003,7 +778,7 @@ namespace preisachVector_ns{
             
             iend = static_cast<unsigned int>(this->size()* (i+1.0)/numThreads );
             // cout << istart << ":" << iend << endl;
-            threads_update.push_back(std::thread(&preisachVector<DIM>::addForward_thread, this, newVal, bPilot, istart, iend));
+            threads_update.push_back(std::thread(&preisachVector<DIM>::addForward_thread, this, newVal, istart, iend));
             istart = iend ;
         }
 
@@ -1014,29 +789,16 @@ namespace preisachVector_ns{
         }
 
         // cout << "calc the stuff" << endl;        
-        if(!bPilot){
-            InputValue = newVal_in;
-            // calcInputValue();
-            calcOutputValue();
-            calcMat();
-            calcMatDiff();
 
-            return OutputValue;
-        }
-        else{
+        InputValue = newVal_in;
+        // calcInputValue();
+        calcOutputValue();
+        calcMat();
+        calcMatDiff();
+
+        return OutputValue;
+
 
-            if(E->getIsInverse()){
-                pilotHBMuMuDiff.first = newVal_in;
-            }
-            else{
-                pilotHBMuMuDiff.second = newVal_in;
-            }
-            // calcInputValue(true);
-            calcMat(true);
-            calcMatDiff(true);
-            
-            return calcOutputValue(true);
-        }
 
         // cout << "leave addForward" << endl;
         
@@ -1045,7 +807,7 @@ namespace preisachVector_ns{
     }
 
     template<unsigned int DIM>
-    void preisachVector<DIM>::addForward_thread(const vec_t<double, DIM> newVal, bool bPilot, const unsigned int istart, const unsigned int iend){
+    void preisachVector<DIM>::addForward_thread(const vec_t<double, DIM> newVal, const unsigned int istart, const unsigned int iend){
         unsigned int i;
         double in_axis;
         for(i = istart; i != iend; ++i){
@@ -1059,28 +821,22 @@ namespace preisachVector_ns{
                 in_axis += (*this)[i]->getXYZW().z * newVal.third();
             }            
 
-            if(!bPilot){
-                (*this)[i]->addForward(in_axis);
-            }
-            else {
-                // saved in pilots
-                UNUSED((*this)[i]->pilotForward(in_axis));
-            }
+            (*this)[i]->addForward(in_axis);
 
         }
     }
 
     template<unsigned int DIM>
-    vec_t<double, DIM> preisachVector<DIM>::addInverse(const std::initializer_list<double> B, bool bPilot){
+    vec_t<double, DIM> preisachVector<DIM>::addInverse(const std::initializer_list<double> B){
         vec_t<double, DIM> tD(B);
-        return addInverse(tD, bPilot);
+        return addInverse(tD);
     }
     template<unsigned int DIM>
     vec_t<double, DIM> preisachVector<DIM>::addInverse(double B_x, double B_y, double B_z){
-        return addInverse(vec_t<double, DIM>{B_x, B_y, B_z}, false);
+        return addInverse(vec_t<double, DIM>{B_x, B_y, B_z});
     }
     template<unsigned int DIM>
-    vec_t<double, DIM> preisachVector<DIM>::addInverse(const vec_t<double, DIM>& newVal, bool bPilot){
+    vec_t<double, DIM> preisachVector<DIM>::addInverse(const vec_t<double, DIM>& newVal){
         cout << __PRETTY_FUNCTION__ << endl;
         cout << "Inverse mode not correctly implemented!!!!!" << endl;
 
@@ -1095,7 +851,7 @@ namespace preisachVector_ns{
         for(i = 0; i!= numThreads; ++i){
             
             iend = static_cast<unsigned int>(this->size()* (i+1.0)/numThreads );
-            threads_update.push_back(std::thread(&preisachVector<DIM>::addForward_thread, this, newVal, bPilot, istart, iend));
+            threads_update.push_back(std::thread(&preisachVector<DIM>::addForward_thread, this, newVal, istart, iend));
             istart = iend ;
         }
 
@@ -1103,24 +859,20 @@ namespace preisachVector_ns{
             threads_update[i].join();
         }
 
-        if(!bPilot){
-            OutputValue = newVal;
-            calcInputValue();
-            calcMat();
-            calcMatDiff();
 
-            return InputValue;
-        }
-        else{
-            calcMat(true);
-            calcMatDiff(true);
-            return calcInputValue(true);
-        }
+        OutputValue = newVal;
+        calcInputValue();
+        calcMat();
+        calcMatDiff();
+
+        return InputValue;
+
+
     }
 
 
     template<unsigned int DIM>
-    void preisachVector<DIM>::addInverse_thread(const vec_t<double, DIM> newVal, bool bPilot, const unsigned int istart, const unsigned int iend){
+    void preisachVector<DIM>::addInverse_thread(const vec_t<double, DIM> newVal, const unsigned int istart, const unsigned int iend){
         unsigned int i;
         double in_axis;
         for(i = istart; i != iend; ++i){
@@ -1134,37 +886,30 @@ namespace preisachVector_ns{
                 in_axis += (*this)[i]->getXYZW().z * newVal.third();
             }            
 
-            if(!bPilot){
-                (*this)[i]->addInverse(in_axis);
-            }
-            else {
-                // saved in pilots
-                UNUSED((*this)[i]->pilotInverse(in_axis));
-            }
+
+            (*this)[i]->addInverse(in_axis);
+
 
         }
     }
 
     template<unsigned int DIM>
-    pVpilots_t<DIM> preisachVector<DIM>::pilotForward(const vec_t<double, DIM>& HorB){
-        vec_t<double, DIM> BorH = addForward(HorB, true);
-
-        calcMat(pilotHBMuMuDiff.third, true);
-
-        calcMatDiff(pilotHBMuMuDiff.fourth, true);
-
+    void preisachVector<DIM>::updatePilot(){
+        (*pilot) = *this;
+    }
 
+    template<unsigned int DIM>
+    shared_ptr<preisachVector<DIM>> preisachVector<DIM>::pilotForward(const vec_t<double, DIM>& HorB){
         
-        pilotHBMuMuDiff.first = CHECK_E_FW_INV(HorB, BorH);
-        pilotHBMuMuDiff.second = CHECK_E_FW_INV(BorH, HorB);
-        
-        
+        updatePilot();
+        vec_t<double, DIM> BorH = pilot->addForward(HorB);
+
         //cout << "current H " << preisachVector::H.first << "current mat " << preisachVector::mat.xx << "current B" << preisachVector::B.first << endl;
         //cout << " in pilot function " << H.first << " mat: " << mat.xx << endl;
-        return pilotHBMuMuDiff;
+        return pilot;
     }   
     template<unsigned int DIM>
-    pVpilots_t<DIM> preisachVector<DIM>::pilotForward(const std::initializer_list<double> H){
+    shared_ptr<preisachVector<DIM>> preisachVector<DIM>::pilotForward(const std::initializer_list<double> H){
         vec_t<double, DIM> o(H);
         return pilotForward(o);
     }
@@ -1179,23 +924,15 @@ namespace preisachVector_ns{
     // }
 
     template<unsigned int DIM>
-    pVpilots_t<DIM> preisachVector<DIM>::pilotInverse(const std::initializer_list<double> BorH){
-
-        vec_t<double, DIM> HorB = addInverse(BorH, true);
-
-        calcMat(pilotHBMuMuDiff.third, true);
-        calcMatDiff(pilotHBMuMuDiff.fourth, true);
+    shared_ptr<preisachVector<DIM>> preisachVector<DIM>::pilotInverse(const std::initializer_list<double> BorH){
 
+        updatePilot();
+        vec_t<double, DIM> HorB = pilot->addInverse(BorH);
 
-        pilotHBMuMuDiff.first = CHECK_E_FW_INV(HorB, BorH);
-        pilotHBMuMuDiff.second = CHECK_E_FW_INV(HorB, BorH);
-        
-        //cout << "current H " << preisachVector::H.first << "current mat " << preisachVector::mat.xx << "current B" << preisachVector::B.first << endl;
-        //cout << " in pilot function " << H.first << " mat: " << mat.xx << endl;
-        return pilotHBMuMuDiff;
+        return pilot;
     }
     template<unsigned int DIM>
-    pVpilots_t<DIM> preisachVector<DIM>::pilotInverse(const vec_t<double, DIM>& B){
+    shared_ptr<preisachVector<DIM>> preisachVector<DIM>::pilotInverse(const vec_t<double, DIM>& B){
         vec_t<double, DIM> o(B);
         return pilotInverse(o);
     }
diff --git a/C++/preisach/preisachVector/preisachVector.h b/C++/preisach/preisachVector/preisachVector.h
index 3be79a4ff31f6daf07caa50da73a215b6d46c4e7..6351f608a0262c1b4ce0029347cb65b681c585cc 100644
--- a/C++/preisach/preisachVector/preisachVector.h
+++ b/C++/preisach/preisachVector/preisachVector.h
@@ -88,7 +88,7 @@ namespace preisachVector_ns{
         public:
 
             preisachVector(const shared_ptr<EverettMatrix>&, const shared_ptr<distributions>&, const bool interpolate=true, 
-                        const bool differential=true, const std::vector<double> rotationalParameters=std::vector<double>());
+                        const bool differential=true, const std::vector<double> rotationalParameters=std::vector<double>(), bool isPilot=false);
             
             preisachVector(const preisachVector&);
             preisachVector(const shared_ptr<preisachVector>&);
@@ -103,50 +103,44 @@ namespace preisachVector_ns{
 
             void demagnetise();
 
-            const vec_t<double, DIM> calcInputValue(const bool bPilot=false);
-            void calcInputValue(vec_t<double, DIM>&, const bool bPilot=false) const;
+            const vec_t<double, DIM> calcInputValue();
+            const vec_t<double, DIM> calcOutputValue();
 
-            const vec_t<double, DIM> calcOutputValue(const bool bPilot=false);
-            void calcOutputValue(vec_t<double, DIM>&, const bool bPilot=false) const;
+            void updatePilot();
 
-            const mat_t<double, DIM, DIM> calcMat(const bool bPilot=false);
-            void calcMat(mat_t<double, DIM, DIM>&, const bool bPilot=false) const;
-
-            const mat_t<double, DIM, DIM> calcMatDiff(const bool bPilot=false);
-            void calcMatDiff(mat_t<double, DIM, DIM>&, const bool bPilot=false) const;
+            const mat_t<double, DIM, DIM> calcMat();
             
-            const mat_t<double, DIM, DIM> calcMatDiffDiff(const bool bPilot=false);
-            void calcMatDiffDiff(mat_t<double, DIM, DIM>&, const bool bPilot=false) const;
 
+            const mat_t<double, DIM, DIM> calcMatDiff();
+            
+            const mat_t<double, DIM, DIM> calcMatDiffDiff();
+            
             /* forward model     */
             vec_t<double, DIM> addForward(const std::initializer_list<double> H);
             vec_t<double, DIM> addForward(const vec_t<double, DIM> H);
             vec_t<double, DIM> addForward(double H_x, double H_y=0, double H_z=0);
-            vec_t<double, DIM> addForward(const std::initializer_list<double> H, bool bPilot);
-            vec_t<double, DIM> addForward(const vec_t<double, DIM> H, bool bPilot);
-            void addForward_thread(const vec_t<double, DIM> newVal_in, bool bPilot, const unsigned int istart, const unsigned int iend);
+            void addForward_thread(const vec_t<double, DIM> newVal_in, const unsigned int istart, const unsigned int iend);
 
-            pVpilots_t<DIM> pilotForward(const std::initializer_list<double> H);
-            pVpilots_t<DIM> pilotForward(const vec_t<double, DIM>& H);
+            shared_ptr<preisachVector> pilotForward(const std::initializer_list<double> H);
+            shared_ptr<preisachVector> pilotForward(const vec_t<double, DIM>& H);
 
 
             /* inverse model     */
-            vec_t<double, DIM> addInverse(const vec_t<double, DIM>& B, bool bPilot=false);
-            vec_t<double, DIM> addInverse(const std::initializer_list<double>, bool bPilot=false);
+            vec_t<double, DIM> addInverse(const vec_t<double, DIM>& B);
+            vec_t<double, DIM> addInverse(const std::initializer_list<double>);
             vec_t<double, DIM> addInverse(double B_x, double B_y=0, double B_z=0);
-            void addInverse_thread(const vec_t<double, DIM> newVal_in, bool bPilot, const unsigned int istart, const unsigned int iend);
+            void addInverse_thread(const vec_t<double, DIM> newVal_in, const unsigned int istart, const unsigned int iend);
 
             // vec_t<double, DIM> addB_BS(const vec_t<double, DIM>& B, const bool bPilot=false, const double eps = 1e-2);
             // vec_t<double, DIM> addB_LIN(const vec_t<double, DIM>& B, const bool bPilot=false, const double eps = 1e-2);
 
-            pVpilots_t<DIM> pilotInverse(const std::initializer_list<double> B);
-            pVpilots_t<DIM> pilotInverse(const vec_t<double, DIM>& B);
+            shared_ptr<preisachVector> pilotInverse(const std::initializer_list<double> B);
+            shared_ptr<preisachVector> pilotInverse(const vec_t<double, DIM>& B);
 
 
 
 
-            pVpilots_t<DIM> getPilots()const {return pilotHBMuMuDiff;};
-            pVpilots_t<DIM>& getPilots(){return pilotHBMuMuDiff;}; 
+            shared_ptr<preisachVector> getPilot()const {return pilot;};
 
 
             pVpilots_t<DIM> getPrevDiffValInOutMatMatDiff()const {return prevDiffValInOutMatMatDiff;};
@@ -246,10 +240,10 @@ namespace preisachVector_ns{
             vec_t<double, DIM> addB(const std::initializer_list<double> B){return CHECK_E_FW_INV(addInverse(B), addForward(B));}
             vec_t<double, DIM> addB(const vec_t<double, DIM> B){return CHECK_E_FW_INV(addInverse(B), addForward(B));}
             vec_t<double, DIM> addB(double B_x, double B_y=0, double B_z=0){return CHECK_E_FW_INV(addInverse(B_x, B_y, B_z), addForward(B_x, B_y, B_z));}
-            pVpilots_t<DIM> pilotH(const std::initializer_list<double> H){return CHECK_E_FW_INV(pilotForward(H), pilotInverse(H));}
-            pVpilots_t<DIM> pilotH(const vec_t<double, DIM>& H){return CHECK_E_FW_INV(pilotForward(H), pilotInverse(H));}
-            pVpilots_t<DIM> pilotB(const std::initializer_list<double> B){return CHECK_E_FW_INV(pilotInverse(B), pilotForward(B));}
-            pVpilots_t<DIM> pilotB(const vec_t<double, DIM>& B){return CHECK_E_FW_INV(pilotInverse(B), pilotForward(B));}
+            shared_ptr<preisachVector> pilotH(const std::initializer_list<double> H){return CHECK_E_FW_INV(pilotForward(H), pilotInverse(H));}
+            shared_ptr<preisachVector> pilotH(const vec_t<double, DIM>& H){return CHECK_E_FW_INV(pilotForward(H), pilotInverse(H));}
+            shared_ptr<preisachVector> pilotB(const std::initializer_list<double> B){return CHECK_E_FW_INV(pilotInverse(B), pilotForward(B));}
+            shared_ptr<preisachVector> pilotB(const vec_t<double, DIM>& B){return CHECK_E_FW_INV(pilotInverse(B), pilotForward(B));}
             [[nodiscard]] double getNormH()const {return CHECK_E_FW_INV(getNormInputValue(), getNormOutputValue());}
             [[nodiscard]] double getNormB()const {return CHECK_E_FW_INV(getNormOutputValue(), getNormInputValue());}
             void setNormB( double b) {return CHECK_E_FW_INV(setNormOutputValue(b), setNormInputValue(b));}
@@ -265,6 +259,8 @@ namespace preisachVector_ns{
             shared_ptr<EverettMatrix> E = nullptr;
             shared_ptr<distributions> dist = nullptr;
             std::vector<double> rotationalParameters;
+            bool isPilot = false;
+
             double normInputValue = 1;
             double normOutputValue = 1;
             double normEnergy = std::nan("not set yet");
@@ -279,9 +275,9 @@ namespace preisachVector_ns{
             mat_t<double, DIM, DIM> mat;
             mat_t<double, DIM, DIM> mat_diff;
 
+            
 
-
-            pVpilots_t<DIM> pilotHBMuMuDiff;  
+            shared_ptr<preisachVector> pilot;  
             pVpilots_t<DIM> prevDiffValInOutMatMatDiff;  
 
             void calcNormHB();
diff --git a/python/myPackage.py b/python/myPackage.py
index d4cc9fb1644abb95dd5149c4bf608aad4f8c51d9..297031c514bcd338f724f7a40e9c147cd8e61b1c 100644
--- a/python/myPackage.py
+++ b/python/myPackage.py
@@ -128,7 +128,7 @@ class L2Draw:
 
 def myNewtonSolver(gfu,a,f,eps=1e-13,Nmax=10, callback=lambda gfu, it: None,
     printrates=True, printError = True, converge_func=lambda *args : False, 
-    dampfactor = 1, inverse="pardiso", pre = None, itToClacJacobi=1, a4Linearisation=None):
+    dampfactor = 1, inverse="pardiso", pre = None, itToClacJacobi=1):
     """Newton solver for nonlinear problems
 
     :param gfu: solution vector
@@ -145,7 +145,6 @@ def myNewtonSolver(gfu,a,f,eps=1e-13,Nmax=10, callback=lambda gfu, it: None,
     :param pre: optional preconditioner to solve the system in each iteration, defaults to None
     :return: convergence, it
     """
-    a4Linearisation = check_input(a4Linearisation, a)
 
 
 
@@ -170,14 +169,14 @@ def myNewtonSolver(gfu,a,f,eps=1e-13,Nmax=10, callback=lambda gfu, it: None,
 
 
         if it%itToClacJacobi == 0:
-            a4Linearisation.AssembleLinearization(gfu.vec)
+            a.AssembleLinearization(gfu.vec)
         
         f.Assemble()
 
         if True or type(pre) == type(None):
-            du.vec.data = a4Linearisation.mat.Inverse(fes.FreeDofs(), inverse=inverse) * (res.vec - f.vec)
+            du.vec.data = a.mat.Inverse(fes.FreeDofs(), inverse=inverse) * (res.vec - f.vec)
         else:
-            solvers.BVP(bf = a4Linearisation, lf = f, gf = du, pre = pre)
+            solvers.BVP(bf = a, lf = f, gf = du, pre = pre)
         du.vec.data *= dampfactor
         gfu.vec.data -= du.vec
 
@@ -1260,7 +1259,7 @@ def plotPV(pV, nFig = 0, ax = None, field="HB", numbers=False, label=True, file=
 
 
 
-def plotEVF(ev, nFig = 0, title=None, show=True, vec= None, only3D=False):
+def plotEVF(ev, nFig = 0, title=None, show=True, vec= None, only3D=False, only2D=False):
     import matplotlib.pyplot as plt
     from mpl_toolkits.mplot3d import Axes3D
     import numpy as np  
@@ -1280,35 +1279,37 @@ def plotEVF(ev, nFig = 0, title=None, show=True, vec= None, only3D=False):
             if np.isnan(Z[i, j]):
                 print(X[i, j], Y[i, j], "\t->\t", Z[i, j])
 
-    fig = plt.figure(nFig)#, figsize=plt.figaspect(2))
-    plt.clf()    
-    ax = plt.axes(projection='3d')
-    #ax = fig.add_subplot(2, 1, 2, projection='3d')
+    if not only2D:
+        fig = plt.figure(nFig)#, figsize=plt.figaspect(2))
+        plt.clf()    
+        ax = plt.axes(projection='3d')
+        #ax = fig.add_subplot(2, 1, 2, projection='3d')
 
-    # plt.ion()
-    ax.plot_surface( X, Y, Z)#, c=1, vmin=-1, vmax=1,s=50)
-    # ax.scatter( X, Y, Z)#, c=1, vmin=-1, vmax=1,s=50)
+        # plt.ion()
+        ax.plot_surface( X, Y, Z)#, c=1, vmin=-1, vmax=1,s=50)
+        # ax.scatter( X, Y, Z)#, c=1, vmin=-1, vmax=1,s=50)
 
 
-    #pB = ax.scatter( spherePoints[0,:], y, z, s= 10, c='r')
-    # ax.set_zlim3d(-1, 1)
-    # plt.xlim(-1,1)
-    # plt.ylim(-1,1)
+        #pB = ax.scatter( spherePoints[0,:], y, z, s= 10, c='r')
+        # ax.set_zlim3d(-1, 1)
+        # plt.xlim(-1,1)
+        # plt.ylim(-1,1)
 
-    #plt.xticks(vec)
-    #plt.yticks(vec)
-    Z = np.flipud(Z)
-    
-    if type(title) == type(None):
-        plt.title(ev.name + " , dim=" + str(ev.dim))
-    else:
-        plt.title(title)
-    ax.set_xlabel(r"$\beta$")
-    ax.set_ylabel(r"$\alpha$")
-    ax.set_zlabel(r"$E(\alpha, \beta)$")
+        #plt.xticks(vec)
+        #plt.yticks(vec)
+        
+        if type(title) == type(None):
+            plt.title(ev.name + " , dim=" + str(ev.dim))
+        else:
+            plt.title(title)
+        ax.set_xlabel(r"$\beta$")
+        ax.set_ylabel(r"$\alpha$")
+        ax.set_zlabel(r"$E(\alpha, \beta)$")
 
 
     if not only3D:
+        Z = np.flipud(Z)
+
         fig =plt.figure()
         ax = fig.add_subplot(111)
         ax.matshow(Z)