clf
close all
% Fixed geometry for coax
a = 2;
b = 5;
Vab = 10;
e0 = 8.85e-12;


% Potential between cylinders
V = @(r) (Vab.*(1./r - 1./b))./(1./a - 1./b)

% Electrostatic energy
We = (Vab^2*pi*e0)./(a*(1/a-1/b))

% Charge for a fixed potential, geometry and dielectric permittivity
Q = (Vab/a)*((2*pi*e0)/(1/a-1/b))

% Capacitance per length ( C/L )
C_L = @(b) (2*pi*e0)./(a*(1/a - 1./b))

% C/L values from COMSOL
% a = 2
% b = [5 7 9]
C_L_COMSOL_Q  = [9.1673E-11 7.6984E-11 7.0645E-11] 
C_L_COMSOL_We = [9.2722E-11 7.7887E-11 7.1531E-11]

figure
grid on
hold on
i = 1
for k = [5 7 9]
    plot(k, C_L(k), 'bo-')
    plot(k, C_L_COMSOL_Q(i), 'r*-')
    plot(k, C_L_COMSOL_We(i), 'ks-')
    i = i+1;
end

legend('Theoretical', 'COMSOL_Q', 'COMSOL_{We}');

%% Uppgift 11
% Observe the effects of mesh refinement on error.
C0 = 90.6145e-12; % C/L given to 6 correct digits.
% [Normal Fine Finer Extra-Fine Extremely-Fine] mesh

C_L_Q = [8.603265698925808E-11	8.658189215226473E-11	8.746681753844463E-11	8.845806943873758E-11	8.925183852867329E-11];
C_L_We= [9.089053703142881E-11	9.08093184419686E-11	9.07391837070597E-11	9.066997083867336E-11	9.063654599870883E-11];


err_Q = abs(C_L_Q - C0)
err_We = abs(C_L_We - C0)

figure
plot(err_Q)
hold on
plot(err_We)
legend('err_Q', 'err_{We}')
xlabel('Mesh refinement')
ylabel('Error [pf/m]')