%   MAIN - MoM solver - metal structures with the triangular mesh
%   ECE 539 April 2010 C ECE WPI

clear all
%--------------------------------------------------------------------------
%   Step 1
%   Create or load geometry structure: P, t
%   Ground plane/TL parameters - coupled TLs
gp_x            = 75e-3;    %   ground plane width (m)
gp_y            = 20e-3;    %   ground plane height (m)
L               = 60e-3;    %   Length of the coupled-TL section (m)
W               = 5e-3;     %   TL width (m)
s               = 1e-3;     %   gap
h               = 1e-3;     %   ground plane height

mesh;
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%   Step 2
%   Create basis functions and the associated geometry matrices
tic
mom1;
toc
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%   Step 3
%   Find the excitation/termination edges
ports;
%--------------------------------------------------------------------------

%--------------------------------------------------------------------------
%   Step 4
%   Frequency sweep
%   Frequency data
fstart = 0.2e9; fstop = 1.1e9; steps = 50; % Hz
f      = fstart+[0:steps-1]*(fstop-fstart)/max(steps-1,1); 

resistance  = zeros(1, steps);  %   port impedance - resistance
reactance   = zeros(1, steps);  %   port impedance - reactance
S11         = zeros(1, steps);  %   port S11
S21         = zeros(1, steps);  %   S21 for active/terminated port
RL          = zeros(1, steps);  %   Return Loss
IL          = zeros(1, steps);  %   Insertion Loss

for m = 1:steps
    disp(strcat('Percentage done:', strcat(num2str(m/steps*100),' %')));   
    %   driving the first port and terminating the second one (if any)
    Z   = mom2(geom, const, f(m));
    for n = 1:length(Index2)
        Resistance  = 50*length(Index2);    %   resistances in parallel
        term        = Index2(n);    
        Z(term,term)= Z(term,term) + Resistance*geom.EdgeLength(term)^2;
    end
    I   = Z\V;                  %   electric current
    Vg  = 1.00;                 %   source(generator) voltage fixed at 1 V
    %   The rest of code follows Ref. D. M. Pozar, Microwave Engineering, 
    %   Wiley, 2005, 3rd edition, Chapter 4 
    %   total current through excitation edge(s)    
    Current         = sum(port_dir'.*I(Index1).*geom.EdgeLength(Index1)');   
    InputImpedance  = Vg/Current;                                       
    resistance(m)   = real(InputImpedance);         %   input resistance
    reactance(m)    = imag(InputImpedance);         %   input reactance
    S11(m)          = (InputImpedance - 50)/(InputImpedance + 50);  
    RL(m)           = 20*log10(abs(S11(m)));                  
    VSWR(m)         = (1 + abs(S11(m)))/(1 - abs(S11(m)));    
    FeedPower(m)    = 1/2*real(Current*conj(Vg));        
    %   S21 (two-port network only; see Pozar 2005, pp. 175-176)
    if ~isempty(Index2)
        %   total current through termination edge(s)
        Current         = sum(term_dir'.*I(Index2).*geom.EdgeLength(Index2)');    
        %   outcoming voltage at port 2 (terminated in a 50 Ohm load)                                                                   
        V2_minus        = Current*50;
        %   incident wave voltage at port 1; obtained using two equations: 
        %   V1_plus + V1_minus = Vg and S11 =V1_minus/V1_plus
        %   S21 for the two-port network        
        V1_plus         = Vg/(1 + S11(m));                                                                                                              
        S21(m)          = V2_minus/V1_plus;                                 
        IL(m)           = 20*log10(abs(S21(m)));
    end
end
%--------------------------------------------------------------------------