% C3_2.m
% Solution to problem 2 of control 3, 2015/04/17
% J.M.Soler

% Set system parameters
clear all
cm = 1.e-2;                                % one centimeter in m
kc = 9.e9;                                 % Coulomb constant in N*m^2/C^2
rq = [ -4.30, -1.50, +1.29, +4.70 ]*cm;    % charge positions
rp = [ -5.10, -2.70, +3.00, +7.00 ]*cm;    % probe positions
vp = [ -1.30, +3.33, -0.77, +3.00 ];       % probe potentials

% Find charges by solving vmat*q=(kc/(rp-rq))*q=vp
for ip = 1:numel(rp)
    for iq = 1:numel(rq)
          vmat(ip,iq) = kc./abs(rp(ip)-rq(iq));
    end
end
q = vmat\vp';
q = q'

% Check solution
for ip = 1:4
    v(ip) = kc*sum(q./abs(rp(ip)-rq));
end
v

% Plot potential
xmax = 0.98*cm;                   % mesh range
dx = xmax/20;                     % mesh interval
x = (-xmax:dx:+xmax);             % x mesh values
y = x;                            % y mesh values
[x,y] = meshgrid(x,y);            % 2-D mesh
v = 0*x;                          % potential at mesh points
for iq = 1:4
    v = v + kc*q(iq)./sqrt((x-rq(iq)).^2+y.^2);
end
surf(x/cm,y/cm,v)
xlabel('x (cm)')
ylabel('y (cm)')
zlabel('V (V)')
