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

%   FD computer program - direct solution

clear all;
close all;

a = 0.01;
b = 0.03;
Va = 1;
Vb = -1;
EPS0=8.8542*10^(-12);
N = [2 3 5 7 9 10 12 25];
%N = 10;
%m = 1;
for m = 1 : length(N)
    d = a/N(m);
    %number of inner nodes
    N1 = N(m) + 1;
    %number of outer node
    N2 = b/a *N(m) + 1;   
    %K = N2^2; % total nodes
  
    Vstart = zeros(N2,N2);
    Vstart(1,:) = Vb;
    Vstart(:,1) = Vb;
    Vstart(N2,:) = Vb;
    Vstart(:,N2) = Vb;
    
    lim1 = (N2-N1)/2 + 1;
    lim2 = (N2+N1)/2;
    Vstart(lim1:lim2,lim1:lim2)=Va;
    
    % A,C matrix
    [A,C]= mACfd(Vstart,N2);
    V = inv(A)*C;  
   
   %V2D
    for i = 1:N2
    V2D(i,:) = V((i-1)*N2+1:i*N2);  
    end;
    
   [x,y]= meshgrid(0:d:b);
   [Ex,Ey] = gradient(-V2D,d,d);
   
   sigmaOut = zeros(1,N2-1);
   sigmaIn = zeros(1,N1-1);
   
for i = 1:N2-1
sigmaOut(i) = EPS0/2/d*(3/2*V2D(1,i)-2*V2D(2,i)+1/2*V2D(3,i)+...
    3/2*V2D(1,i+1)-2*V2D(2,i+1)+1/2*V2D(3,i+1));
end;
k = (N2-N1)/2 + 1;

for i = k:(N2 + N1)/2 -1
sigmaIn(i-k+1)=EPS0/2/d*(3/2*V2D(k,i)-2*V2D(k-1,i)+1/2*V2D(k-2,i)+...
    3/2*V2D(k,i+1)-2*V2D(k-1,i+1) + 1/2*V2D(k-2,i+1));
end;

Qouter(m) = 4*d*sum(sigmaOut);
Qinner(m) = 4*d*sum(sigmaIn);


Cout(m)= Qouter(m)/(Vb - Va);
Cin(m) = Qinner(m)/(Va - Vb);

figure(4*m - 3);
 quiver (x,y,Ex,Ey); xlabel('x [m]'); ylabel('y [m]');
 title(['Electric field intensity vector at each node for N = '...
     ,num2str(N(m)),', direct FD method']);axis equal;  
     
   figure(4*m - 2);
   surf(x,y,V2D); shading interp; colorbar; 
   xlabel('x [m]'); ylabel('y [m]'); zlabel('V [V]');
   title(['Potential distribution for N = ', num2str(N(m))...
       ', direct FD method']);

   figure(4*m - 1);
   dinner = a/(length(sigmaIn)-1);
   innerCond = 0:dinner:a;
   plot(innerCond, sigmaIn);
   xlabel('x [m]'); ylabel('\rho_{in} [C/m^2]');
   
   figure(4*m);
   douter = b/(length(sigmaOut)-1);
   outerCond = 0:douter:b;
   plot(outerCond,sigmaOut);
   xlabel('x [m]'); ylabel('\rho_{out} [C/m^2]');
   
error(m) = (Qinner(m)+Qouter(m))/Qouter(m)*100;
clear Vstart V V2D sigmaIn sigmaOut;
end;