#include "vs.h"
#define PI 3.141592654
int main() {
   double a_= 1.0; // radius of the circle
	// A. Bode's Integration Formula
   double weight[25] = {14.0/45.0, 64.0/45.0, 24.0/45.0, 64.0/45.0, 28.0/45.0,
                                   64.0/45.0, 24.0/45.0, 64.0/45.0, 28.0/45.0,
                                   64.0/45.0, 24.0/45.0, 64.0/45.0, 28.0/45.0,
                                   64.0/45.0, 24.0/45.0, 64.0/45.0, 28.0/45.0,
                                   64.0/45.0, 24.0/45.0, 64.0/45.0, 28.0/45.0,
                                   64.0/45.0, 24.0/45.0, 64.0/45.0, 14.0/45.0};
   Quadrature qp(weight, 0.0, a_, 25);
   J d_r(a_/24.0);      // per normalized length

   // B. Define Basis Functions
   H2 r((double*)0, qp),   // radius variable
      phi = INTEGRABLE_VECTOR_OF_TANGENT_OF_TANGENT_BUNDLE("int, int, Quadrature",
            3/*vector size*/, 1/*spatial dim.*/, qp);
   phi[0] = cos((PI/2.0/a_)*r); phi[1] = cos((3.0*PI/2.0/a_)*r); phi[2] = cos(5.0*PI/2.0/a_*r);
   H0 d2_phi = INTEGRABLE_VECTOR("int, Quadrature", 3, qp);
   for(int i = 0; i < 3; i++) d2_phi[i] = dd(phi)(i)[0][0]; // degenerated hessian matrix

   // C. Weighted Residual Formulation and the solution of the Generalized Eigenvalue Problem
   C0 A = -( ( d(phi)(0)%((H0)phi) + d2_phi%((H0)phi)*((H0)r) ) | d_r ),
      B = ( ((H0)phi)%((H0)phi) * ((H0)r) ) | d_r;
   C0 L = MATRIX("int, int", 3, 3);
   Cholesky Ch(B);
   for(i = 0; i < 3; i++)
   for(int j = 0; j < 3; j++)
   	if(i >= j) L[i][j] = Ch.rep_ptr()[0][i*(1+i)/2+j];
      else L[i][j] = 0.0;
   C0 L_inv = L.inverse();
   C0 C = L_inv*((A+~A)/2.0*(~L_inv));
   C0 lambda = Eigen(C).Eigenvalues();
   C0 v = Eigen(C).Eigenvectors();
   cout << lambda << endl;
   cout << v << endl;
	return 0;
}