From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Sun, 11 Sep 2011 12:53:36 +0000 (+0000)
Subject: [doc] erster Commit
X-Git-Url: https://git.leopard-lacewing.eu/?a=commitdiff_plain;h=0cbae27ad5fc4e802b0dbd1fb74051c44a180f87;p=bacc.git

[doc] erster Commit
[doc] Kapitel und Sektionen Struktur angedeutet
[doc] einige Integrale geschrieben

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

diff --git a/.cproject b/.cproject
index 9eadfa7..c2da3d3 100644
--- a/.cproject
+++ b/.cproject
@@ -45,7 +45,7 @@
 						</toolChain>
 					</folderInfo>
 					<sourceEntries>
-						<entry excluding="src/bem3d" flags="VALUE_WORKSPACE_PATH|RESOLVED" kind="sourcePath" name=""/>
+						<entry excluding="src/maik|src/bem3d" flags="VALUE_WORKSPACE_PATH|RESOLVED" kind="sourcePath" name=""/>
 					</sourceEntries>
 				</configuration>
 			</storageModule>
@@ -508,7 +508,7 @@
 						</toolChain>
 					</folderInfo>
 					<sourceEntries>
-						<entry excluding="src/bem3d" flags="VALUE_WORKSPACE_PATH|RESOLVED" kind="sourcePath" name=""/>
+						<entry excluding="src/maik|src/bem3d" flags="VALUE_WORKSPACE_PATH|RESOLVED" kind="sourcePath" name=""/>
 					</sourceEntries>
 				</configuration>
 			</storageModule>
diff --git a/doc/doc1.dvi b/doc/doc1.dvi
new file mode 100644
index 0000000..48044d5
Binary files /dev/null and b/doc/doc1.dvi differ
diff --git a/doc/doc1.pdf b/doc/doc1.pdf
new file mode 100644
index 0000000..8d6f89e
Binary files /dev/null and b/doc/doc1.pdf differ
diff --git a/doc/doc1.tex b/doc/doc1.tex
new file mode 100644
index 0000000..d8baa51
--- /dev/null
+++ b/doc/doc1.tex
@@ -0,0 +1,85 @@
+\documentclass[a4paper,10pt,fleqn]{report}
+
+\usepackage{fullpage}	%Seiten etwas Größer
+\usepackage{amsmath}	%Mathematische Symbole
+%\usepackage{moreverb}
+\usepackage{graphicx,psfrag,subfig}	%Grafiken einbinden/Texte ersetzen/Bilder nebeneinander
+%\usepackage{ifthen}
+%\usepackage{showkeys}
+%\usepackage[final,numbered]{mcode}
+\usepackage{color}	%Farben benutzen und Definieren
+\usepackage{hyperref}	%Links im Inhaltsverzeichnis
+\hypersetup{linkbordercolor={1 1 1},citebordercolor={1 1 1},urlbordercolor={1 1 1}}
+
+
+\usepackage[ngerman]{babel}	%Sprachpacket für Überschriften
+\usepackage[utf8]{inputenc}	%Eingabekodierung
+\usepackage{fixltx2e}	%Deutschsprach Bugs
+
+
+\def\todo#1{\textcolor{red}{#1}}
+\def\Matlab{{\sc Matlab}}
+
+
+\date{29.08.2011}
+\author{P. Schaefer}
+
+\begin{document}
+\tableofcontents
+\clearpage
+\chapter{Einleitung}
+\section{Problem}
+Zwei Elemente $E = v + [(0,l_1) \times (0,l_2) \times \{0\}]$ und $\hat E = \hat v + [(0,\hat l_1) \times (0,\hat l_2) \times \{0\}]$
+
+\noindent
+Berechnet werden soll:
+\begin{eqnarray*}
+ -\frac{1}{4\pi} \int_E \int_{\hat E} \frac{1}{|x-y|} ds_y ds_x
+\end{eqnarray*}
+
+\section{Gitter}
+Alle Gitter bestehen nur aus achsenorientierten Rechtecken.
+Entscheidend für die Tests sind Gitter in einer Ebene und Ganze Körper...
+
+\chapter{Analytisch}
+
+\section{einfach Inetgral}
+\begin{eqnarray*}
+ g(p;y;x;\lambda) &:=& \int \{(x-y)^2 + \lambda^2 \}^p dy
+\end{eqnarray*}
+
+
+\section{doppel Integral}
+\begin{eqnarray*}
+ G(p;y_1;y_2;x_1;x_2;\lambda) &:=& \int \int \{(x_1-y_1)^2 + (x_2-y_2)^2 + \lambda^2\}^p dy_2 dy_1
+\end{eqnarray*}
+
+
+
+\section{Parallele Elemente}
+Für parallele Elemente lässt sich zeigen:
+\begin{eqnarray*}
+F(x_1,x_2,y_1,y_2,\delta_1,\delta_2,\delta_3) &:=& \int \int \int \int \left\{(x_1-y_1-\delta_1)^2+(x_2-y_2-\delta_2)^2+\delta_3^2\right\}^{-1/2}dy_2 dy_1 dx_2 dx_1
+\end{eqnarray*}
+
+\section{Ortogonale Elemente}
+Für ortogonale Elemente lässt sich zeigen:
+\begin{eqnarray*}
+F(x_1,x_2,y_1,y_2,\delta_1,\delta_2,\delta_3) &:=& \int \int \int \int \left\{(x_1-\delta_1)^2+(x_2-y_2-\delta_2)^2+(-y_3-\delta_3)^2\right\}^{-1/2}dy_3 dy_2 dx_2 dx_1
+\end{eqnarray*}
+
+\chapter{Quadratur}
+\section{Was ist das}
+gauss-Quadratur
+\section{Umsetzung für BLA}
+
+
+
+\chapter{Schätzer und Kreterien}
+\section{mögliche Laagen der Elemente}
+so so oder so...
+\section{Fehlerschätzer}
+\section{möglichkeiten für Kreterien}
+
+ 
+\end{document}
diff --git a/src/mex_Fort.cpp b/src/mex_Fort.cpp
index 0f4b141..4741cb0 100644
--- a/src/mex_Fort.cpp
+++ b/src/mex_Fort.cpp
@@ -21,22 +21,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
 
 	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_ort(*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_perp(*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_Fpar.cpp b/src/mex_Fpar.cpp
index 9103e4d..369f683 100644
--- a/src/mex_Fpar.cpp
+++ b/src/mex_Fpar.cpp
@@ -21,22 +21,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
 
 	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]) = F_parY1Y2_2(*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_FparQY2X2.cpp b/src/mex_FparQY2X2.cpp
deleted file mode 100755
index 883158a..0000000
--- a/src/mex_FparQY2X2.cpp
+++ /dev/null
@@ -1,44 +0,0 @@
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.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 (x1,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_parQY2X2(*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_FparY1Y2_2.cpp b/src/mex_FparY1Y2_2.cpp
deleted file mode 100755
index 7c76461..0000000
--- a/src/mex_FparY1Y2_2.cpp
+++ /dev/null
@@ -1,44 +0,0 @@
-#include <iostream>
-#include <cmath>
-#include <cassert>
-#include "mex.h"
-#include <stdlib.h>
-
-#include "slpRectangle.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 (x1,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_parY1Y2_2(*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
index 24379f5..1c946f8 100644
--- a/src/mex_G00.cpp
+++ b/src/mex_G00.cpp
@@ -23,14 +23,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 	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]) = doubleQuad(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
-
+	*mxGetPr(plhs[0]) = G_AY1Y2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
+	*mxGetPr(plhs[1]) = G_QY1Y2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]));
 
 }
 
diff --git a/src/mex_build_A.cpp b/src/mex_build_A.cpp
index 0cdec5f..5cd019f 100644
--- a/src/mex_build_A.cpp
+++ b/src/mex_build_A.cpp
@@ -233,7 +233,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 					//		y1[0], y0[1], y3[1], d[0], d[1], d[2]);
 					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
 					tmp
-							= apply0Int4(F_par, fabs(xa[rxa]), fabs(xb[rxb]),
+							= calcParIntA( fabs(xa[rxa]), fabs(xb[rxb]),
 									fabs(ya[rxa]), fabs(yb[rxb]), d[rxa],
 									d[rxb], d[rx]);
 					//	 printf("%d%d|%.2f\n",j,k,tmp);
@@ -242,7 +242,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 					//		y1[0], y0[1], y3[1], d[0], d[1], d[2]);
 					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
 					tmp
-							= apply0Int4(F_par, fabs(xa[rxa]), fabs(xb[rxb]),
+							= calcParIntA( fabs(xa[rxa]), fabs(xb[rxb]),
 									fabs(yb[rxa]), fabs(ya[rxb]), d[rxa],
 									d[rxb], d[rx]);
 					//	 printf("%d%d|%.2f\n",j,k,tmp);
@@ -255,26 +255,25 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
 				if (rxa == rya) {
 					tmp
-							= apply0Int4(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),
+							= calcOrtIntA( fabs(xb[rxb]), fabs(xa[rxa]),
 									fabs(ya[rya]), fabs(yb[ryb]), d[rxb],
 									d[rxa], d[rx]);
 				} else if (rxa == ryb) {
 					tmp
-							= apply0Int4(F_ort, fabs(xb[rxb]), fabs(xa[rxa]),
+							= calcOrtIntA( fabs(xb[rxb]), fabs(xa[rxa]),
 									fabs(yb[ryb]), fabs(ya[rya]), d[rxb],
 									d[rxa], d[rx]);
 				} else if (rxb == rya) {
 					tmp
-							= apply0Int4(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),
+							= calcOrtIntA( fabs(xa[rxa]), fabs(xb[rxb]),
 									fabs(ya[rya]), fabs(yb[ryb]), d[rxa],
 									d[rxb], d[rx]);
 				} else {
 					tmp
-							= apply0Int4(F_ort, fabs(xa[rxa]), fabs(xb[rxb]),
+							= calcOrtIntA( fabs(xa[rxa]), fabs(xb[rxb]),
 									fabs(yb[ryb]), fabs(ya[rya]), d[rxa],
 									d[rxb], d[rx]);
 				}
-				A[(k * em) + j] = tmp;
 
 			}
 			A[(k * em) + j] = 1. / (4 * M_PI) * tmp;
diff --git a/src/mex_build_As1.cpp b/src/mex_build_As1.cpp
index 2356c72..cc1b7e5 100644
--- a/src/mex_build_As1.cpp
+++ b/src/mex_build_As1.cpp
@@ -87,7 +87,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
 	//LageInformationen
 	int rx, rxa, rxb, ry, rya, ryb;
-	int rtxa, rtxb, rtx;
 
 	//Ausrechnen
 	for (j = 0; j < em; ++j) {
@@ -121,17 +120,17 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 		// Lageeigenschaften des Flaechenstuecks
 		//if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])
 
-		rtx = dimOfVec(xn);
+		rx = dimOfVec(xn);
 
-		rtxa = dimOfVec(xa);
+		rxa = dimOfVec(xa);
 
-		rtxb = dimOfVec(xb);
+		rxb = dimOfVec(xb);
 
 
 		//kleinste Ecke finden und fuer \delta verwenden
 
-		if (xa[rtxa] > 0) {
-			if (xb[rtxb] > 0) {
+		if (xa[rxa] > 0) {
+			if (xb[rxb] > 0) {
 				for (i = 0; i < 3; ++i)
 					dt[i] = -x0[i];
 			} else {
@@ -139,7 +138,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 					dt[i] = -x3[i];
 			}
 		} else {
-			if (xb[rtxb] > 0) {
+			if (xb[rxb] > 0) {
 				for (i = 0; i < 3; ++i)
 					dt[i] = -x1[i];
 			} else {
@@ -149,20 +148,11 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 			}
 		}
 
-		//for (i=0;i<3;++i)
-		//	dt[i] = 0;
-
-
-		// printf("(%d n- %f | %f | %f)\n",j,xn[0],xn[1],xn[2]);
+		xda = xa[rxa];
+		xdb = xb[rxb];
 
 		for (k = 0; k < em; ++k) {
 
-			xda = xa[rtxa];
-			xdb = xb[rtxb];
-
-			rx = rtx;
-			rxa = rtxa;
-			rxb = rtxb;
 
 
 			y0[0] = C[(int) E[k] - 1];
@@ -221,34 +211,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 			yda = ya[rya];
 			ydb = yb[ryb];
 
-
-				if(xda*xda+xdb*xdb>yda*yda+ydb*ydb){
-
-				/*	printf("+(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
-					printf("[%.2f,%.2f@%.2f,%.2f]",xda,xb[rxb],ya[rya],yb[ryb]);
-					printf("{%d%d} %.1f %.1f %.1f",j,k,d[0],d[1],d[2]);
-					printf("\n");*/
-
-
-					dtmp = xda; xda = yda; yda = dtmp;
-					dtmp = xdb; xdb = ydb; ydb = dtmp;
-
-					itmp = rxa; rxa = rya; rya = itmp;
-					itmp = rxb; rxb = ryb; ryb = itmp;
-					itmp = rx; rx = ry; ry = itmp;
-
-					for (i = 0; i<3;++i)
-						d[i] = -d[i];
-
-
-
-			/*		printf("-(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
-					printf("[%.2f,%.2f@%.2f,%.2f]",xda,xb[rxb],ya[rya],yb[ryb]);
-					printf("{%d%d} %.1f %.1f %.1f",j,k,d[0],d[1],d[2]);
-					printf("\n");*/
-				}
-
-
 			//printf("(%d,%d)",rx,ry);
 			if (rx == ry) {
 				//printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
@@ -260,29 +222,17 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 					//		y1[0], y0[1], y3[1], d[0], d[1], d[2]);
 					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
 
-			//		if(xda*xda+xdb*xdb<=yda*yda+ydb*ydb){
-
 					tmp
-							= apply0Int4(F_par, fabs(xda), fabs(xdb),
+							= calcParIntO1( fabs(xda), fabs(xdb),
 									fabs(yda), fabs(ydb), d[rxa],
 									d[rxb], d[rx]);
 
-			//		}else{
-
-			//		tmp
-			//				= apply0Int4(F_par, fabs(yda), fabs(ydb),
-			//						fabs(xda), fabs(xdb), -d[rxa],
-			//						-d[rxb], -d[rx]);
-			//		}
-			//		if(fabs(tmp-tmp2)>10e-16)
-			//			printf("wtf");
-					//	 printf("%d%d|%.2f\n",j,k,tmp);
 				} else {
 					//printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],
 					//		y1[0], y0[1], y3[1], d[0], d[1], d[2]);
 					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
 					tmp
-							= apply0Int4(F_par, fabs(xda), fabs(xdb),
+							= calcParIntO1( fabs(xda), fabs(xdb),
 									fabs(ydb), fabs(yda), d[rxa],
 									d[rxb], d[rx]);
 					//	 printf("%d%d|%.2f\n",j,k,tmp);
@@ -295,22 +245,22 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
 				if (rxa == rya) {
 					tmp
-							= apply0Int4(F_ort, fabs(xdb), fabs(xda),
+							= calcOrtIntO1( fabs(xdb), fabs(xda),
 									fabs(yda), fabs(ydb), d[rxb],
 									d[rxa], d[rx]);
 				} else if (rxa == ryb) {
 					tmp
-							= apply0Int4(F_ort, fabs(xdb), fabs(xda),
+							= calcOrtIntO1( fabs(xdb), fabs(xda),
 									fabs(ydb), fabs(yda), d[rxb],
 									d[rxa], d[rx]);
 				} else if (rxb == rya) {
 					tmp
-							= apply0Int4(F_ort, fabs(xda), fabs(xdb),
+							= calcOrtIntO1( fabs(xda), fabs(xdb),
 									fabs(yda), fabs(ydb), d[rxa],
 									d[rxb], d[rx]);
 				} else {
 					tmp
-							= apply0Int4(F_ort, fabs(xda), fabs(xdb),
+							= calcOrtIntO1( fabs(xda), fabs(xdb),
 									fabs(ydb), fabs(yda), d[rxa],
 									d[rxb], d[rx]);
 				}
diff --git a/src/mex_build_As2.cpp b/src/mex_build_As2.cpp
index f254131..26b1cb4 100644
--- a/src/mex_build_As2.cpp
+++ b/src/mex_build_As2.cpp
@@ -1,7 +1,7 @@
 /*
  * voll Analytisch
  * Element mit groesserer Diagonale nach aussen gedreht
- *
+ * Quadratur ueber E_j bei dist(E_j,E_k) >= eta*dia(E_j) | dia(E_j) <= dia(E_k)
  *
  */
 
@@ -17,7 +17,8 @@
 
 #include "slpRectangle.hpp"
 
-#define EPS 0.02
+//Quadraturverwendung!!!
+#define eta 0.8
 
 using namespace std;
 
@@ -87,7 +88,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
 	//LageInformationen
 	int rx, rxa, rxb, ry, rya, ryb;
-	int rtxa, rtxb, rtx;
 
 	//Ausrechnen
 	for (j = 0; j < em; ++j) {
@@ -121,17 +121,17 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 		// Lageeigenschaften des Flaechenstuecks
 		//if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])
 
-		rtx = dimOfVec(xn);
+		rx = dimOfVec(xn);
 
-		rtxa = dimOfVec(xa);
+		rxa = dimOfVec(xa);
 
-		rtxb = dimOfVec(xb);
+		rxb = dimOfVec(xb);
 
 
 		//kleinste Ecke finden und fuer \delta verwenden
 
-		if (xa[rtxa] > 0) {
-			if (xb[rtxb] > 0) {
+		if (xa[rxa] > 0) {
+			if (xb[rxb] > 0) {
 				for (i = 0; i < 3; ++i)
 					dt[i] = -x0[i];
 			} else {
@@ -139,7 +139,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 					dt[i] = -x3[i];
 			}
 		} else {
-			if (xb[rtxb] > 0) {
+			if (xb[rxb] > 0) {
 				for (i = 0; i < 3; ++i)
 					dt[i] = -x1[i];
 			} else {
@@ -149,20 +149,11 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 			}
 		}
 
-		//for (i=0;i<3;++i)
-		//	dt[i] = 0;
-
-
-		// printf("(%d n- %f | %f | %f)\n",j,xn[0],xn[1],xn[2]);
+		xda = xa[rxa];
+		xdb = xb[rxb];
 
 		for (k = 0; k < em; ++k) {
 
-			xda = xa[rtxa];
-			xdb = xb[rtxb];
-
-			rx = rtx;
-			rxa = rtxa;
-			rxb = rtxb;
 
 
 			y0[0] = C[(int) E[k] - 1];
@@ -221,25 +212,6 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 			yda = ya[rya];
 			ydb = yb[ryb];
 
-				//groesseres Element nach aussen drehen
-				if(xda*xda+xdb*xdb>yda*yda+ydb*ydb){
-
-					dtmp = xda; xda = yda; yda = dtmp;
-					dtmp = xdb; xdb = ydb; ydb = dtmp;
-
-					itmp = rxa; rxa = rya; rya = itmp;
-					itmp = rxb; rxb = ryb; ryb = itmp;
-					itmp = rx; rx = ry; ry = itmp;
-
-					for (i = 0; i<3;++i)
-						d[i] = -d[i];
-
-				}
-
-				//if(xda*xda+xdb*xdb>=d[0]*d[0]+d[1]*d[1]+d[2]*d[2]){
-				//	printf(".");
-				//}
-
 			//printf("(%d,%d)",rx,ry);
 			if (rx == ry) {
 				//printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
@@ -251,31 +223,19 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 					//		y1[0], y0[1], y3[1], d[0], d[1], d[2]);
 					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
 
-			//		if(xda*xda+xdb*xdb<d[0]*d[0]+d[1]*d[1]+d[2]*d[2]){
-
 					tmp
-							= apply0Int4(F_par, fabs(xda), fabs(xdb),
+							= calcParIntO2( fabs(xda), fabs(xdb),
 									fabs(yda), fabs(ydb), d[rxa],
-									d[rxb], d[rx]);
+									d[rxb], d[rx],eta);
 
-			//		}else{
-
-					tmp2
-							= apply0Int4(F_parY1Y2_2, fabs(xda), fabs(xdb),
-									fabs(yda), fabs(ydb), d[rxa],
-									d[rxb], d[rx]);
-			//		}
-					if(fabs(tmp-tmp2)>10e-8)
-						return;
-					//	 printf("%d%d|%.2f\n",j,k,tmp);
 				} else {
 					//printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],
 					//		y1[0], y0[1], y3[1], d[0], d[1], d[2]);
 					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
 					tmp
-							= apply0Int4(F_par, fabs(xda), fabs(xdb),
+							= calcParIntO2( fabs(xda), fabs(xdb),
 									fabs(ydb), fabs(yda), d[rxa],
-									d[rxb], d[rx]);
+									d[rxb], d[rx],eta);
 					//	 printf("%d%d|%.2f\n",j,k,tmp);
 				}
 
@@ -286,24 +246,24 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 
 				if (rxa == rya) {
 					tmp
-							= apply0Int4(F_ort, fabs(xdb), fabs(xda),
+							= calcOrtIntO2( fabs(xdb), fabs(xda),
 									fabs(yda), fabs(ydb), d[rxb],
-									d[rxa], d[rx]);
+									d[rxa], d[rx],eta);
 				} else if (rxa == ryb) {
 					tmp
-							= apply0Int4(F_ort, fabs(xdb), fabs(xda),
+							= calcOrtIntO2( fabs(xdb), fabs(xda),
 									fabs(ydb), fabs(yda), d[rxb],
-									d[rxa], d[rx]);
+									d[rxa], d[rx],eta);
 				} else if (rxb == rya) {
 					tmp
-							= apply0Int4(F_ort, fabs(xda), fabs(xdb),
+							= calcOrtIntO2( fabs(xda), fabs(xdb),
 									fabs(yda), fabs(ydb), d[rxa],
-									d[rxb], d[rx]);
+									d[rxb], d[rx],eta);
 				} else {
 					tmp
-							= apply0Int4(F_ort, fabs(xda), fabs(xdb),
+							= calcOrtIntO2( fabs(xda), fabs(xdb),
 									fabs(ydb), fabs(yda), d[rxa],
-									d[rxb], d[rx]);
+									d[rxb], d[rx],eta);
 				}
 
 			}
diff --git a/src/mex_build_As3.cpp b/src/mex_build_As3.cpp
new file mode 100644
index 0000000..6a6da3a
--- /dev/null
+++ b/src/mex_build_As3.cpp
@@ -0,0 +1,288 @@
+/*
+ * voll Analytisch
+ * Element mit groesserer Diagonale nach aussen gedreht
+ * Quadratur ueber E_j bei dist(E_j,E_k) >= eta*dia(E_j) | dia(E_j) <= dia(E_k)
+ *
+ */
+
+
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include "mex.h"
+#include <stdlib.h>
+
+//#include "tbb/parallel_for.h"
+
+
+#include "slpRectangle.hpp"
+
+//Quadraturverwendung!!!
+#define eta 0.8
+
+using namespace std;
+
+inline int dimOfVec(double * vec)
+{
+	if (vec[2] != 0)
+		return 2;
+	else if (vec[1] != 0)
+		return 1;
+	else
+		return 0;
+}
+
+
+void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
+
+	int i, j, k;
+
+	//sicherheitsabfragen zu Datengroessen
+	if (nrhs != 2)
+		mexErrMsgTxt("expected (coordinates(Nx3),elements(Mx4))");
+	if (nlhs > 1)
+		mexErrMsgTxt("has only one output argument");
+
+	int cm = mxGetM(prhs[0]);
+	int cn = mxGetN(prhs[0]);
+
+	if (cn != 3)
+		mexErrMsgTxt("expected coordinates (Nx3)");
+	int em = mxGetM(prhs[1]);
+	int en = mxGetN(prhs[1]);
+	if (en != 4)
+		mexErrMsgTxt("expected elements (Mx4)");
+	//Vorbereitung der Daten
+
+	plhs[0] = mxCreateDoubleMatrix(em, em, mxREAL);
+	double * A = mxGetPr(plhs[0]);
+	double * C = mxGetPr(prhs[0]);
+	double * E = mxGetPr(prhs[1]);
+
+	double * x0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * x1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * x2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * x3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * xn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * xa = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * xb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+
+	double * y0 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * y1 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * y2 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * y3 = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * yn = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * ya = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * yb = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+
+	double * d = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+	double * dt = mxGetPr(mxCreateDoubleMatrix(3, 1, mxREAL));
+
+	double tmp,tmp2;
+
+	int itmp;
+	double dtmp;
+	double * vtmp;
+
+	long double xda,xdb,yda,ydb;
+
+	//LageInformationen
+	int rx, rxa, rxb, ry, rya, ryb;
+
+	//Ausrechnen
+	for (j = 0; j < em; ++j) {
+		x0[0] = C[(int) E[j] - 1];
+		x0[1] = C[cm + (int) E[j] - 1];
+		x0[2] = C[2 * cm + (int) E[j] - 1];
+
+		x1[0] = C[(int) E[em + j] - 1];
+		x1[1] = C[cm + (int) E[em + j] - 1];
+		x1[2] = C[2 * cm + (int) E[em + j] - 1];
+
+		x2[0] = C[(int) E[2 * em + j] - 1];
+		x2[1] = C[cm + (int) E[2 * em + j] - 1];
+		x2[2] = C[2 * cm + (int) E[2 * em + j] - 1];
+
+		x3[0] = C[(int) E[3 * em + j] - 1];
+		x3[1] = C[cm + (int) E[3 * em + j] - 1];
+		x3[2] = C[2 * cm + (int) E[3 * em + j] - 1];
+
+		for (i = 0; i < 3; ++i)
+			xa[i] = x1[i] - x0[i];
+
+		for (i = 0; i < 3; ++i)
+			xb[i] = x3[i] - x0[i];
+
+		//Kreuzprodukt
+		xn[0] = (xa[1]) * (xb[2]) - (xa[2]) * (xb[1]);
+		xn[1] = (xa[2]) * (xb[0]) - (xa[0]) * (xb[2]);
+		xn[2] = (xa[0]) * (xb[1]) - (xa[1]) * (xb[0]);
+
+		// Lageeigenschaften des Flaechenstuecks
+		//if(xn[0]*xn[0]+xn[1]*xn[1]<xn[2]*xn[2])
+
+		rx = dimOfVec(xn);
+
+		rxa = dimOfVec(xa);
+
+		rxb = dimOfVec(xb);
+
+
+		//kleinste Ecke finden und fuer \delta verwenden
+
+		if (xa[rxa] > 0) {
+			if (xb[rxb] > 0) {
+				for (i = 0; i < 3; ++i)
+					dt[i] = -x0[i];
+			} else {
+				for (i = 0; i < 3; ++i)
+					dt[i] = -x3[i];
+			}
+		} else {
+			if (xb[rxb] > 0) {
+				for (i = 0; i < 3; ++i)
+					dt[i] = -x1[i];
+			} else {
+				for (i = 0; i < 3; ++i)
+					dt[i] = -x2[i];
+				//printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");
+			}
+		}
+
+		xda = xa[rxa];
+		xdb = xb[rxb];
+
+		for (k = 0; k < em; ++k) {
+
+
+
+			y0[0] = C[(int) E[k] - 1];
+			y0[1] = C[cm + (int) E[k] - 1];
+			y0[2] = C[2 * cm + (int) E[k] - 1];
+
+			y1[0] = C[(int) E[em + k] - 1];
+			y1[1] = C[cm + (int) E[em + k] - 1];
+			y1[2] = C[2 * cm + (int) E[em + k] - 1];
+
+			y2[0] = C[(int) E[2 * em + k] - 1];
+			y2[1] = C[cm + (int) E[2 * em + k] - 1];
+			y2[2] = C[2 * cm + (int) E[2 * em + k] - 1];
+
+			y3[0] = C[(int) E[3 * em + k] - 1];
+			y3[1] = C[cm + (int) E[3 * em + k] - 1];
+			y3[2] = C[2 * cm + (int) E[3 * em + k] - 1];
+
+			for (i = 0; i < 3; ++i)
+				ya[i] = y1[i] - y0[i];
+
+			for (i = 0; i < 3; ++i)
+				yb[i] = y3[i] - y0[i];
+
+			yn[0] = (ya[1]) * (yb[2]) - (ya[2]) * (yb[1]);
+			yn[1] = (ya[2]) * (yb[0]) - (ya[0]) * (yb[2]);
+			yn[2] = (ya[0]) * (yb[1]) - (ya[1]) * (yb[0]);
+
+			ry = dimOfVec(yn);
+
+			rya = dimOfVec(ya);
+
+			ryb = dimOfVec(yb);
+
+			//kleinste Ecke finden und fuer \delta verwenden
+
+			if (ya[rya] > 0) {
+				if (yb[ryb] > 0) {
+					for (i = 0; i < 3; ++i)
+						d[i] = dt[i] + y0[i];
+				} else {
+					for (i = 0; i < 3; ++i)
+						d[i] = dt[i] + y3[i];
+				}
+			} else {
+				if (yb[ryb] > 0) {
+					for (i = 0; i < 3; ++i)
+						d[i] = dt[i] + y1[i];
+				} else {
+					for (i = 0; i < 3; ++i)
+						d[i] = dt[i] + y2[i];
+					//printf("Error! Daten existieren nicht. Zugriff auf vierten punkt");
+				}
+			}
+
+			yda = ya[rya];
+			ydb = yb[ryb];
+
+			//printf("(%d,%d)",rx,ry);
+			if (rx == ry) {
+				//printf("(%d,%d@%d,%d)",rxa,rxb,rya,ryb);
+				//printf("[%.2f,%.2f@%.2f,%.2f]",xda,xdb,yda,ydb);
+				//printf("%d%d %.1f %.1f %.1f|||",j,k,dt[0],dt[1],dt[2]);
+				//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
+				if (rxa == rya) {
+					//printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],
+					//		y1[0], y0[1], y3[1], d[0], d[1], d[2]);
+					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
+
+					tmp
+							= calcParIntO3( fabs(xda), fabs(xdb),
+									fabs(yda), fabs(ydb), d[rxa],
+									d[rxb], d[rx],eta);
+
+				} else {
+					//printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x0[0], x1[0], x0[1], x3[1], y0[0],
+					//		y1[0], y0[1], y3[1], d[0], d[1], d[2]);
+					//printf("%d%d %.1f %.1f %.1f\n",j,k,d[0],d[1],d[2]);
+					tmp
+							= calcParIntO3( fabs(xda), fabs(xdb),
+									fabs(ydb), fabs(yda), d[rxa],
+									d[rxb], d[rx],eta);
+					//	 printf("%d%d|%.2f\n",j,k,tmp);
+				}
+
+			} else {
+				//printf("(%d,%d@%d,%d)", rxa, rxb, rya, ryb);
+				//printf("(%d,%d)", rx, ry);
+				//printf("\n");
+
+				if (rxa == rya) {
+					tmp
+							= calcOrtIntO3( fabs(xdb), fabs(xda),
+									fabs(yda), fabs(ydb), d[rxb],
+									d[rxa], d[rx],eta);
+				} else if (rxa == ryb) {
+					tmp
+							= calcOrtIntO3( fabs(xdb), fabs(xda),
+									fabs(ydb), fabs(yda), d[rxb],
+									d[rxa], d[rx],eta);
+				} else if (rxb == rya) {
+					tmp
+							= calcOrtIntO3( fabs(xda), fabs(xdb),
+									fabs(yda), fabs(ydb), d[rxa],
+									d[rxb], d[rx],eta);
+				} else {
+					tmp
+							= calcOrtIntO3( fabs(xda), fabs(xdb),
+									fabs(ydb), fabs(yda), d[rxa],
+									d[rxb], d[rx],eta);
+				}
+
+			}
+			A[(k * em) + j] = 1. / (4 * M_PI) * tmp;
+
+		}
+
+	}
+	//printf("\n");
+	//Rueckgabe (eventuell zurueck schreiben)
+	//   mxFree(x0);
+	//mxFree(x1);
+	//mxFree(x3);
+	//mxFree(y0);
+	//mxFree(y1);
+	//mxFree(y3);
+	//mxFree(xn);
+	//mxFree(yn);
+	//mxFree(d);
+
+	return;
+}
diff --git a/src/mex_calcOrtInt.cpp b/src/mex_calcOrtInt.cpp
new file mode 100755
index 0000000..f4b35ae
--- /dev/null
+++ b/src/mex_calcOrtInt.cpp
@@ -0,0 +1,31 @@
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include "mex.h"
+#include <stdlib.h>
+
+#include "slpRectangle.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 (x1,x2,y1,y2,d1,d2,d3)");
+	if (nlhs > 2)
+		mexErrMsgTxt("has maximal two output arguments");
+
+	plhs[0] = mxCreateDoubleMatrix(1, 3, mxREAL);
+
+	*mxGetPr(plhs[0]) = calcortIntA(*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) = calcParIntQX1X2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+//	*(mxGetPr(plhs[0])+2) = calcParIntQY1X1(*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_calcParInt.cpp b/src/mex_calcParInt.cpp
new file mode 100755
index 0000000..08b26a9
--- /dev/null
+++ b/src/mex_calcParInt.cpp
@@ -0,0 +1,32 @@
+#include <iostream>
+#include <cmath>
+#include <cassert>
+#include "mex.h"
+#include <stdlib.h>
+
+#include "slpRectangle.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 (x1,x2,y1,y2,d1,d2,d3)");
+	if (nlhs > 2)
+		mexErrMsgTxt("has maximal two output arguments");
+
+
+	plhs[0] = mxCreateDoubleMatrix(1, 3, mxREAL);
+
+	*mxGetPr(plhs[0]) = calcParIntA(*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) = calcParIntQX1X2(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]),*mxGetPr(prhs[4]),*mxGetPr(prhs[5]),*mxGetPr(prhs[6]));
+	*(mxGetPr(plhs[0])+2) = calcParIntQY1X1(*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_g0.cpp b/src/mex_g0.cpp
index 7dfb788..4353069 100644
--- a/src/mex_g0.cpp
+++ b/src/mex_g0.cpp
@@ -23,12 +23,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 	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]));
-
-	*mxGetPr(plhs[0]) = g0(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
-	*mxGetPr(plhs[1]) = *mxGetPr(prhs[1])!=0?singleQuad(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3])):0;
-
+	*mxGetPr(plhs[0]) = g_AY(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
+	*mxGetPr(plhs[0]) = g_QY(*mxGetPr(prhs[0]),*mxGetPr(prhs[1]),*mxGetPr(prhs[2]),*mxGetPr(prhs[3]));
 
 }
 
diff --git a/src/mexit.m b/src/mexit.m
index d103b15..f868508 100644
--- a/src/mexit.m
+++ b/src/mexit.m
@@ -4,8 +4,6 @@ mex mex_G00.cpp slpRectangle.cpp
 mex mex_Fpar.cpp slpRectangle.cpp
 mex mex_Fort.cpp slpRectangle.cpp
 
-mex mex_FparY1Y2_2.cpp slpRectangle.cpp
-
 mex mex_build_A.cpp slpRectangle.cpp
 mex mex_build_As1.cpp slpRectangle.cpp
 mex mex_build_As2.cpp slpRectangle.cpp
diff --git a/src/slpRectangle.cpp b/src/slpRectangle.cpp
index f4a0e45..8fad411 100644
--- a/src/slpRectangle.cpp
+++ b/src/slpRectangle.cpp
@@ -3,7 +3,7 @@
 #include <cassert>
 #include <stdlib.h>
 
-#include <mex.h>
+//#include <mex.h>
 
 #include "slpRectangle.hpp"
 
@@ -14,9 +14,14 @@ using namespace std;
 double gauss_w5[] = {0.1185,0.2393,0.2844,0.2393,0.1185};
 double gauss_n5[] = {0.0469,0.2308,0.5000,0.7692,0.9531};
 
+
+
 int sign(double);
 //double arsinh(double);
 
+
+/* ============================== */
+
 int inline sign(double x) {
 	return x > 0 ? 1 : (x < 0 ? -1 : 0);
 }
@@ -26,7 +31,11 @@ int inline sign(double x) {
  }
  */
 
-double singleQuad(double p, double y, double x, double l) {
+double inline max(double x, double y){
+	return x<y?y:x;
+}
+
+double g_QY(double p, double y, double x, double l) {
 
 	double sol = 0;
 
@@ -37,21 +46,8 @@ double singleQuad(double p, double y, double x, double l) {
 	return sol;
 }
 
-
-double doubleQuad(double p, double y1, double y2, double x1, double x2, double l) {
-
-	double sol = 0;
-	for(int i = 0;i<5;++i){
-		for(int j=0;j<5;++j){
-			sol += pow((x1-y1*gauss_n5[i])*(x1-y1*gauss_n5[i])
-					+(x2-y2*gauss_n5[j])*(x2-y2*gauss_n5[j])+l*l,p)*y1*gauss_w5[i]*y2*gauss_w5[j];
-		}
-	}
-	return sol;
-}
-
 //y-x muss != 0 sein
-double g0(double p, double y, double x, double l) {
+double g_AY(double p, double y, double x, double l) {
 	//printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);
 
 	double sol = 0;
@@ -71,7 +67,7 @@ double g0(double p, double y, double x, double l) {
 			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);
+					* g_AY(p - 1, y, x, l) / (2 * p + 1);
 	} else {
 		if (p == -0.5)
 			sol = sign(y - x) * log(fabs(y - x));
@@ -82,7 +78,34 @@ double g0(double p, double y, double x, double l) {
 	return sol;
 }
 
-double G00(double p, double y1, double y2, double x1, double x2, double l) {
+double G_QY1Y2(double p, double y1, double y2, double x1, double x2, double l) {
+
+	double sol = 0;
+	for(int i = 0;i<5;++i){
+		for(int j=0;j<5;++j){
+			sol += pow((x1-y1*gauss_n5[i])*(x1-y1*gauss_n5[i])
+					+(x2-y2*gauss_n5[j])*(x2-y2*gauss_n5[j])+l*l,p)*y1*gauss_w5[i]*y2*gauss_w5[j];
+		}
+	}
+	return sol;
+}
+
+
+double G_AY2X2(double y1, double y2, double x1, double x2, double d3) {
+
+	//	 printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);
+
+	double hlp = ((x1 - y1 ) * (x1 - y1 ) + (x2 - y2) * (x2 - y2) + d3 * d3);
+
+	double sol = sqrt(hlp);
+
+	if ((x2 - y2) != 0)
+		sol += (x2 - y2) * log(sqrt(hlp)-(x2 - y2));
+
+	return sol;
+}
+
+double G_AY1Y2(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 pt = p + 1.5;
 	double sol = 0;
@@ -99,84 +122,45 @@ double G00(double p, double y1, double y2, double x1, double x2, double l) {
 
 	} else if ((pt > 0) && ((int) pt == pt)) {
 		if (l != 0)
-			sol = 2 * p * l * l * G00(p - 1, y1, y2, x1, x2, l);
+			sol = 2 * p * l * l * G_AY1Y2(p - 1, y1, y2, x1, x2, l);
 		if ((y1 - x1) != 0)
-			sol += (y1 - x1) * g0(p, y2, x2,
+			sol += (y1 - x1) * g_AY(p, y2, x2,
 					sqrt((y1 - x1) * (y1 - x1) + l * l));
 		if ((y2 - x2) != 0)
-			sol += (y2 - x2) * g0(p, y1, x1,
+			sol += (y2 - x2) * g_AY(p, y1, x1,
 					sqrt((y2 - x2) * (y2 - x2) + l * l));
 		sol /= 2 * p + 2;
 	} else {
 		sol = NAN;
-		cout << "warning in G00: no case for p="<< p <<" defined. Possible: [-1.5,-0.5,0.5,...]\n";
-		//mexErrMsgTxt("no case for p defined\n");
+		cout << "warning in G_AY1Y2: no case for p="<< p <<" defined. Possible: [-1.5,-0.5,0.5,...]\n";
 	}
 
 	return sol;
 }
 
-double F_parY1Y2(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 *= G00(-0.5, x1, x2, y1 + d1, y2 + d2, d3);
-
-	if ((x1 - y1 - d1) != 0)
-		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) * 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 F_parY1X2(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 *= G00(-0.5, x1, -y2 - d2, y1 + d1, -x2, d3);
-
-	if ((x1 - y1 - d1) != 0)
-		sol -= (x1 - y1 - d1) * g0(0.5, x1, y1 + d1,
-				sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));
+double Gs_AX2Y1Y2(double p, double y1, double y2, double x1, double x2, double l) {
+	double sol = 0;
 
-	if ((x2 - y2 - d2) != 0)
-		sol -= (x2 - y2 - d2) * g0(0.5, -y2 - d2, -x2,
-				sqrt((x1 - y1 - d1) * (x1 - y1 - d1) + d3 * d3));
+	//sol += y2*G_AY1Y2(y1,x2,)
 
-	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 F_parX1Y2(double x1, double x2, double y1, double y2, double d1, double d2,
+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 *= G00(-0.5, -y1 - d1, x2, -x1, y2 + d2, d3);
+		sol *= G_AY1Y2(-0.5, x1, x2, y1 + d1, y2 + d2, d3);
 
 	if ((x1 - y1 - d1) != 0)
-		sol -= (x1 - y1 - d1) * g0(0.5, -y1 - d1, -x1,
+		sol -= (x1 - y1 - d1) * g_AY(0.5, x1, y1 + d1,
 				sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));
 
 	if ((x2 - y2 - d2) != 0)
-		sol -= (x2 - y2 - d2) * g0(0.5, x2, y2 + d2,
+		sol -= (x2 - y2 - d2) * g_AY(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
@@ -185,71 +169,34 @@ double F_parX1Y2(double x1, double x2, double y1, double y2, double d1, double d
 	return sol;
 }
 
-double F_parY1Y2_2(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 *= doubleQuad(-0.5, x1, x2, y1 + d1, y2 + d2, d3);
-
-	if ((x1 - y1 - d1) != 0)
-		sol -= (x1 - y1 - d1) * singleQuad(0.5, x1, y1 + d1,
-				sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + d3 * d3));
 
-	if ((x2 - y2 - d2) != 0)
-		sol -= (x2 - y2 - d2) * singleQuad(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 F_parQY2X2(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 hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2
-			- d2) + d3 * d3);
-
-	double sol = sqrt(hlp);
-
-	if ((x2 - y2 - d2) != 0)
-		sol += (x2 - y2 - d2) * log(sqrt(hlp)-(x2 - y2 - d2));
-
-	return sol;
-}
 
 double F_ort(double x1, double x2, double y2, double y3, double d1, double d2,
 		double d3) {
 
 	//	 printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y2,y3,d1,d2,d3);
-	double sol = -G00(0.5, y3, x1, -d3, d1, x2 - y2 - d2);
+	double sol = -G_AY1Y2(0.5, y3, x1, -d3, d1, x2 - y2 - d2);
 
 	if ((x1 - d1) * (x2 - y2 - d2) != 0)
-		sol -= (x1 - d1) * (x2 - y2 - d2) * G00(-0.5, x2, y3, y2 + d2, -d3,
+		sol -= (x1 - d1) * (x2 - y2 - d2) * G_AY1Y2(-0.5, x2, y3, y2 + d2, -d3,
 				x1 - d1);
 
 	if ((x1 - d1) != 0)
-		sol += (x1 - d1) * g0(0.5, y3, -d3,
+		sol += (x1 - d1) * g_AY(0.5, y3, -d3,
 				sqrt((x1 - d1) * (x1 - d1) + (x2 - y2 - d2) * (x2 - y2 - d2)));
 
 	if ((y3 + d3) * (x2 - y2 - d2) != 0)
-		sol -= (y3 + d3) * (x2 - y2 - d2) * G00(-0.5, x1, x2, d1, y2 + d2,
+		sol -= (y3 + d3) * (x2 - y2 - d2) * G_AY1Y2(-0.5, x1, x2, d1, y2 + d2,
 				-y3 - d3);
 
 	if ((y3 + d3) != 0)
-		sol += (y3 + d3) * g0(0.5, x1, d1,
+		sol += (y3 + d3) * g_AY(0.5, x1, d1,
 				sqrt((x2 - y2 - d2) * (x2 - y2 - d2) + (y3 + d3) * (y3 + d3)));
 
 	return sol / 2.;
 }
 
-double applyInt4(
+/*double applyInt4(
 		double(*f)(double, double, double, double, double, double, double),
 		double a, double b, double c, double d, double r, double t, double u,
 		double v, double d1, double d2, double d3) {
@@ -261,18 +208,8 @@ double applyInt4(
 			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,
-	 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 apply0Int4(
 		double(*f)(double, double, double, double, double, double, double),
@@ -299,102 +236,187 @@ double apply0Int2(
 
 }
 
-double slpADLO(double y1, double y2, double x1, double x2, double a) {
-	double G3 = 0;
-	double gL = 0;
-	double gK = 0;
-	double tmp;
+// Berechnet das eigentliche Integral für parallele Elemente voll Analytisch
+double calcParIntA(double b, double d, double t, double v, double d1, double d2, double d3) {
+	return apply0Int4(F_par, b, d, t, v, d1, d2, d3);
 
-	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);
+}
+
+double calcParIntQX1X2(double b, double d, double t, double v, double d1,
+		double d2, double d3) {
+	//Fall 2
+	double sol = 0;
+
+	for (int i = 0; i < 5; ++i) {
+		for (int j = 0; j < 5; ++j) {
+			sol += G_AY1Y2(-0.5, t + d1, v + d2, b * gauss_n5[i],
+					d * gauss_n5[j], d3) * b * d * gauss_w5[i] * gauss_w5[j];
+			sol -= G_AY1Y2(-0.5, d1, v + d2, b * gauss_n5[i], d * gauss_n5[j],
+					d3) * b * d * gauss_w5[i] * gauss_w5[j];
+			sol -= G_AY1Y2(-0.5, t + d1, d2, b * gauss_n5[i], d * gauss_n5[j],
+					d3) * b * d * gauss_w5[i] * gauss_w5[j];
+			sol += G_AY1Y2(-0.5, d1, d2, b * gauss_n5[i], d * gauss_n5[j], d3)
+					* b * d * gauss_w5[i] * gauss_w5[j];
+		}
 	}
-	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);
+
+	return sol;
+
+}
+
+double calcParIntQY1X1(double b, double d, double t, double v, double d1,
+		double d2, double d3) {
+	//Fall 3
+	double sol = 0;
+
+	for (int i = 0; i < 5; ++i) {
+		for (int j = 0; j < 5; ++j) {
+			sol += G_AY2X2( t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i],d , d3) * t * b * gauss_w5[i] * gauss_w5[j];
+			sol -= G_AY2X2( t * gauss_n5[j] + d1, v + d2, b * gauss_n5[i], 0, d3) * t * b * gauss_w5[i] * gauss_w5[j];
+			sol -= G_AY2X2( t * gauss_n5[j] + d1, d2, b * gauss_n5[i], d, d3) * t * b * gauss_w5[i] * gauss_w5[j];
+			sol += G_AY2X2( t * gauss_n5[j] + d1, d2, b * gauss_n5[i], 0, d3) * t * b * gauss_w5[i] * gauss_w5[j];
+		}
 	}
-	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 sol;
+
+}
+
+double calcParIntQY1(double b, double d, double t, double v, double d1, double d2, double d3) {
+	//Fall 4
+	return apply0Int4(F_par, b, d, t, v, d1, d2, d3); ///ACHTUNG noch Falsch
+
+}
+
+double calcOrtIntA(double b, double d, double t, double v, double d1, double d2, double d3) {
+	return apply0Int4(F_ort, b, d, t, v, d1, d2, d3);
+
+}
+
+double calcOrtIntQY1Y2(double b, double d, double t, double v, double d1, double d2, double d3) {
+	return 0;
+
+}
+
+
+double calcParIntO1(double b, double d, double t, double v, double d1, double d2, double d3) {
+
+	//Elmente anordnen
+	if(b*b+d*d>t*t+v*v){
+		double tmp = 0;
+
+		tmp = b; b = t; t = tmp;
+		tmp = d; d = v; v = tmp;
+
+		d1 = -d1;
+		d2 = -d2;
+		d3 = -d3;
 	}
 
-	return (y1 - x1) * gL + (y2 - x2) * gK - a * a * G3;
+	return calcParIntA( b, d, t, v, d1, d2, d3);
 }
 
-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:
-			cout << "p must be either 0, -1/2, -1 or -3/2. (" << p << ")\n";
-			return NAN;
-		}
+double calcOrtIntO1(double b, double d, double t, double v, double d1, double d2, double d3) {
+
+	//Elmente anordnen
+	if(b*b+d*d>t*t+v*v){
+		double tmp = 0;
+
+		tmp = b; b = t; t = tmp;
+		tmp = d; d = v; v = tmp;
+
+		d1 = -d1;
+		d2 = -d2;
+		d3 = -d3;
+	}
+
+	return calcOrtIntA( b, d, t, v, d1, d2, d3);
+}
+
+
+double calcParIntO2(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {
+
+	//Elmente anordnen
+	if(b*b+d*d>t*t+v*v){
+		double tmp = 0;
+
+		tmp = b; b = t; t = tmp;
+		tmp = d; d = v; v = tmp;
+
+		d1 = -d1;
+		d2 = -d2;
+		d3 = -d3;
+	}
+
+	if((b*b+d*d) *eta <= d1*d1+d2*d2+d3*d3){
+		return calcParIntQX1X2( b, d, t, v, d1, d2, d3);
+	}else{
+		return calcParIntA( b, d, t, v, d1, d2, d3);
+	}
+
+
+}
+
+double calcOrtIntO2(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {
+
+	//Elmente anordnen
+	if(b*b+d*d>t*t+v*v){
+		double tmp = 0;
+
+		tmp = b; b = t; t = tmp;
+		tmp = d; d = v; v = tmp;
+
+		d1 = -d1;
+		d2 = -d2;
+		d3 = -d3;
+	}
+
+	if((b*b+d*d) *eta <= d1*d1+d2*d2+d3*d3){
+		return 0;
+	}else{
+		return calcOrtIntA( b, d, t, v, d1, d2, d3);
 	}
 }
 
-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;
+
+double calcParIntO3(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {
+
+	//Achsen vertauschen
+	if(b*b+t*t>d*d+v*v){
+		double tmp = 0;
+
+		tmp = b; b = d; d = tmp;
+		tmp = t; t = v; v = tmp;
+
+		tmp = d1; d1 = d2; d2 = tmp;
+	}
+
+	if(max(b,t) *eta <= d1){
+		return calcParIntQY1X1(b, d, t, v, d1, d2, d3);
+	}else{
+		return calcParIntA( b, d, t, v, d1, d2, d3);
+	}
+
+
 }
 
-double FLO_perp(double x1, double x2, double y2, double y3, double delta1,
-		double delta2, double a) {
-	double xd1 = x1 - delta1;
-	double xd2 = x2 - delta2;
-	double yd2 = y2 - delta2;
-	double yd3 = y3 - a;
-	double tmp2 = x2 - y2 - delta2;
-	double rval = 0;
-
-	rval = xd1 * compute_g0(y3, a, sqrt(xd1 * xd1 + tmp2 * tmp2), -0.5) + yd3
-			* compute_g0(x1, delta1, sqrt(tmp2 * tmp2 + yd3 * yd3), -0.5) - xd1
-			* tmp2 * slpADLO(x2, y3, yd2, a, xd1) - yd3 * tmp2 * slpADLO(x1,
-			x2, delta1, yd2, -y3 - a) - slpADLO(y3, x1, -a, delta1, tmp2);
-	return rval;
+double calcOrtIntO3(double b, double d, double t, double v, double d1, double d2, double d3,double eta) {
+
+	//Elmente anordnen
+	if(b*b+d*d>t*t+v*v){
+		double tmp = 0;
+
+		tmp = b; b = t; t = tmp;
+		tmp = d; d = v; v = tmp;
+
+		d1 = -d1;
+		d2 = -d2;
+		d3 = -d3;
+	}
+
+	if(max(b,t) *eta > d1){
+		return 0;
+	}else{
+		return calcOrtIntA( b, d, t, v, d1, d2, d3);
+	}
 }
diff --git a/src/slpRectangle.hpp b/src/slpRectangle.hpp
index 5133e06..f81dfd5 100644
--- a/src/slpRectangle.hpp
+++ b/src/slpRectangle.hpp
@@ -1,66 +1,60 @@
 #ifndef HILBERT3D_LAPLACE_SLPRECTANGLE_HPP_GUARD_
 #define HILBERT3D_LAPLACE_SLPRECTANGLE_HPP_GUARD_
 
-int sign(double);
-double arsinh(double);
+//int sign(double);
+//double arsinh(double);
 
-// sol = g0(p,y,x,l); Quadratur !!! 0-y
-double singleQuad(double, double, double, double);
-
-// sol = G00(p,y1,y2,x1,x2,l); Quadratur !!! 0-y1 & 0-y2
-double doubleQuad(double, double, double, double, double, 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);
-#define F_par F_parY1Y2
-// sol = F_ort(x1,x2,y2,y3,d1,d2,d3);
-double F_ort(double, double, double, double, double, double, double);
-
-
+double g_AY(double, double, double, double);
+// sol = g0(p,y,x,l);
+double g_QY(double, double, double, double);
 
-// sol = F_par(x1,x2,y1,y2,d1,d2,d3);
-double F_parY1Y2(double, double, double, double, double, double, double);
-// sol = F_par(x1,x2,y1,y2,d1,d2,d3);
-double F_parY1X2(double, double, double, double, double, double, double);
-// sol = F_par(x1,x2,y1,y2,d1,d2,d3);
-double F_parX1Y2(double, double, double, double, double, double, double);
 
-// sol = F_par(x1,x2,y1,y2,d1,d2,d3); Quadratur!!!
-double F_parY1Y2_2(double, double, double, double, double, double, double);
-// sol = F_par(x1,x2,y1,y2,d1,d2,d3); Quadratur!!!
-double F_parQY2X2(double, double, double, double, double, double, double);
+// sol = G00(p,y1,y2,x1,x2,l);
+double G_AY1Y2(double, double, double, double, double, double);
+// sol = G00(p,y1,y2,x1,x2,l);
+double G_AY2X2(double, double, double, double, double, double);
+// sol = G00(p,y1,y2,x1,x2,l);
+double G_QY1Y2(double, double, double, double, double, double);
 
-// funktionen von Maik
-double slpADLO(double, double, double, double, double);
-double compute_g0(double, double, double, double);
-double FLO_plane(double, double, double, double, double, double, double);
-double FLO_perp(double, 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,y2,y3,d1,d2,d3);
+double F_ort(double, double, double, double, double, double, double);
 
-// sol = quadInt((F_par/F_ort),a,b,c,d,r,t,u,v,d1,d2,d3);
-double applyInt4(
-		double(*)(double, double, double, double, double, double, double),
-		double, double, double, double, double, double, double, double, double,
-		double, double);
 
-// sol = quad0Int((F_par/F_ort),b,d,t,v,d1,d2,d3);
+// sol = quad0Int4((F_par/F_ort),b,d,t,v,d1,d2,d3);
 double apply0Int4(
 		double(*)(double, double, double, double, double, double, double),
 		double, double, double, double, double, double, double);
 
-double apply0Int2(
-		double(*)(double, double, double, double, double, double, double),
-		double, double, double, double, double, double, double);
+// sol = calcParInt.(b,d,t,v,d1,d2,d3);
+double calcParIntA(double, double, double, double, double, double, double);
+double calcParIntQX1X2(double, double, double, double, double, double, double);
+double calcParIntQY1X1(double, double, double, double, double, double, double);
 
+double calcParIntQX1(double, double, double, double, double, double, double);
 
-double calcParIntA(double, double, double, double, double, double, double);
 
+// sol = calcOrtInt.(b,d,t,v,d1,d2,d3);
+double calcOrtIntA(double, double, double, double, double, double, double);
+double calcOrtIntQX1X2(double, double, double, double, double, double, double);
 
+// Elemente Vertauschen
+// sol = calc...Int01(b,d,t,v,d1,d2,d3);
+double calcParIntO1(double, double, double, double, double, double, double);
+double calcOrtIntO1(double, double, double, double, double, double, double);
 
+// Quadratur ueber kleineres Element
+// sol = calc...Int02(b,d,t,v,d1,d2,d3,eta);
+double calcParIntO2(double, double, double, double, double, double, double, double);
+double calcOrtIntO2(double, double, double, double, double, double, double, double);
 
+// Quadratur kuerzere Achse/Seite
+// sol = calc...Int03(b,d,t,v,d1,d2,d3,eta);
+double calcParIntO3(double, double, double, double, double, double, double, double);
+double calcOrtIntO3(double, double, double, double, double, double, double, double);
 
 #endif
diff --git a/src/test.cpp b/src/test.cpp
deleted file mode 100644
index 04b69a7..0000000
--- a/src/test.cpp
+++ /dev/null
@@ -1,79 +0,0 @@
-/*
- * test.cpp
- *
- *  Created on: 06.04.2011
- *      Author: treecity
- */
-
-#include <iostream>
-#include <cmath>
-
-using namespace std;
-
-
-	int static sign(double x)
-	{
-		return x > 0 ? 1 : (x < 0 ? -1 : 0);
-	}
-
-	/***
-	 * ! Diese Funktionen benutzen nur die ersten beiden Dimensionen !
-	 */
-
-	double g(double p, double y, double x, double lambda){
-		double sol = NAN;
-
-		if(lambda==0){	//Lemma B1
-			if(p!=-0.5)
-				sol = (y-x)*pow(fabs(y-x),2*p)/(2*p+1); //(p<.5) & x außerhalb ?
-			else
-				sol = sign(y-x)*log(fabs(y-x));	//x außerhalb ?
-		}else{	//Lemma B2 (wird eigententlich nicht gebraucht)
-
-
-		}
-	return sol;
-	}
-
-	double G(double p, double * y, double * x, double lambda) {
-		double sol =NAN;
-		if(p==-1.5){
-			if(lambda==0){
-
-			}else{
-
-			}
-
-		}else if(p==-0.5){
-
-		}
-
-
-		return sol;
-	}
-
-	double F(double* x, double * y, double delta1, double delta2,double delta3){
-
-		return NAN;
-	}
-
-int main(){
-
-	double p=-2,y=2,x=3.5,l=0;
-
-
-	cout << "g(" << p << "," << y << "," << x << "," << l << ")=" << g(p,y,x,l) << "\n";
-		p = -0.5;
-	cout << "g(" << p << "," << y << "," << x << "," << l << ")=" << g(p,y,x,l) << "\n";
-		x = 1.5;
-	cout << "g(" << p << "," << y << "," << x << "," << l << ")=" << g(p,y,x,l) << "\n";
-		p = 0;
-	cout << "g(" << p << "," << y << "," << x << "," << l << ")=" << g(p,y,x,l) << "\n";
-		l = -2;
-	cout << "g(" << p << "," << y << "," << x << "," << l << ")=" << g(p,y,x,l) << "\n";
-
-
-
-
-	return 0;
-}
diff --git a/src/test_Fpar.m b/src/test_Fpar.m
index 920441c..0f394ef 100644
--- a/src/test_Fpar.m
+++ b/src/test_Fpar.m
@@ -32,6 +32,13 @@ function [ret] = test_Fpar(a,b,c,d,r,t,u,v,l1,l2,l3)
     
     sol3 = surfQuad(@(x1,y1)(int2(x1,y1,c,d,u,v)),a,b,r,t,5);
     
-    ret = [quad_sol sol sol2 sol3];
+    sol4 =...
+      surfQuad(@(x1,y1)mex_FparQY2X2(x1,d,y1,v,l1,l2,l3),a,b,r,t,5)...
+      -surfQuad(@(x1,y1)mex_FparQY2X2(x1,c,y1,v,l1,l2,l3),a,b,r,t,5)...
+      -surfQuad(@(x1,y1)mex_FparQY2X2(x1,d,y1,u,l1,l2,l3),a,b,r,t,5)...
+      +surfQuad(@(x1,y1)mex_FparQY2X2(x1,c,y1,u,l1,l2,l3),a,b,r,t,5);
+
+    
+    ret = [quad_sol sol sol2 sol3 sol4];
 end
 
diff --git a/src/test_solveError_Grid3.m b/src/test_solveError_Grid3.m
new file mode 100755
index 0000000..df116a6
--- /dev/null
+++ b/src/test_solveError_Grid3.m
@@ -0,0 +1,132 @@
+
+mex mex_build_As3.cpp slpRectangle.cpp
+
+dataIso =[];
+dataAniso =[];
+maxSize = 500;
+fileName = './test_save';
+Netz = 'exmpl_2DLShape'; %Energienorm sollte gegen 8.28466 konvergieren
+% Netz = 'exmpl_3DFichCube'; %Energienorm sollte gegen 162265 konvergieren
+% Netz = 'exmpl_2DQuad';
+
+
+%% Anisotrope Verfeinerung
+
+if(~exist('d1'))
+  
+  disp 'Normal'
+  load(Netz)
+
+  ref = 0;
+
+while size(elements,1)<maxSize
+  tic
+  ref = ref+1; 
+  
+
+
+  %Netz Isotrop verfeinern (4 Teile)
+  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
+  
+  %Matrix mit feinem Netz berechnen
+  A_fine = mex_build_A(coordinates_fine,elements_fine);
+  
+  %rechte Seite Aufbauen
+  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
+  
+  %Lösungsvektor bestimmen
+  x_fine = A_fine\b;
+  
+  %Energienorm^2 bestimmen
+  xe_fine = x_fine'*A_fine*x_fine;
+
+  % Welche Elemente sollen verfeinert werden
+  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);
+  
+  % Wie soll verfeinert werden
+  marked = mark(x_fine(f2s)',ind,0.5,0.5);
+
+  %Daten zur Auswertung zwischenspeichern (testweise?)
+  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];
+  
+  %Sicherheit falls keine Elemente markiert sind
+  if(sum(marked~=1)==0)
+    disp 'Abgebrochen weil nicht mehr verfeinert wird!'
+    break;
+  end
+  
+  %Verfeinern
+  %[coordinates, elements] = refineQuad(coordinates,elements,marked);
+  file = ['save' int2str(ref)];
+  load(file);
+  toc
+end
+d1 = dataAniso
+end
+
+%% Anisotrope Verfeinerung
+disp 'Optimiert'
+load(Netz)
+
+ref = 0;
+while size(elements,1)<maxSize
+  tic
+  ref = ref+1; 
+  
+
+
+  %Netz Isotrop verfeinern (4 Teile)
+  [coordinates_fine,elements_fine, f2s]=refineQuad(coordinates,elements,2*ones(1,size(elements,1)));
+  
+  %Matrix mit feinem Netz berechnen
+  A_fine = mex_build_As3(coordinates_fine,elements_fine);
+  
+  %rechte Seite Aufbauen
+  b = sqrt(sum(quadNorm(coordinates_fine,elements_fine,'w').^2,2));
+  
+  %Lösungsvektor bestimmen
+  x_fine = A_fine\b;
+  
+  %Energienorm^2 bestimmen
+  xe_fine = x_fine'*A_fine*x_fine;
+
+  % Welche Elemente sollen verfeinert werden
+  ind = computeEstSlpMuTilde(x_fine,coordinates,elements,f2s);
+  
+  % Wie soll verfeinert werden
+  marked = mark(x_fine(f2s)',ind,0.5,0.5);
+
+  %Daten zur Auswertung zwischenspeichern (testweise?)
+  dataAniso(ref,:) = [size(elements,1) sqrt(sum(ind)) xe_fine];
+  
+  %Sicherheit falls keine Elemente markiert sind
+%   if(sum(marked~=1)==0)
+%     disp 'Abgebrochen weil nicht mehr verfeinert wird!'
+%     break;
+%   end
+  
+  %Verfeinern
+  %[coordinates, elements] = refineQuad(coordinates,elements,marked);
+  file = ['save' int2str(ref)];
+  load(file);
+  toc
+end
+d2 = dataAniso
+
+clear ref b A_fine coordinates_fine elements_fine marked ind f2s x_fine
+
+
+%% Alles Zeichnen
+figure (1)
+loglog(d1(:,1),d1(:,2),d2(:,1),d2(:,2));
+
+legend('Normal','optimiert')
+xlabel('Elemente')
+ylabel('Fehler')
+
+figure (2)
+loglog(d1(:,1),d1(:,3),d2(:,1),d2(:,3));
+
+legend('Normal','Optimiert')
+xlabel('Elemente')
+ylabel('EnergieNorm')