Changes and additions to CurrentCode.m and GetCurrents.m

CondCode/CurrentCode.m

@@ -3,7 +3,7 @@
 clearvars -GLOBAL
 close all
 % set(0,'DefaultFigureWindowStyle','docked')
-global C fc
+global im fig fc map C 
 
 addpath ../geom2d/geom2d
 
@@ -22,42 +22,38 @@
 Acond = 1;
 Bcond = 10;
 
-SimType = 'e';
+SimType = 'c';
 
 Max = 5;
 ncircs = 20;
 doPlot = 1; % set to 1 if you would like to draw the plot
 
-if SimType == 'c'
+if SimType == 'c' || SimType == 'r'
     n = 20;
-    nSims = 1;
+    nSims = 5;
     Res = zeros(n, nSims);
-
     for k = 1:n
         Max(k) = k;
         fc = k;
         for i = 1:nSims
             % V = 1 --> R = 1/I
             Res(k,i) = 1/GetCurrents(ncircs,Max(k),nx,ny,...
-                Acond,Bcond,doPlot,SimType);
+                Acond,Bcond,doPlot,SimType,i);
         end
     end
 else %if SimType == 'e'
-
     n = 20;
-    nSims = 30;
+    nSims = 1;
     Res = zeros(n, nSims);
-
     for k = 1:n
         Max(k) = k * 3;
         fc = k;
         for i = 1:nSims
             % V = 1 --> R = 1/I
             Res(k,i) = 1 / GetCurrents(Max(k), 10, nx, ny,...
-                Acond, Bcond, doPlot, SimType);
+                Acond, Bcond, doPlot, SimType,i);
         end
     end
-
 end
 
 if doPlot
@@ -88,4 +84,4 @@
 
 subplot(4, 1, 4), plot(Max, StdCurr./ AveCurr * 100, 'r');
 xlabel('Num');
-ylabel('Variance (% of mean)');
+ylabel('Variance (% of mean)');

CondCode/GetCurrents.m

@@ -1,5 +1,5 @@
-function [ Curr ] = GetCurrents(ngeo, Max, nx, ny, Acond, Bcond, plot, SimType)
-global im fig fc map
+function [ Curr ] = GetCurrents(ngeo, Max, nx, ny, Acond, Bcond, plot, SimType,FigNum)
+global im fig fc map C
 
 if SimType == 'c'
     pgeo = rand(ngeo, 2);
@@ -8,16 +8,36 @@
     rgeo = rand(ngeo,1)*Max;
     rcircs2 = rgeo.^2;
 
-    cMap = zeros(nx,ny);
+    cMap = zeros(ny,nx);
 
     for i = 1:nx
         for j = 1:ny
-            cMap(i,j) = Acond;
+            cMap(j,i) = Acond;
             for p = 1:ngeo
                 dx = (pgeo(p, 1) - i);
                 dy = (pgeo(p, 2) - j);
-                if dx * dx + dy * dy < rcircs2(p)
-                    cMap(i, j) = Bcond;
+                if (dx^2 + dy^2) < rcircs2(p)
+                    cMap(j, i) = Bcond;
+                end
+            end
+        end
+    end
+elseif SimType == 'r'
+    pgeo = rand(ngeo,2);
+    pgeo(:,1) = pgeo(:,1)*nx;
+    pgeo(:,2) = pgeo(:,2)*ny;
+    dgeo = ones(ngeo,1)*Max;
+
+    cMap = zeros(ny,nx);
+
+    for i = 1:nx
+        for j = 1:ny
+            cMap(j,i) = Acond;
+            for p = 1:ngeo
+                dx = abs(pgeo(p, 1) - i);
+                dy = abs(pgeo(p, 2) - j);
+                if (dx < dgeo(p)) && (dy < dgeo(p))
+                    cMap(j, i) = Bcond;
                 end
             end
         end
@@ -29,23 +49,22 @@
     pgeo(:, 3) = pgeo(:, 2) * 180;
     rgeo = ones(ngeo, 1) * Max;
 
-    cMap = zeros(nx, ny);
+    cMap = zeros(ny, nx);
 
     for i = 1:nx
         for j = 1:ny
-            cMap(i, j) = Acond;
+            cMap(j,i) = Acond;
             for p = 1:ngeo
-                point = [i, j];
+                point = [j,i];
                 el = [pgeo(p, 1), pgeo(p, 2),rgeo(p), rgeo(p) / 4, pgeo(p, 3)];
                 if isPointInEllipse(point, el)
-                    cMap(i, j)= Bcond;
+                    cMap(j,i)= Bcond;
                 end
             end
         end
     end
 end
 
-
 G = sparse(nx*ny);
 B = zeros(1,nx*ny);
 
@@ -65,9 +84,9 @@
             nxp = j + (i) * ny;
             nyp = j + 1 + (i - 1) * ny;
 
-            rxm = (cMap(i, j) + cMap(i - 1, j)) / 2.0;
-            rxp = (cMap(i, j) + cMap(i + 1, j)) / 2.0;
-            ryp = (cMap(i, j) + cMap(i, j + 1)) / 2.0;
+            rxm = (cMap(j,i) + cMap(j,i-1)) / 2.0;
+            rxp = (cMap(j,i) + cMap(j,i+1)) / 2.0;
+            ryp = (cMap(j,i) + cMap(j+1,i)) / 2.0;
 
             G(n, n) = -(rxm+rxp+ryp);
             G(n, nxm) = rxm;
@@ -79,9 +98,9 @@
             nxp = j + (i) * ny;
             nym = j - 1 + (i - 1) * ny;
 
-            rxm = (cMap(i, j) + cMap(i - 1, j)) / 2.0;
-            rxp = (cMap(i, j) + cMap(i + 1, j)) / 2.0;
-            rym = (cMap(i, j) + cMap(i, j - 1)) / 2.0;
+            rxm = (cMap(j,i) + cMap(j,i-1)) / 2.0;
+            rxp = (cMap(j,i) + cMap(j,i+1)) / 2.0;
+            rym = (cMap(j,i) + cMap(j-1,i)) / 2.0;
 
             G(n, n) = -(rxm + rxp + rym);
             G(n, nxm) = rxm;
@@ -93,10 +112,10 @@
             nym = j-1 + (i-1)*ny;
             nyp = j+1 + (i-1)*ny;
 
-            rxm = (cMap(i,j) + cMap(i-1,j))/2.0;
-            rxp = (cMap(i,j) + cMap(i+1,j))/2.0;
-            rym = (cMap(i,j) + cMap(i,j-1))/2.0;
-            ryp = (cMap(i,j) + cMap(i,j+1))/2.0;
+            rxm = (cMap(j,i) + cMap(j,i-1))/2.0;
+            rxp = (cMap(j,i) + cMap(j,i+1))/2.0;
+            rym = (cMap(j,i) + cMap(j-1,i))/2.0;
+            ryp = (cMap(j,i) + cMap(j+1,i))/2.0;
 
             G(n,n) = -(rxm+rxp+rym+ryp);
             G(n,nxm) = rxm;
@@ -110,33 +129,9 @@
 
 V = G\B';
 
-Vmap = zeros(nx,ny);
-for i = 1:nx
-    for j = 1:ny
-        n = j + (i - 1) * ny;
+Vmap = reshape(V,[ny,nx]);
 
-        Vmap(i, j) = V(n);
-    end
-end
-
-for i = 1:nx
-    for j = 1:ny
-        if i == 1
-            Ex(i, j) = (Vmap(i + 1, j) - Vmap(i, j));
-        elseif i == nx
-            Ex(i, j) = (Vmap(i, j) - Vmap(i - 1, j));
-        else
-            Ex(i, j) = (Vmap(i + 1, j) - Vmap(i - 1, j)) * 0.5;
-        end
-        if j == 1
-            Ey(i, j) = (Vmap(i, j + 1) - Vmap(i, j));
-        elseif j == ny
-            Ey(i, j) = (Vmap(i, j) - Vmap(i, j - 1));
-        else
-            Ey(i, j) = (Vmap(i, j + 1) - Vmap(i, j - 1)) * 0.5;
-        end
-    end
-end
+[Ex,Ey] = gradient(Vmap,1,1);
 
 Ex = -Ex;
 Ey = -Ey;
@@ -148,18 +143,18 @@
     if fc == 1
         fig = figure();
     end
-
-    subplot(2, 2, 1), H = surf(cMap');
+    figure(FigNum)
+    subplot(2, 2, 1), H = surf(cMap);
     set(H, 'linestyle', 'none');
     view(0, 90)
-    subplot(2, 2, 2), H = surf(Vmap');
+    subplot(2, 2, 2), H = surf(Vmap);
     set(H, 'linestyle', 'none');
     view(0, 90)
-    subplot(2, 2, 3), quiver(Ex', Ey');
+    subplot(2, 2, 3), quiver(Ex, Ey);
     axis([0 nx 0 ny]);
-    subplot(2, 2, 4), quiver(eFlowx', eFlowy');
+    subplot(2, 2, 4), quiver(eFlowx, eFlowy);
     axis([0 nx 0 ny]);
-
+    
     if fc == 1
         f = getframe(fig);
         [im, map] = rgb2ind(f.cdata, 256, 'nodither');
@@ -170,8 +165,8 @@
     end
 end
 
-C0 = sum(eFlowx(1, :));
-Cnx = sum(eFlowx(nx, :));
+C0 = sum(eFlowx(:,1));
+Cnx = sum(eFlowx(:,ny));
 Curr = (C0 + Cnx) * 0.5;
 
-end
+end