program CirAnl; { Performs Nodal Circuit Analysis for Freq Resp }

{             Jeff Crawford         December 10, 1984              }
{                                                                  }
{                          Version  1.0                            }
{                                                                  }
{                           REFERENCES                             }
{                                                                  }
{   [1]  "Unify Two-Port Calculations", D. J. Daruvala, Electronic }
{        Design, January 4, & January 18, 1974 - Two Articles      }
{   [2]  "High-Frequency Amplifiers, 2nd Edition", Ralph S. Carson,}
{        John Wiley & Sons, 1982                                   }


type     matrix    = array[1..25,1..25] of real;
         vector    = array[1..100] of real;
         XYData    = record
                       XX,YY: real;
                     end;
         stat      = string[1];


var      L,Cap         : matrix; { Inductance, Capacitance Terms Sqmodel }
         Fyi,Fyr       : matrix; { Admittance Matrices; Real, Imag       }
         F,F_low,F_high: real;   { Frequency Information                 }
         F_inc         : real;   { Increment of Frequency                }
         Nodes         : integer;{ Number of Nodes in Circuit; Passed from}
                                 { Sqmodel Main Program                  }
         omega         : real;   { Radian Frequency                      }
         N             : integer;{ Parameter in Node Reduction Routine   }
                                 { N Starts at Nodes and is Reduced by 1 }
                                 { Until Arriving at N=2 for a Two-Port  }
         No_turns      : integer;{ # Turns in Inductor - Passed from Main}
         A,B           : matrix; { Equiv of Fyr,Fyi Matrices; Used in    }
                                 { Matrix Reduction Routine to a 2 x 2   }
         i,j           : integer;{ Used in Looping                       }
         Freq,Mag      : vector; { Arrays to Hold Information of Z v.s. F}
         Ans,Ans1,Ans2 : stat;   { Controls Function Routing          }
         Num           : integer; { Number of Data Points Store on Disk  }


Procedure Write_File( X,Y: vector);

var      j              : integer;
         write_file     : file of XYData;
         read_file      : file of XYData;
         file_rec       : XYData;
         file_name      : string[10];
         file_size      : integer;
         OK             : Boolean;
         Ans1,Ans2      : string[1]; { Controls Program Routing      }

label 1,2;

begin
2:   clrscr;
     gotoXY(10,10);
     write('Write Data to What File  ? ');
     readln(file_name);
     assign(write_file,file_name);
     {$I-} Reset(write_file) {$I+};
     OK:= (IOresult = 0);     { Determines if File Already Exists }
     if OK then
         begin
         clrscr;
         gotoXY(10,10);
         write('File Already Exists . . . Overwrite ?  Y / N ');
         readln(Ans2);
         Ans2:= upcase(Ans2);
         if Ans2 <> 'Y' then goto 1;
         end;
    assign(write_file,file_name);
    rewrite(write_file);
    with file_rec do
         begin
         for j:= 1 to Num do   { Writes Data to Disk File }
              begin
              XX:= X[j];
              YY:= X[j];
              write(write_file,file_rec);
              end;
         end;
    flush(write_file);
    close(write_file);
1:  if Ans2 = 'N' then
         begin
         clrscr;
         gotoXY(10,10);
         write('Write Data to Different File   Y / N ?  ');
         readln(Ans1);
         Ans1:= upcase(Ans1);
         if Ans1 = 'Y' then goto 2;
         end;
end;


Procedure Swap(var X1,Y1: real); { Used to Swap Rows & Columns to   }
                             { Change Input and/or Output Terminals }
                             { for Calculation; Declared "var" so   }
                             { Switch is Made in Major Matrices     }

    var   XX: real;

    begin
         XX:= X1;
         X1:= Y1;
         Y1:= XX;
    end;



Procedure Save_C(V: real; i,j: integer);  { Saves Capacitors }

    var  temp: real;

    begin
         temp:= omega*V;
         Fyi[i,i]:= Fyi[i,i] + temp;
         Fyi[j,j]:= Fyi[j,j] + temp;
         Fyi[i,j]:= Fyi[i,j] - temp;
         Fyi[j,i]:= Fyi[j,i] - temp;
    end;

Procedure Save_L(V: real; i,j: integer);

    var  temp: real;

    begin
         temp:= -1.0/(omega*V);
         Fyi[i,i]:= Fyi[i,i] + temp;
         Fyi[j,j]:= Fyi[j,j] + temp;
         Fyi[i,j]:= Fyi[i,j] - temp;
         Fyi[j,i]:= Fyi[j,i] - temp;
    end;



Procedure Fill_Caps(var N: integer);
    var  i : integer;

    begin
         for i:= 2 to N-1 do
              begin
                   Save_C(cap[i-1,i-1],i,0);
              end;
         for i:= 2 to N-2 do
              begin
                   Save_C(cap[i-1,i], i, i+1);
              end;
    end;


Procedure Fill_Inductors(var N: integer);

    var  i              : integer;
         Ind_val        : real;

    begin
         Save_L(0.5*L[1,1],1,2);
         Save_L(0.5*L[N-2,N-2], N-1,N);
         for i:= 2 to N-2 do
              begin
                   Ind_val:= 1/(0.5*(L[i-1,i-1] + L[i,i])) + 1.0/L[i-1,i];
                                { i,j are Coupled Turn Designators }
                   Save_L(1.0/Ind_val,i,i+1);
              end;
         for i:= 2 to N-1 do
              begin
                   for j:= i+2 to N-1 do
                        begin
                             Save_L(L[i-1,j-1], i,j);
                        end;
               end;
    end;


Procedure Reduce_Mat(In_node,Out_node:integer; N:integer; var A,B:matrix);

{ Reduces Beginning N*N Floating Admittance Matrix to an Equivalent Two - }
{ Port Representation by Proper Matrix Reduction Techniques Discussed in  }
{ the References Sited.  After the Necessary Swapping to Place Input Node }
{ in Columns/Row 1 and Output Node in 2, Reference Node in 3, The Starting}
{ Matrix is Reduced Down to a Two-Port Representation.                    }

var      Den      : real; { Magnitude of Yrr                              }
         r        : integer; { Row & Column Being Reduced                 }
         i,j      : integer;
         LL       : integer; { Largest Diagonal Element Row - Better      }
         XX1      : real; { Real Part of Common Term Subtracted From      }
                          { Yps in Reduction Process                      }
         YY1      : real; { Imag Part of Common Term Subtracted From      }
                          { Yps in Reduction Process                      }
         S1,S     : real; { Tests for Largest Diagonal Element; Done to   }
                          { Retain as Much Accuracy as Possible in the    }
                          { Reduction Process                             }


    begin
         if (In_node <> 1) or (Out_node <> 2) then   { Swap In & Out }
         begin
              for j:= 1 to N do
                   begin
                   Swap(A[1,j], A[In_node,j]);
                   Swap(B[1,j], B[In_node,j]);
                   Swap(A[2,j], A[Out_node,j]);
                   Swap(B[2,j], B[Out_node,j]);
                   end;
              for i:= 1 to N do
                   begin
                   Swap(A[i,1], A[i,In_node]);
                   Swap(B[i,1], B[i,In_node]);
                   Swap(A[i,2], A[i,Out_node]);
                   Swap(B[i,2], B[i,Out_node]);
                   end;
         end;
    while N > 2 do { Finds Largest Diagonal Element in Rows Remaining     }
                   { to be Reduced to Enhance Accuracy and Performs the   }
                   { Necessary Swapping                                   }
         begin
         LL:= 3;
         S:= abs( A[LL,LL] ) + abs( B[LL,LL] );
         if N > 3 then
              begin
              for i:= 4 to N do
                   begin
                   S1:= abs( A[i,i] ) + abs( B[i,i] );
                   if S1 > S then
                        begin
                        S:= S1;
                        LL:= i;
                        end;
                   end;
              end;
         if LL <> N then
              begin
              for j:= 1 to N do
                   begin
                   Swap( A[LL,j],A[N,j]);
                   Swap( B[LL,j],B[N,j]);
                   end;
              for i:= 1 to N do
                   begin
                   Swap( A[i,LL], A[i,N]);
                   Swap( B[i,LL],B[i,N]);
                   end;
              end;
         Den:= sqr(A[N,N]) + sqr(B[N,N]); { Beginning of Main Reduction }
         r:= N;
         N:= N - 1;
         for i:= 1 to N do
              begin
              if (A[i,r] <> 0) or (B[i,r] <> 0) then
                   begin
                   XX1:= (A[i,r]*A[r,r] + B[r,r]*B[i,r])/Den;
                   YY1:= (A[r,r]*B[i,r] - A[i,r]*B[r,r])/Den;
                   for j:= 1 to N do
                        begin
                        if (A[r,j]<> 0) or (B[r,j]<> 0) then
                             begin
                             A[i,j]:= A[i,j] - A[r,j]*XX1 + B[r,j]*YY1;
                             B[i,j]:= B[i,j] - A[r,j]*YY1 - B[r,j]*XX1;
                             end;
                        end;
                   end;
              end;
         end;
end;

Procedure Out_put(Ans2: stat; var FYr,FYi:matrix );

begin
    write('  ',F:6:3,'  ',Fyr[1,1]:6:4,'  ',Fyi[1,1]:6:4,'  ',Fyr[1,2]:6:4);
    write('  ',Fyi[1,2]:6:4,'  ',Fyr[2,1]:6:4,'  ',Fyi[2,1]:6:4);
    writeln('  ',Fyr[2,2]:6:4,'  ',Fyi[2,2]:6:4);
    if Ans2 = 'Y' then
    begin
    write(LST,'  ',F:6:3,'  ',Fyr[1,1]:6:4,'  ',Fyi[1,1]:6:4,'  ',Fyr[1,2]:6:4);
    write(LST,'  ',Fyi[1,2]:6:4,'  ',Fyr[2,1]:6:4,'  ',Fyi[2,1]:6:4);
    writeln(LST,'  ',Fyr[2,2]:6:4,'  ',Fyi[2,2]:6:4);
    end;
end;


{ //////////////////////////////////////////////////////////////////// }
{ //////////////////////         MAIN              /////////////////// }
{ //////////////////////        BLOCK              /////////////////// }
{ //////////////////////////////////////////////////////////////////// }

begin
No_turns:= 4;
L[1,2]:= 3.08175E-11;
L[1,3]:= 2.54812E-11;
L[1,4]:= 2.27640E-11;
L[2,3]:= 4.78563E-11;
L[2,4]:= 3.85507E-11;
L[3,4]:= 6.57417E-11;
L[1,1]:= 7.7183E-11;
L[2,2]:= 1.2684E-10;
L[3,3]:= 1.8074E-10;
L[4,4]:= 2.378E-10;
cap[1,2]:= 6.21531E-15;
cap[2,3]:= 8.56065E-15;
cap[3,4]:= 1.0906E-14;
cap[2,0]:= 1.41149E-15;
cap[3,0]:= 1.85961E-15;
cap[1,1]:= 4.02368E-15;
cap[4,4]:= 9.6381E-15;
for i:= 1 to 25 do
    begin
    for j:= 1 to 25 do
         begin
         Fyr[i,j]:= 0.0;
         Fyi[i,j]:= 0.0;
         A[i,j]:= 0.0;
         B[i,j]:= 0.0;
         end;
    end;
Ans2:= 'N';
clrscr;
gotoXY(10,10);
write('Enter Lower, Upper, & Freq-Increment,  GHz  ');
readln(F_low, F_high, F_inc);
F:= F_low;
writeln;
write('         Is a Hard-Copy Desired  Y / N ');
readln(Ans2);
Ans2:= upcase(Ans2);
Nodes:= No_turns + 2;
Num:= 0;
clrscr;
gotoXY(1,2);
write('                         Two - Port Admittance Parameters ');
writeln;
write('  Freq           Y11             Y12    ');
writeln('         Y21             Y22');
write('  GHz        Real    Imag     Real Imag  ');
writeln('    Real  Imag      Real  Imag');
writeln;
if Ans2 = 'Y' then
    begin
    writeln(LST,'                   Two - Port Admittance Parameters ');
    writeln(LST);
    write(LST,'  Freq           Y11              Y12    ');
    writeln(LST,'         Y21             Y22');
    write(LST,'  GHz        Real    Imag      Real Imag  ');
    writeln(LST,'    Real  Imag      Real  Imag');
    writeln(LST);
    end;
while F<= F_high do
    begin
    Num:= Num + 1;
    Omega:= F*2*Pi*1.0E9;
    N:= Nodes;
    for i:= 1 to N do
         begin
         for j:= 1 to N do
              begin
              Fyr[i,j]:= 0.0;
              Fyi[i,j]:= 0.0;
              end;
         end;
    Fill_Caps(N);
    Fill_inductors(N);
    Reduce_Mat(1,Nodes,N,Fyr,Fyi);
    Freq[Num]:= F;
    Mag[Num]:= -1.0/Fyi[1,1];
    Out_put(Ans2,Fyr,Fyi);
    F:= F + F_inc;
    end;
    write('                                              ');
    writeln('          Any Key to Continue ');
    repeat
        delay(100);
    until keypressed;
    clrscr;
    gotoXY(10,10);
    write('Store Data to Disk ?  Y/N ');
    readln(Ans);
    Ans:= upcase(Ans);
    if Ans = 'Y' then Write_file(Freq,Mag);
end.