% -------------------------------------------------------------------
% Compute capacitance per unit length of a coax pair of rectangles.
% Coordinate system centered in middle.
% -------------------------------------------------------------------
function cap = mycapacitor(a, b, c, d, n, tol, rel)
% Arguments:
%    a   =  width of inner conductor
%    b   = height of inner conductor
%    c   =  width of outer conductor
%    d   = height of outer conductor
%    n   = number of points in the x-direction (horizontal)
%    tol = relative tolerance for capacitance
%    rel = relaxation parameter
%          (a good choice is 2-c/n, where c is about pi)
% Returns:
%    cap = capacitance per unit length [pF/m]

% Make grids
h  = 0.5*c/n;                   % Grid spacing
na = round(0.5*a/h);            % Number of segments on 'a'
x  = linspace(0, 0.5*c, n+1);   % Grid points along x-axis.
m  = round(0.5*d/h);            % Number of segments on 'd'
mb = round(0.5*b/h);            % Number of segments on 'b'
y  = linspace(0, 0.5*d, m+1);   % Grid points on y-axis.

mode = "GS"; % Select solver.

% Initialize potential and mask array.
f       = zeros(length(x), length(y));      % 2D-array with solution.
mask    = ones(length(x), length(y))*rel;   % 2D-array with relaxation.

% Initialize inside inner conductor.
for i = 1:na+1
    for j = 1:mb+1
        mask(i,j)   = 0;
        f(i,j)      = 1;
    end
end

switch(mode)
    case "J"
        %Jacobi iteration (DOESN'T CURRENTLY WORK! IGNORE!)
        oldcap = 0;
        for iter = 1:100000 % MAX number of iterations
            f = jacobi(f ,mask,n,m);      % Perform single Jacobi iteration.
            cap = gauss(n,m,h,f);   % Compute capacitance from potential.
            if (abs(cap-oldcap)/cap < tol)  % Stop if change in capacitance sufficiently small.
                break
            else
                oldcap = cap;               % Continue until converged.
            end
        end
        str = sprintf('Jacobi. Number of iterations until convergence : %4i', iter); 
        disp(str)
    case "GS"
        % Gauss-Seidel iteration
        oldcap = 0;
        for iter = 1:10000 % MAX number of iterations
            f = GS(f,mask,n,m);     % Perform single gauss-seidel iteration.
            cap = gauss(n,m,h,f);   % Compute capacitance from potential.
            if (abs(cap-oldcap)/cap < tol)  % Stop if change in capacitance sufficiently small.
                break
            else
                oldcap = cap;               % Continue until converged.
            end
        end
        str = sprintf('Gauss-Seidel. Number of iterations until convergence : %4i', iter); 
        disp(str)
end   
        %% Plot potential
        [X,Y] = meshgrid([-fliplr(x(2:end)) x], [-fliplr(y(2:end)) y]);
        F = [flip(f(:,2:end),2) f];  % Flip around lines of symmetry
        F = [flip(F(2:end,:),1); F];
        F = rot90(F); % Rotate matrix anticlockwise to align it with our plot.
        % surf(X,Y,F, 'FaceColor', 'interp'); % Alternate plot. Hides coax.
        % geometry plot.
        s=pcolor(X,Y,F);
        s.FaceColor = 'interp';
        s.EdgeColor = 'none';
        view(2);
        colorbar
        hold on
        rectangle('Position',[-a/2 -b/2 a b],'EdgeColor','r','LineWidth',2) % Draw coax. geometry.
        rectangle('Position',[-c/2 -d/2 c d],'EdgeColor','r','LineWidth',2)
        grid on
        axis square
        xlim([-c/2 c/2]);
        ylim([-d/2 d/2]);
        hold off
    


end

