function Connector_Pins
% ERIC LUNDQUIST
% FDFD Code
% Function: To simulate variances in pin fillings
% 2009

close all
clear all
clc

wire_diameter = 1 ; % wire diameter
wire_space = wire_diameter*2 ; % distance between wire centers
distance = 6*wire_diameter+0.5*wire_space ; % distance between grid center and edge (this equation was calculated to minimize error)

% 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
e0 = 8.854e-12; % Permittivity Constant
u = 3e8; % Speed of light

Cap(1)=0;
Cap(2)=0;
e=0;
e(1) = 1;
e(2) = 5; % Relative permittivity of the insulating material
eS = 5 ; % Dielectric constant of connector surface (i.e. 2.5 for plastic, 1 for air)
eP = 2.5 ; % Dielectric constant of pin filling (i.e. plastic)

global m
for m=1:2   
    
    % Wire Characteristics
    er = e(m);
    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)*e0;           % permittivity vector
    BC = zeros(r,c) ;
    
    % Pin Map
    % 1 2
    % 3 4
    
    % ---CONDUCTOR---
    % A = potential (voltage)
    % BC = 1 if boundary condition (metal), 0 if not
    % E = permittivity
    for i=1:rows
        for j=1:columns       
            if (i-(ht-wire_separation))^2+(j-(lt+wire_separation))^2<=round((a/gs)^2) % Pin 1
                A(i,j)=vd;
                BC(i,j)=1;
                E(i,j)=eP*e0;
            elseif (i-(ht-wire_separation))^2+(j-(lt-wire_separation))^2<=round((a/gs)^2) % Pin 2
                A(i,j)=vd;
                BC(i,j)=1;
                E(i,j)=eP*e0;
            elseif (i-(ht+wire_separation))^2+(j-(lt+wire_separation))^2<=round((a/gs)^2) % Pin 3
                A(i,j)=0;
                BC(i,j)=0;
                E(i,j)=eP*e0;
            elseif (i-(ht+wire_separation))^2+(j-(lt-wire_separation))^2<=round((a/gs)^2) % Pin 4
                A(i,j)=0;
                BC(i,j)=0;
                E(i,j)=eP*e0;
            elseif i==5 % grounded shell
                A(i,j)=0;
                BC(i,j)=1;
            elseif j==5 % grounded shell
                A(i,j)=0;
                BC(i,j)=1;
            elseif i==rows-5 % grounded shell
                A(i,j)=0;
                BC(i,j)=1;
            elseif j==columns-5 % grounded shell
                A(i,j)=0;
                BC(i,j)=1;
            end
        end
    end
    %BC(60:80,130:140)=0; %damage
    if m==2
        figure % Visualization of BC Definitions
        imagesc(fliplr(BC));
        axis equal
        colormap(winter)
        title('Visualization of BC Definitions')
        colorbar;
    end
    
    % ---WIRE INSULATION---
    if m==2
        for i=1:rows
            for j=1:columns
                if (i-(ht-wire_separation))^2+(j-(lt-wire_separation))^2>round((a/gs)^2)
                    E(i,j)=er*e0;
                elseif (i-(ht-wire_separation))^2+(j-(lt+wire_separation))^2>round((a/gs)^2)
                    E(i,j)=er*e0;
                elseif (i-(ht+wire_separation))^2+(j-(lt-wire_separation))^2>round((a/gs)^2)
                    E(i,j)=er*e0;
                elseif (i-(ht+wire_separation))^2+(j-(lt+wire_separation))^2>round((a/gs)^2)
                    E(i,j)=er*e0;
                else
                    E(i,j)=eS*e0; % dielectric constant of connector surface
                end
            end
        end
        figure
        imagesc(fliplr(E))
        axis equal
        title('Visualization of Permittivity (Wire Insulation)')
        colorbar
        colormap(winter)
    end
    
    
    % ---ANALYSIS---
    A = FDFD(A,BC,E);
    if m==2
        figure % Visualization of FDFD Fields
        imagesc(fliplr(A))
        axis equal
        title('Visualization of FDFD Fields')
        colorbar;
    end
    
    % Calculation of Charge
    left = c*.15 ;
    right = c*.85 ;
    up = r*.85 ;
    down = r*.15 ;
    
    Joutr = round(right);
    top = round(up);
    bt = round(down);
    Jole = round(left);
    
    if m==2
        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(fliplr(B))
        axis equal
        colorbar
        colormap(winter)
        title('Visualization of Charge Boundary')
    end
    
    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;
    Cap(m)=C ;
    
end

% ---DISPLAY DATA---
Capacitance = Cap(2)           % Capacitance of center wire
VoP = u*sqrt(Cap(1)/Cap(2));   % Velocity of Propagation (VoP)
VoP_c = VoP/u ;                % VoP as fraction of speed of light (c)

% Calculation of Impedance
ZoS = 1/(u*sqrt( Cap(1)*Cap(2) )) % Simulated Characteristic Impedance (Zo) Calculation

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)
global m
status = ['Simulation ', num2str(m) , ' of 2 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])