[src] Mex und Matlab machen das gleiche

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

diff --git a/src/G00.m b/src/G00.m
index 1112d9f840b1f502d54438aa0e21473a8131e17b..b725d433e1e7c070af210ef6e871c037b4001e7e 100644 (file)
--- a/src/G00.m
+++ b/src/G00.m
@@ -1,6 +1,6 @@
 function sol = G00(p,y1,y2,x1,x2,l)
-
-%   sol = nan;
+% fprintf('%.1f | %.1f %.1f | %.1f %.1f | %.1f +',p,x1,x2,y1,y2,l);
+   sol = 0;
     if(p==-1.5)
        if(l==0)
           sol = -sqrt((y1-x1)^2+(y2-x2)^2)/((y1-x1)*(y2-x2));
@@ -10,7 +10,9 @@ function sol = G00(p,y1,y2,x1,x2,l)
               /(((y1-x1)^2+l^2)*((y2-x2)^2+l^2))+1);
        end
     elseif(p==-0.5)
-       sol = 2*p*l^2*G00(p-1,y1,y2,x1,x2,l); 
+       if(l~=0)
+          sol = 2*p*l^2*G00(p-1,y1,y2,x1,x2,l);
+       end
        if((y1-x1)~=0)
           sol = sol + (y1-x1)*g0(p,y2,x2,sqrt((y1-x1)^2+l^2));
        end

diff --git a/src/build_A.cpp b/src/build_A.cpp
index 3fb7295b71fd001d0016de8c0d17fcde8e648258..58964728d19453af22925f3ce23d9573fe7c9d74 100644 (file)
--- a/src/build_A.cpp
+++ b/src/build_A.cpp
@@ -133,7 +133,7 @@ for (int j=0;j<em;++j){
                 if(rx=2){\r
                         //printf("%d%d %.1f %.1f | %.1f %.1f || %.1f %.1f | %.1f %.1f || %.1f %.1f %.1f\n",j,k,x1[0],x1[1],x2[0],x3[1],y1[0],y1[1],y2[0],y3[1],d[0],d[1],d[2]);\r
                         tmp = quadInt(F_par,x1[0],x1[1],x2[0],x3[1],y1[0],y1[1],y2[0],y3[1],d[0],d[1],d[2]);\r
-                        printf("%d%d|%.2f\n",j,k,tmp);\r
+               //       printf("%d%d|%.2f\n",j,k,tmp);\r
                         A[(k*em)+j] = 1./(4*my_PI)*tmp;\r
 \r
                 }else\r
@@ -185,13 +185,16 @@ for (int j=0;j<em;++j){
     }\r
     \r
     double g0(double p, double y, double x, double l){\r
+    //printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);\r
+    \r
                double sol = 0;\r
 \r
                if(l!=0){\r
-            if(p==0.5)\r
+            if(p==0.5){\r
                 sol = (y-x)/2*sqrt((y-x)*(y-x)+\r
-                        l*l)+l*l/(2*arsinh((y-x)/abs(l)));\r
-            else if(p==0)\r
+                        l*l)+l*l/2*arsinh((y-x)/abs(l));\r
+      //          printf("%.2f |",sol);\r
+            }else if(p==0)\r
                 sol = y-x;\r
             else if(p==-0.5)\r
                 sol = arsinh((y-x)/abs(l));\r
@@ -208,11 +211,14 @@ for (int j=0;j<em;++j){
                        else\r
                                sol = (y-x)*pow(abs(y-x),2*p)/(2*p+1);\r
                }\r
+        \r
         return sol;\r
        }\r
     \r
     double G00(double p, double y1,double y2, double x1,double x2, double l) {\r
-               double sol =0;\r
+//     printf("%.1f | %.1f %.1f | %.1f %.1f | %.1f +",p,x1,x2,y1,y2,l);\r
+        \r
+        double sol =0;\r
                if(p==-1.5){\r
                        if(l==0){\r
                 sol = -sqrt((y1-x1)*(y1-x1)+(y2-x2)*(y2-x2))\r
@@ -242,7 +248,7 @@ for (int j=0;j<em;++j){
     double F_par(double x1,double x2,double y1,double y2,double d1,double d2,double d3)\r
     {\r
 \r
-                printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
+       //       printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);\r
         double sol = (x1-y1-d1)*(x2-y2-d2);\r
     \r
         if(sol!=0)\r

diff --git a/src/g0.m b/src/g0.m
index be704ebda9f43ed5865786fc1522a4c678150cce..15714196bb4fae7a5d8a50d4917d297a22dfbbaa 100644 (file)
--- a/src/g0.m
+++ b/src/g0.m
@@ -1,21 +1,22 @@
-function sol = g0(p, y, x, lambda)
-
+function sol = g0(p, y, x, l)
+ %fprintf('%.1f | %.1f | %.1f | %.1f +',p,x,y,l);
  %sol = nan;
      yx = y-x;
-     if(lambda~=0)   
+     if(l~=0)   
         if(p==0.5)
-            sol = yx/2*((yx^2 + lambda^2)^0.5)+...
-                lambda^2/2*arsinh(yx/abs(lambda));
+            sol = yx/2*((yx^2 + l^2)^0.5)+...
+                l^2/2*arsinh(yx/abs(l));
+      %      fprintf('%.2f |',sol);
         elseif(p==0)
             sol = yx;
         elseif(p==-0.5)
-            sol = arsinh(yx/abs(lambda));
+            sol = arsinh(yx/abs(l));
         elseif(p==-1)
-            sol = atan(yx/abs(lambda));
+            sol = atan(yx/abs(l));
         elseif(p==-1.5)
-            sol =  yx/(lambda^2)*((yx^2 + lambda^2)^-0.5);
+            sol =  yx/(l^2)*((yx^2 + l^2)^-0.5);
         else
-            sol = yx*((yx^2 + lambda^2)^p)+2*p*lambda^2*g0(p-1,y,x,lambda);
+            sol = yx*((yx^2 + l^2)^p)+2*p*l^2*g0(p-1,y,x,l);
         end
      else
         if(p==-0.5)
@@ -24,7 +25,7 @@ function sol = g0(p, y, x, lambda)
             sol = (yx).*abs(yx).^(2*p)/(2*p+1);
         end
      end
-
+    
 end
 
 function sol = arsinh(x)    %inline

diff --git a/src/test_solve.m b/src/test_solve.m
index b40f953d0e1d9f711cfd994e7562906be75a61cc..422b32360b3c8efd89fca3e73b1b88deef7284bf 100644 (file)
--- a/src/test_solve.m
+++ b/src/test_solve.m
@@ -24,7 +24,7 @@ for j = 1:N
         tmp = intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3))...
             ,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2));
 
-        fprintf('%d%d|%.2f\n',j,k,tmp);
+%        fprintf('%d%d|%.2f\n',j,k,tmp);
         A1(j,k) = 1/ (4*pi) *tmp;
 %           intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3))...
 %             ,ej(1,1),ej(1,2),ej(2,1), ej(3,2),ek(1,1), ek(1,2),ek(2,1), ek(3,2));