From 72796aa0c7e30bbdbcbf1f8d37268bc979780ace Mon Sep 17 00:00:00 2001
From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Wed, 11 May 2011 14:22:07 +0000
Subject: [PATCH] =?utf8?q?[src]=20tests=20f=C3=BCr=20g0=20G00=20und=20F=5F?=
 =?utf8?q?par=20hinzugef=C3=BCgt=20[src]=20kleine=20Fehler=20behoben...?=
MIME-Version: 1.0
Content-Type: text/plain; charset=utf8
Content-Transfer-Encoding: 8bit

git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@21 26120e32-c555-405d-b3e1-1f783fb42516
---
 src/SLPrecangle.cpp   |  40 +++++---
 src/SLPrecangle.hpp   |   6 ++
 src/build_A.cpp       | 225 +-----------------------------------------
 src/build_A.mexa64    | Bin 17226 -> 17398 bytes
 src/build_A2.m        |  12 +--
 src/mex_Fpar.cpp      |  44 +++++++++
 src/mex_G00.cpp       |  38 +++++++
 src/mex_g0.cpp        |   2 +-
 src/refineQuad.m      |   4 +
 src/surfDoubleQuad.m  |  45 +++++++++
 src/surfQuad.m        |  19 ++++
 src/test_Fpar.m       |  19 ++++
 src/test_G00.m        |  23 +++++
 src/test_functions.m  |  60 +++++++++++
 src/test_g0.m         |   4 +
 src/test_solveError.m |  17 +++-
 src/test_solveS.m     |  19 +---
 17 files changed, 310 insertions(+), 267 deletions(-)
 create mode 100644 src/mex_Fpar.cpp
 create mode 100644 src/mex_G00.cpp
 create mode 100644 src/surfDoubleQuad.m
 create mode 100644 src/surfQuad.m
 create mode 100644 src/test_Fpar.m
 create mode 100644 src/test_G00.m
 create mode 100644 src/test_functions.m

diff --git a/src/SLPrecangle.cpp b/src/SLPrecangle.cpp
index 5a414b9..d5c9e11 100644
--- a/src/SLPrecangle.cpp
+++ b/src/SLPrecangle.cpp
@@ -12,13 +12,8 @@
 using namespace std;
 
 int sign(double);
-double arsinh(double);
+//double arsinh(double);
 
-// sol = quadInt((F_par/F_ort),s1,s2,k1,k2,t1,t2,l1,l2,d1,d2,d3);
-double quadInt(
-		double(*)(double, double, double, double, double, double, double),
-		double, double, double, double, double, double, double, double, double,
-		double, double);
 /*
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 	//sicherheitsabfragen zu Datengroessen
@@ -166,11 +161,11 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 int inline sign(double x) {
 	return x > 0 ? 1 : (x < 0 ? -1 : 0);
 }
-
+/*
 double inline arsinh(double x) {
 	return log(x + sqrt(x * x + 1));
 }
-
+*/
 //y-x muss != 0 sein
 double g0(double p, double y, double x, double l) {
 	//printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);
@@ -187,7 +182,7 @@ double g0(double p, double y, double x, double l) {
 		else if (p == -0.5)
 			sol = asinh((y - x) / fabs(l));
 		else if (p == -1)
-			sol = atan((y - x) / fabs(l));
+			sol = atan((y - x) / fabs(l))/2;
 		else if (p == -1.5)
 			sol = (y - x) / ((l * l) * sqrt((y - x) * (y - x) + l * l));
 		else
@@ -229,7 +224,9 @@ double G00(double p, double y1, double y2, double x1, double x2, double l) {
 					sqrt((y2 - x2) * (y2 - x2) + l * l));
 		sol /= 2 * p + 2;
 	} else {
-		mexErrMsgTxt("no case for p defined");
+		sol = NAN;
+		printf("warning in G00: no case for p defined\n");
+		//mexErrMsgTxt("no case for p defined\n");
 	}
 
 	return sol;
@@ -242,14 +239,14 @@ double F_par(double x1, double x2, double y1, double y2, double d1, double d2,
 	double sol = (x1 - y1 - d1) * (x2 - y2 - d2);
 
 	if (sol != 0)
-		sol *= slpADLO(x1, x2, y1 + d1, y2 + d2, d3);
+		sol *= G00(-0.5,x1, x2, y1 + d1, y2 + d2, d3);
 
 	if ((x1 - y1 - d1) != 0)
-		sol -= (x1 - y1 - d1) * compute_g0(0.5, x1, y1 + d1,
+		sol -= (x1 - y1 - d1) * g0(0.5, x1, y1 + d1,
 				sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));
 
 	if ((x2 - y2 - d2) != 0)
-		sol -= (x2 - y2 - d2) * compute_g0(0.5, x2, y2 + d2,
+		sol -= (x2 - y2 - d2) * g0(0.5, x2, y2 + d2,
 				sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));
 
 	double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2
@@ -258,11 +255,16 @@ double F_par(double x1, double x2, double y1, double y2, double d1, double d2,
 	return sol;
 }
 
-double inline quadInt(
+double quadInt(
 		double(*f)(double, double, double, double, double, double, double),
-		double s1, double s2, double k1, double k2, double t1, double t2,
-		double l1, double l2, double d1, double d2, double d3) {
+		double a, double b, double c, double d, double r, double t,
+		double u, double v, double d1, double d2, double d3) {
 
+   return f(b,d,t,v,d1,d2,d3)-f(b,d,t,u,d1,d2,d3)-f(b,d,r,v,d1,d2,d3)+f(b,d,r,u,d1,d2,d3)
+    -f(b,c,t,v,d1,d2,d3)+f(b,c,t,u,d1,d2,d3)+f(b,c,r,v,d1,d2,d3)-f(b,c,r,u,d1,d2,d3)
+    -f(a,d,t,v,d1,d2,d3)+f(a,d,t,u,d1,d2,d3)+f(a,d,r,v,d1,d2,d3)-f(a,d,r,u,d1,d2,d3)
+    +f(a,c,t,v,d1,d2,d3)-f(a,c,t,u,d1,d2,d3)-f(a,c,r,v,d1,d2,d3)+f(a,c,r,u,d1,d2,d3);
+/*
 	return f(k1, k2, l1, l2, d1, d2, d3) - f(k1, k2, l1, t2, d1, d2, d3) - f(
 			k1, k2, t1, l2, d1, d2, d3) + f(k1, k2, t1, t2, d1, d2, d3) - f(k1,
 			s2, l1, l2, d1, d2, d3) + f(k1, s2, l1, t2, d1, d2, d3) + f(k1, s2,
@@ -271,8 +273,14 @@ double inline quadInt(
 			d1, d2, d3) - f(s1, l2, t1, t2, d1, d2, d3) + f(s1, s2, l1, l2, d1,
 			d2, d3) - f(s1, 0, l1, t2, d1, d2, d3) - f(s1, s2, t1, l2, d1, d2,
 			d3) + f(s1, s2, t1, t2, d1, d2, d3);
+
+			*/
 }
 
+
+
+
+
 double slpADLO(double y1, double y2, double x1, double x2, double a) {
 	double G3 = 0;
 	double gL = 0;
diff --git a/src/SLPrecangle.hpp b/src/SLPrecangle.hpp
index 6f80150..0c91548 100644
--- a/src/SLPrecangle.hpp
+++ b/src/SLPrecangle.hpp
@@ -18,4 +18,10 @@ double slpADLO(double, double, double, double, double);
 double compute_g0(double, double, double, double);
 double FLO_plane(double, double, double, double, double, double, double);
 
+// sol = quadInt((F_par/F_ort),s1,s2,k1,k2,t1,t2,l1,l2,d1,d2,d3);
+double quadInt(
+		double(*)(double, double, double, double, double, double, double),
+		double, double, double, double, double, double, double, double, double,
+		double, double);
+
 #endif
diff --git a/src/build_A.cpp b/src/build_A.cpp
index f789985..03c5054 100644
--- a/src/build_A.cpp
+++ b/src/build_A.cpp
@@ -4,35 +4,13 @@
 #include "mex.h"
 #include <stdlib.h>
 
+#include "SLPrecangle.hpp"
+
 #define my_PI 3.141592653589793
 #define EPS 0.02
 
 using namespace std;
 
-int sign(double);
-double arsinh(double);
-
-// sol = g0(p,y,x,l);
-double g0(double, double, double, double);
-// sol = G00(p,y1,y2,x1,x2,l);
-double G00(double, double, double, double, double, double);
-// sol = F_par(x1,x2,y1,y2,d1,d2,d3);
-double F_par(double, double, double, double, double, double, double);
-// sol = F_ort(x1,x2,y1,y2,d1,d2,d3);
-//double F_ort(double, double, double, double, double, double, double);
-
-
-double slpADLO(double, double, double, double, double);
-double compute_g0(double, double, double, double);
-double FLO_plane(double, double, double, double, double, double, double);
-
-
-// sol = quadInt((F_par/F_ort),s1,s2,k1,k2,t1,t2,l1,l2,d1,d2,d3);
-double quadInt(
-		double(*)(double, double, double, double, double, double, double),
-		double, double, double, double, double, double, double, double, double,
-		double, double);
-
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 	//sicherheitsabfragen zu Datengroessen
 	if (nrhs != 2)
@@ -138,8 +116,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 				if (rx = 2) {
 					//printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x1[0],x1[1],x2[0],x3[1],y1[0],y1[1],y2[0],y3[1],d[0],d[1],d[2]);
 					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
-					tmp = quadInt(FLO_plane, x1[0], x1[1], x2[0], x3[1], y1[0],
-							y1[1], y2[0], y3[1], d[0], d[1], d[2]);
+					tmp = quadInt(F_par, x1[0], x2[0], x1[1], x3[1], y1[0],
+							y2[0], y1[1], y3[1], d[0], d[1], d[2]);
 					//	 printf("%d%d|%.2f\n",j,k,tmp);
 					A[(k * em) + j] = 1. / (4 * my_PI) * tmp;
 
@@ -175,198 +153,3 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
 	return;
 }
-
-int inline sign(double x) {
-	return x > 0 ? 1 : (x < 0 ? -1 : 0);
-}
-
-double inline arsinh(double x) {
-	return log(x + sqrt(x * x + 1));
-}
-
-//y-x muss != 0 sein
-double g0(double p, double y, double x, double l) {
-	//printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);
-
-	double sol = 0;
-
-	if (l != 0) {
-		if (p == 0.5) {
-			sol = (y - x) / 2 * sqrt((y - x) * (y - x) + l * l) + l * l / 2
-					* arsinh((y - x) / fabs(l));
-			//          printf("%.2f |",sol);
-		} else if (p == 0)
-			sol = y - x;
-		else if (p == -0.5)
-			sol = asinh((y - x) / fabs(l));
-		else if (p == -1)
-			sol = atan((y - x) / fabs(l));
-		else if (p == -1.5)
-			sol = (y - x) / ((l * l) * sqrt((y - x) * (y - x) + l * l));
-		else
-			sol = (y - x) * pow((y - x) * (y - x) + l * l, p) + 2 * p * l * l
-					* g0(p - 1, y, x, l) / (2 * p + 1);
-	} else {
-		if (p == -0.5)
-			sol = sign(y - x) * log(fabs(y - x));
-		else
-			sol = (y - x) * pow(fabs(y - x), 2 * p) / (2 * p + 1);
-	}
-
-	return sol;
-}
-
-double G00(double p, double y1, double y2, double x1, double x2, double l) {
-	//	printf("%.1f | %.1f %.1f | %.1f %.1f | %.1f +",p,x1,x2,y1,y2,l);
-
-	double sol = 0;
-	if (p == -1.5) {
-		if (l == 0) {
-			sol = -sqrt((y1 - x1) * (y1 - x1) + (y2 - x2) * (y2 - x2)) / ((y1
-					- x1) * (y2 - x2));
-		} else {
-			sol = sign((y1 - x1) * (y2 - x2)) / (2 * fabs(l)) * acos(
-					-2 * (y1 - x1) * (y1 - x1) * (y2 - x2) * (y2 - x2) / (((y1
-							- x1) * (y1 - x1) + l * l) * ((y2 - x2) * (y2 - x2)
-							+ l * l)) + 1);
-		}
-
-	} else if (p == -0.5) {
-		if (l != 0)
-			sol = 2 * p * l * l * G00(p - 1, y1, y2, x1, x2, l);
-		if ((y1 - x1) != 0)
-			sol += (y1 - x1) * g0(p, y2, x2,
-					sqrt((y1 - x1) * (y1 - x1) + l * l));
-		if ((y2 - x2) != 0)
-			sol += (y2 - x2) * g0(p, y1, x1,
-					sqrt((y2 - x2) * (y2 - x2) + l * l));
-		sol /= 2 * p + 2;
-	} else {
-		mexErrMsgTxt("no case for p defined");
-	}
-
-	return sol;
-}
-
-double F_par(double x1, double x2, double y1, double y2, double d1, double d2,
-		double d3) {
-
-	//	 printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);
-	double sol = (x1 - y1 - d1) * (x2 - y2 - d2);
-
-	if (sol != 0)
-		sol *= slpADLO( x1, x2, y1 + d1, y2 + d2, d3);
-
-	if ((x1 - y1 - d1) != 0)
-		sol -= (x1 - y1 - d1) * compute_g0(0.5, x1, y1 + d1,
-				sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));
-
-	if ((x2 - y2 - d2) != 0)
-		sol -= (x2 - y2 - d2) * compute_g0(0.5, x2, y2 + d2,
-				sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));
-
-	double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2
-			- d2) + d3 * d3);
-	sol += 1. / 3 * hlp * sqrt(hlp);
-	return sol;
-}
-
-double inline quadInt(
-		double(*f)(double, double, double, double, double, double, double),
-		double s1, double s2, double k1, double k2, double t1, double t2,
-		double l1, double l2, double d1, double d2, double d3) {
-
-	return f(k1, k2, l1, l2, d1, d2, d3) - f(k1, k2, l1, t2, d1, d2, d3) - f(
-			k1, k2, t1, l2, d1, d2, d3) + f(k1, k2, t1, t2, d1, d2, d3) - f(k1,
-			s2, l1, l2, d1, d2, d3) + f(k1, s2, l1, t2, d1, d2, d3) + f(k1, s2,
-			t1, l2, d1, d2, d3) - f(k1, s2, t1, t2, d1, d2, d3) - f(s1, k2, l1,
-			l2, d1, d2, d3) + f(s1, k2, l1, t2, d1, d2, d3) + f(s1, k2, t1, l2,
-			d1, d2, d3) - f(s1, l2, t1, t2, d1, d2, d3) + f(s1, s2, l1, l2, d1,
-			d2, d3) - f(s1, 0, l1, t2, d1, d2, d3) - f(s1, s2, t1, l2, d1, d2,
-			d3) + f(s1, s2, t1, t2, d1, d2, d3);
-}
-
-double slpADLO(double y1, double y2, double x1, double x2, double a) {
-	double G3 = 0;
-	double gL = 0;
-	double gK = 0;
-	double tmp;
-
-	tmp = y1 - x1;
-	if (fabs(tmp) >= EPS * y1) {
-		tmp = sqrt(y1 * y1 + x1 * x1 + a * a - 2 * y1 * x1);
-		gL = compute_g0(-0.5, y2, x2, tmp);
-	}
-	tmp = y2 - x2;
-	if (fabs(tmp) >= EPS * y2) {
-		tmp = sqrt(y2 * y2 + x2 * x2 + a * a - 2 * y2 * x2);
-		gK = compute_g0(-0.5, y1, x1, tmp);
-	}
-	if (fabs(a * a) > EPS) {
-		if ((y1 - x1) * (y2 - x2) * a >= 0)
-			tmp = 1.;
-		else
-			tmp = -1.;
-
-		G3 = tmp * acos(
-				(-2. * (y1 - x1) * (y1 - x1) * (y2 - x2) * (y2 - x2)) / (((y1
-						- x1) * (y1 - x1) + a * a) * ((y2 - x2) * (y2 - x2) + a
-						* a)) + 1.) / (2. * a);
-	}
-
-	return (y1 - x1) * gL + (y2 - x2) * gK - a * a * G3;
-}
-
-double compute_g0(double p, double y, double x, double a) {
-	int sp = (int) 2 * (p - EPS); // MK (p-EPS) instead of (p+EPS)
-	//printf("\n compute_g0, p = %lf, sp = %d\n",p,sp);
-	assert(
-			p == 0 || (p == -0.5) || (p == -1) || (p == -1.5) || (fabs(a)
-					<= EPS));
-	if (fabs(a) <= EPS) {
-		// printf("\n a < eps\n");
-		switch (sp) {
-		case 0:
-			return y - x;
-		case -1:
-			return log(fabs(x - y)) * (y - x) / fabs(y - x);
-		case -2:
-			return -(y - x) / fabs(y - x) / fabs(y - x);
-		case -3:
-			return -0.5 * (y - x) / fabs(y - x) / fabs(y - x) / fabs(y - x);
-		}
-	} else {
-		//  printf("\n a > eps\n");
-		switch (sp) {
-		case 0:
-			return y - x;
-		case -1:
-			return asinh((y - x) / fabs(a));
-		case -2:
-			return atan((y - x) / fabs(a));
-		case -3:
-			return (y - x) * pow((x * x + y * y + a * a - 2 * x * y), -0.5)
-					/ (a * a);
-		default:
-			printf("p must be either 0, -1/2, -1 or -3/2.");
-			return NAN;
-		}
-	}
-}
-
-double FLO_plane(double x1, double x2, double y1, double y2, double delta1, double delta2, double a){
-  double yd1  = y1+delta1;
-  double yd2  = y2+delta2;
-  double tmp1 = x1-y1-delta1;
-  double tmp2 = x2-y2-delta2;
-  double tmp3 = sqrt( tmp2*tmp2 + a*a );
-  double tmp4 = sqrt( tmp1*tmp1 + a*a );
-  double tmp5 = pow(tmp1*tmp1 + tmp2*tmp2 + a*a,3./2.);
-  double rval = 0;
-
-  rval = tmp1*tmp2*slpADLO(x1,x2,yd1,yd2,a)
-    - tmp1*compute_g0(-0.5,x1,yd1,tmp3)
-    - tmp2*compute_g0(-0.5,x2,yd2,tmp4)
-    + tmp5/3.;
-  return rval;
-}
diff --git a/src/build_A.mexa64 b/src/build_A.mexa64
index 89a686873c4c08e9b820971886b0aede28b6eec9..ee6b1c6a38f6067822ca43508fe35aba8cbbd3fa 100755
GIT binary patch
literal 17398
zcmeHOeRP}EdB3uv#7S`UiZ(PrL$na6fmm_AAbbSuI7(!&^Fb_-f*nP%<&PfQa`gcw
zlwsnyjZoaWE41#$O17S{tnKNtj&a;>ZS5o_4rFI3Q0hWY@o7mX*J&403LcbtzvsF4
z%CBBaXFJ>e%=Phm@ALiK=YGF>?_=JYEk!n)(8(_D6a*ccFEJI0aCEgQAgd4?#5_Eg
zimRB{F+-YFT5-KZ>0AgmBgn>F{CAZxh}_Kuof~$EFIBmgD*~G2jw-pMN>1knF(B!@
zSb|~>a{3=tRkh|T6rawCm(q@Ni6$ynD1tSwXD%UhuH|+@j{N+eChh0F%3h(k45<3)
zTvdnb_E8w4aBbb@7pw31%?Gdj$zLAWAANdyaM>%bjSZ|hMZS<Y@#***6##nNcu4%W
z?1vVvfBY-YGorx5E*6ME?XcpBBfmTk|E>b~?-jt4UvB(c(=s63RwUp@@So4_Hwxf?
zTmb*Y0{9aJ@CNY9@o!C|K;-l1Yyo^t0sLJB@Vg7(hYH}&6u{RY&@S=zG{`w=EuhQA
zW#ZEdBznU<MGl2pI-;GScp@B2ghC<|YU~b&n%g?V?QQ!bB6RnIp?f1OZSh1TR@EMk
z$0Kpk5$WEN>}*W5MLVhNjm2u?Ep^=q(b2s%lBkv0He^+?NH`JM98K<PkJN?}v9|8e
zE#>8sz9S|=ySMF1Y-o$dL;J$<$hzvzwnWux?}Pd9n;FlpL1TAYLc~86ONem1t+Q2x
z6X8zL6@5swM_WX=G0Gn9h_!VlnniotzK(8BJj#VgwiHe%Ts+a#c+*V)JU25-w$#$t
z7>X0Ln($qTR@9}5;BtZ+CEvbeTYFQe(t{|3Z(b*Ay#cYcrg~FVXpLu$XT2;30-?2@
zb=kI%)7~6GbL$+8Y6t%HR7j;F9?3M0<==|RJLl85{y6?O(jY>jtO(DeJLT9W9JM<o
zl3I&0^<y{wlUKE2#&?JvCi)+M*Yv~<e})l1N7S3>3za<Bl)ODF->v3lJxx64A}&ld
z7JR)LL_@3GOxOIwA7fCATkzfGl27NFU&F(SpXOfkx-@FSo9E@TCcNtMY<SFszf8I!
zpUOTB3pEupzSxA<?L-8_gx6~W;g^{3dWaF;ZNlri5x&%fS4%UYMY##TkTLSH!Gt&a
zUtz+RnD~7r{2~*6hY9aA;R7Z-^%0%+n(!`-MBZw`o9D+a6aFJ6eqZnT;*1j=Ec11G
zoV!v8U+>A{6VvN&`-EFeUv~$d&Se!S5!_0W6XVlJ*R3I#dg;VjhHoU8nrdR0;pGHV
z6HOdr_$q>_X(o;^d^y3?btVQFUO@0%f_oS)BAA+NqKo18uLn#`HL;iBiv&{>P3&O!
z*9247ny6s-7X(vNO_VeII>FRL6K;lICYYLL!eIDCf~iR+1YnAbFYRyu>FXO#6iuJy
zung$r@7Z`1G~r8~tE}6V+~@1vc;Pz0iOYQHjh_W2a~o;tU-}o<x<yU@b%&8b-=m<^
z-0(}^$Om>`@1)Il>sZ|FT(+5lK5X-)91WmXPJi?ySnov<^w_x<9CxB<d!VxZ^oeHY
zvV(Au_&55XQuyIO;wrE=!a8}n$k)Gd6)2~t$En>;pBNY*p33_x@2m8~Ro~<1#~@#s
z`mj3njxY66@B2mYxqsvBOVIo;{Be4^viC#VpAu#lVV(oVm%csoH1d@9IV##v?47c+
zR1NVAFwY~!nTM#5>MKRfQ{KbAG*GX8IKzA!F*kzeyHta8@^I?COcTKB)SMN*)TzuM
zcp*;BOLkHKt5XOC)I7>0?;uNlpL*R7sB2}eV~eLU4U&GuN_QywMoI6n($6veRAwpX
zeF2!ZvK&_|@Be|l+mMIl8O8DwW;tQi(4)vhAXle;XeEag`72C*!Ajnw$opBjci2i_
zuIRfZeW>@*!@~JF4@Oqn`<!h!<vo`6AGW=mdTJfVwRX^#>ZcKa%!lk$Umc#*Ainfd
zgt+KjS>zk6PH#Rr<xAHd@uh6O-jQNos&;@VcJu*F^n<vm*w=g9@TH!jLF!AhO7Pk2
z-eb9ZM!KK-MeU(9TTFSsX&dzo`qKVEU;murkTHCXFIT7QUL~_VzLbklcK>r(Y-(pQ
zP{O1+RcWOAs9LGMc08qGJ7WrI?x!Y(a#^ST3P4snX~knQcd^oG)SaBV=ydHgd};5I
zi%z2g&KN#G{vH%|1F9dx|Bo&@gRc<60PxdhXX#4-f=mT)0@M)@&LUFOk)ROOG%CPs
zkUo|g^QE~X0wyI|<5ebw0V2R$@Fgfk?CPP^Rf!49Pp%c=|84-a<k^*y)7U{iXkF}(
z_HKz)7wjwNDh%|U)Ki~^t2OEGQLj1WJwts30a-auH*m_%ps8!5s)`jV+Y3@vM%sJE
z?ma{;gHG5#2R_fVd3DM!EDZq^d^R_Db;`?1xV52YmX0fCskm=Og>J-3r|UinIDB*(
z_ZWn3py13jA7LoMGD3j%Fg>>cy~Bt(>!HU&RzpH#RDvTa*rW?JGwG5{xoyY2X(9CN
zqQEH$Ds3UlXyM=q#uM~9{CI_+YqDyc!52}cPEcVqJMTZ;?DUN6zZOr&Ti667uLN!w
zP(`ux&yPEwe|d-|Q0M03la8SxYD}kLmtCj62yHV2H8$X?s6P!rFZvTY7h%xd`)Qy9
zX#<fsHyK!T5-4n50V<}0)aDfsEb?JeCnx`u1mj$*4zhGDRLD_hT#2Cs{e~Vy3t;>R
zV5r`B@EW(sRAagKHG02^>N~DQMKD9{m7`|P@yU-8g(?Px8avYKJ!7lxJyo3Uqdsf%
zelzVqbF2Tb^Yiy$mh9z;2h&*}&2z=+ewz7itp(Ng0VajVPmUoNSY|65D)&@|PM@H4
z6YJoQ$d#4v(e$=>G0MLF_1~u~z45)JpezSP^6tW>!G}nAGkDaZJGuO4@n{%F=$Af#
zUE2P+&VzR$;;5rHL#^{Hkj^K*Os&!P1=$iLS6@qVR5JYIsQk;U3eAodruid)4zY8T
zSV+t1M;2tt-~?c@Y&ecVD~mQe>%eJ{qCKM|p7Fm+d`3x`VgV&clY$?DL9!YoJ_rRy
zCN@Y<{bXbD)7X6OcO{(=_L7I*vlC0qHF=BZHzjYFp_7^GO`f2K?itCv`upkWiJznY
z5hH1Ke48ZKCRgW5o(C2fKVnxWC0~am{By)WnRwc)p_?@9Ol}8%!RS#KPorAUwQ~e!
zhaZ{G0eMGi^@47S)5IN8zn`kgjp(PiK$r%aLEOP%3KB_WG9gyt@(vLVkqm(5qWOWL
zddHjC(F150$Ah%NQsyOeNGhQjgEZs}Klp>(-29-1Tl_2w28gsrHp<sjFQL$tL_yxa
zHR%PbYLtJbkxfxOkDcs8E;Sd7I=^(n`TQH$aDubs3y?xoM>56i_KjGff+Lf)<O2E5
zoec%8lB$9*7>=N<MeP^~Q#3Ly?{Jeaxt3&0eh;(M&e$yb9>_N;sAE%9T@9$XD}su-
zI%s@2{t=SSCl<kfXJ0X%z?FRhxc&Df>(B^r(nbBV{#6u6K@gK1b<`lm1-f!t^%{1}
zD0v<V5Cd)|>fAgMMkY#C&nA6-$f|BLtTc{fqxdY&2OPyW;Leyiiqw`tVA%!7BCa;M
zLc^2)5yk*E0^I3-iu`qi1<QZH@iKH~HvWzxcO&f-N^C5Q-hsHHngMjvJv1hdPxQ!M
zD96bG$T^>=;{a`wgj$GrN^BtPFG_5It`2x*P*4pDuR(D#C~AhIj;wZIHw>WU8306v
z`2DI+2dOz~k4<7!-2VoO4adOXr2u(};vJ$k`ZA7Buv7yQby=!g&RE!quHkv(k@=ag
zUYwprY{)$IU5+otLguyiWru-qzzYj5dg$nr1Eip&+++qS3Y)pbWQLm&{*sxfG6P|g
z89g!sjusdo*I<0lS!%<{D>2WX#$j;?!$)L7Z*!a+tH5XK$(Y`Ekk6=*n%+?lR|AZ?
zB2eY(phzXh5i@jz`jY^>Q_zkcQID3d{~;kRP<^LRDULUf;!4;>RBr=CwjUG*B2=FS
z=nrs9$d;r~I`)G-^EGVtp@df&p7hC*$MG)Czfct*Bd}CBavxrpX&^E}p)hJ$|K+5=
z>_4S`c!~)Q&P#Uzc`uNmev={a6gKpv$q>~mJKPM%-$6&_<dif7Z!Ly6O5DS!CUA{L
z`Am0XhhREq++PoCq=;N9NkKnygg@ZT=(*_e^GDIbzSIZQ7`jp}>ZDkw@rVrJ(J_Wc
z5+QLt3ZIRn8;_19cto(fDt!|L*&Yx;cbonR7`;bPHZuY=wVdL_cz1b^<l!VUw4E|)
z&76aIIN?<u&Vf9ftP?en?Fv8H`k_3WtaBD;PaaOz$$k!ls&w{s^QJ6+^{vd@Tdjeb
zah0n-%C7dwoo0Tu?5&x(HNOut4x5?felfphcG$$sTZ6(54^teJ!-4!hSi&ordC|<P
zb|zwqlj>mQWStgHqZXa8p0?2SwS(yObYicsMS}<EBpB$I%~a3hS#9G3{W5L=wS}&y
z{ukf_Fv*V51LgxG*d5crQENJ8LN~fH%`q<j6x=IWF^>it9YJXGvl9Qt;!@zSroDm*
z<3;CzQHb|W{V@^2<FQ4(Qy)6}(pbhO9%uY?;$wZo&c27`!ba0%y?1KjlL~~1dYX9Z
z9IxUHVP^YzIn0TR$2ztnbtT^ciS`~7>vEKkcV7hwETiPrUvv3>OwoJbE!Z8i*rQqY
zpktxZnq_}Y)#Gle$3m_LMP_n)@6>c+Y2Pqj(cyibvyXo1f<0Z8XL3MjJDwT+HTJ@v
zosGW${A3JM7cayaD49JMQlPP@FokvuL23&MzflsztDnpn=pgb;h+xv2nY<bcT`?Hs
zLieELa%=%+67LH6#M67ef~G0NTs81nl;n1e`U3(d&lx2?NU5k2chIqn5uVXlW7cO{
zT|vTsfcA`MP{zO9I75bA#}<S<&+Z=f!YlB?m?4)7RKjo^n1gCP#=B-L!MssHFXPo%
z1sToj;?Y|2E%0j3ARKh`F#?`7tDwjg3lcN?c-SMX^bJxntYFD;?!?55a|jGKfZ;fY
zw~?7&@S==0e#SZSoJKlX;oShA%^mN8-T9e8)%0lLARj5bwEc0c1~EnlEyv%YL~>^*
zO9(uqK~hRWjxRBWw*gASIB0yBqx~ljO&ZZ$FYZK$HFf)d;~t1+p3lf$f{M_FgqJls
z0PLYtjy#aS4u_ZZ)4?-Hdg@2Aon<e1h<ft&+?B+ZhR=_@{h>n>ZF+g@Gpsbt?2g7k
zhBGs}qim#xKzRYS@XtMw_nWh}`M0q9a_!Y;GPR|IkF(f|75fo#WC#qj`wVbLr^ikV
zU7QOS5Opf^#=lc&f?28&p-41PC1j^W%n+E&Y1QaVJ({H=Aco`LFtcRN@Mcu*{2Iu#
z+!2$w<6Fe=4KSd(qd=%$Kz2!=JPkq>zX9Uem_|4?z@C!t^0XDu4bi``Y>``J4*rVY
zQ*@N@dy1=o!h4GCfHM0*p<e@JqLEJTEAXnuF}?s?G<#teyz~_S&Siini}9OB8b4D^
zO()-`A3ren4-Ak<W}1x9uSDp;Dn+`xB8`bile@Gr8jUsK_v%C>Ub?M&?F}m<?U9a1
zXChu&3+M*X8jic8o$U|f5phS8iLPYA9gejmiAQABn6$WA%?wHFbF)J2fXaJEPj{4u
zp1bRzaICYfv&G%k>E2piew({9>TblZ+V19P%-!W~ir^RZNE0@QUG9!#JmKCKaYxz`
zt&y0!e5Jc=^{O@Wbc45S?W#4Ntc~%zWO4ilGQXhmDGe+7)MH9ySmhOdPCvfu?bcr&
z`lqS4?v(iA3Pl&#IZ)Da8r`D{wyW*8)3+~hS=@o}cOy^yW_p^Io+HR<b$uSY-wVil
zaJXE8-mwYi`whsSz#e-q@@8~2I#PZSUG@m_kK;HvjGPkvrcKGV|6XD1Hf&d4Hm`U9
zw3WcoZ{7cL4!W?Oi^@et)n|%I=I-j9cX)2U<8yO9Z@0VWfg3a=;;q5|Ni5MMo8@nj
zEGvL&$NyQ-ufP<Pr8oB4?Y}L`>LI=-@&En*n4W$S$Mr1Vr@;5*3?KUj@?IPTk5gT=
zfA;U$NEzw*8UBm#M!;vbcMtSzuxIU&o;i@Gw|PIy(eo45L-ZBE{~WO_Gy7vNWoL-K
z8T3iSIF%!RFe~2#d>vw+FD?+jlrBZ!QUoqV;8FxGMc`5d{*Mum+o|bk{XLxi{!OdZ
zbPemAmTpRK%1mfDV}UPP;6GSk`eu)kmS1dvueHGRjUc7VWS0Hz7jNo0jV_n@@BcU*
zrP2k(j}teiQXFt8;cW#c`f-jD?H4KWms-g2a+lK*mD84pQf8H#b0wDL*ZV@+`cYD!
z3yN}u(D;BdthRw-OtK684V|`!?GO7utQ^OAo?&elTT+In6dmovu=W#MT!se}9q*|b
z)^Si~|KBqGedK?qP=6PxzmK=_zn}v1n#zBr@^@7}R|V<{m0z#&n^nF=<@c$)S>+F^
z{4*;5lFFY|x%O*oRn=|oQvbeWXCk?(2FJ5l)V<EL-m|7`Lz0o}Zji7Lp7_HZiSRz;
ziI~h=b*VF&h<I8$lOF!kwydp5FoYvrtMD{E+zCmUCt{NFU?di&?~<8K0U;%&L9oUm
z?O~EvrLOjb@U-FkW8uNq#dzYh)&^R%3EvhAPoy=}91C|uLaj}Z)sXNsCZe%8?91$v
zjUYm2xTCER_M=dPhlh1bJ^SL|Zj5&D$<x|k|4aC4$kCY4{khQn2d!QWTMY_2YvMQh
zsFv5`T<4{VQ;!GLeX{Z!P@*xe<@G$EbA4Zf>=KWb*W<GSIMS`<^}M2UJ&%BelbVJ!
zD10c<nAP%noa&s<KH$;*Yd)R(L8CcA!+M_3In|Gn)qjZ<q6r0>L$th}M|8edF_Jv_
zpzE*DB94OAO*lQT>3mEvYX7yoj^Cpec|A|+T*r^(shzCyp9PL$s^@Wi|3l~64p=ym
zujX_BfLZ>SGNAM9f~17BJU*e|)Q@M5yqZV!J*6f1K-iqLoW{{xZ?n8!w{(u*x^t4Y
zuaRFtf#h}jjW3rJo$n>0DOv6R6=-JpLidQY3akCUx5yt^Avtw^&|+V=pVs#sC9nI3
zzGtQH8M!s1%B}Llz>^6rKTzfdB*Z`&XM9gUmgLhae+p=`{1GL8M9EwEt?(JhQ{U6}
z^?Ii7ftlwI@<7Y$cA$GdX8ZcS!O$w=#&`5G%RaTd&R++`EU)(^86`iEm6Qc7r}JN0
z<n{X3tqf~Et;PyZfIv81|DJL;v-T)?>-e+E{}$ZjzqOuP{zld|^B`VII>!f7S=&rd
zkFP6*@0pP(YdMW3SdSfzqdc=gLdVpGjow$$OGxd%t``k=U3V!VpMw=L-(w+K?VDw+
zWjnmp>!F!Ow`}zKmW#Iwy>8{==Lo$X<>DPeuQR#$xk9fWx%hcPuM4^O;_SSii=Us}
z_vYef-_s+y9W#7(-Zcr>vCG%<sR<4K5~1hGT>K)T=f7OMGkae$7jFnXU*+Om+54Zl
z_{%ZBxXmao+r>wOo<Gb)JKh%Sc_A0SSm^Pdi@zd!zR1OYOz82Mi(ex2xXZ;`79S90
zf1e#s8b@>T+`FTGZ^Ml|JuhqhG%jtpk!Otq4I5jYICz0CLOY&S>!a550Q1lIugAwH
z75<Q#PxQF!R(f{mgPKxQ$p3?^$0ZI~+UpC9H-!6or9frRFn)GBegk;&e?YCTT51G^
zd~y3%#jnrPE0w@o$ZcYY7*O@n_5ED|J$5jVU40(6>Uq+l=NiSY&td$IAG$pX{~|Rg
z62G%YUdi~`<7g}JHuML*p0E%e4NA{A1ssWglR)00@XubuXu<FNQG7_@-&F@={%rzy
z+80p04yZ25zcm2=H^AF)YPR<CZ!`XKF=%PW9{?}lr?vmTr1ZbyCb|&(8w2t;75|`$
z8x!$(tAL(K;3*E)@%s+zpFNMg4?MoSo|THQ$mGM%FM$6z@a|b!gt(~yo<5V#r@smK
zeD`VT=2<@e0|oGh3*i5u0R9B<%k%jEW&!>S1@J{!KlAxHpYgNr11~PXzq|nc_5%3Y
z0(j9FZf`eTITnqvL_CpfZuT_d?r~M!_IpD$)w}BWCi3PS=<X^LsVN$2X^-v;w}<%j
z9SSGA1y0so?GgMY?74aU`db7MLTydm1l=m+0lq2O(ebdh<K4EI&1|{1vep~Y5Sh{C
z87q9$XDj;7Z(~=N2yLy|zNxY%6siv8pm*;|)NBZa+M@hUVOPzLSfnxB+0q{ISfpl=
zcCX(O>I%o2WWrVG)#X&JWFiu3DJP_$yV-47%iFhX+2ySZ)m3h)@#6CK?lqcg_bu`E
zuFB0d+gTGRYw@#mR_iUFN`{-NI}<y$NK{jkR~L{YKuz&zs5RW#)Q&6QA#Mj96THNb
zm&Kv;)-5&LLtX7*{Io5T2t_&}9^v;Lq0PZ<m9^DXBE(_f+u^wO&G|NYce|p(d|%pn
z6+O@0^8$Cqm7ajSYfbmC<vn(JSv=PYE^DjH;dzwNMe|%`bh|uPLGBg#TpF&S=OXAn
adM-xak}rIVzQDbCy7Hc@3K!mU5&s3vbZ-&>

literal 17226
zcmeHOeRP}Ed4FX^iIZR?rvyTPMr#sg+C*^zfk2mH$5A3PJLJO#DJT(+EywP$E!UF3
zByBM@h6tnL71EP->@@4TF6-8kKU$i#rLmKgB>i-u?MY$hG&?1&JTn?XDP9ZJe!u75
zE5DMZZfAe&pK*Qs-urw%_qosgcwfEmgZ_q1We$gBiqra#MQC8I!qjNOk=43^s>WJp
z&B66LYq9XUrb#nek6*7)rnW4vAZTMY{s&eHh}tU+Q(I1jU!-dzR|YiYj%vB1TF%rS
zGN72cT!Gds)cilAyINQ;(|o2TFU!^oi8g63SFlx3H)yn}jodcK(a-;BF@A2>_DaQN
zK=;qox;xUhD`1S_+PtON>i^LfzWAfgU2nYXc{((@<gun@ANfg^zEGTeCjOop00UlJ
z6#gsD1M}BD^qFS`QR3pX=2;!aVa*dqy`TvHjuQBXOW^637ypH22@u{wB;e2EznI+@
zO5hKbz<;9z{;?8x5Ag59e_?q7h+_UcSpv_zDaL<$3H;6y`28jD110bc2y~_Oy5s~{
za^*&}7Fd^ATjwd%=jUruB(k?V))R>*TKW=^h!u&n_P0dZJ9}EXI`>7b$Q}1Zwnz7N
z#uL%Lx~`UZJQ}yUqy3wbJ*|n(SP$F&zP`rz-lqP9)!n~2nrKwj7F2b8(UwGXV=TF+
zE85tS=<Do{+)`bw=v(`&_(%E@R#$AV)e`UQ>9AT7Ej?Cm?0&1IH6}*3_I36o+O5dW
zEjtqHI%Dz3o|bs@ruv@FMBQrteZ}w_1us2-;r>pj@7&Yf?~BKz5zVz)5*imzw6%WV
z0|0zC3rnuGx3x7AC$yUQ-b4rb(nh$NaI50mlkDtji`4oMg_fIdvKswCYjZ>WhPuca
z-x}Xq)eZ(DH~MbMje}bDW?5WYXJJ;m@Nbqvw#sBCa~`X|Wwjrg%Xxh*{?~I7p|C2$
z^~g3gw~6Dpvry_rw3)|V{3jQ8V#Rk^TW$2Wfj9KT4WGh{pJjz@^!ZvIY+Bx#lRu!>
zWwT7IT1qa<I={**wZZpSE7;V9_)*QzwbZ`PoU!5U>+HA<uTxm;owMOD6O6uS`*Q&5
zrBw*!HoS>D2_75Xto6h%vf-)A6t4|$`bJun4X?La5vjJ}&D0}doegjIzs82I5Dn@Q
zu;DMa;kVlGZW}&m!*k3{xz~oTG)UAPHoSfP>b2qDW8)7DzFMAfBVjUM;^o3(%L)vh
zDnFTB`{6(LTG{1YxVo3rphdWYk`w2%D9iT{X4XxN3%rvshiYO};O&GtL=#5^zJoA_
zX5z5G8wfLXCPoCljqq&30|Kuh%psfT75M#xIaCw(3j9989HNP>0xu@a)S9RfcmZJ!
z)kL+x^9XZ@CcFZd5$4cLcm)3a^?*4f6Bb~`B{1w-3uItuG*Olvm#~bO(mb&K9B5V`
zb)mLtM{-YKaJ`ps;<CW-`mdvvxs6)V*S)yZYc-^opF#zFUk0V&hF=E8{@@G@PC5d&
zo{M|kOExm-qmDqzl?J^wd&McRegaL<`z}21>Os?SPi^?@$#(Y=?$G34KLC|hGaN`P
z2K$rC!46L^52V-c1?6n^Jjd<qN#!y5YVWSStF{^L1|E8qBUhVxr#|)DK<aygzb}K|
z>Gii<jTAWX$85HC@Eym;iTM~YM}Y|pe=u_hb;|z~oA11n%??gEMe08CJT5%@%QJB{
zQbYaZ9lD7z7zXwSHvSg!he7AwZVi8pOsTgrJAtWB%~}>noz8p`>`+oJ;vdXB2(Uhd
zpisXeag0U&+f2dZ%(KiQF?~9-Rk1KF3R&J1(&@}f#Zpy-<ps?$Pq8o=3bj0=S$_St
zxOu?z#pcW7n&p(RoGjGPt;v56a((JMh2$-o{D_dBFC?$h<h`QYf25G^*7VJaet7UX
z#@KRy@)69JVgFN((UgB=_#lyv7g7fwX55C;3~3;h?v*H}IoPL%-o}-q6c|28#3lEN
zvcRGG;f<%JAiw|zQjWmjSUKd@0V3Ur3^6PR70ZV5z~Bi_Aa!sNC5A;eSRGFP$P8A`
za9VOFT@Bu0aUeDHQM4Ur0*By+c$uDc0`i_f>kIY6O+TS89LYecl2~W+Q#nFvTRBjc
z%}}b=Gd#o|r-rV`S#S#GPw8hGK)>$mh>As?so}KZv01!S>-IFAn!4n!+~x@k4>1^*
z+@2b^<Ou+AkWQn$6NvCR{QvlpJM<zE9soaYcUOH6KuD+nPJ%iH!Z^Cil!IFI*HZ&#
z4^@iPIoP2efT_f2Tx?SqBmw3^--A-bF$|@ZwOH+Bv>;f9|2qInP|?0x@q4z?6Qfbw
zGCnRU)Ce!Mi!~mVXsLtC;8w%%bIh;P{!`2<#An4ElL0BGfVT8W)s<^hE=^K(7(Hf&
z(|?M?rII#13tgB_hx$~rRgefk!FhWk)TjJV#S~DgX6U$bhKjqURhW^jbywcS_`_4T
z=Px}(?_!juN0S%^T!jkIZlULor~e3sK=kl3U(`_O99wWiN1eLRGn+2Ols}Y`U-O}7
z2g9c&*xE$P81T?Z!4rB3e!NKNy*ah+(09>hPSkkX-FKgDcl*ZnEydM!4hM(i)xeDc
zswsDW_JsS{7mjiLb#FW|={i=%A$5D4xl9U|1%?|BLz4sOI%b$!B*MrrWE_#025GJe
zFbSqg|CzvWGbJ{7hDO~F@1G?aH-NR$e+>Nd&}ckrI#cOC#AYL`RyH0G8)yx6h(3hT
zU75ZYZ6u?bupTlF^(vH}d_?&Nz@l7JECI!$C#*;XmD1!HnCYe~(WAKRI)qCzO{b|A
zwSOS90;dQXx+s&>Qy}s;+Py-l-{CAmsZ3ERF?5i<$+tDCrOp?XGPa5cKB@)J6&1Ag
zi?LAsYSJR(MMaD))vsx(_iCv#MWyolH7%9FNsZAMEh=Sf<@z-(_>2<t55VVw$lE%l
z2WyC<n)#@<b?i^t${Q~wSf>S#78T6zTAsyuT2Lbk{4DHUBa=*7{4ox+49k(y7W3k-
zr2bk<9WE+m%rzh2bkCV;S3-U!cY%ZXJ=3vQrmD3RT`rI^b_xeg3BHA+d!SWzdWfdg
z+zF@EYzU`Seh1SkUBhXWgTb`Qs$g2>Nw^Wi9!{%)45rn%1ohMmr&ZK~X%&Ysmy9rj
z8^mfCayf9{g(a$7^1X7_iOH+9KV0#!Oc}!1p#PMke(-ep@DR5&hxk0)eCk#?P`RJD
z{y%XxmV*n9Mni`X^YY;|4^p?LdE#<E@yzeACOveDM{peKYj@S|u8o{M$rJhx+;#*o
zzWxj@jpuM_9>-<Zc`Ufsg|Mvbo9%x19;~wnkw3#Vxd}-3W8dT&JoGiS{!%X7NIAAV
z&F9g%7qqGn#y5!decKp>NB7Ug@;{cjc!@I{YcIH6t+IG|8ey^tVevFW$8IWXhP^6?
za(f6?ydtwYPt8AynxH?Ec}7~EimwZks0i`g&<qP9dKNkZUD%0&V-wZ(-rNU)eU(W!
zyT@Mm$0zu}7+@oNTOsrhByX6eqg8QKo1x2aZueNG=k;uM;%A6286{2|sd<JZ^?{{%
zS~5Y&r|pspD7h_pC-_Q6mLV<{*-t|oNPf^ocsKh!Ht`oqKSLXq0d6K<$JRsc(5T8v
z^gH7aC2Id;Vgx~}^y7lAL4hkIi-~6!#5~jnOBvpwuq6RMMgsUlk_vJ`>^%#;A{ZoT
zw;GJU(-|beo(uB6%}GC4^<aDxep}4CA7Yl@U#U+AXWXAY>3;TQ+|GitVkxAUI+-tB
zqzmuI5jHe7*+>`YyCfnSMkTv~FnC;(vT(q63Nsp+xuPH><VMO?#Dr(}7;eP__Y|8k
zlBZBtxeFbyjG|+e-56V+s}Yj!$1aEe?xAvAfm=BjxP5men-D`dS;@Q%{{#&VO9;mv
zCTNIpfv&tnAN`F4w;})qhyiDbBsJ4x_t?Y|m9ZR(Ff(4K&8XHmmW$#6EO`>em*LJi
zdlWgAAz;M?*CFsJS9t!He}plJI|WI)pQC;WmT`X&bTxoC^NZIc^t_d*cv6Is9f&Ks
z8AO`y=3G26algt!H9xKeoBOdQ3D6cr;1b}g_$6)`_mw3!L031t@-V0#hF7Mwbll^5
z4NRGCwYyZTRu(y_?$$XS;&3z`oy2-`_seMRItqC|1H?rKsm(EZKv`-4sqvWVTi#sQ
zhSc!Aynk+{2Hit!5N%0oi7)0t<{Q6fw?jNA2K}&5$%~0T9iW1W1vWEK!4MWV^Yj0h
zHY39be`#iiHUnXs88b73uGjI9imt)<?s1Ocqz6w@27%4Ai%RS!(+P|f-%H;v;?g_?
zxr#B7=LwuMeFfva@CU$Q0AEB0zU%(P_c09zr~a7e#N|QH;M6<qBkrL`F=`WoLdmv|
zB(B2wc_^-j7|O%`sfjiX!c5qI5&B$7dS$^Y7@-T{M$kf-e@zuT#SP@MFHZ?Qyj^+2
zpe(GD_T88#yWug|`MdxGzp_boPmVq0`sy4~6|1n1>%MGcU*@AP?xD{EHCaD6l~tYK
zPc+jCCq9=4LdS{BgOgaaADxN+2jD0B&>_jW43tRsWW?c~S!9ElIYQy&2uT<?rz^_g
zYi1BSNS+Q3O!_mE*YJ7SY|A<ZvsPxM;&QzM%p@MVaWlur)-VYXQ9ZolXsMOy972td
z$wu#~_zk2Ga>kR~A=kGB;XBjUsy)kFkGMR3A8h#N$z|Gg%r>S?ZlYWTugpY|3`QqB
zuGLVUc}O-0Y$Gx!N9hVKdep#MDSK8$v&}O`%(cvT)=tAjNi-p0qA3$j7AElxEGn#E
z$#ZVszfU_iu6yBeJ%uO2nOE@KiY{E7b`EM}%b*p6RBt`5gJ5?*!XFPw&%g`WJ^h^N
zt|}1wc$RS8j27j_Co2fvVn8G0S}T~4O~ZN6xMT6?6qJeC2$b|ol1L63Ah}+`J!+<V
zLS+d$!nID;M(%aH!KU|K=#ub?he7Tw)Ds@dCCk0iWnAR!mc$+`N2F6k^L)Vd{W}Ck
zwI*gzSjbpvs}SulSi-`B{8EuR2FmjQ@GqGtOY4js?X`>w_z`MmH(^4rp73!7`*F>F
znC!>Ez=bv_iO$QmQAB@q|2$HsGM^IWP>$+}G7`Jk3AH<N@ckDvXZ4^9m5wMX2j2B@
zJm1ayP}W^!1u(m4ng<uX0qr8gm%)JUo&i>80eZqARcGL42+HFRImVfZv%4Y$S-3wv
z5L<PS^$8putlTw@3-9<VyP>`^d0wef(MLAsM<0gNN<&_Ca(x<y0&JseK#4-wLxn-L
zjWU4TEc8A$nmPUlj)r)y76|6M&lOTCg5-pS(i@le*V8j=S~-K|S`HnVA>5!){}eSL
zJYnS~AZ`ji*F`M&@Kp?fr?jE3|7zM0`;$BEcwFCtj?A^%5WFogG;<w>54qH1!c~$b
z@7PrjYg9y+D*jnEC7y+F!=a~8di`A$UMy`F_{eABXD`DKBI&htK;b8uj{?fv4~qUl
z#N2J+p_FUUJg{)Pxf;LM%m?6J0(i0<zxfQ~XPT*O@^$|BgA?b-2!9vKWNCrFl5y=<
zCEDK`ZB0bmyj88SSYI1{e@{f?Ra^RRykSMOE7~3HNyMuf0o|bad%fMsc*46U>Wy|L
zI--5v>J{FVt5>bz)f?;cuDo&88XtH&TH@YVPuByuM7^<OqBoiFw)E{yQlphKYP0DT
zqthf~-7Ds;o|v~4Kf`<5VWZdE7R8$a(KhW9tHCp8yQ?E#y<MTs>-wUBk1PB|T~F!S
z8dZ3!M%7af?!V$IoBwU!Q(CUPrdi5QHEi(_HOqrq&Qt@Y;kd^6!3F?0I^bQOY_<dS
zCEO5jsbBWnY<3j&1>8qxP#>SdEjIGuRm}e?)JgvAfqE4Vq#dYP_!~TnW8Zem(eH6w
zbJ?8o5ztlu$KU49zmm<$I~BE;d+Hu3tC+oGaL(hi)2>g>dcx^^Y(BUt4cqvqU)$l!
z8Y~=^;cp#uI_Ej__!1z8-^gaiFzIr99|7NRS&omkpT_^2Z)USE(O09-`BMj(NdFf8
zpLr{rU1`_jtP&-p|0{Vg)lzx#59j3H27Ni!lyBtGzm}t43;c1!jIqes-<xxh*mujj
z5qLKO??&L=2)rAC{}&NZ$CGT<ydP&eZSV#*HMdrlU#JRi@u|Y|35$U*DKH0fZ)P!c
zzPo3K7Z$)v3t+x0$a0yga_|4*p|zAlOI7{mAG0yGc(!N3{jQWM+`Y2kwqD92o-t7H
z3{Q$VKjBfQltsGc(T63o$}6=NQ*VR;W6oQ6gd*iCUE>9HDFz?ZhV}8vI->dTkWPxR
zhvSL(KdK$a;|qa}T^#iUp3-!TlfcGL9M=RM(R3V_^61*!|5;_;t2FP&7aIPw4#*dE
z{j9E^*Y(eI{VQGnPS>+_;I7p5^}4=U*PC>Gm#*7&eV?w4Uz_XdZu3?(?@9I~lB*hU
zJJT2Q-sD^BTeEUqQjj;@pkT}L#UJQSwCq8h=u>rvY4yYsQQzL4q)*<oUD?@Y351)P
z4$If}Ko2BUo#<1P`=WhuzLzX)8i;5q1A?_L+SNjN-RkX1SiVlYSZw+5+AyxTtLX$S
z)`pjoEnl=F(%#q79gTFfLDoQ)uQd_ti^IOEKHLf-bhdPNw!(f4YH;z1Zlz~W9New3
zZn^g=9I*c-d=2C{7fc?X(eu;Pyx*`FJzaAs-Xr2RC}!l%{5Ex!<}~v_r%6tJ9a@~*
zM&8VSQ=9KKXqP-j-ptDy;HcZkoAt!hX8iz*6x}@oS^=~;SB<=xpQe_367U%R4WFr-
zLE~CrV6#4$n*C!b^j~2us|^jV9Y)@)AEv%nGg6*DnEq?D6-UG9CeEy9raqz>jsHg8
z#P4GT@@9QBwTU0)IZlQ0j{?V-n)Ta!|6yul2P{(PtG#>zfL(q-8!&ZlW70xK-n@&P
zE|k~ni1|)w(PFzb_1?N2!-xBJ`)0o~HGYlGE5^P-eg+N7oAEokR8dU5l|)-HejC)^
zqG6XW^&OGnFSP$mf&A9xiqq8msn=f2_!->4pg|UsALctN^PN%Q_!Y{J7RVo3=>^8J
z4y}|*z7tR_buE-X4Yb|<5iNg2%NOz&!rzBHbI<s1_A~Pxn0@`A2S(nE!%u;;+c)18
zjIJU#-p5x}?rP*s{Sqj4d2=4gX!((xq-q#BQ~$C+-t2$<)%pBJjgd3(1RBJd{txK)
z><6^`e>Z-wL6-h&cC=J6^6xhYU6Yr^)O`2f*cO7md`T;O%Z^0b$Qd+YGj|M*b!MG{
zw(5f#pEKQtf75@{FHZMj_Se;bfU=-~SZH5y6~8*+t=R|d9Q?}0>~Hyar)Bo7eEclS
z>__={mu2>ueEe+7>>v5~IhNTM^6}-lbw3}E=;-+6<7a-SN4gU$d~V&f2{>`OH|wbl
z4gLzttdsfp%Pq71<>TFW5NOA!wiBy_S+DF=C%!~8>rp;_0oE6<9VKn2^&ZQtA9kV>
zcLQc!$j2|V%zV$sUzxkF$j4u0nfaNIUu2nimygGn1_eds`^<Q99?dH9-5vAZQ6vtW
zmkxZ9R~QFQ8%L2i_<=9OI3CvfBS0zl;OelZ{WtUD!y11;uP0{S^#e~mUUP#cTKLUi
zasM9{J(bo-!Fc_Z;62uX>$L)1eNON*$MK(lr`>y(=JSsMEEcyPXnu45yh4k;f|~su
z(fu<0eM9pn|4cE;vuZS*V4&R#JSd~cGjh~d3Vvq(^I5h7<7oB)5yEAS)^lm0!pO6G
zG`T-h|MBY-UY@a|$?x%qU#3$?p5X%@68$stuwC#AtiuIy<{KCm{3*=O2etlJmym8*
z^6Vb<lbZjKj*k#=IaET==YVH_3+LAn(LZyY`3ms%@4}4!XG`FpD}f&e-fM(y^)E``
z-zcGn>vgg3U9Sbcn4dS7z}J+(M@!)E2mXCFzfJSW68v8%fq%9H{yD+V{9gBL3I3mz
z!2hNM-hq9!SbVJ3maZ<_C&5;0Um~7JwzvCQ@ojKj)1BKR4fQ*k<V)d=dC(noHd0$G
zvbQU?r==?*cgK;IWWR+w-`=h$e&h7rymsv^773BgwthmlTIvS7E!o}ufU)D>vQf-z
z+Fsk}j~Iw%Oncgj+=`18^G$ARZ?6^E+;HcH+J;D^K9Yytxg*iAE)waC$vuB$=h{t?
z-j=>LRq$c)>T3F$OhhAltBJJuL3C$M-JP2@?eI56nrb&R`0>T`&NYT>=PmKB-r9`~
zcZw!ZZrofg3aBE|7LP?bT6)^L@L6?4VrTk+2LtK@Y>3>tso~B@Z&wR`R#nA{M0+3}
zm1hW%jiD{IjrDa_L|m6IuJI|g)LYbt)HOBsZ;%T=?=JE^cZm<IwI2QbGQVwB--fF%
ztMjek1L)FUYxAq@d}aL9J6}QmXXp7ed=s9J;0NOQ81v<M=`Yhue2dO+(DPN{gY<mF
Fe*wrhE(!nu

diff --git a/src/build_A2.m b/src/build_A2.m
index f6765a9..9caa83a 100644
--- a/src/build_A2.m
+++ b/src/build_A2.m
@@ -4,11 +4,11 @@ function A = build_A2(coordinates,elements)
     A1 = zeros(N);
 
     % untere schranke s t obere schranke k l
-    intF = @(f,s1,s2,k1,k2,t1,t2,l1,l2) ...
-        f(k1,k2,l1,l2)-f(k1,k2,l1,t2)-f(k1,k2,t1,l2)+f(k1,k2,t1,t2)...
-        -f(k1,s2,l1,l2)+f(k1,s2,l1,t2)+f(k1,s2,t1,l2)-f(k1,s2,t1,t2)...
-        -f(s1,k2,l1,l2)+f(s1,k2,l1,t2)+f(s1,k2,t1,l2)-f(s1,l2,t1,t2)...
-        +f(s1,s2,l1,l2)-f(s1,0,l1,t2)-f(s1,s2,t1,l2)+f(s1,s2,t1,t2);
+    intF = @(f,a,b,c,d,r,t,u,v)...
+        f(b,d,t,v)-f(b,d,t,u)-f(b,d,r,v)+f(b,d,r,u)...
+        -f(b,c,t,v)+f(b,c,t,u)+f(b,c,r,v)-f(b,c,r,u)...
+        -f(a,d,t,v)+f(a,d,t,u)+f(a,d,r,v)-f(a,d,r,u)...
+        +f(a,c,t,v)-f(a,c,t,u)-f(a,c,r,v)+f(a,c,r,u);
     tic
     for j = 1:N
         for k = 1:N
@@ -17,7 +17,7 @@ function A = build_A2(coordinates,elements)
             d = ek(1,:) - ej(1,:);
 %             d = ones(1,3);
             A1(j,k) = intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3))...
-                ,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2));
+                ,ej(1,1),ej(2,1),ej(1,2), ej(3,2),ek(1,1), ek(2,1),ek(1,2), ek(3,2));
         end
     end
     
diff --git a/src/mex_Fpar.cpp b/src/mex_Fpar.cpp
new file mode 100644
index 0000000..9af5fac
--- /dev/null
+++ b/src/mex_Fpar.cpp
@@ -0,0 +1,44 @@
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include "mex.h"
+#include <stdlib.h>
+
+#include "SLPrecangle.hpp"
+
+#define my_PI 3.141592653589793
+#define EPS 0.02
+
+using namespace std;
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
+	//sicherheitsabfragen zu Datengroessen
+	if (nrhs != 7)
+		mexErrMsgTxt("expected (1,x2,y1,y2,d1,d2,d3)");
+	if (nlhs > 2)
+		mexErrMsgTxt("has maximal two output arguments");
+
+
+	plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
+	plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
+
+//	plhs[0] = mxCreateDoubleMatrix(1, 2, mxREAL);
+
+
+	//printf("%f,%f,%f,%f",*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
+
+	//sol = G00(p,y1,y2,x1,x2,l);
+
+	*mxGetPr(plhs[0]) = F_par(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+	*mxGetPr(plhs[1]) = FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+
+//	printf("%f",FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6])));
+
+//	*(mxGetPr(plhs[0])+1) = FLO_plane(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+
+
+}
+
+
+
diff --git a/src/mex_G00.cpp b/src/mex_G00.cpp
new file mode 100644
index 0000000..0bc77f3
--- /dev/null
+++ b/src/mex_G00.cpp
@@ -0,0 +1,38 @@
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include "mex.h"
+#include <stdlib.h>
+
+#include "SLPrecangle.hpp"
+
+#define my_PI 3.141592653589793
+#define EPS 0.02
+
+using namespace std;
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
+	//sicherheitsabfragen zu Datengroessen
+	if (nrhs != 6)
+		mexErrMsgTxt("expected (p, y1, y2, x1, x2, l)");
+	if (nlhs > 2)
+		mexErrMsgTxt("has maximal two output arguments");
+
+
+	plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
+	plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL);
+
+
+	//printf("%f,%f,%f,%f",*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
+
+	//sol = G00(p,y1,y2,x1,x2,l);
+
+	*mxGetPr(plhs[0]) = G00(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
+	*mxGetPr(plhs[1]) = slpADLO(*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
+
+
+}
+
+
+
diff --git a/src/mex_g0.cpp b/src/mex_g0.cpp
index b31afe4..a93de87 100644
--- a/src/mex_g0.cpp
+++ b/src/mex_g0.cpp
@@ -17,7 +17,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 	if (nrhs != 4)
 		mexErrMsgTxt("expected (p, y, x, l)");
 	if (nlhs > 2)
-		mexErrMsgTxt("has maximal two output argument");
+		mexErrMsgTxt("has maximal two output arguments");
 
 
 	plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
diff --git a/src/refineQuad.m b/src/refineQuad.m
index 8d31fc9..bef2205 100644
--- a/src/refineQuad.m
+++ b/src/refineQuad.m
@@ -1,5 +1,9 @@
 function [coordinates,elements,fa2so] = refineQuad(coordinates,elements,type)
 
+if([1 1] == size(type))
+    type = repmat(type, size(elements,1),1);
+end
+
 c_loop = size(elements,1);
 
 fa2so = repmat([1:c_loop]',1,4);
diff --git a/src/surfDoubleQuad.m b/src/surfDoubleQuad.m
new file mode 100644
index 0000000..49a9e3e
--- /dev/null
+++ b/src/surfDoubleQuad.m
@@ -0,0 +1,45 @@
+function in = surfDoubleQuad(f,a,b,c,d,r,t,u,v,w)
+%Vierfachintegral ueber f(x1,x2,y1,y2) auf [a b]x[c d]x[r t]x[u v] mit 
+%der Simpsonregel und w Auswertungspunkten
+    s = (b-a)/w;
+    in = 0;
+    for i = a:s:b-s
+        in = in + fx2(f,i,c,d,r,t,u,v,w) + ...
+                4*fx2(f,(2*i+s)/2,c,d,r,t,u,v,w) + ...
+                  fx2(f,i+s,c,d,r,t,u,v,w);
+    end
+    in = in *s/6;
+end
+
+function in = fx2(f,x1,c,d,r,t,u,v,w)
+    s = (d-c)/w;
+    in = 0;
+    for i = c:s:d-s
+        in = in + fy1(f,x1,i,r,t,u,v,w) + ...
+                4*fy1(f,x1,(2*i+s)/2,r,t,u,v,w) + ...
+                  fy1(f,x1,i+s,r,t,u,v,w);
+    end
+    in = in *s/6;
+end
+
+function in = fy1(f,x1,x2,r,t,u,v,w)
+    s = (t-r)/w;
+    in = 0;
+    for i = r:s:t-s
+        in = in + fy2(f,x1,x2,i,u,v,w) + ...
+                4*fy2(f,x1,x2,(2*i+s)/2,u,v,w) + ...
+                  fy2(f,x1,x2,i+s,u,v,w);
+    end
+    in = in *s/6;
+end
+
+function in = fy2(f,x1,x2,y1,u,v,w)
+    s = (v-u)/w;
+    in = 0;
+    for i = u:s:v-s
+        in = in + f(x1,x2,y1,i) + ...
+                4*f(x1,x2,y1,(2*i+s)/2) + ...
+                  f(x1,x2,y1,i+s);
+    end
+    in = in *s/6;
+end
\ No newline at end of file
diff --git a/src/surfQuad.m b/src/surfQuad.m
new file mode 100644
index 0000000..b089675
--- /dev/null
+++ b/src/surfQuad.m
@@ -0,0 +1,19 @@
+function in = surfQuad(f,a,b,c,d,v)
+%Doppelintegral ueber f(x,y) auf [a b]x[c d] mit der Simpsonregel und v
+%Auswertungspunkten
+    s = (b-a)/v;
+    in = 0;
+    for i = a:s:b-s
+        in = in + fy(f,i,c,d,v) + 4*fy(f,(2*i+s)/2,c,d,v) + fy(f,i+s,c,d,v);
+    end
+    in = in *s/6;
+end
+
+function in = fy(f,x,c,d,v)
+    s = (d-c)/v;
+    in = 0;
+    for i = c:s:d-s
+        in = in + f(x,i) + 4*f(x,(2*i+s)/2) + f(x,i+s);
+    end
+    in = in *s/6;
+end
\ No newline at end of file
diff --git a/src/test_Fpar.m b/src/test_Fpar.m
new file mode 100644
index 0000000..27c3e6c
--- /dev/null
+++ b/src/test_Fpar.m
@@ -0,0 +1,19 @@
+function [ret] = test_Fpar(a,b,c,d,r,t,u,v,l1,l2,l3)
+
+
+    quad_sol = surfDoubleQuad(@(x1,x2,y1,y2)1/sqrt((x1-y1-l1).^2+(x2-y2-l2).^2+l3.^2)...
+                ,a,b,c,d,r,t,u,v,4);
+
+    intF = @(f,a,b,c,d,r,t,u,v)...
+        f(b,d,t,v)-f(b,d,t,u)-f(b,d,r,v)+f(b,d,r,u)...
+        -f(b,c,t,v)+f(b,c,t,u)+f(b,c,r,v)-f(b,c,r,u)...
+        -f(a,d,t,v)+f(a,d,t,u)+f(a,d,r,v)-f(a,d,r,u)...
+        +f(a,c,t,v)-f(a,c,t,u)-f(a,c,r,v)+f(a,c,r,u);
+            
+    sol = intF(@(x1,x2,y1,y2)mex_Fpar(x1,x2,y1,y2,l1,l2,l3),a,b,c,d,r,t,u,v);
+
+
+
+    ret = [quad_sol sol];
+end
+
diff --git a/src/test_G00.m b/src/test_G00.m
new file mode 100644
index 0000000..c23354f
--- /dev/null
+++ b/src/test_G00.m
@@ -0,0 +1,23 @@
+function [ret] = test_G00(p, x1, x2,l,a,b,c,d)
+
+
+    quad_sol = surfQuad(@(y1,y2)((x1-y1).^2+(x2-y2).^2+l.^2).^p,a,b,c,d,4);
+
+
+    [my(1) mk(1)] = mex_G00(p,a,c,x1,x2,l);
+    [my(2) mk(2)] = mex_G00(p,b,d,x1,x2,l);
+    [my(3) mk(3)] = mex_G00(p,b,c,x1,x2,l);
+    [my(4) mk(4)] = mex_G00(p,a,d,x1,x2,l);
+
+
+    my_sol = my(1)+my(2)-my(3)-my(4);
+    if(p==-0.5)
+        mk_sol = mk(1)+mk(2)-mk(3)-mk(4);
+    else
+        mk_sol = nan;
+    end
+
+
+    ret = [quad_sol my_sol mk_sol];
+end
+
diff --git a/src/test_functions.m b/src/test_functions.m
new file mode 100644
index 0000000..00ba036
--- /dev/null
+++ b/src/test_functions.m
@@ -0,0 +1,60 @@
+function done=test_functions(eps)
+
+    done = 0;
+%% g0 testen  
+p_g0 = [0.5 0 -0.5 -1 -1.5];
+
+    for p = p_g0
+        s = test_g0(p,1.4,2,1,4);
+        if(abs(s(1)-s(2))>eps)
+            p
+            s
+            return;
+        end
+    end
+    
+    % x muss ausserhalb des Intervalls liegen
+    for p = p_g0
+        s = test_g0(p,1.4,0,2,4);
+        if(abs(s(1)-s(2))>eps)
+            p
+            s
+            return;
+        end
+    end
+
+%% G00 testen
+p_G00 = [-0.5 -1.5];
+
+    for p = p_G00
+        s = test_G00(p,1.4,0.7,1,2,4,1,3);
+        if(abs(s(1)-s(2))>eps)
+            p
+            s
+            return;
+        end
+    end
+    
+    for p = p_G00
+        s = test_G00(p,1.1,0.3,0,2,4.5,1,2.5);
+        if(abs(s(1)-s(2))>eps)
+            p
+            s
+            return;
+        end
+    end
+    
+%% F_par testen
+    
+        
+        s = test_Fpar(1,3,2,3,4,5,1,3,1,3,1);
+        if(abs(s(1)-s(2))>eps)
+            p
+            s
+            return;
+        end
+        
+        
+        
+    done = 1;
+end
\ No newline at end of file
diff --git a/src/test_g0.m b/src/test_g0.m
index b08aaf0..1bd8d15 100644
--- a/src/test_g0.m
+++ b/src/test_g0.m
@@ -3,6 +3,10 @@ function [ret] = test_g0(p, x,l,a,b)
 
 quad_sol = quad(@(y)((x-y).^2+l.^2).^p,a,b);
 
+X = 0:0.2:5;
+Y = ((x-X).^2+l.^2).^p;
+
+plot(X,Y)
 
 [my(1) mk(1)] = mex_g0(p,a,x,l);
 [my(2) mk(2)] = mex_g0(p,b,x,l);
diff --git a/src/test_solveError.m b/src/test_solveError.m
index 7424977..9a1c633 100644
--- a/src/test_solveError.m
+++ b/src/test_solveError.m
@@ -2,29 +2,36 @@
 
 load exmpl_2DLShape
 
+
+[coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
+[coordinates,elements]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
+
+
 A = build_A(coordinates,elements);
+sum(sum (A-A'))
 b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2));
 x = A\b
 xe = x'*A*x;
 
 [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,ones(1,size(elements,1)));
 A_fine = build_A(coordinates_fine,elements_fine);
+sum(sum(A_fine-A_fine'))
 b = sqrt(sum(quad_norm(coordinates_fine,elements_fine,'w').^2,2));
 x_fine = A_fine\b;
 xe_fine = x_fine'*A_fine*x_fine;
 
 x_interpol = zeros(size(elements_fine,1),1);
 for i = 1:size(elements,1)
-    x_interpol(f2s(i,:)) = x(i)/4;
+    x_interpol(f2s(i,:)) = x(i);
 end
 xe_interpol = x_interpol'*A_fine*x_interpol;
 
 [x_fine x_interpol x_fine-x_interpol]
 
-l = x_fine(f2s(1,:));
-T4 = [1 1 1 1; 1 -1 1 -1; 1 1 -1 -1;1 -1 -1 1];
-l
-T4*l*1/4
+% l = x_fine(f2s(1,:));
+% T4 = [1 1 1 1; 1 -1 1 -1; 1 1 -1 -1;1 -1 -1 1];
+% l
+% T4*l*1/4
 
 [xe xe_fine xe_interpol]
 
diff --git a/src/test_solveS.m b/src/test_solveS.m
index 1de19e0..215c0c8 100644
--- a/src/test_solveS.m
+++ b/src/test_solveS.m
@@ -14,24 +14,7 @@ intF = @(f,s1,s2,k1,k2,t1,t2,l1,l2) ...
     -f(s1,k2,l1,l2)+f(s1,k2,l1,t2)+f(s1,k2,t1,l2)-f(s1,l2,t1,t2)...
     +f(s1,s2,l1,l2)-f(s1,0,l1,t2)-f(s1,s2,t1,l2)+f(s1,s2,t1,t2);
 tic
-for j = 1:N
-    for k = 1:N
-        ej = coordinates(elements(j,:)',:);
-        ek = coordinates(elements(k,:)',:);
-        d = ek(1,:) - ej(1,:);
-%         fprintf('%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n'...
-%             ,j,k,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2)...
-%             ,d(1),d(2),d(3));
-        tmp = intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3))...
-            ,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2));
-
-%        fprintf('%d%d|%.2f\n',j,k,tmp);
-        A1(j,k) = 1/ (4*pi) *tmp;
-%           intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3))...
-%             ,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2));
-
-    end
-end
+A2 = build_A2(coordinates,elements);
 t(1) = toc;
 disp ' '
 tic
-- 
2.47.3