%% Lab B

close all
format shorteng

%% My parameters

p = [9 3 0 8 1 4 7 9 5 9];

p(4)
p(6)
p(9)
p(10)

w = 44e-3;
h = 30e-3;
x1 = 12e-3;
y1 = 9e-3;
r1 = 4e-3;
d1 = 2e-3;

x2 = (28 + p(4)/2)*10^-3
y2 = (13 + p(6)/2)*10^-3
r2 = (5 + p(9)/10)*10^-3
d2 = (3-p(9)/10)*10^-3

er1 = 3+p(10)/2 % Relative permittivity 
er2 = 8-p(10)/2

%
%% Capacitance coefficients per unit length (for air-filled transmission line)
% From COMSOL
% These are actually c/l, ie. capacitance coeff per unit length.
% Derived in COMSOL using: \frac{q_i }{l_z }=\frac{c_i }{l_z }=\int_L \hat{\;n} \cdot \overrightarrow{\;D} \;\textrm{dl} on the conductor/dielectric interface, for all three conductors (including outer shield).
% See lab report for details.

cl =[6.975614733066389E-11	-1.1335964260201521E-11	-5.842018307083462E-11;
   -1.1339293016111567E-11	 7.370702897590696E-11	-6.23677359601164E-11;
   -5.846036581240197E-11	-6.239672443221213E-11	 1.2085709024422877E-10]

% We use the expressions derived analytically (with outer shield set at V=0). \frac{Q}{L} is a free variable that we can set so that the scaling is reasonable.

QoverL = 10^-10;
v1 = QoverL*(cl(1,2)+cl(2,2))/(cl(1,1)*cl(2,2)-cl(1,2)*cl(2,1))
v2 = -QoverL*(cl(1,1)+cl(2,1))/(cl(1,1)*cl(2,2)-cl(1,2)*cl(2,1))

% These potential values give use back our set charge, QoverL, in COMSOL. This ensures that we have done things correctly.
% With this we can compute potential difference, and then capacitance per unit length, based on the relationship C=\frac{Q}{U}

CL = (cl(1,1)*cl(2,2) - cl(1,2)*cl(2,1))/(cl(1,1) + cl(1,2) + cl(2,1) + cl(2,2))

%
%% Capacitance coefficients per unit length (for partially water-filled transmission line)
% Same procedure as before.

cl =[
5.551685082490112E-10	-9.840296233661254E-12	-5.453282120162526E-10;
-9.843674915312611E-12	2.1873870627708876E-10	-2.0889503136234942E-10;
-5.457317602019538E-10	-2.0898944429498694E-10	7.547212044825389E-10;
]
QoverL = 10^-10;
v1 = QoverL*(cl(1,2)+cl(2,2))/(cl(1,1)*cl(2,2)-cl(1,2)*cl(2,1))
v2 = -QoverL*(cl(1,1)+cl(2,1))/(cl(1,1)*cl(2,2)-cl(1,2)*cl(2,1))
CL = (cl(1,1)*cl(2,2) - cl(1,2)*cl(2,1))/(cl(1,1) + cl(1,2) + cl(2,1) + cl(2,2))

%% MISC functions
% Magnitude and normal-relative angle of E-field on interface.
% Args:
% E1: magnitude of E-field in region 1.
% eps1: relative permittivity in region 1
% eps2: relative permittivity in region 2
% alpha1: normal-relative incident angle in region 1 in degrees.
% Returns:
% E2: mag of E-field in region 2.
% alpha2: normal-relative angle of E-field in region 2.

finde2 = @(E1, eps1, eps2, alpha1) ([E1*sqrt(sin(alpha1)^2+eps1/eps2*cos(alpha1)^2), atan(eps2/eps1*tan(alpha1)) ])
finde2(1, 1, 80, deg2rad(3))