% This code calculates the 'discrete' charges over two parallel plates -
% forming a capacitor.

%Due to symmetry the individual charges have been reduced to six unique
%charges

% sigma_1, sigma_2, sigma_6, sigma_17, sigma_18 and sigma_22 they are
% embedded in potential equations, however they may be found using
% simultaneous equations.

%Defining geometric variables:

h = 1;              %Distance between centroids [m]
d = 10;             %Plate Separation [m]
E_0 = 8.85E-12;     %epsilon_0 [F/m]
k = ((h^2)/(4*pi*E_0)); %constant which is removed for simplification

vt = 1;
vb = -1;

% Creating the large potential equation:

potential_1_sig_1 = (k*((4*log(1 + sqrt(2))/h) + (1/(3*(sqrt(2)*h))) + (2/(3*h))));
potential_1_sig_2 = (k*((2/(sqrt(13)*h) + (2/(sqrt(10)*h)) + (3/h))));
potential_1_sig_6 = (k*((2/(sqrt(5)*h)) + (3/(2*sqrt(2)*h))));
potential_1_sig_17 = (k*((1/(sqrt(18*h^2 + d^2))) + (2/(sqrt(9*h^2 + d^2))) + (1/d)));
potential_1_sig_18 = (k*((2/(sqrt(13*h^2 + d^2))) + (2/(sqrt(10*h^2 + d^2))) + (2/(sqrt(4*h^2 + d^2)))+ (2/(sqrt(h^2 + d^2))))); 
potential_1_sig_22 = (k*((1/(sqrt(8*h^2 + d^2))) + (2/(sqrt(5*h^2 + d^2))) + (1/(sqrt(2*h^2 + d^2)))));
                
potential_2_sig_1 = k*((1/(sqrt(13)*h)) + (1/(sqrt(10)*h)) + (3/(2*h)));
potential_2_sig_2 = k*((4*log(1 + sqrt(2))/h) + (1/(sqrt(10)*h)) + (2/(sqrt(5)*h)) + (3/(2*sqrt(2)*h)) + (4/(3*h)));
potential_2_sig_6 = k*((1/(sqrt(5)*h)) + (1/(sqrt(2)*h)) + (3/(2*h)));
potential_2_sig_17 = k*((1/(sqrt(13*h^2 + d^2))) + (1/(sqrt(10*h^2 + d^2))) + (1/(sqrt(4*h^2 + d^2))) + (1/(sqrt(h^2 + d^2))));
potential_2_sig_18 = k*((1/(sqrt(10*h^2 + d^2))) + (1/(sqrt(9*h^2 + d^2))) + (1/(sqrt(8*h^2 + d^2))) + (2/(sqrt(5*h^2 + d^2))) + (1/(sqrt(2*h^2 + d^2))) + (1/(sqrt(h^2 + d^2)))+ (1/d));
potential_2_sig_22 = k*((1/(sqrt(5*h^2 + d^2))) + (1/(sqrt(4*h^2 + d^2))) + (1/(sqrt(2*h^2 + d^2))) + (1/(sqrt(h^2 + d^2))));
                
potential_6_sig_1 = k*((2/(sqrt(5)*h)) + (3/(2*sqrt(2)*h)));
potential_6_sig_2 = k*((2/(sqrt(5)*h)) + (2/(sqrt(2)*h)) + (3/h));
potential_6_sig_6 = k*((4*log(1 + sqrt(2))/h) + (1/(sqrt(2)*h)) + (2/h));
potential_6_sig_17 = k*((1/(sqrt(8*h^2 + d^2))) + (2/(sqrt(5*h^2 + d^2))) + (1/(sqrt(2*h^2 + d^2))));
potential_6_sig_18 = k*((2/(sqrt(5*h^2 + d^2))) + (2/(sqrt(4*h^2 + d^2))) + (2/(sqrt(2*h^2 + d^2))) + (2/(sqrt(h^2 + d^2))));
potential_6_sig_22 = k*((1/(sqrt(2*h^2 + d^2))) + (2/(sqrt(h^2 + d^2))) + (1/d));
                
potential_17_sig_1 = k*((1/(sqrt(18*h^2 + d^2))) + (2/(sqrt(9*h^2 + d^2))) + (1/d));
potential_17_sig_2 = k*((2/(sqrt(13*h^2 + d^2))) + (2/(sqrt(10*h^2 + d^2))) + (2/(sqrt(4*h^2 + d^2))) + (2/(sqrt(h^2 + d^2))));
potential_17_sig_6 = k*((1/(sqrt(8*h^2 + d^2))) + (2/(sqrt(5*h^2 + d^2))) + (1/(sqrt(2*h^2 + d^2))));
potential_17_sig_17 = k*((4*log(1 + sqrt(2))/h) + (1/(3*sqrt(2)*h)) + (2/(3*h)));
potential_17_sig_18 = k*((2/(sqrt(13)*h)) + (2/(sqrt(10)*h)) + (3/h));
potential_17_sig_22 = k*((2/(sqrt(5)*h)) + (3/(2*sqrt(2)*h)));
                
potential_18_sig_1 = k*((1/(sqrt(13*h^2 + d^2))) + (1/(sqrt(10*h^2 + d^2))) + (1/(sqrt(4*h^2 + d^2))) + (1/(sqrt(h^2 + d^2))));
potential_18_sig_2 = k*((1/(sqrt(10*h^2 + d^2))) + (1/(sqrt(9*h^2 + d^2))) + (1/(sqrt(8*h^2 + d^2))) + (2/(sqrt(5*h^2 + d^2))) + (1/(sqrt(2*h^2 + d^2))) + (1/(sqrt(h^2 + d^2))) + (1/d));
potential_18_sig_6 = k*((1/(sqrt(5*h^2 + d^2))) + (1/(sqrt(4*h^2 + d^2))) + (1/(sqrt(2*h^2 + d^2))) + (1/(sqrt(h^2 + d^2))));
potential_18_sig_17 = k*((1/(sqrt(13)*h)) + (1/(sqrt(10)*h)) + (3/(2*h)));
potential_18_sig_18 = k*((4*log(1 + sqrt(2))/h) + (1/(sqrt(10)*h)) + (2/(sqrt(5)*h)) + (3/(2*(sqrt(2)*h))) + (4/(3*h)));
potential_18_sig_22 = k*((1/(sqrt(5)*h)) + (1/(sqrt(2)*h)) + (3/(2*h)));

potential_22_sig_1 = k*((1/(sqrt(8*h^2 + d^2))) + (2/(sqrt(5*h^2 + d^2))) + (1/(sqrt(2*h^2 + d^2))));
potential_22_sig_2 = k*((2/(sqrt(5*h^2 + d^2))) + (2/(sqrt(4*h^2 + d^2))) + (2/(sqrt(2*h^2 + d^2))) + (2/(sqrt(h^2 + d^2))));
potential_22_sig_6 = k*((1/(sqrt(2*h^2 + d^2))) + (2/(sqrt(h^2 + d^2))) + (1/d));
potential_22_sig_17 = k*((2/(sqrt(5)*h)) + (3/(2*sqrt(2)*h)));
potential_22_sig_18 = k*((2/(sqrt(5)*h)) + (2/(sqrt(2)*h)) + (3/h));
potential_22_sig_22 = k*((4*log(1 + sqrt(2))/h) + (1/(sqrt(2)*h)) + (2/h));

charges = [potential_1_sig_1 potential_1_sig_2 potential_1_sig_6 potential_1_sig_17 potential_1_sig_18 potential_1_sig_22;
            potential_2_sig_1 potential_2_sig_2 potential_2_sig_6 potential_2_sig_17 potential_2_sig_18 potential_2_sig_22;
            potential_6_sig_1 potential_6_sig_2 potential_6_sig_6 potential_6_sig_17 potential_6_sig_18 potential_6_sig_22;
            potential_17_sig_1 potential_17_sig_2 potential_17_sig_6 potential_17_sig_17 potential_17_sig_18 potential_17_sig_22;
            potential_18_sig_1 potential_18_sig_2 potential_18_sig_6 potential_18_sig_17 potential_18_sig_18 potential_18_sig_22;
            potential_22_sig_1 potential_22_sig_2 potential_22_sig_6 potential_22_sig_17 potential_22_sig_18 potential_22_sig_22]
 
v = [vt;
    vt;
    vt;
    vb;
    vb;
    vb]

charge = charges\v

qt = (h^2*(4*charge(1)+(8*charge(2))+(4*charge(3))))
qb = (h^2*(4*charge(4)+(8*charge(5))+(4*charge(6))))

Ctotal = (qt/(vt-vb))

charge = charge*1E12

%z = [30 20 20 30; 20 10 10 20; 20 10 10 20; 30 20 20 30];
z = [charge(1) charge(2) charge(2) charge(1); charge(2) charge(3) charge(3) charge(2); charge(2) charge(3) charge(3) charge(2); charge(1) charge(2) charge(2) charge(1)];
x=[0.5:1:3.5];
y=[0.5:1:3.5];

surf(x,y,z)
axis([0 4 0 4 0 16])

shading interp