% chargedSpheres.m
% JMS solution to problem 1 of control 2, Feb 14, 2022

clear all
cm = 1e-2;                                  % one cm in m
nC = 1e-9;                                  % one nanoC in C
k = 8.99e9;                                 % Coulomb constant in N*m^2/C^2
Rq = [8;3;5]*cm;                            % sphere radii in m
rq = [[20;10;0],[-30;15;0],[10;-25;0]]*cm;  % sphere coordinates, by columns, in m
Vq = [50;-60;-50];                          % sphere potentials, in V
dx = 1*cm;                                  % grid interval, in m
xmax = 50*cm;                               % grid range, in m

% Set matrix A such that V=A*Q
nq = numel(Rq);
for iq = 1:nq
  for jq = 1:nq
    if iq==jq
      A(iq,jq) = k/Rq(iq);
    else
      rij = norm(rq(:,iq)-rq(:,jq));
      A(iq,jq) = k/rij;
    end
  end
end

% Find sphere charges
Q = A\Vq;
fprintf('Q = %f, %f, %f nC\n',Q/nC)

% Set (x,y) mesh grid
x = -xmax:dx:xmax;
y = x;
[x,y] = ndgrid(x,y);

% Find potential at mesh points external to the spheres
V = zeros(size(x));
for iq = 1:nq
  rrq = sqrt( (x-rq(1,iq)).^2 + (y-rq(2,iq)).^2 );
  V = V + k*Q(iq)./rrq;
end

% Find potential at mesh points internal to the spheres
for iq = 1:nq
  rrq = sqrt( (x-rq(1,iq)).^2 + (y-rq(2,iq)).^2 );
  V(rrq<=Rq(iq)) = Vq(iq);
end

% Plot potential using surf
figure(11)
surf(x'/cm,y'/cm,V')
axis equal
xlabel('x (cm)')
ylabel('y (cm)')
zlabel('V (V)')
set(gca,'FontSize',15)
grid on

% Plot equipotentail lines
figure(12)
contour(x'/cm,y'/cm,V')
colorbar
axis equal
xlabel('x (cm)')
ylabel('y (cm)')
set(gca,'FontSize',15)
grid on
