% Unit 3 (MATLAB Course @ UAM)
% (c) Alexander Knebe

clear all
close all
clc

% physical constants
e0 = 8.854e-12;

% charge positions
rq = [-4.3e-2, -1.5e-2, 1.29e-2, 4.7e-2];

% probe positions
ra  = -5.1e-2; 
rb  = -2.7e-2;
rc  = 3.0e-2;
rd  = 7.0e-2;

% probe potentials
Va  = -1.3;
Vb  = 3.33;
Vc  = -0.77;
Vd  = 3.0;

lim = 0.98e-2;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% linear system; Vj = (rq,i-rj) qj
A = [ 1./abs(rq(1)-ra) 1./abs(rq(2)-ra) 1./abs(rq(3)-ra) 1./abs(rq(4)-ra);
      1./abs(rq(1)-rb) 1./abs(rq(2)-rb) 1./abs(rq(3)-rb) 1./abs(rq(4)-rb);
      1./abs(rq(1)-rc) 1./abs(rq(2)-rc) 1./abs(rq(3)-rc) 1./abs(rq(4)-rc);
      1./abs(rq(1)-rd) 1./abs(rq(2)-rd) 1./abs(rq(3)-rd) 1./abs(rq(4)-rd)] ./ (4.*pi.*e0)
V = [Va;Vb;Vc;Vd];
q = A\V

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
N = 50;
xm = linspace(-lim,lim,N);
ym = linspace(-lim,lim,N);
[x,y] = meshgrid(xm,ym);
V = zeros(N);
for i=1:length(q)
    V = V + q(i)./sqrt((x-rq(i)).^2+y.^2)./(4.*pi.*e0);
end

figure(1)
surf(x,y,V)
xlabel('x (m)')
ylabel('y (m)')
zlabel('V (V)')
