% Set up matrix A and right-hand side s.

function[A, s] = setAs(a,b,c,d,n,m)
% 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)
%    m    = number of points in the y-direction (vertical)
% Returns:
%    A    = matrix on sparse storage format
%    s    = right-hand side on sparse storage format

hx = 0.5*c/n;           % Grid size in x-direction        
na = round(0.5*a/hx);   % Number of cells for half width of
                        % inner conductor 
hy = 0.5*d/m;           % Grid size in y-direction
m  = round(0.5*d/hy);   % Number of cells for half height of
                        % outer conductor
mb = round(0.5*b/hy);   % Number of cells for half height of
                        % inner conductor
p  = 1;                 % Potential on inner conductor

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%          The upper right corner is discretized
%
%               --------------------+
%                        c/2        |
%                                   |
%               ------------+       | d/2
%                   a/2     |       |
%                           | b/2   |
%                           |       |
%
%                     (Dimensions)
%
%
%             The nodes are numbered like this
%        y
%        ˆ
%        |
% (m-1)hy| (m-1)n+1 (m-1)n+2 (m-1)n+3  ...    mn
%        |     :        :        :            :
%   2hy  |   2n+1     2n+2     2n+3    ...    3n
%    hy  |    n+1      n+2      n+3    ...    2n
%    0   |     1        2        3      ...     n
%           --------------------------------------> x
%              0       hx      2 hx        (n-1)hx
%
%                   (Discretization)

N = n * m; % Total number of unknowns (grid-points)
cx = hx^-2;
cy = hy^-2; % Centimeters?

% Generate a matrix withN = m*n rows (-> nodes on the grid),
% and five columns, one for each nonzero diagonal of A.
% The first column gives contribution from nodes beneath.
% The second column gives contribution from nodes to the left.
% The third column gives self-contribution.
% The fourth column gives contribution from nodes to the right.
% The fifth column gives contribution from nodes above.
% The following lines assume some knowledge of MATLAB.  If you
% feel uncertain, insert the ’keyboard’ command.  This causes
% MATLAB to stop.  Then execute lines by ’dbstep’ and examine
% the result.

keyboard
C = repmat([cy cx -2*(cx+cy) cx cy], N, 1);

idx0R = n:n:N-n;          % Nodes with V=0 to the right
idxNB = na+2:n;           % Nodes with dV/dy = 0 beneath
idxNL = 1+n*(mb+1):n:N;   % Nodes with dV/dx = 0 to the left


idx1C = repmat((1:na+1)', 1, mb+1) + repmat((0:n:n*mb),na+1,1);   % 'x-index + n*(y-index-1)’ for all
idx1C = idx1C(:)';        % nodes on (or inside) the inner
                          % conductor where V = 1
                          % and convert to row vector
C(idx1C,[1 2 4 5]) = 0;
C(idx1C, 3) = 1;
C(idx0R, 4) = 0;
C(idxNB, 5) = 2*cy;
C(idxNL, 4) = 2*cx;
C(idxNL, 2) = 0;


end