function FDFD_Microstrip

% Eric Lundquist
% Microstrip FDFD Simulation
% March 2009

close all
clear all
clc

epsilon = 2.5 ; % dielectric constant of substrate
thickness = 8 ; % substrate thickness
width = 8 ;     % microstrip width

% NORMALIZATION
gs = thickness/50;      % Grid size
w = width/gs;           % number of units in width
t = thickness/gs;       % number of units in substrate thickness

Cap(1)=0;
Cap(2)=0;

global m
for m=1:2                % Start loop for calculation of both capacitances
    
    e(1) = 1;
    e(2) = epsilon;      % Relative permittivity of the insulating material 
    er = e(m);           % Permittivity thus changes for the second loop

    % Constants
    e0 = 8.854e-12;      % Permittivity 
    u = 3e8;             % Speed of light 

    % Wire Characteristics 
    sf = 5 ;                   % scaling factor for width of grid relative to microstrip width
    medium = 1;                % relative permittivity of the surrounding region (i.e. air=1, water=80)
    rows = round(4*t);         % number of rows
    columns = round(sf*w);     % number of columns
    r=rows;
    c=columns;
    ht = round(rows/2);        % vertical Center of the grid 
    lt = round(columns/2);     % horizontal Center of the grid
    vd = 100;                  % potential of the conductor 
    A = zeros(r,c);            % grid vector
    E = ones(r,c);             % permittivity vector
    mloc = round(t) ;          % location of microstrip (row #)
    BC = zeros(r,c) ;
    
    % ---CONDUCTOR---
    A(mloc,columns*2/5:columns*3/5)=vd;
    BC(mloc,columns*2/5:columns*3/5)=1;

    % ---WIRE INSULATION---
    E(:,:)=e0;
    E(1:mloc-1,:)=er*e0;
    BC(1,:)=1;
    if m==2
        figure
        imagesc(flipdim(E,1))
        title('Visualization of Substrate (Permittivity)');
        colorbar
        colormap(winter)
    end
    
    % ---ANALYSIS---
    
    A = FDFD(A,BC,E) ;
    if m==2
    figure
    imagesc(flipdim(A,1));
    title('FDFD Fields')
    colormap('hot')
    colorbar
    end

    % ---CALCULATION OF CHARGE---
    
    % Charge Boundary Location
    left = c*.25 ;
    right = c*.75 ;
    up = r*.75;
    down = mloc*.25 ;
    
    % Charge Boundary Visualization
    if m==2
    B = 0 ;
    B = A ; % copy A to B for visualization
    B(down:up,left:right) = vd-B(down:up,left:right) ;
    figure
    imagesc(flipdim(B,1))
    title('Visualization of Charge Boundaries')
    colormap(winter)
    end
    
    Jole = round(left);
    Joutr = round(right);
    top = round(up);
    bt = round(down);
    
    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

% Calculation of Impedance
Z=1/(u*sqrt( Cap(1)*Cap(2) ));

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

% ---FIND EXPECTED VALUES---
Er = e(2);                                              % Permittivity of substrate
d = thickness;                                           % Substrate thickness
W = width;                                              % Microstrip width
Eeff = (Er+1)/2 + ((Er-1)/2) * 1/sqrt(1+12*d/W);        % Effective Permittivity

if W/d <= 1                                             % Calculation of Zo
    Zo = 60/sqrt(Eeff)*log(8*d/W + 0.25*W/d);
elseif W/d > 1
    Zo = 120*pi/(sqrt(Eeff)*(W/d+1.393+0.667*log(W/d+1.444)));
end

Capacitance = Cap(2);
Analytic_Impedance = Zo

Percent_Error = abs(round((Z-Analytic_Impedance)/Analytic_Impedance*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)
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])