%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% HW 1, Problem 1.2 solution key
% Carlos J. Cela, Jan 2012
%
% Calculates and plots E field distribution in a 2D region due to the 
% presence of point electrical charges.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% First, calculate E field
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Define region size
xmin = 0.1;
xmax = 0.2;
ymin = 0.1;
ymax = 0.2;

% Define permittivity of space
epsilon0 = 8.85e-12;

% Define position and magnitude of electric charges
chargeX = [ 0.07  0.15  0.19  0.09];
chargeY = [ 0.21  0.25  0.08  0.08];
charge  = [-1e-9  2e-9 -5e-9  3e-9];

% Scale to 

% Initialize arrays to store E field vector components. The
% region will be divided in a number of points per axis for plotting.
Ex(1:51,1:51) = 0;
Ey(1:51,1:51) = 0; 

% Iterate through each charge on outer loop
% x,y are 'plotting' coordinates
% i,j are 'E field calculation' coordinates
for q = 1:4
  x = 0;
 
  % Now iterate in the region of interest, dividing it in
  % discrete points, to be able to plot
  for i = xmin:(xmax-xmin)/50:xmax
    x = x+1;
    y = 0;
    for j = ymin:(ymax-ymin)/50:ymax
      y = y+1;
    
      % Calculate vector components in the charge-to-point direction
      rx = i-chargeX(q);
      ry = j-chargeY(q);

      % Calculate distance r between current point and this charge
      r = sqrt(rx^2+ry^2);

      % Do not count current charge if calculating point on
      % of charge position
      if( r == 0 )
        continue;
      end

      % Divide vector components by its own magnitude to obtain 
      % versor (unit vector)
      rx = rx/r;
      ry = ry/r;

      % Calculate X and Y contributions for this charge and point, adding
      % it as a vector to the result obtained with previous charges.
      E = (charge(q)/(4*pi*epsilon0*r^2));

      Ex(y,x) = Ex(y,x) + E*rx;
      Ey(y,x) = Ey(y,x) + E*ry;

    end
  end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Now, plot
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Calculate and plot E field magnitude
figure(1);
imagesc(sqrt(Ex.^2+Ey.^2));
colormap hot
grid on

% Flip Y-axis of image (by default this 
% type of plot is upside-down)
set(gca,'YDir','normal');

% Set tick marks on axes
set(gca,'XTick',[1 10 20 30 40 50]);
set(gca,'XTickLabel','0.10|0.12|0.14|0.16|0.18|0.20');
set(gca,'YTick',[1 10 20 30 40 50]);
set(gca,'YTickLabel','0.10|0.12|0.14|0.16|0.18|0.20');

% Add title and units, label axes
title('Electric field magnitude [V/m] and direction at each point', 'fontsize', 14);
xlabel('X position [m]', 'fontsize', 14);
ylabel('Y position [m]', 'fontsize', 14);

% Hold plot to overlap quiver on top
hold on 

% Add color bar scale for reference
colorbar

% Overlap E field vectors, scaling arrows by a factor
% of 2 so the direction can be seen clearly.
h = quiver(Ex, Ey,2);
set(h, 'LineWidth',1,'Color' , [.5 .5 .8]);
%axes([1 51 1 51]);
hold off




