From 9dfc460f7825213034351f413616137f8b299217 Mon Sep 17 00:00:00 2001
From: treecity <treecity@26120e32-c555-405d-b3e1-1f783fb42516>
Date: Fri, 22 Apr 2011 09:51:28 +0000
Subject: [PATCH] [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
---
 src/G00.m        |  8 +++++---
 src/build_A.cpp  | 18 ++++++++++++------
 src/g0.m         | 21 +++++++++++----------
 src/test_solve.m |  2 +-
 4 files changed, 29 insertions(+), 20 deletions(-)

diff --git a/src/G00.m b/src/G00.m
index 1112d9f..b725d43 100644
--- 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 3fb7295..5896472 100644
--- a/src/build_A.cpp
+++ b/src/build_A.cpp
@@ -133,7 +133,7 @@ for (int j=0;j<em;++j){
 		 if(rx=2){
 			 //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]);
 			 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]);
-			 printf("%d%d|%.2f\n",j,k,tmp);
+		//	 printf("%d%d|%.2f\n",j,k,tmp);
 			 A[(k*em)+j] = 1./(4*my_PI)*tmp;
 
 		 }else
@@ -185,13 +185,16 @@ for (int j=0;j<em;++j){
     }
     
     double g0(double p, double y, double x, double l){
+    //printf("%.1f | %.1f | %.1f | %.1f +",p,x,y,l);
+    
 		double sol = 0;
 
 		if(l!=0){
-            if(p==0.5)
+            if(p==0.5){
                 sol = (y-x)/2*sqrt((y-x)*(y-x)+
-                        l*l)+l*l/(2*arsinh((y-x)/abs(l)));
-            else if(p==0)
+                        l*l)+l*l/2*arsinh((y-x)/abs(l));
+      //          printf("%.2f |",sol);
+            }else if(p==0)
                 sol = y-x;
             else if(p==-0.5)
                 sol = arsinh((y-x)/abs(l));
@@ -208,11 +211,14 @@ for (int j=0;j<em;++j){
 			else
 				sol = (y-x)*pow(abs(y-x),2*p)/(2*p+1);
 		}
+        
         return sol;
 	}
     
     double G00(double p, double y1,double y2, double x1,double x2, double l) {
-		double sol =0;
+//	printf("%.1f | %.1f %.1f | %.1f %.1f | %.1f +",p,x1,x2,y1,y2,l);
+        
+        double sol =0;
 		if(p==-1.5){
 			if(l==0){
                 sol = -sqrt((y1-x1)*(y1-x1)+(y2-x2)*(y2-x2))
@@ -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)
     {
 
-		 printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);
+	//	 printf("%.1f %.1f | %.1f %.1f | %.1f %.1f %.1f",x1,x2,y1,y2,d1,d2,d3);
         double sol = (x1-y1-d1)*(x2-y2-d2);
     
         if(sol!=0)
diff --git a/src/g0.m b/src/g0.m
index be704eb..1571419 100644
--- 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 b40f953..422b323 100644
--- 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));
-- 
2.47.3