#include "vs.h"
#if defined(__TRANSCENDENTAL)
H1 cosh(const H1& a) { return (exp(a)+exp(-a))/2.0; }
#endif
int main() {
	const double f_ = 1.0;
   double weight[13] = {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, 14.0/45.0};
   Quadrature qp(weight, 0.0, 1.0, 13);
   J d_l(1.0/12.0);      // per normalized length
	H1 X = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE( "int, int, Quadrature", 2, 2, qp),
		phi = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE( "int, int, Quadrature", 2, 2, qp),
      x, y; x &= X[0]; y &= X[1];
   ((H0)x) = 1.0; d(y) = 0.0;
   H0 v_bar =  f_/4.0 *(((H0)x).pow(2)+((H0)y).pow(2));
#if defined(__TRANSCENDENTAL)
	const double PI = 3.141592654;
	phi[0] = cos(PI/2.0*x)*cosh(PI/2.0*y) + cos(PI/2.0*y)*cosh(PI/2.0*x);
   phi[1] = cos(3.0*PI/2.0*x)*cosh(3.0*PI/2.0*y) + cos(3.0*PI/2.0*y)*cosh(3.0*PI/2.0*x);
#else
   phi[0] = x.pow(4) - 6.0*x.pow(2)*y.pow(2) + y.pow(4);
   phi[1] = x.pow(8) -28.0*x.pow(6)*y.pow(2)+70.0*x.pow(4)*y.pow(4)
                     -28.0*x.pow(2)*y.pow(6)+y.pow(8);
#endif
   H0 dphi_dx = d(phi)(0);
   C0 M = dphi_dx % ((H0)phi) | d_l,
      b = dphi_dx * v_bar | d_l;
   C0 c = b / M;
	cout << c << endl;
   C0 k = -((- v_bar + ((H0)phi)*c) | d_l);
   cout << k << endl;
	return 0;
}