[src] neuer Input deshalb Zwischenspeicherung

[src] funzt schon so einiges

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

diff --git a/src/mex_FparQY2X2.cpp b/src/mex_FparQY2X2.cpp
new file mode 100755 (executable)
index 0000000..883158a
--- /dev/null
+++ b/src/mex_FparQY2X2.cpp
@@ -0,0 +1,44 @@
+#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_FparY1X2.cpp b/src/mex_FparY1X2.cpp
deleted file mode 100755 (executable)
index 81453f1..0000000
--- a/src/mex_FparY1X2.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_parX1Y2(*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
new file mode 100755 (executable)
index 0000000..7c76461
--- /dev/null
+++ b/src/mex_FparY1Y2_2.cpp
@@ -0,0 +1,44 @@
+#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_build_As2.cpp b/src/mex_build_As2.cpp
index c60250927060d90ee1056e4c8f380485c5d4ebba..f2541312c6c4d639912c0ae685797748ce4f1baf 100644 (file)
--- a/src/mex_build_As2.cpp
+++ b/src/mex_build_As2.cpp
@@ -261,7 +261,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) {
                        //              }else{\r
 \r
                                        tmp2\r
-                                                       = apply0Int2(F_parY1Y2_2, fabs(xda), fabs(xdb),\r
+                                                       = apply0Int4(F_parY1Y2_2, fabs(xda), fabs(xdb),\r
                                                                        fabs(yda), fabs(ydb), d[rxa],\r
                                                                        d[rxb], d[rx]);\r
                        //              }\r

diff --git a/src/mexit.m b/src/mexit.m
index f868508cb42165a61125148f8ae48db23eb6b7c1..d103b15e0daecd1e6f43c19d1f6e97971536e47c 100644 (file)
--- a/src/mexit.m
+++ b/src/mexit.m
@@ -4,6 +4,8 @@ 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 b92c6a40093524eaa3f984d27ea15ae4c28ff016..f4a0e45b0e9bdf39a3611093aa92c431195171cc 100644 (file)
--- a/src/slpRectangle.cpp
+++ b/src/slpRectangle.cpp
@@ -208,6 +208,22 @@ double F_parY1Y2_2(double x1, double x2, double y1, double y2, double d1, double
        return sol;\r
 }\r
 \r
+double F_parQY2X2(double x1, double x2, double y1, double y2, double d1, double d2,\r
+               double d3) {\r
+\r
+       //       printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
+\r
+       double hlp = ((x1 - y1 - d1) * (x1 - y1 - d1) + (x2 - y2 - d2) * (x2 - y2\r
+                       - d2) + d3 * d3);\r
+\r
+       double sol = sqrt(hlp);\r
+\r
+       if ((x2 - y2 - d2) != 0)\r
+               sol += (x2 - y2 - d2) * log(sqrt(hlp)-(x2 - y2 - d2));\r
+\r
+       return sol;\r
+}\r
+\r
 double F_ort(double x1, double x2, double y2, double y3, double d1, double d2,\r
                double d3) {\r
 \r

diff --git a/src/slpRectangle.hpp b/src/slpRectangle.hpp
index 5b3a715ec35acf6f12ba7b65ac6a9ad6326c8b64..5133e06f2bf6aaae6260ae080c1d5740ab868572 100644 (file)
--- a/src/slpRectangle.hpp
+++ b/src/slpRectangle.hpp
@@ -31,14 +31,15 @@ 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);
 
-/*
 // 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 = quadInt((F_par/F_ort),a,b,c,d,r,t,u,v,d1,d2,d3);
 double applyInt4(
@@ -56,4 +57,10 @@ double apply0Int2(
                double, double, double, double, double, double, double);
 
 
+double calcParIntA(double, double, double, double, double, double, double);
+
+
+
+
+
 #endif

diff --git a/src/surfQuad.m b/src/surfQuad.m
index 36b22893016b223642de19cf68c47f89b7835d16..b5967c87cc0c97b2024b96403453bbef3a033b2a 100644 (file)
--- a/src/surfQuad.m
+++ b/src/surfQuad.m
@@ -12,7 +12,9 @@ function in = surfQuad(f,a,b,c,d,v)
     \r
     in = 0;\r
     for i = 1:length(yn)\r
-        in = in + yw(i) * xw*f(xn,yn(i));\r
+      for j = 1:length(xn)\r
+        in = in + yw(i) * xw(j)*f(xn(j),yn(i));\r
+      end\r
     end\r
 \r
 end
\ No newline at end of file

diff --git a/src/test_Fpar.m b/src/test_Fpar.m
index e1f70a7ae7575105dfbca182c3832852ed1498f9..920441c9e7382328de7d5be93894be001f108fb7 100644 (file)
--- a/src/test_Fpar.m
+++ b/src/test_Fpar.m
@@ -1,5 +1,6 @@
 function [ret] = test_Fpar(a,b,c,d,r,t,u,v,l1,l2,l3)
-
+%
+% [ret] = test_Fpar(a,b,c,d,r,t,u,v,l1,l2,l3)
 
     quad_sol = surfDoubleQuad(@(x1,x2,y1,y2)1/sqrt((x1-y1-l1).^2+(x2-y2-l2).^2+l3.^2)...
                 ,a,b,c,d,r,t,u,v,5);
@@ -11,23 +12,26 @@ function [ret] = test_Fpar(a,b,c,d,r,t,u,v,l1,l2,l3)
         +f(a,c,t,v)-f(a,c,t,u)-f(a,c,r,v)+f(a,c,r,u);
       
       
-    intF2 = @(f,a,b,c,d,r,t,u,v)...
-        f(b,d,t,v)-f(a,d,t,v)-f(b,c,t,v)+f(a,c,t,v);
-      
       
-     [dm sol] = mex_Fpar(b,d,t,v,l1,l2,l3)
-
-     [dm sol1] =  mex_Fpar(b,d,0,v,l1,l2,l3)
-     sol = sol - sol1;
-     [dm sol1] =  mex_Fpar(b,d,t,0,l1,l2,l3)
-     sol = sol - sol1;
-     [dm sol1] =  mex_Fpar(b,d,0,0,l1,l2,l3)
-    q_sol = sol + sol1;
       
+      [ mex_Fpar(b,d,t,v,l1,l2,l3) mex_FparY1Y2_2(b,d,t,v,l1,l2,l3);...
+        mex_Fpar(b,d,t,u,l1,l2,l3) mex_FparY1Y2_2(b,d,t,u,l1,l2,l3);...
+        mex_Fpar(b,d,r,v,l1,l2,l3) mex_FparY1Y2_2(b,d,r,v,l1,l2,l3)];
       
+      [i j] = mex_G00(-0.5, b, d, t + l1, v + l2, l3);
             
     sol = intF(@(x1,x2,y1,y2)mex_Fpar(x1,x2,y1,y2,l1,l2,l3),a,b,c,d,r,t,u,v);
-
-    ret = [quad_sol sol q_sol sol/q_sol];
+    
+   sol2 = intF(@(x1,x2,y1,y2)mex_FparY1Y2_2(x1,x2,y1,y2,l1,l2,l3),a,b,c,d,r,t,u,v);
+   
+    int2 = @(x1,y1,c,d,u,v)...
+      mex_FparQY2X2(x1,d,y1,v,l1,l2,l3)...
+      -mex_FparQY2X2(x1,c,y1,v,l1,l2,l3)...
+      -mex_FparQY2X2(x1,d,y1,u,l1,l2,l3)...
+      +mex_FparQY2X2(x1,c,y1,u,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];
 end
 

diff --git a/src/test_G00.m b/src/test_G00.m
index 66e7f7af4a382be85e4a3ccd5691010951c0c990..097f8c4ff2eabb8a02043773a75ee91a998ea19d 100644 (file)
--- a/src/test_G00.m
+++ b/src/test_G00.m
@@ -1,16 +1,18 @@
 function [ret] = test_G00(p, x1, x2,l,a,b,c,d)
-
+%
+% [ret] = test_G00(p, x1, x2,l,a,b,c,d)
 
     quad_sol = surfQuad(@(y1,y2)((x1-y1).^2+(x2-y2).^2+l.^2).^p,a,b,c,d,4);
 
 
-    [my(1) ] = mex_G00(p,a,c,x1,x2,l);
-    [my(2) mq_sol] = mex_G00(p,b,d,x1,x2,l);
-    [my(3) ] = mex_G00(p,b,c,x1,x2,l);
-    [my(4) ] = mex_G00(p,a,d,x1,x2,l);
+    [my(1) mq(1)] = mex_G00(p,a,c,x1,x2,l)
+    [my(2) mq(2)] = mex_G00(p,b,d,x1,x2,l)
+    [my(3) mq(3)] = mex_G00(p,b,c,x1,x2,l)
+    [my(4) mq(4)] = mex_G00(p,a,d,x1,x2,l)
 
 
     my_sol = my(1)+my(2)-my(3)-my(4);
+    mq_sol = mq(1)+mq(2)-mq(3)-mq(4);
 
     ret = [quad_sol my_sol mq_sol];
 end

diff --git a/src/test_g0.m b/src/test_g0.m
index dbd75dcfa00c73969f87edac4842109a2c282989..4163f9590b524776f6fc9dbdb1a3b8f5f5325811 100644 (file)
--- a/src/test_g0.m
+++ b/src/test_g0.m
@@ -6,17 +6,18 @@ function [ret] = test_g0(p, x,l,a,b)
 f = @(y)((x-y).^2+l.^2).^p;
 quad_sol = we*f(no);
 
-cumsum(f(no))'
+% cumsum(f(no))'
 
 % X = 0:0.2:5;
 % Y = ((x-X).^2+l.^2).^p;
 % 
 % plot(X,Y)
 
-my(1) = mex_g0(p,a,x,l);
-[my(2) mq_sol]= mex_g0(p,b,x,l);
+[my(1) mq(1)] = mex_g0(p,a,x,l);
+[my(2) mq(2)]= mex_g0(p,b,x,l);
 
 my_sol = my(2) - my(1);
+mq_sol = mq(2) - mq(1);