#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>
#include <iomanip.h>
#include "vs.h"
#define PI 3.141592654
ofstream ofs("temp.cpp", ios::out | ios::trunc);
H0 cosh(const H0& a) { return (exp(a)+exp(-a))/2.0; }
int main() {
   double f_0 = -1.0,
	// A. 2-D 17x17-point Bode's Integration Formula
          bode[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},
	       weight[17][17];
   for(int i = 0; i < 17; i++)
   for(int j = 0; j < 17; j++) weight[i][j] = bode[i] * bode[j];
   Quadrature qp(weight[0], 0.0, 1.0, 17, 0.0, 1.0, 17);
   J d_a(pow( (1.0/16.0), 2.0));      // per normalized area
   double x0, x1;
   H1 z(2, (double*)0, qp);
   H0 integrand(qp), zai, eta; zai &= ((H0)z)[0]; eta &= ((H0)z)[1];
   for(i = 0; i < 11; i++) {
   	for(int j = 0; j < 11; j++) {
   		x0 = 0.1 * (double)j; x1 = 0.1 * (double)i;
      	integrand = 0.0;
      	for(int n = 1; n <= 20; n++) {
        		integrand += (f_0/(n*PI*sinh(n*PI))* sin(n*PI*x0)*sin(n*PI*zai)*
            	(cosh(n * PI *(1.0-x1-eta))-cosh(n * PI *(1.0-sqrt((x1-eta)*(x1-eta))))));
         }
      	C0 integal = integrand | d_a;
      	ofs << "u(" << x0 << ", " << x1 << ") = " << ((double)integal) << ", ";
   	}
      ofs << endl;
   }
   ofs.close();
	return 0;
}
/*
int main() {
   double x0, x1, f_0 = 1.0,
	// A. 2-D 6x6 semi open (left) formula
          semi_open[11] = {0.0, 23.0/12.0, 7.0/12.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 13.0/12.0, 5.0/12.0},
	       w_q1[11][11], w_q2[11][11], w_q3[11][11], w_q4[11][11];
   for(int i = 0; i < 11; i++)
   for(int j = 0; j < 11; j++) {
   	w_q1[i][j] = semi_open[i] * semi_open[j];
      w_q2[i][j] = semi_open[10-i] * semi_open[j];
      w_q3[i][j] = semi_open[10-i] * semi_open[10-j];
      w_q4[i][j] = semi_open[i] * semi_open[10-j];
   }

   for(int i = 0; i < 11; i++)
   for(int j = 0; j < 11; j++) {
   	x0 = 0.1 * (double)j; x1 = 0.1 * (double)i;
   	Quadrature q1(w_q1[0], x0, 1.0, 11, x1, 1.0, 11),
                 q2(w_q2[0], 0.0, x0, 11, x1, 1.0, 11),
                 q3(w_q3[0], 0.0, x0, 11, 0.0, x1, 11),
                 q4(w_q4[0], x0, 1.0, 11, 0.0, x1, 11);
      H0 z[4],
      	z0(2, (double*)0, q1), z1(2, (double*)0, q2),
         z2(2, (double*)0, q3), z3(2, (double*)0, q4);
      z[0] &= z0; z[1] &= z1; z[2] &= z2; z[3] &= z3;
      double a[4];
      a[0] = (1.0-x0)/10.0*(1.0-x1)/10.0,
      a[1] = x0/10.0*(1.0-x1)/10.0,
      a[2] = x0/10.0*x1/10.0,
      a[3] = (1.0-x0)/10.0*x1/10.0;
      C0 integal = 0.0;
      for(int k = 0; k < 4; k++) {
   		J d_a(a[k]);
      	if(a[k] != 0.0) {
      		H0 r = sqrt((x0-z[k][0]).pow(2)+(x1-z[k][1]).pow(2));
         	integal += 1.0/(2.0*PI) * (log(1.0/r)) | d_a;
         }
      }
      cout << "u(" << x0 << ", " << x1 << ") = " << ((double)integal) << endl;
   }
   return 0;
}*/