simax = 16;
romax = 6;
K = simax*romax;
V = 100;

nbc = simax*3;
[X Y] = meshgrid(1:simax,1:romax);
mdl = (Y-1)*simax +X;
mdl(2:end-2,1:end) = 0;
bc = zeros(nbc,2);
bc(:,1) = mdl(mdl~=0);

for t = 0:simax-1  %Apply voltage boundary conditions
    t = 2+t*3;
    bc(t,2) = 100;
end

%Finite Difference Method contsnts in Cylindrical
delro = 0.25;
delsi = pi/8;
ronot = 0;

n=1;
for k = 1:K
    if k <= n*simax
        ll = n*simax-simax+2;   %lower limit
        ul = n*simax-1;         %upper limit
        lb = ll-1;              %lower boundary
        ub = n*simax;           %upper boundary
        ronot = n*delro;
        %inner node ~closer to origin
        Ao = (1-delro/2/ronot)/delro^2;
        %outer node ~towards edge
        Bo = (1+delro/2/ronot)/delro^2;
        %counterclockwise node and clockwise node
        Co = (delro/ronot/delsi)^2/delro^2;
        %original node ~the center node
        Do = -(2+2*(delro/ronot/delsi)^2)/delro^2;
        for k = ll:ul
        A(k,k) = Do;
        if k-simax > 0
            A(k,k-simax)=Ao;
        end
        if k+simax < K
            A(k,k+simax)=Bo;
        end
        if k+1 < K
            A(k,k+1) = Co;
        end
        if k-1 < 0
            A(k,k-1)=Co;
        end
        end
        %Boundary
        for k = ub
        Ao = 1-delro/2/ronot;
        Bo = 1+delro/2/ronot;
        Co = (delro/ronot/delsi)^2;
        Do = -(2+2*(delro/ronot/delsi)^2);
        A(k,k)=Do;
            if k-simax > 0
                A(k,k-simax)=Ao;
            end
            if k+simax < K
                A(k,k+simax)=Bo;
            end
            if k-simax+1 > 0
                A(k,k-simax+1) = Co;
            end
            if k-1 > 0
                A(k,k-1) = Co;
            end
        end
        for k = lb
        Ao = 1-delro/2/ronot;
        Bo = 1+delro/2/ronot;
        Co = (delro/ronot/delsi)^2;
        Do = -(2+2*(delro/ronot/delsi)^2);
        A(k,k)=Do;
            if k-simax > 0
                A(k,k-simax) = Ao;
            end
            if k+simax < K
                A(k,k+simax) = Bo;
            end
            if k+1 > 0
                A(k,k+1) = Co;
            end
            if k+simax-1 < K
                A(k,k+simax-1) = Co;
            end
        end 
    else
        n = n+1;
    end
end

%Boundary Conditions
for k = 1:simax
    A(k,k) = 0;
end
for k = 33:48
    A(k,k) = 0;
end
%B-Matrix
B(17:32) = -V;
[I J] = size(bc);

if J == 3 % bc in the form of (i, j, v)
    bc = [(bc(:,1)-1)*romax+bc(:,2) bc(:,3)]; % bc changed into (k, v)
end

for i = 1:I % boundary condition enforcement
    m = bc(i,1);
    A(m,:) = 0; A(m,m) = 1; B(m) = bc(i,2);
end
x = A\B;
theta(1:96,1) = 0;
rho(1:96,1) = 0;
for g = 1:6
    up = g*simax;
    low = g*simax+1-simax;
    Ro = g*delro;
    theta(low:up)=1:16;
    rho(low:up) = Ro;
end
theta = theta*pi/8;
%x(:,3)=x;
%x(:,1)=theta;
%x(:,2)=rho;
[X Y Z] = pol2cart(theta,rho,x);
Phi = reshape(x,simax,romax);
%figure;image(Phi'); axis image; set(gcf,'color','w');
figure;contour(Phi',romax*simax); grid on; set(gcf,'color','w');

%Calculating the impedance
v=3*10^8;
eo=8.85*10^-12;
qo = 0;
for r=0:simax
qo = qo + eo*((x(5*simax+r)-100)/2/delro)*delro;

end
Co = abs(qo/V);

Zo = 1/v/sqrt(Co^2)