[src] Funktionen laufen jetzt

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

diff --git a/src/F_par.m b/src/F_par.m
index 95b78423ac8711c146ebe2c1f59e569f1581991c..9a57744ded82266b977b2b225bd8410534ab6220 100644 (file)
--- a/src/F_par.m
+++ b/src/F_par.m
@@ -2,10 +2,27 @@ function sol = F_par(x1,x2,y1,y2,d1,d2,d3)
 
 %     xyd=x-y-delta(1:2);
     
-    sol = (x1-y1-d1)*(x2-y2-d2)*G00(-0.5,x1,x2,y1+d1,y2+d2,d3)...
-        -(x1-y1-d1)*g0(0.5,x1,y1+d1,sqrt((x2-y2-d2)^2+d3^2))...
-        -(x2-y2-d2)*g0(0.5,x2,y2+d2,sqrt((x1-y1-d1)^2+d3^2))...
-        +1/3*((x1-y1-d1)^2+(x2-y2-d2)^2+d3^2)^(1.5);
+    sol = (x1-y1-d1)*(x2-y2-d2);
+    
+    if(sol~=0)
+        sol = sol * G00(-0.5,x1,x2,y1+d1,y2+d2,d3);
+    end
+    
+    if((x1-y1-d1)~=0)
+        sol = sol - (x1-y1-d1)*g0(0.5,x1,y1+d1,sqrt((x2-y2-d2)^2+d3^2));
+    end
+    
+    if((x2-y2-d2)~=0)
+        sol = sol -(x2-y2-d2)*g0(0.5,x2,y2+d2,sqrt((x1-y1-d1)^2+d3^2));
+    end
+    
+    sol = sol + 1/3*((x1-y1-d1)^2+(x2-y2-d2)^2+d3^2)^(1.5);
+
+
+%     sol = (x1-y1-d1)*(x2-y2-d2)*G00(-0.5,x1,x2,y1+d1,y2+d2,d3)...
+%         -(x1-y1-d1)*g0(0.5,x1,y1+d1,sqrt((x2-y2-d2)^2+d3^2))...
+%         -(x2-y2-d2)*g0(0.5,x2,y2+d2,sqrt((x1-y1-d1)^2+d3^2))...
+%         +1/3*((x1-y1-d1)^2+(x2-y2-d2)^2+d3^2)^(1.5);
 
 %     sol = xyd(1)*xyd(2)*G00(-0.5,x,y+delta(1:2),delta(3))...
 %         -xyd(1)*g0(0.5,x(1),y(1)+delta(1),norm([xyd(2) delta(3)]))...

diff --git a/src/G00.m b/src/G00.m
index f181612c5ad9b59219d8087ce03b3759a9877652..1112d9f840b1f502d54438aa0e21473a8131e17b 100644 (file)
--- a/src/G00.m
+++ b/src/G00.m
@@ -3,16 +3,21 @@ function sol = G00(p,y1,y2,x1,x2,l)
 %   sol = nan;
     if(p==-1.5)
        if(l==0)
-           sol = -sqrt((y1-x1)^2+(y2-x2)^2)/((y1-x1)*(y2-x2));
+          sol = -sqrt((y1-x1)^2+(y2-x2)^2)/((y1-x1)*(y2-x2));
        else
           sol = sign((y1-x1)*(y2-x2))/(2*abs(l)) ...
               * acos(-2*(y1-x1)^2*(y2-x2)^2 ...
               /(((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) ...
-           + (y1-x1)*g0(p,y2,x2,sqrt((y1-x1)^2+l^2)) ...
-           + (y2-x2)*g0(p,y1,x1,sqrt((y2-x2)^2+l^2)))/(2*p+2);
+       sol = 2*p*l^2*G00(p-1,y1,y2,x1,x2,l); 
+       if((y1-x1)~=0)
+          sol = sol + (y1-x1)*g0(p,y2,x2,sqrt((y1-x1)^2+l^2));
+       end
+       if((y2-x2)~=0)
+          sol = sol + (y2-x2)*g0(p,y1,x1,sqrt((y2-x2)^2+l^2));
+       end
+       sol = sol / (2*p+2);
     else
         error 'no case for p defined'
     end

diff --git a/src/test_solve.m b/src/test_solve.m
index 5f2e0e326a10a984aff916d5c341541cc709c30e..99c5f2327dc2d18e4a4f028af5f83ae7cf49f5a8 100644 (file)
--- a/src/test_solve.m
+++ b/src/test_solve.m
@@ -3,7 +3,6 @@ N = size(elements);
 
 A = zeros(N);
 
-
 % untere schranke s t obere schranke k l
 intF = @(f,s1,s2,k1,k2,t1,t2,l1,l2) ...
     f(k1,k2,l1,l2)-f(k1,k2,l1,t2)-f(k1,k2,t1,l2)+f(k1,k2,t1,t2)...
@@ -11,18 +10,12 @@ intF = @(f,s1,s2,k1,k2,t1,t2,l1,l2) ...
     -f(s1,k2,l1,l2)+f(s1,k2,l1,t2)+f(s1,k2,t1,l2)-f(s1,l2,t1,t2)...
     +f(s1,s2,l1,l2)-f(s1,0,l1,t2)-f(s1,s2,t1,l2)+f(s1,s2,t1,t2);
 
-% intF(@(s1,s2,k1,k2,t1,t2,l1,l2) F_par(s1,s2,k1,k2,t1,t2,l1,l2,d(1),d(2),d(3))...
-%     ,s1,s2,k1,k2,t1,t2,l1,l2);
 for j = 1:N
-    
-    
     for k = 1:N
-        
-        
         ej = coordinates(elements(j,:)',:);
         ek = coordinates(elements(k,:)',:);
-        d = ek(1,:) - ej(1,:)
-             
+        d = ek(1,:) - ej(1,:);
+
         A(j,k) = 1/ (4*pi) *...
           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));
@@ -32,6 +25,4 @@ end
 
 A
 b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2))
-
-
-x = A\b;
\ No newline at end of file
+x = A\b
\ No newline at end of file