#include "vs.h"
#define PI 3.141592654
int main() {
	// A. Bode's Integration Formula--exact up to 5-th order
	double w[17] = {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, 14.0/45.0};
   Quadrature qp(w, 0.0, 1.0, 17);
   ;      // the normalized length of the integration segments

   // B. Define Basis Functions
   H1 x(qp),
      phi = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature",
                                      4/*vector size*/, 1/*spatial dim.*/, qp);
   phi[0] = x*(1.0-x);
   for(int i = 1; i < 4; i++) phi[i] = phi[i-1]*x;
   // C. Variational Formulation
   C0 M = ( d(phi)*(~d(phi)) ) | J(1.0/16.0);            // LSH bilinear form
   C0 b = ( ((H0)phi) * cos(PI*((H0)x)) ) | J(1.0/16.0); // RHS forced vector

   // D. Solutions--the Ritz Coefficients
   C0 c = b / M;
   // coefficients of semi-analytical Rayleigh-Ritz approximation, p. 270
   //double semi[4] = {0.20220, -0.28910, -0.34610, 0.23070};
   for(i = 0; i < 4; i++) cout << c[i] << endl;
   //	assert(fabs( ( ((double)c[i]) - semi[i] )/semi[i] ) < 1.e-2);
   return 0;
}