% chargedSpheres.m
% Solution to problem 1 of control 2. 2018/02/12
% J.M.Soler, Feb.2018

% Set data
clear
cm = 0.01;                      % one cm in m
dat = [ ...
  0.5   1.2  -0.1  12    7 
 -0.3   0.8  -0.4  15  -12 
  0.4  -0.4   0.1  18   15 
 -0.8   0.1   1.1  20  -15 
 -0.9  -0.7   0.5  10   10 ];
r = dat(:,1:3)';                % sphere center coordinates (m)
R = dat(:,4)*cm;                % sphere radii (m)
V = dat(:,5);                   % sphere potentials (V)
n = size(dat,1);                % number of spheres
k = 8.99e9;                     % Coulomb constant (N*m^2/C^2)

% Find matrix U such that V=U*q
for i = 1:n
    for j = 1:n
        if i==j
            U(i,j) = k/R(i);
        else
            rij = norm(r(:,i)-r(:,j));
            U(i,j) = k/rij;
        end
    end
end

% Solve V=U*q to find q, in C
q = (U\V)'   % note: the transpose is purely aesthetical

% Find the electric field at each sphere
E = zeros(3,n);
for i = 1:n
    for j = 1:n
        if i~=j
            rij = r(:,i)-r(:,j);
            E(:,i) = E(:,i) + k*q(j)*rij/norm(rij)^3; % electric field (V/m)
            f(:,i) = E(:,i)*q(i);                     % electric force (N)
        end
    end
end
E
f

% Plot spheres and fields
[xs,ys,zs] = sphere();      % points on a unit sphere
figure(1)
clf
hold on
for i = 1:n
    x = xs*R(i)+r(1,i);     % points on a sphere of radius R(i) at r(:,i)
    y = ys*R(i)+r(2,i);
    z = zs*R(i)+r(3,i);
    col = (V(i)-min(V))/(max(V)-min(V));  % color of sphere
    col = col*ones(size(x));              % color of each point
    surf(x,y,z,col)                       % plot sphere
end
quiver3(r(1,:),r(2,:),r(3,:),E(1,:),E(2,:),E(3,:))  % plot electric field
colorbar
hold off
xlabel('x (m)') 
ylabel('y (m)') 
zlabel('z (m)')
axis equal
grid on
