From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Wed, 11 May 2011 14:22:07 +0000 (+0000)
Subject: [src] tests für g0 G00 und F_par hinzugefügt
X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=72796aa0c7e30bbdbcbf1f8d37268bc979780ace;p=bacc.git

[src] tests für g0 G00 und F_par hinzugefügt
[src] kleine Fehler behoben...

git-svn-id: https://drops.fb12.tu-berlin.de/svn/bacc/trunk@21 26120e32-c555-405d-b3e1-1f783fb42516
---

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 89a6868..ee6b1c6 100755
Binary files a/src/build_A.mexa64 and b/src/build_A.mexa64 differ
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