]> git.leopard-lacewing.eu Git - bacc.git/commitdiff
[src] dist Funktionen korrigiert
authorPeter Schaefer <peter.schaefer@tuwien.ac.at>
Thu, 21 Mar 2013 19:37:10 +0000 (20:37 +0100)
committerPeter Schaefer <peter.schaefer@tuwien.ac.at>
Thu, 21 Mar 2013 19:37:10 +0000 (20:37 +0100)
[src] laengenparser vereinfacht

src/mex_build_V.cpp
src/slpRectangle.cpp

index c36156951512f23a52f55cf7fbe62426bd73fe7f..3e2dbfebea4a03cc663e0762e3270eff3b58265a 100644 (file)
@@ -26,7 +26,7 @@
 /***************************************************************************/\r
 \r
 #define DEBUG 1\r
-#define PARALLEL\r
+//#define PARALLEL\r
 \r
 //#include <iostream>\r
 #include <cmath>\r
@@ -51,20 +51,47 @@ int inline dimOfVec(double* vec) {
        else if (vec[0] != 0)\r
                return 0;\r
        else {\r
-//             mexErrMsgTxt("length of Site is zero");\r
+               mexErrMsgTxt("length of Site is zero");\r
                return -1;\r
        }\r
 }\r
 \r
+//// Gibt von [0 1 2] die Fehlende Zahl zurueck\r
+//int inline dimOfThird(int a, int b) {\r
+//     return ((-(a + b) % 3) + 3) % 3;\r
+//}\r
+\r
 // Gibt von [0 1 2] die Fehlende Zahl zurueck\r
 int inline dimOfThird(int a, int b) {\r
-       return ((-(a + b) % 3) + 3) % 3;\r
+       switch (a + b) {\r
+       case 1:\r
+               return 2;\r
+       case 2:\r
+               return 1;\r
+       case 3:\r
+               return 0;\r
+       }\r
 }\r
 \r
 void inline add(double* a, double* b) {\r
        for (int i = 0; i < 3; ++i)\r
                a[i] += b[i];\r
 }\r
+\r
+bool inline iseq(double* a, double* b) {\r
+       for (int i = 0; i < 3; ++i)\r
+               if (a[i] != b[i])\r
+                       return false;\r
+       return true;\r
+}\r
+\r
+bool inline iseq(double* a, double b) {\r
+       for (int i = 0; i < 3; ++i)\r
+               if (a[i] != b)\r
+                       return false;\r
+       return true;\r
+}\r
+\r
 void inline add(double* a, double* b, double* c) {\r
        for (int i = 0; i < 3; ++i)\r
                a[i] = b[i] + c[i];\r
@@ -102,19 +129,19 @@ void inline getSCorner(double* d, double x[4][3]) {
        }\r
 }\r
 \r
-double inline distT(double b, double d, double t, double v, double d1,\r
+double inline distT_par(double b, double d, double t, double v, double d1,\r
                double d2, double d3, double * dummy) {\r
        double dis1 = 0, dis2 = 0;\r
        if (d1 < 0)\r
                dis1 = -d1 - t;\r
-       else\r
+       else if(d1 > 0)\r
                dis1 = d1 - b;\r
        if (dis1 < 0)\r
                dis1 = 0;\r
 \r
        if (d2 < 0)\r
                dis2 = -d2 - v;\r
-       else\r
+       else if(d2 > 0)\r
                dis2 = d2 - d;\r
        if (dis2 < 0)\r
                dis2 = 0;\r
@@ -123,9 +150,35 @@ double inline distT(double b, double d, double t, double v, double d1,
 //             mexErrMsgTxt("expected d3 == 0");\r
 \r
 //     return dis1 * dis1 + dis2 * dis2 + d3 * d3;\r
-       return (dis1 * dis1 + dis2 * dis2 + d3 * d3) * (4 * M_PI);\r
+       return sqrt(dis1 * dis1 + dis2 * dis2 + d3 * d3) * (4 * M_PI);\r
 }\r
+double inline distT_ort(double b, double d, double t, double v, double d1,\r
+               double d2, double d3, double * dummy) {\r
+       double dis1 = 0, dis2 = 0, dis3 = 0;\r
+       if (d1 < 0)\r
+               dis1 = -d1;\r
+       else\r
+               dis1 = d1 - b;\r
+       if (dis1 < 0)\r
+               dis1 = 0;\r
 \r
+       if (d2 < 0)\r
+               dis2 = -d2 - t;\r
+       else\r
+               dis2 = d2 - d;\r
+       if (dis2 < 0)\r
+               dis2 = 0;\r
+\r
+       if (d3 < 0)\r
+               dis3 = -d3 - v;\r
+       else\r
+               dis3 = d3;\r
+       if (dis3 < 0)\r
+               dis3 = 0;\r
+\r
+//     return dis1 * dis1 + dis2 * dis2 + dis3 * dis3;\r
+       return sqrt(dis1 * dis1 + dis2 * dis2 + dis3 * dis3) * (4 * M_PI);\r
+}\r
 void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {\r
        /*\r
         mexPrintf("%d",omp_get_max_threads());\r
@@ -203,8 +256,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
                ctypeO = cOrtO1;\r
                break;\r
        case 0:\r
-               ctypeP = distT;\r
-               ctypeO = distT;\r
+               ctypeP = distT_par;\r
+               ctypeO = distT_ort;\r
                break;\r
        }\r
 \r
@@ -257,143 +310,165 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
 #pragma omp parallel for schedule(static,1) private(i,j,k,tmp,x,y,xa,xb,ya,yb,d,dt,rx,rxa,rxb,ry,rya,ryb,lastRow,firstRow) shared(C,E,ctypeO,ctypeP,em,targetSize,actualNumberOfThreads)\r
 \r
 #endif\r
-               for (i = 0; i < actualNumberOfThreads; ++i) {\r
-                       firstRow = (int) sqrt(2 * i * targetSize);\r
-                       if (i == actualNumberOfThreads - 1) {\r
-                               lastRow = em - 1;\r
-                       } else {\r
-                               lastRow = (int) sqrt(2 * (i + 1) * targetSize) - 1;\r
-                       }\r
+       for (i = 0; i < actualNumberOfThreads; ++i) {\r
+               firstRow = (int) sqrt(2 * i * targetSize);\r
+               if (i == actualNumberOfThreads - 1) {\r
+                       lastRow = em - 1;\r
+               } else {\r
+                       lastRow = (int) sqrt(2 * (i + 1) * targetSize) - 1;\r
+               }\r
 \r
 #if DEBUG >= 4\r
 #pragma omp critical\r
-                       mexPrintf("thread Nr.: (%d) | firstRow/lastRow: [%d / %d] \n",\r
+               mexPrintf("thread Nr.: (%d) | firstRow/lastRow: [%d / %d] \n",\r
 \r
 #ifdef PARALLEL\r
-                                       omp_get_thread_num()\r
+                               omp_get_thread_num()\r
 #else\r
-                                       i\r
+                               i\r
 #endif\r
-                                       , firstRow, lastRow);\r
+                               , firstRow, lastRow);\r
 #pragma omp end critical\r
 #endif\r
-                       {\r
-                               for (j = firstRow; j <= lastRow; ++j) {\r
-                                       x[0][0] = C[(int) E[j] - 1];\r
-                                       x[0][1] = C[cm + (int) E[j] - 1];\r
-                                       x[0][2] = C[2 * cm + (int) E[j] - 1];\r
+               {\r
+                       for (j = firstRow; j <= lastRow; ++j) {\r
+                               x[0][0] = C[(int) E[j] - 1];\r
+                               x[0][1] = C[cm + (int) E[j] - 1];\r
+                               x[0][2] = C[2 * cm + (int) E[j] - 1];\r
 \r
-                                       x[1][0] = C[(int) E[em + j] - 1];\r
-                                       x[1][1] = C[cm + (int) E[em + j] - 1];\r
-                                       x[1][2] = C[2 * cm + (int) E[em + j] - 1];\r
+                               x[1][0] = C[(int) E[em + j] - 1];\r
+                               x[1][1] = C[cm + (int) E[em + j] - 1];\r
+                               x[1][2] = C[2 * cm + (int) E[em + j] - 1];\r
 \r
-                                       x[2][0] = C[(int) E[2 * em + j] - 1];\r
-                                       x[2][1] = C[cm + (int) E[2 * em + j] - 1];\r
-                                       x[2][2] = C[2 * cm + (int) E[2 * em + j] - 1];\r
+                               x[2][0] = C[(int) E[2 * em + j] - 1];\r
+                               x[2][1] = C[cm + (int) E[2 * em + j] - 1];\r
+                               x[2][2] = C[2 * cm + (int) E[2 * em + j] - 1];\r
 \r
-                                       x[3][0] = C[(int) E[3 * em + j] - 1];\r
-                                       x[3][1] = C[cm + (int) E[3 * em + j] - 1];\r
-                                       x[3][2] = C[2 * cm + (int) E[3 * em + j] - 1];\r
+                               x[3][0] = C[(int) E[3 * em + j] - 1];\r
+                               x[3][1] = C[cm + (int) E[3 * em + j] - 1];\r
+                               x[3][2] = C[2 * cm + (int) E[3 * em + j] - 1];\r
 \r
-                                       // Seitenvektoren aufbauen\r
-                                       sub(xa, x[1], x[0]);\r
-                                       sub(xb, x[3], x[0]);\r
+                               // Seitenvektoren aufbauen\r
+                               sub(xa, x[1], x[0]);\r
+                               sub(xb, x[3], x[0]);\r
 \r
-                                       // Lageeigenschaften des Flaechenstuecks\r
-                                       rxa = dimOfVec(xa);\r
-                                       rxb = dimOfVec(xb);\r
-                                       rx = dimOfThird(rxa, rxb);\r
+                               // Lageeigenschaften des Flaechenstuecks\r
+                               rxa = dimOfVec(xa);\r
+                               rxb = dimOfVec(xb);\r
+                               rx = dimOfThird(rxa, rxb);\r
+\r
+                               //kleinste Ecke finden und fuer \delta verwenden\r
+//                                     getSCorner(dt, x);\r
+\r
+                               if (xa[rxa] <= 0 || xb[rxb] <= 0)\r
+                                       mexWarnMsgTxt("X Element hat Seitenlängen <=0!");\r
+\r
+                               if(rxa == rxb)\r
+                                       mexWarnMsgTxt("X Laenge ist nur in einer Richtung");\r
+\r
+//                                     if (!iseq(dt, x[0]))\r
+//                                             mexWarnMsgTxt("x[0] ist nicht die kleinste Ecke");\r
+\r
+                               for (k = 0; k <= j; ++k) {\r
+\r
+                                       y[0][0] = C[(int) E[k] - 1];\r
+                                       y[0][1] = C[cm + (int) E[k] - 1];\r
+                                       y[0][2] = C[2 * cm + (int) E[k] - 1];\r
+\r
+                                       y[1][0] = C[(int) E[em + k] - 1];\r
+                                       y[1][1] = C[cm + (int) E[em + k] - 1];\r
+                                       y[1][2] = C[2 * cm + (int) E[em + k] - 1];\r
 \r
-                                       //kleinste Ecke finden und fuer \delta verwenden\r
-                                       getSCorner(dt, x);\r
+                                       y[2][0] = C[(int) E[2 * em + k] - 1];\r
+                                       y[2][1] = C[cm + (int) E[2 * em + k] - 1];\r
+                                       y[2][2] = C[2 * cm + (int) E[2 * em + k] - 1];\r
 \r
-                                       for (k = 0; k <= j; ++k) {\r
+                                       y[3][0] = C[(int) E[3 * em + k] - 1];\r
+                                       y[3][1] = C[cm + (int) E[3 * em + k] - 1];\r
+                                       y[3][2] = C[2 * cm + (int) E[3 * em + k] - 1];\r
 \r
-                                               y[0][0] = C[(int) E[k] - 1];\r
-                                               y[0][1] = C[cm + (int) E[k] - 1];\r
-                                               y[0][2] = C[2 * cm + (int) E[k] - 1];\r
+                                       // Seiten Vektoren aufbauen\r
+                                       sub(ya, y[1], y[0]);\r
+                                       sub(yb, y[3], y[0]);\r
 \r
-                                               y[1][0] = C[(int) E[em + k] - 1];\r
-                                               y[1][1] = C[cm + (int) E[em + k] - 1];\r
-                                               y[1][2] = C[2 * cm + (int) E[em + k] - 1];\r
+                                       // Lageeigenschaften des Flaechenstuecks\r
+                                       rya = dimOfVec(ya);\r
+                                       ryb = dimOfVec(yb);\r
+                                       ry = dimOfThird(rya, ryb);\r
+\r
+                                       //kleinste Ecke finden\r
+//                                             getSCorner(d, y);\r
 \r
-                                               y[2][0] = C[(int) E[2 * em + k] - 1];\r
-                                               y[2][1] = C[cm + (int) E[2 * em + k] - 1];\r
-                                               y[2][2] = C[2 * cm + (int) E[2 * em + k] - 1];\r
+                                       if (ya[rya] <= 0 || yb[ryb] <= 0)\r
+                                               mexWarnMsgTxt("Y Element hat Seitenlängen <=0!");\r
 \r
-                                               y[3][0] = C[(int) E[3 * em + k] - 1];\r
-                                               y[3][1] = C[cm + (int) E[3 * em + k] - 1];\r
-                                               y[3][2] = C[2 * cm + (int) E[3 * em + k] - 1];\r
+                                       if(rya == ryb)\r
+                                               mexWarnMsgTxt("Y Laenge ist nur in einer Richtung");\r
 \r
-                                               // Seiten Vektoren aufbauen\r
-                                               sub(ya, y[1], y[0]);\r
-                                               sub(yb, y[3], y[0]);\r
+//                                     if (!iseq(d, y[0]))\r
+//                                             mexWarnMsgTxt("y[0] ist nicht die kleinste Ecke");\r
 \r
-                                               // Lageeigenschaften des Flaechenstuecks\r
-                                               rya = dimOfVec(ya);\r
-                                               ryb = dimOfVec(yb);\r
-                                               ry = dimOfThird(rya, ryb);\r
+                                       //\delta berechnen\r
+//                                             sub(d, dt);\r
 \r
-                                               //kleinste Ecke finden\r
-                                               getSCorner(d, y);\r
+                                       sub(d, y[0], x[0]);\r
+\r
+                                       if(j==1 && k==0){\r
+                                               mexPrintf("[%d %d](%f,%f)(%f,%f)+(%f,%f,%f)", j, k,\r
+                                                               (xa[rxa]), (xb[rxb]), (ya[rya]),\r
+                                                               (yb[ryb]), d[rxa], d[rxb], d[rx]);\r
+                                               mexPrintf("{%d,%d;%d,%d}\n",rxa,rxb,rya,rxb);\r
+                                       }\r
 \r
-                                               //\delta berechnen\r
-                                               sub(d, dt);\r
 #if DEBUG >= 5\r
 #pragma omp critical\r
-                                               mexPrintf("[%d %d](%f,%f)(%f,%f)+(%f,%f,%f)\n", j, k,\r
-                                                               fabs(xa[rxa]), fabs(xb[rxb]), fabs(ya[rxa]),\r
-                                                               fabs(yb[rxb]), d[rxa], d[rxb], d[rx]);\r
+                                       mexPrintf("[%d %d](%f,%f)(%f,%f)+(%f,%f,%f)\n", j, k,\r
+                                                       (xa[rxa]), (xb[rxb]), (ya[rya]),\r
+                                                       (yb[ryb]), d[rxa], d[rxb], d[rx]);\r
 #pragma omp end critical\r
 #endif\r
-                                               if (rx == ry) {\r
-                                                       if (rxa == rya) {\r
-                                                               tmp = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]),\r
-                                                                               fabs(ya[rxa]), fabs(yb[rxb]), d[rxa],\r
-                                                                               d[rxb], d[rx], zeta);\r
-\r
-                                                       } else {\r
-                                                               tmp = ctypeP(fabs(xa[rxa]), fabs(xb[rxb]),\r
-                                                                               fabs(yb[rxa]), fabs(ya[rxb]), d[rxa],\r
-                                                                               d[rxb], d[rx], zeta);\r
-                                                       }\r
+                                       if (rx == ry) {\r
+                                               if (rxa == rya) {\r
+                                                       tmp = ctypeP((xa[rxa]), (xb[rxb]), (ya[rxa]),\r
+                                                                       (yb[rxb]), d[rxa], d[rxb], d[rx], zeta);\r
 \r
                                                } else {\r
+                                                       tmp = ctypeP((xa[rxa]), (xb[rxb]), (yb[rxa]),\r
+                                                                       (ya[rxb]), d[rxa], d[rxb], d[rx], zeta);\r
+                                               }\r
 \r
-                                                       if (rxa == rya) {\r
-                                                               tmp = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]),\r
-                                                                               fabs(ya[rya]), fabs(yb[ryb]), d[rxb],\r
-                                                                               d[rxa], d[rx], zeta);\r
-                                                       } else if (rxa == ryb) {\r
-                                                               tmp = ctypeO(fabs(xb[rxb]), fabs(xa[rxa]),\r
-                                                                               fabs(yb[ryb]), fabs(ya[rya]), d[rxb],\r
-                                                                               d[rxa], d[rx], zeta);\r
-                                                       } else if (rxb == rya) {\r
-                                                               tmp = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]),\r
-                                                                               fabs(ya[rya]), fabs(yb[ryb]), d[rxa],\r
-                                                                               d[rxb], d[rx], zeta);\r
-                                                       } else {\r
-                                                               tmp = ctypeO(fabs(xa[rxa]), fabs(xb[rxb]),\r
-                                                                               fabs(yb[ryb]), fabs(ya[rya]), d[rxa],\r
-                                                                               d[rxb], d[rx], zeta);\r
-                                                       }\r
-\r
+                                       } else {\r
+\r
+                                               if (rxa == rya) {\r
+                                                       tmp = ctypeO((xb[rxb]), (xa[rxa]), (ya[rya]),\r
+                                                                       (yb[ryb]), d[rxb], d[rxa], d[rx], zeta);\r
+                                               } else if (rxa == ryb) {\r
+                                                       tmp = ctypeO((xb[rxb]), (xa[rxa]), (yb[ryb]),\r
+                                                                       (ya[rya]), d[rxb], d[rxa], d[rx], zeta);\r
+                                               } else if (rxb == rya) {\r
+                                                       tmp = ctypeO((xa[rxa]), (xb[rxb]), (ya[rya]),\r
+                                                                       (yb[ryb]), d[rxa], d[rxb], d[rx], zeta);\r
+                                               } else {\r
+                                                       tmp = ctypeO((xa[rxa]), (xb[rxb]), (yb[ryb]),\r
+                                                                       (ya[rya]), d[rxa], d[rxb], d[rx], zeta);\r
                                                }\r
 \r
+                                       }\r
+\r
 #pragma omp critical\r
 \r
-                                               A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r
-                                               if (k != j)\r
-                                                       A[(j * em) + k] = 1. / (4 * M_PI) * tmp;\r
+                                       A[(k * em) + j] = 1. / (4 * M_PI) * tmp;\r
+                                       if (k != j)\r
+                                               A[(j * em) + k] = 1. / (4 * M_PI) * tmp;\r
 \r
 #pragma omp end critical\r
 \r
-                                       }\r
                                }\r
                        }\r
                }\r
-#ifdef PARALLEL\r
        }\r
+#ifdef PARALLEL\r
+}\r
 #endif\r
 \r
        return;\r
index 9828f44cb2ecea1cc35737d5fa4ba3fd41de59c1..04a0c8ab5f7fde33b18bac963f2dd1214e998abe 100644 (file)
@@ -97,7 +97,7 @@ void inline switch_dim(double& b, double& d, double& t, double& v, double& d1,
 /*\r
  * gibt den Abstand zum Quadrat zurueck\r
  */\r
-double inline dist2(double b, double d, double t, double v, double d1,\r
+double inline dist2_par(double b, double d, double t, double v, double d1,\r
                double d2, double d3) {\r
        double dis1 = 0, dis2 = 0;\r
        if (d1 < 0)\r
@@ -116,14 +116,38 @@ double inline dist2(double b, double d, double t, double v, double d1,
 \r
        return dis1 * dis1 + dis2 * dis2 + d3 * d3;\r
 }\r
+\r
 /*\r
- * gibt den Abstand zurueck\r
+ * gibt den Abstand zum Quadrat zurueck\r
  */\r
-double inline dist(double b, double d, double t, double v, double d1, double d2,\r
-               double d3) {\r
-       return sqrt(dist2(b, d, t, v, d1, d2, d3));\r
+double inline dist2_ort(double b, double d, double t, double v, double d1,\r
+               double d2, double d3) {\r
+       double dis1 = 0, dis2 = 0, dis3 = 0;\r
+       if (d1 < 0)\r
+               dis1 = -d1;\r
+       else\r
+               dis1 = d1 - b;\r
+       if (dis1 < 0)\r
+               dis1 = 0;\r
+\r
+       if (d2 < 0)\r
+               dis2 = -d2 - t;\r
+       else\r
+               dis2 = d2 - d;\r
+       if (dis2 < 0)\r
+               dis2 = 0;\r
+\r
+       if (d3 < 0)\r
+               dis3 = -d3 - v;\r
+       else\r
+               dis3 = d3;\r
+       if (dis3 < 0)\r
+               dis3 = 0;\r
+\r
+       return dis1 * dis1 + dis2 * dis2 + dis3 * dis3;\r
 }\r
 \r
+\r
 /*\r
  * gibt den Abstand in einer Richtung zum Quadrat zurueck\r
  */\r
@@ -138,12 +162,7 @@ double inline dist_s2(double b, double t, double d1) {
 \r
        return dis1 * dis1;\r
 }\r
-/*\r
- * gibt den Abstand in einer Richtung zurueck\r
- */\r
-double inline dist_s(double b, double t, double d1) {\r
-       return sqrt(dist_s2(b, t, d1));\r
-}\r
+\r
 \r
 double inline f_par(double x1, double x2, double y1, double y2, double d1,\r
                double d2, double d3) {\r
@@ -288,7 +307,7 @@ double inline G_AY1Y2(double p, double y1, double y2, double x1, double x2,
 double inline Gs_AX2Y1Y2(double p, double y1, double y2, double x1, double x2,\r
                double l) {\r
        double sol = 0;\r
-       mexPrintf("warning: Gs_AX2Y1Y2 nicht implementiert!");\r
+       mexWarnMsgTxt("Gs_AX2Y1Y2 nicht implementiert!");\r
        return sol;\r
 }\r
 \r
@@ -320,7 +339,7 @@ double inline F_par(double x1, double x2, double y1, double y2, double d1,
 double inline F_ort(double x1, double x2, double y2, double y3, double d1,\r
                double d2, double d3) {\r
 \r
-       //       printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y2,y3,d1,d2,d3);\r
+//     mexPrintf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f \n",x1,x2,y2,y3,d1,d2,d3);\r
        double sol = -G_AY1Y2(0.5, y3, x1, -d3, d1, x2 - y2 - d2);\r
 \r
        if ((x1 - d1) * (x2 - y2 - d2) != 0)\r
@@ -348,19 +367,16 @@ double inline F_ort(double x1, double x2, double y2, double y3, double d1,
        return sol / 2.;\r
 }\r
 \r
-\r
-double inline FY1Y2_par(double x1, double x2, double y1, double y2, double d1, double d2,\r
-               double d3) {\r
+double inline FY1Y2_par(double x1, double x2, double y1, double y2, double d1,\r
+               double d2, double d3) {\r
        return G_AY1Y2(-0.5, y1 + d1, y2 + d2, x1, x2, d3);\r
 }\r
 \r
-double inline FY2X2_par(double y1, double x1, double y2, double x2, double d1, double d2,\r
-               double d3) {\r
+double inline FY2X2_par(double y1, double x1, double y2, double x2, double d1,\r
+               double d2, double d3) {\r
        return G_AY2X2(y1 + d1, y2 + d2, x1, x2, d3);\r
 }\r
 \r
-\r
-\r
 /*\r
  * Quadratur über zwei Integrale und Grenzen in Stammfunktion über zwei Integrale\r
  * b,d - Grenzen für Quadratur\r
@@ -439,11 +455,15 @@ double inline intA4(double b, double d, double t, double v, double d1,
 \r
 double cParO1(double b, double d, double t, double v, double d1, double d2,\r
                double d3, double* zeta) {\r
+//     if (d3 && fabs(d3)!=2)\r
+//             mexPrintf("%f \n", d3);\r
        return intA4<F_par>(b, d, t, v, d1, d2, d3);\r
 }\r
 \r
 double cOrtO1(double b, double d, double t, double v, double d1, double d2,\r
                double d3, double* zeta) {\r
+//     if (!d3)\r
+//             mexPrintf("%f \n", d3);\r
        return intA4<F_ort>(b, d, t, v, d1, d2, d3);\r
 }\r
 \r
@@ -457,13 +477,12 @@ double cParO2(double b, double d, double t, double v, double d1, double d2,
        //kurze Achse nach vorn\r
 //     switch_dim(b, d, t, v, d1, d2);\r
 \r
-       if (zeta[0] * zeta[0] * (t * t + v * v) < dist2(b, d, t, v, d1, d2, d3))\r
+       if (zeta[0] * zeta[0] * (t * t + v * v) < dist2_par(b, d, t, v, d1, d2, d3))\r
                return intQ4<f_par>(b, d, t, v, d1, d2, d3);\r
 \r
        return intA4<F_par>(b, d, t, v, d1, d2, d3);\r
 }\r
 \r
-\r
 double cOrtO2(double b, double d, double t, double v, double d1, double d2,\r
                double d3, double* zeta) {\r
        double tmp = 0, dis2 = 0;\r
@@ -481,23 +500,15 @@ double cOrtO2(double b, double d, double t, double v, double d1, double d2,
                tmp = b;\r
                b = v;\r
                v = tmp;\r
-               d1 = -d1;\r
-               d3 = -d3;\r
+               tmp = -d1;\r
+               d1 = -d3;\r
+               d3 = tmp;\r
        }\r
 \r
-       if (d2 < 0)\r
-               tmp = -d2 - t;\r
-       else\r
-               tmp = d2 - d;\r
-       if (tmp < 0)\r
-               tmp = 0;\r
-\r
-       dis2 = tmp * tmp + d1 * d1 + d3 * d3;\r
-\r
        //kurze Achse nach vorn\r
 //     switch_dim(b, d, t, v, d1, d2);\r
 \r
-       if (zeta[0] * zeta[0] * (t * t + v * v) < dis2)\r
+       if (zeta[0] * zeta[0] * (t * t + v * v) < dist2_ort(b,d,t,v,d1,d2,d3))\r
                return intQ4<f_ort>(b, d, t, v, d1, d2, d3);\r
 \r
        return intA4<F_ort>(b, d, t, v, d1, d2, d3);\r
@@ -512,10 +523,10 @@ double cParO3(double b, double d, double t, double v, double d1, double d2,
        //kurze Achse nach vorn\r
 //     switch_dim(b, d, t, v, d1, d2);\r
 \r
-       if (zeta[1] * zeta[1] * (b * b + d * d) < dist2(b, d, t, v, d1, d2, d3))\r
+       if (zeta[1] * zeta[1] * (b * b + d * d) < dist2_par(b, d, t, v, d1, d2, d3))\r
                return intQ2A2<FY1Y2_par>(b, d, t, v, d1, d2, d3);\r
 \r
-       if (zeta[0] * zeta[0] * (t * t + v * v) < dist2(b, d, t, v, d1, d2, d3))\r
+       if (zeta[0] * zeta[0] * (t * t + v * v) < dist2_par(b, d, t, v, d1, d2, d3))\r
                return intQ4<f_par>(b, d, t, v, d1, d2, d3);\r
 \r
        return intA4<F_par>(b, d, t, v, d1, d2, d3);\r
@@ -532,10 +543,10 @@ double cParO4(double b, double d, double t, double v, double d1, double d2,
        //kurze Achse nach vorn\r
        switch_dim(b, d, t, v, d1, d2);\r
 \r
-       if ((t * t + v * v) < zeta[0] * dist2(b, d, t, v, d1, d2, d3))\r
+       if ((t * t + v * v) < zeta[0] * dist2_par(b, d, t, v, d1, d2, d3))\r
                return intQ4<f_par>(b, d, t, v, d1, d2, d3);\r
 \r
-//     if ((b * b + d * d) < zeta[0] * dist2(b, d, t, v, d1, d2, d3))\r
+//     if ((b * b + d * d) < zeta[0] * dist2_par(b, d, t, v, d1, d2, d3))\r
 //             return calcParIntQX1X2(b, d, t, v, d1, d2, d3);\r
 \r
        return intA4<F_par>(b, d, t, v, d1, d2, d3);\r