from impedance of shielde microstrip transmission line by ujwol palanchoke
calculates the impedance and phase velocity in shielded microstrip transmission line.

impedance_of_microstrip_txline_using_finite_difference_method.m
clear
clc

% created by ujwol palanchoke

display('for the given problem, lx=120, ly=60,h=15, w=30')
display('the program will evaluate the impedance and phatse velocity in microstrip transmission line')
display('enter n=1 to simulate the above dimension, enter n= integer greater than 1 to simulate other dimensions')
display('you can simulate for different strip voltages, for that enter voltage')

V_m= input('Enter Voltage in strip=');   %voltage in metallic strip
n=input('Enter the integer number to scale the dimensions (note only lx and ly could be changed)=');   % change the dimensions


%define geometry
Lx = n*120; Ly = n*60; h=15;  w=30;  %all dimensions are in mm

% define permitivity
e_1=8.854e-12;  %permitivity of air
e_2=3.5*e_1;    %permitivity of dielectric. you can change premitivity accourding to dielectric used

%define differential lengths and number of nodes in each direction,
%discretize the domain
dx = 5; %differential x length of cells
Nx = (Lx/dx)+1; %number of x nodes
dy = 5; %differential y length of cells
Ny = (Ly/dy)+1; %Number of Y nodes

Ny_d=(h/dy)+1;  %number of Y nodes in dielectics
Ny_a=Ny-Ny_d;   %number of Y nodes in air
Nx_m=w/dx;      %number of x node in metal strip
Nx_nm=(Lx-w)/10; % number of nodes in non-metal , 


% WE Assume the symmetricity, ie the distance between metallic wall and one
% end of strip is equal in both sides

%Initialize grid
V(:,:)=zeros(Ny,Nx); %create 13*25 matrix which carries the grid voltages

%initialize the voltage at metal strip. The v(:,:) forms the matrix
%which depicts the geometry for the given problem. here v(1,1) is the
%top-left corner of the transmission line with shielding. following this
%convention, the metal strip lies in (Ny_a+1,?) where? is the grid points which
%represents the actual position of metal strip


for j=(Nx_nm+1):(Nx_nm+Nx_m)    % we write j for x-axis in problem, x axis 10-16
    V(Nx_nm+1,j)=V_m;           % the voltage at metal strip is the input value
end
% caluclate the voltage at the grid points with dielectric
for iteration=1:200 % for good approximatation we compute the value for 200 times
    for i=2:Ny-1
        for j=2:Nx-1
            V_temp=V(i,j);
            e_com=2*(e_1+e_2);
            if (Ny==Ny_a+1)        %at the air dielectric boundary, 10th row in the matrix. Nx_a=9,
                V(i,j)=((e_1/e_com)*V(i+1,j))+((1/4)*V(i,j-1))+((e_2/e_com)*V(i-1,j))+((1/4)*V(i,j+1));
            else
                V(i,j)=(1/4)*(V(i+1,j)+V(i,j-1)+V(i-1,j)+V(i,j+1));
            end
          for j=(Nx_nm+1):(Nx_nm+Nx_m)    % we write j for x-axis in problem
                V(Nx_nm+1,j)=V_m;           % the voltage at metal strip is the input value
          end
        end
    end
end
V;
imagesc(V) % image for the voltage inside the shielded microstrip.
xlabel('x')
ylabel('y')
title('voltage distribution with dielectric')


% now we calculate the total charge inclosed with dielectric in botttom
% layer

% for this condition e_2=3.5*e_1, and e_1=e_0(permitivity of free space)

% evaluating in top face (L1 in lecture) and bottom face (L2 in lecture)
L_1=0;  %initilize the sum
L_2=0;  %initilize the sum
for j=2:Nx-1
    L_2=L_2+(e_2*(V(Ny,j)-V(Ny-1,j)));        %bottom face
    L_1=L_1+(e_1*(V(1,j)-V(2,j)));          %top face
end
% evaluating in right face L3 in lecture. To do this we have to evaluate in
% dielectric, boundary and air.
% initialize the values 
L_3_1=0;        %first summation in lecture, in dielectric
L_3_2=0;        %second term in lecture, in boundary
L_3_3=0;        %third summation in lecture, in air

for i=2:Ny_a   %Nx_a=9, and row 1-9 of the matrix corresponds to air
    L_3_3=L_3_3+ e_1*(V(i,Nx)-V(i,Nx-1));   %voltage at 25- at 24
end
for i=Ny_a+2:Ny    % row 11-13 corresponds to dielectric
    L_3_1=L_3_1+e_2*(V(i,Nx)-V(i,Nx-1));
end
L_3_2=L_3_2+((e_1+e_2)/2)*(V(Ny_a+1,Nx)-V(Ny_a+1,Nx-1));% row 10 corresponds to boundary

% evaluate the total sum in right face
L_3=L_3_1+L_3_2+L_3_3;


% evaluating in left face L4 in lecture. To do this we have to evaluate in
% dielectric, boundary and air.
% initialize the values 
L_4_1=0;        %first summation in lecture, in dielectric
L_4_2=0;        %second term in lecture, in boundary
L_4_3=0;        %third summation in lecture, in air

for i=2:Ny_a   %Nx_a=9, and row 1-9 of the matrix corresponds to air
    L_4_3=L_4_3+ e_1*(V(i,1)-V(i,2));
end
for i=Ny_a+2:Ny    % row 11-13 corresponds to dielectric
    L_4_1=L_4_1+e_2*(V(i,1)-V(i,2));
end
L_4_2=L_4_2+((e_1+e_2)/2)*(V((Ny_a+1),1)-V((Ny_a+1),2));% row 10 corresponds to boundary

% evaluate the total sum at left face

L_4=L_4_1+L_4_2+L_4_3;

%evaluate the total charge inclosed

Q_ed=-1*(L_1+L_2+L_3+L_4)

%evaluate the capacitance with dielectric

C_ed=Q_ed/V_m      %V_m is the potential of the strip

%**************************************************************************
%**************************************************************************

% now evaluating the voltage and charge enclosed without dielectric
 % for this we repalce the dielectric with air, ie e_2=e_1
 
 e_2=e_1;       % repalce dielectric with air
 
 % same code can be used to evaluate the charge withou dielectric
 
 for j=(Nx_nm+1):(Nx_nm+Nx_m)    % we write j for x-axis in problem, x axis 10-16
    V(Ny_a+1,j)=V_m;           % the voltage at metal strip is the input value
end
% caluclate the voltage at the grid points without dielectric
for iteration=1:200 % for good approximatation we compute the value for 200 times
    for i=2:Ny-1
        for j=2:Nx-1
            V_temp=V(i,j);
            e_com=2*(e_1+e_2);
            if (Ny==Ny_a+1)        %at the air dielectric boundary, 10th row in the matrix. Nx_nm=9,
                V(i,j)=((e_1/e_com)*V(i+1,j))+((1/4)*V(i,j-1))+((e_2/e_com)*V(i-1,j))+((1/4)*V(i,j+1));
            else
                V(i,j)=(1/4)*(V(i+1,j)+V(i,j-1)+V(i-1,j)+V(i,j+1));
            end
          for j=(Nx_nm+1):(Nx_nm+Nx_m)    % we write j for x-axis in problem
                V(Nx_nm+1,j)=V_m;           % the voltage at metal strip is the input value
          end
        end
    end
end
V;
figure
imagesc(V) % image for the voltage inside the shielded microstrip.
xlabel('x')
ylabel('y')
title('voltage distribution without dielectric')

 % evaluating in top face (L1 in lecture) and bottom face (L2 in lecture)
L_1=0;  %initilize the sum
L_2=0;  %initilize the sum
for j=2:Nx-1
    L_2=L_2+(e_2*(V(Ny,j)-V(Ny-1,j)));
    L_1=L_1+(e_1*(V(1,j)-V(2,j)));
end
% evaluating in right face L3 in lecture. To do this we have to evaluate in
% dielectric, boundary and air.
% initialize the values 
L_3_1=0;        %first summation in lecture, in dielectric
L_3_2=0;        %second term in lecture, in boundary
L_3_3=0;        %third summation in lecture, in air

for i=2:Ny_a   %Ny_a=9, and row 1-9 of the matrix corresponds to air
    L_3_3=L_3_3+ e_1*(V(i,Nx)-V(i,Nx-1));
end
for i=Ny_a+2:Ny    % row 11-13 corresponds to dielectric
    L_3_1=L_3_1+e_2*(V(i,Nx)-V(i,Nx-1));
end
L_3_2=L_3_2+((e_1+e_2)/2)*(V(Ny_a+1,Nx)-V(Ny_a+1,Nx-1));% row 10 corresponds to boundary

% evaluate the total sum in right face
L_3=L_3_1+L_3_2+L_3_3;


% evaluating in left face L4 in lecture. To do this we have to evaluate in
% dielectric, boundary and air.
% initialize the values 
L_4_1=0;        %first summation in lecture, in dielectric
L_4_2=0;        %second term in lecture, in boundary
L_4_3=0;        %third summation in lecture, in air

for i=2:Ny_a   %Ny_a=9, and row 1-9 of the matrix corresponds to air
    L_4_3=L_4_3+ e_1*(V(i,1)-V(i,2));
end
for i=Ny_a+2:Ny    % row 11-13 corresponds to dielectric
    L_4_1=L_4_1+e_2*(V(i,1)-V(i,2));
end
L_4_2=L_4_2+((e_1+e_2)/2)*(V(Ny_a+1,1)-V(Ny_a+1,2));% row 10 corresponds to boundary

% evaluate the total sum at left face

L_4=L_4_1+L_4_2+L_4_3;

%evaluate the total charge enclosed

Q_eo=-1*(L_1+L_2+L_3+L_4)

C_eo=Q_eo/V_m      %V_m is the potential of the strip

% evaluating the impedance Z=1/(velocity of light*(sqrt(C_ed*C_eo))
V_l=3e8;        % speed of light
den=sqrt(C_ed*C_eo);
Z_imp=1/(V_l*den)

%evaluating phase velocity
m_t=sqrt(C_eo/C_ed);
U_phase=V_l*m_t

Contact us at files@mathworks.com