home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / PASCAL / NUMPAS.ZIP / PARFRACT.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-06-25  |  3.3 KB  |  127 lines
  1. program parfract(input,output);
  2. (* Program evaluates a rational function N(s)/D(s) for complex s. *)
  3.  
  4. type vector = array[0..10] of real;
  5. var coefn,coefd:vector;
  6.     num1,num2,denom1,denom2,z1,z2,magr,magj,magnitude,angle: real;
  7.     mn,md: integer;  (*  degree of numerator and denominator  *)
  8.     q: char;
  9.  
  10. procedure readin(var a1: vector; n1: integer);
  11. (* Reads in polynomial coefficients from the keyboard. *)
  12.  
  13. var i: integer;
  14.  
  15. begin
  16.    for i := 0 to n1 do
  17.       begin
  18.          write('Input the ',i,' coefficent: ');
  19.          readln(a1[i]);
  20.       end;
  21.    writeln('The polynomial coefficients are:');
  22.    writeln;
  23.    for i := 0 to n1 do
  24.       begin
  25.          write(a1[i]:6:3,'  ');
  26.       end;
  27.    writeln;
  28. end;
  29.  
  30. procedure evalpoly(coef: vector; w1,w2: real; var b1,b2: real; n1: integer);
  31. (* Evaluates a polynomial N(s) in the form *)
  32. (*  A(0) + A(1)s + A(2)s^2 + ...+ A(n1)s^n1 for  s = w1 + jw2. *)
  33.  
  34. var a1,a2: vector;
  35.     i: integer;
  36.  
  37. begin
  38.    a1[0] := 1.0;
  39.    a2[0] := 0.0;
  40.    b1 := a1[0]*coef[0];
  41.    b2 := a2[0]*coef[0];
  42.    for i := 1 to n1 do
  43.       begin
  44.          a1[i] := w1*a1[i-1] - w2*a2[i-1];  (* real part of w^i *)
  45.          a2[i] := w1*a2[i-1] + w2*a1[i-1];  (* imag part of w^i *)
  46.          b1 := b1 + a1[i]*coef[i];   (* real part of trunc polynomial *)
  47.          b2 := b2 + a2[i]*coef[i];   (* imag part of trunc polynomial *)
  48.       end;
  49. end;
  50.  
  51. procedure evaltf(a1,a2,b1,b2: real; var d1,d2: real);
  52. (* Evaluates a complex quotient d1 + jd2 = (a1 + ja2)/(b1 + jb2). *)
  53.  
  54. var d: real;
  55.  
  56. begin
  57.    d := b1*b1 + b2*b2;
  58.    d1 := (a1*b1 + a2*b2)/d;
  59.    d2 := (a2*b1 - a1*b2)/d;
  60. end;
  61.  
  62. procedure correct(var d1,d2,w2: real);
  63. (* Corrects partial-fraction calculation for a complex pair of poles. *)
  64.  
  65. var d11: real;
  66.  
  67. begin
  68.    d11 := d2/(2*w2);
  69.    d2 := -d1/(2*w2);
  70.    d1 := d11;
  71. end;
  72.  
  73. procedure rectpol(a1,a2: real; var cmag1,angle1: real);
  74. (* Converts rectangular form a1 + ja2 to polar form. *)
  75. var eta: real;
  76.  
  77. begin
  78.    cmag1 := sqrt(a1*a1 + a2*a2);
  79.    if (a1 = 0.0) and (a2 = 0.0) then
  80.       angle1 := 0.0;
  81.    if (a1 = 0.0) and (a2 > 0.0) then
  82.       angle1 := 90.0;
  83.    if (a1 = 0.0) and (a2 < 0.0) then
  84.       angle1 := -90.0;
  85.    if abs(a1) > 0.0 then
  86.       eta := arctan(a2/a1)*57.296;
  87.    if a1 > 0 then
  88.       angle := eta;
  89.    if (a1 < 0.0) and (a2 >= 0.0) then
  90.       angle1 := eta + 180.0;
  91.    if (a1 < 0.0) and (a2 < 0.0) then
  92.       angle1 := eta - 180.0;
  93. end;
  94.  
  95. begin   (*  main  *)
  96.    write('Input degree of numerator polynomial: ');
  97.    readln(mn);
  98.    writeln;
  99.    write('Input degree of denominator polynomial: ');
  100.    readln(md);
  101.    writeln;
  102.    writeln('Input numerator coefficients as prompted:');
  103.    writeln;
  104.    readin(coefn,mn);
  105.    writeln('Input denominator coefficients as prompted:');
  106.    readin(coefd,md);
  107.    write('Input complex value of z (z1 z2): ');
  108.    readln(z1,z2);
  109.    evalpoly(coefn,z1,z2,num1,num2,mn);
  110.    evalpoly(coefd,z1,z2,denom1,denom2,md);
  111.    evaltf(num1,num2,denom1,denom2,magr,magj);
  112.    writeln('Is this a P-F expansion for a complex root? (y,n):');
  113.    writeln;
  114.    readln(q);
  115.    if q = 'y'
  116.       then correct(magr,magj,z2);
  117.    rectpol(magr,magj,magnitude,angle);
  118.    writeln('The real and imaginary parts are: ');
  119.    writeln(magr:9,', ',magj:9);
  120.    writeln;
  121.    writeln('The magnitude and angle are: ');
  122.    writeln(magnitude:9,', ',angle:5:2,' degrees');
  123. end.
  124.  
  125.  
  126.  
  127.