clear all
close all

% Problem domain, in meters.
X = 1;
Y = 1;
D = 0.025;

% Problem domain, in cells.
NX = X/D;
NY = Y/D;

x_axis = linspace(0,X, NX + 1);
y_axis = linspace(0,Y, NY + 1);

V_last = zeros(NY+1,NX+1);
V_next = zeros(NY+1,NX+1);

% Write boundary conditions into V_last.
V_last(1,:) = -1;
V_last(end,:) = 1;
V_last(:,1) = linspace(-1,1, NY+1);
V_last(:,end) = linspace(-1,1, NY+1);
V_next = V_last;

% Cylinder description
cyl_x = 0.5;
cyl_y = 0.5;
cyl_r = 0.2;
    
ITERATIONS = 1000;

for i=1:ITERATIONS
  for iy=2:NY
    for ix=2:NX
      V_next(iy,ix) = 1/4* ( V_last(iy-1,ix) + V_last(iy+1,ix) + V_last(iy,ix-1) + V_last(iy,ix+1) );
    end
  end
  % Write zero voltages into interior of cylinder.
  for iy=1:NY+1
    for ix=1:NX+1
      if ( (x_axis(ix)-cyl_x)^2 + (y_axis(iy)-cyl_y)^2  < cyl_r^2)
        V_next(iy,ix) = 0;
      end
    end
  end
  % Iterate.
  V_last = V_next;
end

[x,y] = meshgrid(x_axis, y_axis);
[Ex,Ey] = gradient(-V_last,D,D);
contour(x,y,V_last,20);
rectangle('Position', [-cyl_r+cyl_x -cyl_r+cyl_y 2*cyl_r 2*cyl_r], 'Curvature', [1 1]);
hold on
quiver(x,y,Ex,Ey);
hold off
axis equal;
set(gcf,'Color',[1 1 1]);
xlabel ('x, meters');
ylabel ('y, meters');

title ('ECE 311 - HW 7, Problem 3.');

disp ('Conceptual Question: Because E is the -gradient of V, E must be orthogonal');
disp ('to the contours of equal potential. The inner conducting body is all at the');
disp ('same potential, V=0, including its circular boundary. So E must be perpendicular');
disp ('to the conductors surface. An alternate explanation is that the tangential component');
disp ('of is continuous across any media interface (including conductor/air, like here).');
disp ('Inside the conductor the E field is zero, so the tangential component of E outside');
disp ('the cylinder must be zero too (that is, E is perfectly normal to the surface).');