// To compile this code:
// c++ plotMatrixCalcs.c -o plotMatrix
// or with root libraries:
// g++ plotMatrixEigen.c $(root-config --cflags --glibs) -o plotMatrix
// To run the executable:
// ./plotMatrix

int dimension = 100; // Matrix dimension
double R = 18.; // m
double a = 3.; // m
double PI = 3.1415926;
double mu0 = 4*PI*1.e-7; // N/A/A
int nLoops = 10; // Number of current loops

#include <TMinuit.h>
#include <vector>
#include <cstdio>
#include <cmath>
#include </usr/include/eigen3/Eigen/Dense>  // Full path here helps to compile.
//#include <Eigen/Dense>
#include <iostream>




Eigen::MatrixXd GetM()
{
    Eigen::MatrixXd M(dimension, dimension);
    for (int i=0; i<dimension; i++)
    {
        for (int j=0; j<dimension; j++)
        {
            double tElement = 1. / pow( (i-j)*(i-j) + R/a*R/a, 1.5);             
            M(i,j) = tElement;
        }
    }

    return M;
}

Eigen::VectorXd GetU()
{
    Eigen::VectorXd U(dimension);
    for (int i=0; i<dimension; i++)
    {
        double tElement = 1.0;             
        U(i) = tElement;
    }   
    return U;
}

Eigen::VectorXd GetB()
{

    Eigen::VectorXd B(dimension);
    for (int i=0; i<dimension; i++)
    {
        double tElement = 1.0;             
        B(i) = tElement;
    }

    return B;
}


/*
void PlotHisto(Eigen::MatrixXd aMatrix)
{
    int aSize = dimension;
    double xValues[aSize];
    double yValues[aSize];
    for (int i=0; i<aSize; i++)
    {
        xValues[i] = i;
        yValues[i] = aMatrix(i,0);
    }
    TGraph* gMatrix = new TGraph(aSize, xValues, yValues);
    gMatrix->SetTitle("U = M^{-1} B, B = 1.0, 0 < i < 100, 0 < j < 100");
    gMatrix->GetXaxis()->SetTitle("index i");
    gMatrix->GetYaxis()->SetTitle("U ( i ), unnormalized");    
    gMatrix->SetLineWidth(2);
    gMatrix->GetXaxis()->SetLabelSize(0.045);
    gMatrix->GetYaxis()->SetLabelSize(0.045);
    gMatrix->GetXaxis()->SetTitleSize(0.05);        
    gMatrix->GetYaxis()->SetTitleSize(0.05);
    gMatrix->Draw("alp");
}
*/
 
int main()
{
    printf("hello world\n");

//    Eigen::VectorXd U = GetU();
    Eigen::MatrixXd M = GetM();
//    Eigen::MatrixXd B = M*U;
//    PlotHisto(B);
    Eigen::VectorXd B = GetB();
    Eigen::MatrixXd U = ((M).inverse())*B;
//    PlotHisto(U);


    return 0;
}