% Capacitance_problem_2D_pg20
clear
EM_constants
close all
format compact
format long

%% 2D capacitance problem
% Consider coax. transmission line of cross-sectional geometry below.
% Inner conducter, outer waveguide.

%
% +-------------------------+
% |                         |
% |                         |
% |       +---------+       |
% |       |         |       |
% |       |         |b      |d
% |       |         |       |
% |       +---------+       |
% |            a            |
% |                         |
% +-------------------------+
%              c

%% Iterative solutions to the discretized Laplace equation.
% Start with arbitrary initial guess for solution f_i,j at all internal grid-points. (eg.  f_i,j = 0), except boundaries which are kept fixed.
% Iteratively update values of all grid points until a converged solution is reached.

%         *
%         fi,j+1
% 
% 
% *       *       *
% fi-1,j  fi,j    fi+1,j
% 
% 
%         *
%         fi,j-1

%% Jacobi iteration
% At every grid-point, the potential is the average of the potential at the four nearest neighbours.
% Slow convergence.
% f_{i,j}^{\left(n+1\right)} =\frac{1}{4}\left(f_{i-1,j}^{\left(n\right)} +f_{i+1,j}^{\left(n\right)} +f_{i,j-1}^{\left(n\right)} +f_{i,j+1}^{\left(n\right)} \right)
%% Gauss-Seidel iteration
% Modification of the Jacobi iteration: "old" values of f are immediately overwritten by new ones, as soon as they are computed.
% Example below assumes i,j updated in that order.
% f_{i,j}^{\left(n+1\right)} =\frac{1}{4}\left(f_{i-1,j}^{\left(n+1\right)} +f_{i+1,j}^{\left(n\right)} +f_{i,j-1}^{\left(n+1\right)} +f_{i,j+1}^{\left(n\right)} \right)
% A further improvement is "overrelaxation": the correction between iterations is multiplied by a relaxation parameter R.
% R>1 greatly improves convergence, but must be kept below 2 for stability.
% f_{i,j}^{\left(n+1\right)} =f_{i,j}^{\left(n\right)} +R\left(\frac{1}{4}\left(f_{i-1,j}^{\left(n+1\right)} +f_{i+1,j}^{\left(n\right)} +f_{i,j-1}^{\left(n+1\right)} +f_{i,j+1}^{\left(n\right)} \right)-f_{i,j}^{\left(n\right)} \right)
%% The problem at hand
% We generate a grid such that the conducting boundaries fall on grid points. We exploit symmetry such that f_{i+n,j} =f_{i-n,j} and f_{i,j+n} =f_{i,j-n}, that is, we only compute the upper right quadrant.
% Introduce boundary conditions by setting inner conductor f=V=1 and outer conductor f=0.
% Set up array to identify whether grid-point is inside region where Laplace's equation governs potential.
% Iterate with Gauss-Seidel scheme over internal points to solve for potential.
% Capacitance per unit length is C=\frac{Q}{V}=Q. Charge Q=\varepsilon_0 \oint \mathit{\mathbf{E}}\cdot \hat{\mathit{\mathbf{n}}} \;\textrm{dl}=-\varepsilon_0 \oint \frac{\partial \phi }{\partial n}\textrm{dl} (Gauss' Law)
% Iterate until capacitance converges (ie. change within tolerance). Stop iterating.
% Refine grid and extrapolate to zero cell size.

% Coax. dimensions (see above)
a =10*10^(-3);
b =10*10^(-3);
c =20*10^(-3);
d =20*10^(-3);
% Grid resolution (points per half grid)
n =20;

cap = capacitor(a,b,c,d,n,1e-6,1.9)

%

%% Compute for various resolutions, then extrapolate to zero cell size.
i = 1;
na = 10:10:50;
for n = na
    capa(i) = capacitor(a,b,c,d,n,1e-9,1.9);
    i = i + 1;
end
alpha = 1.5; % Given in problem!
h = 0.5.*c./na;
plot(h.^1, capa, '-+')
xlabel('Grid resolution h^{1.5}');ylabel('pF/m');
title('Capacitance per unit length per grid resolution');

A = [ ones(length(na),1)  (h.^alpha)'];
b = capa';

xVec = A\b
hcont = linspace(0, 0.5.*c./10, 1000);
hold on
%plot(hcont.^alpha, xVec(1) + xVec(2)*hcont.^alpha, 'r-')
legend('Computed', 'Extrapolated')

%I0 = xVec(1)

p = polyfit(h,capa,2)
y1 = polyval(p,hcont)
%plot(hcont, y1)
plot(hcont.^1, y1, 'r-')

I0 = y1(1)