#include "vs.h"
int main() {
	C0 c_t_n(2, (double*)0), c_t_n1(2, (double*)0);
	{
	// A. Bode's Integration Formula
   double weight[5] = {14.0/45.0, 64.0/45.0, 24.0/45.0, 64.0/45.0, 14.0/45.0};
   Quadrature qp(weight, 0.0, 1.0, 5);
   J d_l(1.0/4.0);      // per normalized length

   // B. Define Basis Functions
   H1 x(qp),
      phi = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature",
            2/*vector size*/, 1/*spatial dim.*/, qp);
   phi[0] = x; phi[1] = x.pow(2);

   // initial condition
   C0 C = (((H0)phi)%((H0)phi))|d_l,
      b = ((H0)phi)|d_l;
   c_t_n = b / C;
   }
   cout << c_t_n << endl;
	// A. Simpson's Integration Formula
   double weight[3] = {1.0/3.0, 4.0/3.0, 1.0/3.0};
   Quadrature qp(weight, 0.0, 1.0, 3);
   J d_l(1.0/2.0);      // per normalized length

   // B. Define Basis Functions
   H1 x(qp),
      phi = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature",
            2/*vector size*/, 1/*spatial dim.*/, qp);
   phi[0] = x; phi[1] = x.pow(2);

   // B. Variational Formulation 
   C0 C = (((H0)phi)%((H0)phi))|d_l,
      K = (d(phi)(0)%d(phi)(0))|d_l;

   // C. Time Discretization
   double theta_ = 0.5, // central difference
          dt_ = 0.05;   // time step,
	C0 LHS = C + theta_ * dt_ * K,
      RHS = C - (1.0 - theta_) * dt_ * K;

   // C. Solution; time-integration
   //c_n[0] = 4.0; c_n[1] = -10.0/3.0; // initial condition; i.e., u(x, 0) = 1.0
   C0 d_LHS = !LHS; // decomposed once
   for(int i = 0; i < 30; i++) {
   	c_t_n1 = d_LHS*(RHS*c_t_n); // forward/backward substitution
		double iptr;
      if(modf( ((double)(i+1))/2.0, &iptr)==0)cout << c_t_n1 << endl;
      //double iptr;
      //if(modf( ((double)(i+1))/4.0, &iptr)==0)
      	//cout << "t= " << ((i+1)*dt_) << ", u(0.5, 1.0) = " <<
      	//(c_t_n1[0]*0.5+c_t_n1[1]*0.25) << ", " << (c_t_n1[0]+c_t_n1[1]) << endl;
      c_t_n = c_t_n1; // update the solution
   }
	return 0;
}