/***************************************************************************/\r
\r
#define DEBUG 1\r
-#define PARALLEL\r
+//#define PARALLEL\r
\r
//#include <iostream>\r
#include <cmath>\r
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
}\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
// 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
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
#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
/*\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
\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
\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
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
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
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
\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
//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
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
//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
//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