%
%   Book: MATLAB-Based Electromagnetics (Pearson Prentice Hall)
%   Author: Branislav M. Notaros
%   Instructor Resources
%   (c) 2011
%
%   This MATLAB code or any part of it may be used only for 
%   educational purposes associated with the book
%
%
%

%   Electric force due to multiple charges

clear all;
close all;
EPS0 = 8.8542*10^(-12); 

% Input
N = input('Enter the number of charges (positive integer): ');
for i = 1:N
    x(i) = input(['Enter the X coordinate in cm for charge ',...
        int2str(i) ,':  ']);
    y(i) = input(['Enter the Y coordinate in cm for charge ',...
        int2str(i),':  ']);
    z(i) = input(['Enter the Z coordinate in cm for charge ',...
        int2str(i),':  ']);
    Q(i) = input(['Enter the charge in nC for charge ',...
        int2str(i),':  ']);
end
disp('Properties of point of interest');
xp = input('Enter the X coordinate in cm : ');
yp = input('Enter the Y coordinate in cm : ');
zp = input('Enter the Z coordinate in cm : ');
Qp = input('Enter the charge in nC for point of interest: ');

x = x * 10^(-2);
y = y * 10^(-2);
z = z * 10^(-2);
Q = Q * 10^(-9);

xp = xp * 10^(-2);
yp = yp * 10^(-2);
zp = zp * 10^(-2);
Qp = Qp * 10^(-9);

% Compute distance and direction between point of interest and each charge

r = sqrt((xp - x).^2 + (yp - y).^2 + (zp - z).^2);
ux = (xp - x)./r;
uy = (yp - y)./r;
uz = (zp - z)./r;
uVec = [ux; uy; uz];

% Electric force computation

F = (ones(3,1)*(Qp*Q./(4*pi*EPS0*r.^2))).*uVec; 
Ftot = sum(F,2);
Fmag = vectorMag(Ftot);
Fuv = (Ftot/Fmag)';

% Output

fprintf(['Magnitude of resultant force at '...
'point P is %f mN.\n'],Fmag*1000 );
disp('Unit vector of resultant force :');
disp(Fuv);

figure(1);
plot3(0,0,0,'k');
hold on;
for i=1:N
plot3(x(i),y(i),z(i),'o','MarkerSize',10,'MarkerFaceColor','b');
line ([xp , x(i)],[yp,y(i)], [zp,z(i)],'LineStyle',':');
hold on;
vecPlot3D([xp yp zp],[xp yp zp] + F(:,i)',1/Fmag/100,'b',0);
hold on;
end;
vecPlot3D([xp yp zp],[xp yp zp] + Ftot',1/Fmag/100,'r',1);
text(1.4*xp,1.4*yp,1.4*zp,'Total electric force','Color','r');

hold off;
axis equal;
xlabel('x');
ylabel('y');
zlabel('z');
title({'Force from ',int2str(N),'-point charges'});
