[gGF] Parameter aufgesplittet

[src] test_solve hinzugefügt. funktioniert aber noch nicht

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

diff --git a/src/F_par.m b/src/F_par.m
index cb7d43bf05404fac74582416c86f8fe07ad33d05..a1794aa7b82141ecaac1de03f6d049b9a92942c9 100644 (file)
--- a/src/F_par.m
+++ b/src/F_par.m
@@ -1,11 +1,16 @@
-function sol = F_par(x, y, delta)
+function sol = F_par(x1,x2,y1,y2,d1,d2,d3)
 
-    xyd=x-y-delta(1:2);
+%     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,x1+d1,sqrt((x2-y2-d2)^2+d3^2))...
+        -(x2-y2-d2)*g0(0.5,x2,x2+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)]))...
-        -xyd(2)*g0(0.5,x(2),y(2)+delta(2),norm([xyd(1) delta(3)]))...
-        +1/3*(xyd(1)^2+xyd(2)^2+delta(3)^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)]))...
+%         -xyd(2)*g0(0.5,x(2),y(2)+delta(2),norm([xyd(1) delta(3)]))...
+%         +1/3*(xyd(1)^2+xyd(2)^2+delta(3)^2)^(1.5);
 
 end
 

diff --git a/src/G00.m b/src/G00.m
index dd0f40dd6eaa078510b3807a907646e4bc96fdc5..f181612c5ad9b59219d8087ce03b3759a9877652 100644 (file)
--- a/src/G00.m
+++ b/src/G00.m
@@ -1,18 +1,19 @@
-function sol = G00(p,y,x,lambda)
+function sol = G00(p,y1,y2,x1,x2,l)
 
-    sol = nan;
-    yx = y-x;
+%   sol = nan;
     if(p==-1.5)
-       if(lambda==0)
-           sol = -norm(yx)/prod(yx);
+       if(l==0)
+           sol = -sqrt((y1-x1)^2+(y2-x2)^2)/((y1-x1)*(y2-x2));
        else
-          sol = sign(prod(yx))/(2*abs(lambda))*...
-              acos(-2*prod(yx)^2/prod((yx).^2+lambda^2)+1);
+          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*lambda^2*G00(p-1,y,x,lambda)...
-           + sum( yx(1).*g0(p,y(2),x(2),norm([yx(1) lambda])))...
-           + sum( yx(2).*g0(p,y(1),x(1),norm([yx(2) lambda]))))/(2*p+2);
+       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);
+    else
+        error 'no case for p defined'
     end
-
 end

diff --git a/src/g0.m b/src/g0.m
index c8e163295e68b24fce7b79374c14233f58448f08..d576073ea8c9bbcd2f4a51741917fe1edeed8351 100644 (file)
--- a/src/g0.m
+++ b/src/g0.m
@@ -1,6 +1,6 @@
 function sol = g0(p, y, x, lambda)
 
- sol = nan;
+ %sol = nan;
      yx = y-x;
      if(lambda~=0)   
         if(p==0.5)
@@ -23,6 +23,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
new file mode 100644 (file)
index 0000000..4aa6aed
--- /dev/null
+++ b/src/test_solve.m
@@ -0,0 +1,37 @@
+
+N = size(elements);
+
+A = zeros(N);
+
+
+
+% intF = @(f,k1,k2,l1,l2) ...
+%     f(k1,k2,l1,l2)-f(k1,k2,l1,0)-f(k1,k2,0,l2)+f(k1,k2,0,0)...
+%     -f(k1,0,l1,l2)+f(k1,0,l1,0)+f(k1,0,0,l2)-f(k1,0,0,0)...
+%     -f(0,k2,l1,l2)+f(0,k2,l1,0)+f(0,k2,0,l2)-f(0,l2,0,0)...
+%     +f(0,0,l1,l2)-f(0,0,l1,0)-f(0,0,0,l2)+f(0,0,0,0);
+% 
+% intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3)),x(1),x(2),y(1),y(2));
+for j = 1:N
+    
+    
+    for k = j+1:N
+        
+        
+        ej = coordinates(elements(j,:)',:);
+        ek = coordinates(elements(k,:)',:);
+        d = ek(1,:) - ej(1,:)
+        
+        x = [sum(ej(2,:)-ej(1,:)) sum(ej(3,:)-ej(1,:))]
+        y = [sum(ek(2,:)-ek(1,:)) sum(ek(3,:)-ek(1,:))]
+        
+        A(j,k) = 1/ (4*pi) *... intF(@(x1,x2,y1,y2) F_par(x1,x2,y1,y2,d(1),d(2),d(3)),x(1),x(2),y(1),y(2));
+        F_par(x(1),x(2),y(1),y(2),d(1),d(2),d(3));
+    end
+end
+
+A
+b = sqrt(sum(quad_norm(coordinates,elements,'w').^2,2))
+
+
+x = A\b;
\ No newline at end of file