From 440e08e033dbf44482a29d328fe3f87455c4bc03 Mon Sep 17 00:00:00 2001 From: treecity Date: Thu, 1 Sep 2011 13:39:42 +0000 Subject: [PATCH] [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 --- src/{mex_FparY1X2.cpp => mex_FparQY2X2.cpp} | 6 +-- src/mex_FparY1Y2_2.cpp | 44 +++++++++++++++++++++ src/mex_build_As2.cpp | 2 +- src/mexit.m | 2 + src/slpRectangle.cpp | 16 ++++++++ src/slpRectangle.hpp | 11 +++++- src/surfQuad.m | 4 +- src/test_Fpar.m | 32 ++++++++------- src/test_G00.m | 12 +++--- src/test_g0.m | 7 ++-- 10 files changed, 107 insertions(+), 29 deletions(-) rename src/{mex_FparY1X2.cpp => mex_FparQY2X2.cpp} (73%) create mode 100755 src/mex_FparY1Y2_2.cpp diff --git a/src/mex_FparY1X2.cpp b/src/mex_FparQY2X2.cpp similarity index 73% rename from src/mex_FparY1X2.cpp rename to src/mex_FparQY2X2.cpp index 81453f1..883158a 100755 --- a/src/mex_FparY1X2.cpp +++ b/src/mex_FparQY2X2.cpp @@ -21,7 +21,7 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL); - plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL); +// plhs[1] = mxCreateDoubleMatrix(1, 1, mxREAL); // plhs[0] = mxCreateDoubleMatrix(1, 2, mxREAL); @@ -30,8 +30,8 @@ void mexFunction(int nlhs, mxArray *plhs[], int nrhs, const mxArray *prhs[]) { //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])); + *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]))); diff --git a/src/mex_FparY1Y2_2.cpp b/src/mex_FparY1Y2_2.cpp new file mode 100755 index 0000000..7c76461 --- /dev/null +++ b/src/mex_FparY1Y2_2.cpp @@ -0,0 +1,44 @@ +#include +#include +#include +#include "mex.h" +#include + +#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 c602509..f254131 100644 --- 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{ tmp2 - = apply0Int2(F_parY1Y2_2, fabs(xda), fabs(xdb), + = apply0Int4(F_parY1Y2_2, fabs(xda), fabs(xdb), fabs(yda), fabs(ydb), d[rxa], d[rxb], d[rx]); // } diff --git a/src/mexit.m b/src/mexit.m index f868508..d103b15 100644 --- 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 b92c6a4..f4a0e45 100644 --- 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; } +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) { diff --git a/src/slpRectangle.hpp b/src/slpRectangle.hpp index 5b3a715..5133e06 100644 --- 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 36b2289..b5967c8 100644 --- a/src/surfQuad.m +++ b/src/surfQuad.m @@ -12,7 +12,9 @@ function in = surfQuad(f,a,b,c,d,v) in = 0; for i = 1:length(yn) - in = in + yw(i) * xw*f(xn,yn(i)); + for j = 1:length(xn) + in = in + yw(i) * xw(j)*f(xn(j),yn(i)); + end end end \ No newline at end of file diff --git a/src/test_Fpar.m b/src/test_Fpar.m index e1f70a7..920441c 100644 --- 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 66e7f7a..097f8c4 100644 --- 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 dbd75dc..4163f95 100644 --- 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); -- 2.47.3