%==========================================================================
%
%             CHARACTERISTIC IMPEDANCE OF A MICROSTRIP LINE
%
%       This function calculates the characteristic impedance of a
%   microstrip line by using the finite-difference frequency-domain (FDFD)
%   algorithm.  It also computes the velocity of propagation in the
%   microstrip line.
%
%           [Z,vp] = microstrip_impedance(width,height,e_r),
%
%   where "width" is the width of the microstrip line in meters, "height"
%   is the distance from the ground plane to the strip line, and e_r is the
%   relative permittivity of the substrate.
%
%       The function models a region of space around the microstrip line
%   and assumes a potential of 1.0 V on the stripline.  It then computes
%   the potential across the entire domain by using the FDFD algorithm.
%   Next, the total charge on the strip line is computed by evaluating
%   Gauss's law along a closed contour around the conducting strip.  This
%   value is then used to calculate the capacitance of the stripline by
%   using C = Q/V.  The capacitance is compared with that of an identical
%   system where the substrate is free-space.  This gives the
%   characteristic impedance as 
%           
%                    Z_0 = (1/c_0)*1/sqrt(C*C_0)
%
%   This allows us to solve for the characteristic inductance, which can
%   then be used to find the propagation velocity.
%
%==========================================================================
%
%  James Nagel
%  Project 1 - ECE6340 (Numerical Techniques in EM)
%  2/12/2008
%

function [Z,vp] = microstrip_impedance(width,height,e_r)

%  Physical constants.
e_0 = 8.854e-12;        %  Electric permittivity of free space (F/m).
e_1 = e_0*e_r;          %  Electric permittivity of the dielectric (F/m).
c_0 = 2.994e8;          %  Velocity of light in free space (m/s).

%  Solve for an ideal grid spacing.  A good value is one that crams 40 or
%  more grid points beneath the stripline and the ground plane.
h = height/50;

%  Simulation constants
W = 4*width;            %  "Width" of the simulation domain.
H = 5*height;           %  "Height" of the simulation domain.
Nx = round(W/h);        %  Number of x-grid points.
Ny = round(H/h);        %  Number of y-grid points.
V_0 = 1.0;              %  Voltage on the microstrip line (V)

%  Initialize matrices.
phi = zeros(Ny,Nx);       %  Potential function (V/m)
BC  = zeros(Ny,Nx);       %  Boundary condition (boolean)
q_v = zeros(Ny,Nx);       %  Charge distribution (C/m^2).
eps_0 = e_0*ones(Ny,Nx);  %  Permittivity matrix for free space (F/M).
eps_1 = e_0*ones(Ny,Nx);  %  Permittivity matrix with dielectric (F/M).

%  Indices of the microstrip line.
Iw1 = round(Nx/2 - Nx*0.5*width/W);
Iw2 = round(Nx/2 + Nx*0.5*width/W);
Ih1 = round(Ny - Ny*height/H);

%  Insert initial values and boundary conditions.
eps_1(Ny - height/h + 1:Ny,:) = e_1;    %  Substrate permittivity
phi(Ih1,Iw1:Iw2) = V_0;                 %  Microstrip voltage
BC(Ih1,Iw1:Iw2) = 1;                    %  Microstrip Dirichlet BC.
BC(end,:) = 1;                          %  Ground plane Dirichlet BC.

%  Solve for the potentials using the FDFD algorithm.
disp('Calculating potentials.');
phi_1 = FDFD(phi,BC,eps_1,q_v,h);

%  Generate a second phi-matrix with an air-filled dielectric.
disp('Calculating free-space potentials.');
phi_0 = FDFD(phi,BC,eps_0,q_v,h);

%  Contour indices for Gauss's law.
Ix1 = round(Iw1 - width/(2*h));
Ix2 = round(Iw2 + width/(2*h));
Iy1 = round(Ih1 - height/(.5*h));
Iy2 = round(Ih1 + height/(2*h));
Cx = [Ix1:Ix2];
Cy = [Iy1:Iy2];

%  Calculate the charge on the microstrip line by integrating a contour
%  around the strip.
phi = phi_0;
Q1 = sum(phi(Iy1-1,Cx).*eps_0(Iy1-1,Cx) - phi(Iy1,Cx).*eps_0(Iy1,Cx));
Q2 = sum(phi(Iy2+1,Cx).*eps_0(Iy2+1,Cx) - phi(Iy2,Cx).*eps_0(Iy2,Cx));
Q3 = sum(phi(Cy,Ix2+1).*eps_0(Cy,Ix2+1) - phi(Cy,Ix2).*eps_0(Cy,Ix2));
Q4 = sum(phi(Cy,Ix1-1).*eps_0(Cy,Ix1-1) - phi(Cy,Ix1).*eps_0(Cy,Ix1));
Q_0 = Q1 + Q2 + Q3 + Q4;
C_0 = -Q_0/V_0;             %  Capacitance (free-space).

%  Repeat for the case with the dielectric.
phi = phi_1;
Q1 = sum(phi(Iy1-1,Cx).*eps_1(Iy1-1,Cx) - phi(Iy1,Cx).*eps_1(Iy1,Cx));
Q2 = sum(phi(Iy2+1,Cx).*eps_1(Iy2+1,Cx) - phi(Iy2,Cx).*eps_1(Iy2,Cx));
Q3 = sum(phi(Cy,Ix2+1).*eps_1(Cy,Ix2+1) - phi(Cy,Ix2).*eps_1(Cy,Ix2));
Q4 = sum(phi(Cy,Ix1-1).*eps_1(Cy,Ix1-1) - phi(Cy,Ix1).*eps_1(Cy,Ix1));
Q_1 = Q1 + Q2 + Q3 + Q4;
C_1 = -Q_1/V_0;             %  Capacitance (dielectric).

%  Finally, we have the characteristic impedance.
Z_numerical = (1/c_0)*1/sqrt(C_0*C_1);

%  Solve for inductance and propagation velocity.
L = Z_numerical^2*C_1;
vp = 1/sqrt(L*C_1);

%  Compare with an analytically derived impedance.
Z_analytical = 87*(1/sqrt(1.41+e_r))*log(5.98*height/(0.8*width + h));

%  This equation was taken from Wikipedia.  It gives slightly closer
%  results.
% a1 = 60/sqrt(2*(1+e_r));
% a2 = 4*height/width;
% a3 = (14 + 8/e_r)*(4*height)/(11*width);
% a4 = a3;
% a5 = 0.5*pi^2*(1+1/e_r);
% Z_analytical = a1*log(1 + a2*(a3 + sqrt(a4^2 + a5)));

disp(['Z_0 (numerical) = ', num2str(Z_numerical), ' ohms']);
disp(['Z_0 (analytical) = ', num2str(Z_analytical), ' ohms']);
disp(['Propagation velocity = ',num2str(vp/3e8), ' c']);

%  Plot the potential functions.
figure(1)
clf
colormap bone;
imagesc(phi_0);
title('Free-space potential');
figure(2)
clf
colormap bone;
imagesc(phi_1);
title(['Potential with substrate (\epsilon_r = ',num2str(e_r),')']);

Z = Z_numerical;

return;

