program Project1;

uses
  Forms;

type
       complex = record
                re : real;
                im : real;
        end;{for compled number}

        data_array = record
          Ch    : char;
          n1 : integer;
          n2 : integer;
          v1 : real;
          v2: real;
          ree : real;
          imm : real;
          coeff : complex;
        end; {of record definition}



        frequency = record
            start : real;
            stop  : real;
            step  : real;
            no    : integer;
        end;{frequency completed }

       matrixtable = array [0..9,0..10] of complex;{Matrix Table}
       dataTable = array [0..20] of data_array; { data of the matrix table }
       answer = array [0..11] of complex;
var
    matrix : matrixtable;
    freq : frequency;
    outfile : TextFile;
    maxNode,p,i : integer;
    data : datatable;
    ansset:answer;
    ans:complex;
    freq1 : real;
    ofile : string[20];

function add(a,b: complex):complex;   {Addition of complex number}
begin
   result.re:=a.re+b.re;              {addition of real parts}
   result.im:=a.im+b.im;              {addition imaginary parts}
end;


function minus(a,b:complex):complex;  {Subtraction of complex number}
begin
   result.re:=a.re-b.re;                  {subtraction of real parts}
   result.im:=a.im-b.im;                  {subtraction of imaginary parts}
end;



function times(a,b:complex):complex; {Multiplication of complex number}
begin
   result.re:=a.re*b.re-a.im*b.im;       {multiplication to get real part}
   result.im:=a.re*b.im+b.re*a.im;       {multiplication to get imaginary part}
end;



function divide (a,b:complex):complex;    {divide of complex number}
begin
  if(b.re*b.re+b.im*b.im)<> 0 then   {check if it something/0 }
  begin
    result.re:=(a.re*b.re+a.im*b.im)/(b.re*b.re+b.im*b.im);{get the real part}
    result.im:=(b.re*a.im-a.re*b.im)/(b.re*b.re+b.im*b.im);{get imaginary part}
  end
  else
  begin
    result.re:=0;
    result.im:=0;
  end;
end; {procedure of divide finished}



function modulus(a:complex):real;
begin
   result:= sqrt(a.re*a.re + a.im*a.im);
end;


function arc_tan(var a:complex):real;   {to find argument of a complex number}
begin
  if a.re=0 then
  begin
    if a.im>0 then
    arc_tan:=0;
    if a.im<0 then
    arc_tan:=pi;
  end
  else
  if a.im=0 then
     begin
      if a.re>0 then
         arc_tan:=pi/2;
      if a.re<0 then
        arc_tan:=-pi/2;
      end
  else
  arc_tan:=arctan(ans.im/ans.re);
end;


procedure readfile(var freq:frequency ;var  data:dataTable;var  p:integer;var   maxNode:integer);
var
 i : integer;
 infile : Textfile;
 ifile : string[20];
begin
  writeln('Pleas specify the name of the input file:  ');
  readln(ifile); 
  assignfile(infile,ifile);
  reset(infile);
   i:=0;
   repeat
    i:= i+1;
    read(infile,data[i].ch);

    case data[i].ch of { compare value of Ch within the file}
    'L','R','C':begin
                read(infile,data[i].n1,data[i].n2,data[i].v1);
                readln(infile);
                writeln(data[i].ch,' ',data[i].n1,' ',data[i].n2,' ',data[i].v1:4:2);
                end;

        'I' :   begin
                read(infile,data[i].n1,data[i].n2,data[i].v1,data[i].v2);
                readln(infile);
                  data[i].ree := (data[i].v1)*(cos(data[i].v2));
                  data[i].imm := (data[i].v1)*(sin(data[i].v2));

                writeln(data[i].ch,' ',data[i].n1,' ',data[i].n2,' ',data[i].v1:4:2,data[i].v2:4:2);
                end;

        'F' :   begin
                read(infile,freq.no,freq.start,freq.stop);
                readln(infile);
                writeln(data[i].ch,' ',freq.no,' ',freq.start:4:2,' ',freq.stop:4:2);
                end;

        'P' :   begin
                read(infile,p);
                readln(infile);
                writeln(data[i].ch,' ',p);
                end;
       end;

     if data[i].n1 > maxNode then {to find the maxnode}
      maxNode:=data[i].n1;
    if data[i].n2 > maxNode then
      maxNode:=data[i].n2;
  until data[i].ch = 'E';
end;



procedure LoadMatrix(maxNode:integer;var matrix:matrixtable; data:DataTable);
var
i,row,column:integer;
begin
  i:=1;
  for row:=0 to 9 do
     for column:=0 to 10 do
    begin
      matrix[row,column].re:=0;
      matrix[row,column].im:=0;
    end; {set the whole matrix to be zero}

while data[i].ch <> 'E' do
  begin
    with data[i] do
    begin
      if data[i].ch = 'I' then {to add currents together}
      begin
        matrix[n1][maxNode+1].re:=matrix[n1][maxNode+1].re-data[i].coeff.re;
        matrix[n2][maxNode+1].re:=matrix[n2][maxNode+1].re+data[i].coeff.re;
        matrix[n1][maxNode+1].im:=matrix[n1][maxNode+1].im-data[i].coeff.im;
        matrix[n2][maxNode+1].im:=matrix[n2][maxNode+1].im+data[i].coeff.im;
      end
      else
         begin {to add component admittances together}
           matrix[n1][n1].re:=matrix[n1][n1].re+data[i].coeff.re;
           matrix[n2][n2].re:=matrix[n2][n2].re+data[i].coeff.re;
           matrix[n1][n2].re:=matrix[n1][n2].re-data[i].coeff.re;
           matrix[n2][n1].re:=matrix[n2][n1].re-data[i].coeff.re;
           matrix[n1][n1].im:=matrix[n1][n1].im+data[i].coeff.im;
           matrix[n2][n2].im:=matrix[n2][n2].im+data[i].coeff.im;
           matrix[n1][n2].im:=matrix[n1][n2].im-data[i].coeff.im;
           matrix[n2][n1].im:=matrix[n2][n1].im-data[i].coeff.im;
         end;
    end;
  i:=i+1;
  end;
end;



Procedure Pivot(Size, StartCol: integer; var matrix: matrixTable);

type holdRow = array[0..10] of complex;

var row, col: integer;
    hold : holdRow;
begin
  col:=0;
  row:=StartCol;

  for row:=i+1 to maxnode do {find the max coefficient}
   begin
    if modulus(matrix[startcol,col]) > modulus(matrix[col,col]) then
      row := col;
   end;
  for col := i to maxnode+1 do
   begin   {change the row}
   hold[Col].re:=matrix[col,row].re;
   hold[Col].im:=matrix[col,row].im;

   matrix[col,row].re:=matrix[col,startcol].re;
   matrix[col,row].im:=matrix[col,startcol].im;

   matrix[col,StartCol].re:=hold[Col].re;
   matrix[col,StartCol].im:=hold[Col].im;

   end;
end;




procedure elimination(maxNode:integer;var matrix:MatrixTable);
var
row,col,i:integer;
ans:complex;
begin
  for i:=1 to maxnode  do
  begin
    Pivot(maxNode,i,matrix);
    for row:=i+1 to maxNode do
    begin
       ans:=divide(matrix[row,i],matrix[i,i]);
      for col:=i+1 to maxNode+1 do
       matrix[row,col]:= minus(matrix[row,col],times(ans,matrix[i,col]));
    end;
  end;
end;


function BackSubstitution(maxnode:integer; matrix:MatrixTable):answer;
var
i,col,row :integer;
sum:complex;
x:answer;
begin
  for i:=maxnode downto 1 do
  begin
    sum.re:=0;
    sum.im:=0;
    for col:=i+1 to maxnode do
    begin
      sum:=add(sum, times(matrix[i,col],x[col]));
    end;
    x[i]:=divide(minus(matrix[i,maxnode+1],sum),matrix[i,i]);

  end;
  result:=x;
end;




procedure calculation(freq:real;var data:dataTable);
var
  i,row : integer;
  w : real;
begin
  i:=1;
  w:=2*pi*freq;
  for row:=1 to 20 do  {set all data.coeff to be zero}
   begin
    data[row].coeff.re:=0;
    data[row].coeff.im:=0;
   end;
 While (data[i].ch <>'E') do {calculate the admittance until the eof}
  begin
    case data[i].ch of
    'L': begin
            data[i].coeff.re := 0;
            data[i].coeff.im := -1/(w*(data[i].v1));
         end;

    'R': begin
            data[i].coeff.re:= 1/data[i].v1;
            data[i].coeff.im:= 0;
         end;

    'C': begin
             data[i].coeff.re := 0;
             data[i].coeff.im := w*(data[i].v1);
         end;

    'I' : begin
            data[i].coeff.re:=data[i].ree;
            data[i].coeff.im:=data[i].imm;
          end;
     end;
    i:=i+1;
  end;

end;



begin
  maxNode:=0;
  readfile(freq,data,p,maxnode);

  freq.step:=(freq.stop-freq.start)/(freq.no-1);

  freq1:=freq.start;

  writeln('Please sepcify the name of the outfile :');
  readln(ofile);
  assignfile(outfile,'ofile');
  rewrite(outfile);
  writeln(outfile,'frequency',' ':8,'magnitude',' ':8,'phase',' ':8);

  for i:=1 to freq.no do
    begin

      calculation(freq1, data);

      loadmatrix(maxnode,matrix,data);

      elimination(maxNode,matrix);

      ansset:=BackSubstitution(maxNode,matrix);

      ans:=ansset[p];

      write(i,'When frequency ',freq1:6:2, 'Hz');
      writeln('the modulus at ',p,' is ',modulus(ans):6:4);
      writeln('and the phase is ',arc_tan(ans):6:2,'in radians');


      writeln(outfile,freq1:6:6,' ':8,modulus(ans):6:6, ' ':8,arc_tan(ans):6:6);

      freq1:=freq1+freq.step;

    end;
  readln;
  closefile(outfile);
end.