// to execute this code, do
//   root > .L MFunction.C+
//   root > .x TestMFunction.C+

#include "TF2.h"
#include "TMatrixD.h"
#include "MFunction.h"
#include "Riostream.h"
   
void TestMFunction() {

     TF2* g00 = new TF2("g00", "x^2 - y", -1, 1, -1, 1);
     TF2* g01 = new TF2("g01", "2*x*y - 1", -1, 1, -1, 1);
     TF2* g10 = new TF2("g10", "sin(3*x) - y^2", -1, 1, -1, 1);
     TF2* g11 = new TF2("g11", "3*x - 7*y -2", -1, 1, -1, 1);

     MFunction m1;
     m1.SetFunctions(g00, g01, g10, g11);

     cout<<"Simple MFunction example:"<<endl;
     m1.PrintExpFormula();

     cout<<"Now it has been (incorrectly) inverted:"<<endl;
     m1.Invert();
     m1.PrintExpFormula();

     cout<<endl<<endl;

     TF2* f00 = new TF2("f00", "[0] + [1]*x + [2]*y + [3]*x^2 + [4]*x*y + [5]*y^2 + [6]*x^3 + [7]*x^2*y + [8]*x*y^2 + [9]*y^3", -1, 1, -1, 1);
     TF2* f01 = new TF2("f01", "[0] + [1]*x + [2]*y + [3]*x^2 + [4]*x*y + [5]*y^2 + [6]*x^3 + [7]*x^2*y + [8]*x*y^2 + [9]*y^3", -1, 1, -1, 1);
     TF2* f10 = new TF2("f10", "[0] + [1]*x + [2]*y + [3]*x^2 + [4]*x*y + [5]*y^2 + [6]*x^3 + [7]*x^2*y + [8]*x*y^2 + [9]*y^3", -1, 1, -1, 1);
     TF2* f11 = new TF2("f11", "[0] + [1]*x + [2]*y + [3]*x^2 + [4]*x*y + [5]*y^2 + [6]*x^3 + [7]*x^2*y + [8]*x*y^2 + [9]*y^3", -1, 1, -1, 1);

     MFunction m2;
     m2.SetFunctions(f00, f01, f10, f11);

     cout<<"MFunction of 3rd order bi-polynomials:"<<endl;

     m2.PrintExpFormula();

     m2.Invert();

     cout<<"Here it has been (incorrectly) inverted:"<<endl;

     m2.PrintExpFormula();
     
     TMatrixD matrix(2,2);
     m2.Eval(matrix,0.5,0.5);
     matrix.Print();
}