%Microstrip dimensions
a=input('Enter Height:\n');             %height
b=input('Enter Width:\n');              %width
d=input('Enter lower material height:\n'); %upper material depth
c=a-d;                                  %lower material depth
W=input('Enter Conductor width:\n');    %Conductor width
V=input('Enter Conductor Voltage:\n');  %Conductor voltage
v=3*10^8;                               %Speed of Light

h=input('Enter Mesh size:\n');      %Mesh size


e1=input('Enter upper permittivity:\n')*8.85*10^-12;%Permittivity of upper material
e2=input('Enter lower permittivity:\n')*8.85*10^-12;%Permittivity of lower material

imax=(b/h)+1; %# of x nodes
jmax=(a/h)+1; %# of y nodes
K=imax*jmax;
nbc = 2*imax+2*jmax-4;
[X Y] = meshgrid(1:jmax,1:imax);
mdl = (Y-1)*jmax + X;
mdl(2:end-1,2:end-1) = 0;
bc = zeros(nbc,2);
bc(:,1) = mdl(mdl~=0);
bc(end-imax+1:end,2) = 0;
K = imax*jmax;

A = eye(K)*-4;
for k = 1:K % for each row
    if k-1 > 0  % node at left arm
        A(k,k-1) = 1;
    end
    if k+1 < K  % node at right arm
        A(k,k+1) = 1;
    end
    if k-jmax > 0   % node at upper arm
        A(k,k-jmax) = 1;
    end
    if k+jmax < K   % node at lower arm
        A(k,k+jmax) = 1;
    end
end

B = zeros(K,1);

[I J] = size(bc);

if J == 3 % bc in the form of (i, j, v)
    bc = [(bc(:,1)-1)*jmax+bc(:,2) bc(:,3)]; % bc changed into (k, v)
end

for i = 1:I % boundary condition enforcement
    m = bc(i,1);
    A(m,:) = 0; A(m,m) = 1; B(m) = bc(i,2);
end

%Setting Microstrip dimensions & BC
edge = imax*(c/h+1); %calculates node index at far right of interface
mid = edge-(imax/2); %calculates node index at center of conductor
left = mid-W/2/h;
right = mid+W/2/h;

B(left:right)=-V;

x = A\B;
Phi = reshape(x,jmax,imax);
figure;image(Phi'); axis image; set(gcf,'color','w');
figure;contour(Phi',imax*jmax); grid on; set(gcf,'color','w'); 

%Impedance

%Sets node corners of contour
ContA = floor(imax*(c/2/h)+1+(b-W)/4/h);
ContB = ContA+imax*(c+d)/2/h;
ContC = ContB+(W+(b-W)/2)/h;
ContD = ContA+(W+(b-W)/2)/h;
LeftI = ContA+imax*c/2/h;       %node at Left side of interface
RightI = ContD+imax*c/2/h;      %node at Right side of interface

%Left Side Contour from A to B
lsc = 0;
LU = LeftI-imax;
LD = LeftI+imax;
for ls = ContA:imax:ContB
    %upper material
    for ls = ContA:imax:LU
        lsc = lsc+e1*((x(ls-1)-x(ls+1))/2/h)*h;
    end
    %lower material
    for ls = LD:imax:ContB
        lsc = lsc+e2*((x(ls-1)-x(ls+1))/2/h)*h;
    end
end

lsc = lsc + (e1+e2)*((x(LeftI-1)-x(LeftI+1))/2/h)*h/2;  %Left Side Interface

%Bottom Contour from B to C
bsc = 0;
for bs = ContB:ContC
    bsc = bsc+e2*((x(bs+imax)-x(bs-imax))/2/h)*h;
end

%Right Side Contour from C to D
rsc = 0;
RL = RightI-imax;
RU = RightI+imax;
for rs = ContD:imax:ContC
    %lower material
    for rs = RL:imax:ContC
        rsc = rsc+e2*((x(rs+1)-x(rs-1))/2/h)*h;
    end
    %upper material
    for rs = ContD:imax:RU
        rsc = rsc+e1*((x(rs+1)-x(rs-1))/2/h)*h;
    end
end

rsc = rsc +(e1+e2)*((x(RightI+1)-x(RightI-1))/2/h)*h/2; %Right Side Interface

%Top Contour from D to A

tsc = 0;
for ts = ContA:ContD
    tsc = tsc+e1*((x(ts-imax)-x(ts+imax))/2/h)*h;
end

q = -(lsc+bsc+rsc+tsc);

C = abs(q/V);

%Calculate Co
eo = 8.85*10^-12;
lsco = 0;   %Left contour
bsco = 0;   %Bottom contour
rsco = 0;   %Right contour
tsco = 0;   %Top contour

for lso = ContA:imax:ContB
    lsco = lsco + eo*((x(lso-1)-x(lso+1))/2/h)*h;
end

for bso = ContB:ContC
    bsco = bsco + eo*((x(bso+imax)-x(bso-imax))/2/h)*h;
end

for rso = ContD:imax:ContC
    rsco = rsco + eo*((x(rso+1)-x(rso-1))/2/h)*h;
end

for tso = ContA:ContD
    tsco = tsco + eo*((x(tso-imax)-x(tso+imax))/2/h)*h;
end

qo = -(lsco+bsco+rsco+tsco);
Co = abs(qo/V);

%Impedance
Zo = 1/v/sqrt(C*Co)