%
%   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
%
%
%


%   MoM program for a charged parallelepiped


clear all;
close all;

EPS0 = 8.8542*10^(-12);
a = 1;
b = 2;
c = 3;
n = 5;
m = 10;
k = 15;
V0 = 1;

% Local coordinates of the central points in xy plane
[p1,q1,S1] = localCoordinates(n,m,a,b);
% Local coordinates of the central points in xz plane
[p2,q2,S2] = localCoordinates(n,k,a,c);
% Local coordinates of the central points in yz plane
[p3,q3,S3] = localCoordinates(m,k,b,c);


% Make global coordinates
x = [p1,p1,p2,p2,zeros(1,m*k),a*ones(1,m*k)];
y = [q1,q1,zeros(1,n*k),b*ones(1,n*k),p3,p3];
z = [zeros(1,n*m),c*ones(1,n*m),q2,q2,q3,q3];
S = [S1, S1, S2, S2, S3, S3];
 
A = matrixA(EPS0,S,x,y,z);
%matrix B
N = length(x); 
B = V0 * ones(N,1);  
rhos = A\B;  %solving the matrix equation

% Total charge computation
Qtot = totalCharge(S,rhos);
fprintf('\nTotal charge of the parallelepiped is %2.3d pC', Qtot*10^12);
 
% Form 2D representation of charge distribution over one side
 
for i = 1:n
rhos2D1(i,:) = rhos((i-1)*m +1:i*m);
end;
 
figure(1);
surf(rhos2D1);
colormap('cool');
shading interp;
title('Surface charge distribtion');
xlabel('x');
ylabel('y');
zlabel('RhoS');
 
t = 2*m*n;
for i = 1:n
rhos2D2(i,:) = rhos(t+(i-1)*k +1:t+i*k);
end;

figure(2);
surf(rhos2D2);
colormap('cool');
shading interp;
title('Surface charge distribution');
xlabel('x');
ylabel('y');
zlabel('RhoS');
 
   
t = t + 2*n*k;

for i = 1:m
rhos2D3(i,:) = rhos(t+(i-1)*k +1:t+i*k);
end;

figure(3);
surf(rhos2D3);
colormap('cool');
shading interp;
title('Surface charge distribution');
xlabel('x');
ylabel('y');
zlabel('RhoS');

 