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

%==========================================================================
% charge distribution
%==========================================================================
clear all
close all
clc

% physical constants
e0 = 8.854e-12;

% charges
q = 1.0e-12 * [-2.07, 9.08, -16.99, 12.48];

% charge positions
xq = [-4.3e-2, -1.5e-2, 1.29e-2, 4.7e-2];
yq = [      0,       0,       0,      0];

% plotting area
lim   = 0.98e-2;
N     = 50;
xm    = linspace(-lim,lim,N);
ym    = linspace(-lim,lim,N);
[x,y] = meshgrid(xm,ym);

%--------------------------------------------------------------------------
% ACTUAL CHARGE DISTRIBUTION
%--------------------------------------------------------------------------
% obatin potential
V = zeros(N);
for i=1:length(q)
    V = V + q(i)./sqrt((x-xq(i)).^2+(y-yq(i)).^2)./(4.*pi.*e0);
end

% obtain force field
[Fx,Fy] = gradient(V);

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

% plot force field
figure(2)
quiver(x,y,Fx,Fy,3) % scale vectors by 3
axis image
xlabel('x (m)')
ylabel('y (m)')
legend('force field')

%--------------------------------------------------------------------------
% ROTATED CHARGE DISTRIBUTION
%--------------------------------------------------------------------------
phi     = 23.5;
deg2rad = pi/180;
phi_rad = phi*deg2rad;
M       = [cos(phi_rad), -sin(phi_rad); sin(phi_rad), cos(phi_rad)];

% rotate all positions
S  = [xq; yq];
S  = M*S;
xq = S(1,:)
yq = S(2,:)

% obatin potential
V = zeros(N);
for i=1:length(q)
    V = V + q(i)./sqrt((x-xq(i)).^2+(y-yq(i)).^2)./(4.*pi.*e0);
end

% obtain force field
[Fx,Fy] = gradient(V);

% plot potential
figure(3)
surf(x,y,V)
xlabel('x (m)')
ylabel('y (m)')
zlabel('V (V)')

% plot force field
figure(4)
quiver(x,y,Fx,Fy,3) % scale vectors by 3
axis image
xlabel('x (m)')
ylabel('y (m)')
legend('force field')

