#include "vs.h"
#define NODE_NO 3

#if defined(__ARC_LENGTH)
H1 sqrt(const H1& a) {
	H1 b = (H1)a;
	return INTEGRABLE_TANGENT_BUNDLE("const H0&, const H0&", 
										sqrt((H0)b), d(b)/2.0/sqrt((H0)b));
}
#endif

int main() {
	const double r = 1.0; const double PI = 3.141592654;
	double X[NODE_NO][2];
   for(int i = 0; i < NODE_NO; i++) {
   	X[i][0] = r*cos(((double)(i+1))*PI/(2.0*(NODE_NO-1)));
#ifndef __ARC_LENGTH
      X[i][1] = r*sin(((double)(i+1))*PI/(2.0*(NODE_NO-1)));
#endif
   }
	Quadrature qp(1, 4);
	H1 zai(qp),
      N = INTEGRABLE_VECTOR_OF_TANGENT_BUNDLE("int, int, Quadrature", 3, 1, qp),
      x, y;
   N[0] = -zai*(1-zai)/2; N[1] = (1+zai)*(1-zai); N[2] = zai*(1+zai)/2;
   C0 length = 0.0;
   for(i = 0; i < (NODE_NO-1)/2; i++) {
   	x = N[0]*X[i+2][0]+N[1]*X[i+1][0]+N[2]*X[i][0];
#if defined(__ARC_LENGTH)
      y = sqrt(pow(r,2)-x.pow(2));
      length += (sqrt(1.0+(d(y)/d(x)).pow(2)) | J(d(x)));
#else
      y = N[0]*X[i+2][1]+N[1]*X[i+1][1]+N[2]*X[i][1];
   	length += (sqrt(d(x).pow(2)+d(y).pow(2))) | J(1.0);
#endif
   }
   cout << length << endl;     
	return 0;
}