function Shield

% ERIC LUNDQUIST
% FDFD Code
% Shielded Parallel Pair
% 2009

close all
clear all
clc

distance = 12 ; % distance between center and matrix edge (m)
wire_diameter = 1 ; % wire diameter (m)
wire_space = 1.5 ; % distance between wire centers

% Normalization
gs = 0.1  ;                                      % grid size
grid_distance = distance/(wire_diameter) ;       % Normalized distance between wire and grid edge
wire_separation = wire_space/(wire_diameter)/gs; % Distance between wire centers

% Constants
er = 1 ;             % Permittivity of Surrounding Medium (i.e. 1 for air, 80 for water)
e0 = 8.854e-12;      % Permittivity Constant
u = 3e8;             % Speed of light

% Wire Characteristics
a = 1;                      % Radius of the conductor
h=2*grid_distance;          % Horizontal Field Boundary
v=2*grid_distance;          % Vertical Field Boundary
rows = round(h/gs);         % number of rows
columns = round(v/gs);      % number of columns
r=rows;
c=columns;
ht=rows/2;                  % vertical Center of the grid
lt=columns/2;               % horizontal Center of the grid
vd = 100;                   % potential of the conductor
A=zeros(rows,columns);      % Grid vector
E = ones(r,c);              % permittivity vector
BC = zeros(r,c) ;

% ---CONDUCTOR---
% (the addition and subtraction of 'radius' determine distance from
% center wire to surrounding wires
for i=1:rows
    for j=1:columns
        if (i-ht)^2+(j-(lt-wire_separation))^2<round((a/gs)^2) % Left
                A(i,j)=vd;
                BC(i,j)=1;
        elseif (i-ht)^2+(j-(lt+wire_separation))^2<round((a/gs)^2) % Right
                A(i,j)=0;
                BC(i,j)=1;
        elseif (i-ht)^2+(j-lt)^2<round((8*a/gs+5)^2) && (i-ht)^2+(j-lt)^2>round((8*a/gs)^2)% Shield
                A(i,j)=0;
                BC(i,j)=1;
        end
    end
end
BC(1:40,:)=0; % Imposed shield damage
figure % Visualization of BC Definitions
imagesc(BC);
colormap(winter)
title('Visualization of BC Definitions')
colorbar;

% ---WIRE INSULATION---
for i=1:rows
    for j=1:columns
        E(i,j)=e0;
    end
end

% ---ANALYSIS---
A = FDFD(A,BC,E);
figure % Visualization of FDFD Fields
imagesc(A);
title('Visualization of FDFD Fields')
colorbar;

% Quiver Plot
B=A(1:5:end,1:5:end); % reduce resolution of fields vector
[x,y]=meshgrid(1:length(B)); % arrays
[dx,dy]=gradient(B);
figure
contour(x,y,B)
hold on
quiver(x,y,dx,dy)
title('Contour Quiver Plot of Field Gradient')

% Zoomed Quiver Plot
Z=A(1:80,125-39:125+40); % zoom in on damaged section
Z=Z(1:5:end,1:5:end);
[x,y]=meshgrid(1:length(Z)); % arrays
[dx,dy]=gradient(Z);
figure
contour(x,y,Z)
hold on
quiver(x,y,dx,dy)
axis equal
title('Contour Quiver Plot of Field Gradient')

% Calculation of Charge
left = c*.25 ;
right = c*.5 ;
up = r*.75 ;
down = r*.25 ;

Joutr = round(right);
top = round(up);
bt = round(down);
Jole = round(left);

figure
B=A; % generate matrix copy for viewing of charge boundary conditions around conductor
B(bt:top,Jole:Joutr)=B(bt:top,Jole:Joutr)+100;
imagesc(B)
colorbar
colormap(winter)
title('Visualization of Charge Boundary')

rs = sum(A(bt+1:top-1,Joutr).*E(bt+1:top-1,Joutr));        % right column 1
tps= sum(A(top,Jole+1:Joutr-1).*E(top,Jole+1:Joutr-1));    % top side 1
bs = sum(A(bt,Jole+1:Joutr-1).*E(bt,Jole+1:Joutr-1));      % bottom right 1
ls = sum(A(bt+1:top-1,Jole).*E(bt+1:top-1,Jole));          % Left column
outercharge = rs+tps+bs+ls;                                % total outer charge

rsin=sum(A(bt+2:top-2,Joutr-1).*E(bt+2:top-2,Joutr-1)) ;         % right column with dielectric e1
tpsin=sum(A(top-1,Jole+2:Joutr-2).*E(top-1,Jole+2:Joutr-2));     % top side
bsin = sum(A(bt+1,Jole+2:Joutr-2).*E(bt+1,Jole+2:Joutr-2));      % bottom right
c1=2*A(top-1,Joutr-1)*E(top-1,Joutr-1);                          % Corners
c2=2*A(bt+1,Joutr-1)*E(bt+1,Joutr-1);
c3=2*A(top-1,Jole+1)*E(top-1,Jole+1);
c4=2*A(bt+1,Jole+1)*E(bt+1,Jole+1);
lsin = sum(A(bt+2:top-2,Jole+1).*E(bt+2:top-2,Jole+1));
innercharge = lsin+rsin+tpsin+bsin + (c1+c2+c3+c4 );

Q=innercharge-outercharge;

% Calculation of Capacitance
C=Q/vd;

% Calculation of Impedance
Z=1/(u*C);

% ---DISPLAY DATA---
Capacitance = C  ;                  % Capacitance of center wire
%VoP = u*C ;                         % Velocity of Propagation (VoP)
%VoP_c = VoP/u ;                     % VoP as fraction of speed of light (c)
Impedance = Z ;                     % Impedance of wire
final_Capacitance = Capacitance
output_Zo = Impedance
% output_VoP = VoP_c

% % Analytical
% d = wire_diameter ;
% D = distance ;
% N = 376.73 ;
% ZoA = N/pi/2/sqrt(er)*acosh((2*D^2-d^2)/d^2)
% 
% Percent_Error = abs(round((ZoA-Z)/ZoA*100))

close(99)

% FUNCTIONS

function A = FDFD(A,BC,E)
% A   = Initial potentials (V)
% BC  = Defines boundary conditions. If BC(i,j) = 1, it is assumed to
%       be a Dirichlet condition (fixed potential)
%       (On the edges of the domain, any BC not equal to 1 is assumed
%       to be a Neumann condition with d-phi/dn = 0)
% E   = Epsilon (permittivity) matrix (F/m)
% Output A = Final potential of the region after FDFD calculations

% Define Constants
tstart = now ;        %  Initialize time for status bar
[Ny,Nx] = size(A);    %  Simulation domain size.
K = 3*max(Nx,Ny);     %  Total iterations
ww = 1.9;             %  Over-relaxation constant
V_0 = max(max(A));    %  Maximum value of domain potential.
maxError = V_0*1e-5;  %  Maximum error before terminating the loop.

% Generate average permittivity matrix (where each point in e_av is average of the four surrounding points).
e_av = E;
for m = 2:Nx-1
    for n = 2:Ny-1
        e_av(n,m) = (E(n+1,m) + E(n-1,m) + E(n,m+1) + E(n,m-1))/4;
    end
end

%  Find Neumann BCs on top/bottom
Neumann_TOP      = find(BC(1,:) ~= 1);
Neumann_BOTTOM   = find(BC(Ny,:) ~= 1);

% Calculate A
for k = 1:K
    maxR = 0; % initialize max change from previous iteration
    A(1,Neumann_TOP) = A(2,Neumann_TOP); % apply top BC
    for n = 2:Ny-1 % apply left BC
        if (BC(n,1) ~= 1)
            A(n,1) = A(n,2);
        end
        for m = 2:Nx-1
            %  Check for boundaries (BC==1) If BC ~=1, apply Poisson's eqn
            if (BC(n,m) ~= 1)
                R = 0.25*(A(n+1,m)*E(n+1,m) + A(n-1,m)*E(n-1,m) ...
                    + A(n,m+1)*E(n,m+1) + A(n,m-1)*E(n,m-1) ...
                    - 4*A(n,m)*e_av(n,m))/e_av(n,m);
                A(n,m) = A(n,m) + ww*R;
                maxR = max(maxR,abs(ww*R)); %  Update maxR
            end
        end
        if (BC(n,1) ~= 1) % apply right BC
            A(n,Nx) = A(n,Nx-1);
        end
    end
    A(Ny,Neumann_BOTTOM) = A(Ny-1,Neumann_BOTTOM); % apply bottom BC
    if maxR < maxError
        disp(k);  % break the loop if change is minimal
        break;
    end
    if mod(k,10)==0
        status_bar(k,K,tstart)
        pause(0.001);
    end
end
return

% Status Bar Function (displays estimated simulation time remaining)
function status_bar(n,nmax,tstart)
status = ['Simulation is ' , num2str(round(n/nmax*100)) , '% Complete'];
bar = zeros(1,nmax);
bar(1,1:n)=1;
time_status=sprintf('\n\nCalculating remaining time...\n');
if n>(nmax/10)
    tc=now;       %current time
    te=tc-tstart; %elapsed time
    t2=te*nmax/n; %estimated total time
    tr=t2-te;     %estimated remaining time
    tr_str=datestr(tr,'HH:MM:SS');
    time_status = sprintf(['\n\nEstimated time remaining: ',tr_str,'\n']);
end
if n>(nmax*.98)
    time_status = sprintf(['\n\nGenerating Visualization Images...\n']);
end
figure(99);
imagesc(bar)
colormap(winter)
daspect([8,100/nmax,1]) % x/y aspect ratio
axis off
title([status,time_status])