home *** CD-ROM | disk | FTP | other *** search
/ Media Share 9 / MEDIASHARE_09.ISO / pascal / shdk_1.zip / SHCMPLX.PAS < prev    next >
Pascal/Delphi Source File  |  1992-03-23  |  10KB  |  354 lines
  1. {$N+,E+}
  2. unit ShCmplx;
  3. {
  4.                                 ShCmplx
  5.  
  6.                        A Complex Arithmetic Unit
  7.  
  8.                                    by
  9.  
  10.                               Bill Madison
  11.  
  12.                    W. G. Madison and Associates, Ltd.
  13.                           13819 Shavano Downs
  14.                             P.O. Box 780956
  15.                        San Antonio, TX 78278-0956
  16.                              (512)492-2777
  17.                              CIS 73240,342
  18.  
  19.                   Copyright 1991 Madison & Associates
  20.                           All Rights Reserved
  21.  
  22.         This file may  be used and distributed  only in accord-
  23.         ance with the provisions described on the title page of
  24.                   the accompanying documentation file
  25.                               SKYHAWK.DOC
  26. }
  27.  
  28. interface
  29.  
  30. type
  31.   ComplexElement  = extended;
  32.   ComplexBaseType = record
  33.                       Re,
  34.                       Im  : ComplexElement;
  35.                       end;
  36.   Complex         = ^ComplexBaseType;
  37.  
  38. function CmplxError  : integer;
  39.  
  40. function Cmp2Str(A : Complex; Width, Places : byte) : string;
  41. {Converts a complex value to a string of the form (Re + Im i)}
  42.  
  43. function CmpP2Str(A : Complex; Width, Places : byte) : string;
  44. {Converts a complex in polar form to a string with the angle in radians.}
  45.  
  46. function CmpP2StrD(A : Complex; Width, Places : byte) : string;
  47. {Converts a complex in polar form to a string with the angle in degrees.}
  48.  
  49. procedure CAbs(A  : Complex; var Result : ComplexElement);
  50. function CAbsF(A  : Complex)  : ComplexElement;
  51. {Returns the absolute value of a complex number}
  52.  
  53. procedure CConj(A : Complex; Result : Complex);
  54. function CConjF(A : Complex) : Complex;
  55. {Returns the complex conjugate of a complex number}
  56.  
  57. procedure CAdd(A, B : Complex; Result : Complex);
  58. function CAddF(A, B : Complex) : Complex;
  59. {Returns the sum of two complex numbers A + B}
  60.  
  61. procedure RxC(A : ComplexElement; B : Complex; Result : Complex);
  62. function RxCF(A : ComplexElement; B : Complex) : Complex;
  63. {Returns the complex product of a real and a complex}
  64.  
  65. procedure CSub(A, B : Complex; Result : Complex);
  66. function CSubF(A, B : Complex) : Complex;
  67. {Returns the difference between two complex numbers A - B}
  68.  
  69. procedure CMul(A, B : Complex; Result : Complex);
  70. function CMulF(A, B : Complex) : Complex;
  71. {Returns the product of two complex numbers A * B}
  72.  
  73. procedure CDiv(A, B : Complex; Result : Complex);
  74. function CDivF(A, B : Complex) : Complex;
  75. {Returns the quotient of two complex numbers A / B}
  76.  
  77. procedure C2P(A : Complex; Result : Complex);
  78. function C2PF(A : Complex) : Complex;
  79. {Transforms a complex in cartesian form into polar form}
  80. {The magnitude will be stored in Result^.Re and the angle in Result^.Im}
  81.  
  82. procedure P2C(A : Complex; Result : Complex);
  83. function P2CF(A : Complex) : Complex;
  84. {Transforms a complex in polar form into cartesian form}
  85. {The magnitude is stored in A^.Re and the angle in A^.Im}
  86.  
  87. procedure CpPwrR(A : Complex; R : ComplexElement; Result : Complex);
  88. function CpPwrRF(A : Complex; R : ComplexElement) : Complex;
  89. {Raises a complex (in polar form) to a real power}
  90.  
  91. {===========================================================}
  92.  
  93. implementation
  94.  
  95. const
  96.   MaxStack  = 5;
  97.   StackTop  = MaxStack - 1;
  98.   Pi        = 3.14159265358979;
  99.   PiOver2   = Pi / 2.0;
  100.   TwoPi     = Pi * 2.0;
  101.   F180OverPi= 180.0 / Pi;
  102.   SpConj  : byte = 0;
  103.   SpSum   : byte = 0;
  104.   SpDiff  : byte = 0;
  105.   SpRProd : byte = 0;
  106.   SpProd  : byte = 0;
  107.   SpQuot  : byte = 0;
  108.   SpPolar : byte = 0;
  109.   SpCartsn: byte = 0;
  110.   SpPower : byte = 0;
  111.  
  112. var
  113.   Conj,
  114.   Sum,
  115.   Diff,
  116.   RProd,
  117.   Prod,
  118.   Quot,
  119.   Polar,
  120.   Cartsn,
  121.   Power    : array[0..StackTop] of Complex;
  122.   CmpError : integer;
  123.  
  124. function CmplxError  : integer;
  125.   begin
  126.     CmplxError := CmpError;
  127.     CmpError := 0;
  128.     end;
  129.  
  130. function Real2Str(A : ComplexElement; Width, Places : byte) : string;
  131.   var
  132.     T1  : string;
  133.   begin
  134.     Str(A:Width:Places, T1);
  135.     Real2Str := T1;
  136.     end;
  137.  
  138. function Cmp2Str(A : Complex; Width, Places : byte) : string;
  139.   begin
  140.     if A^.Im >= 0.0 then
  141.       Cmp2Str := '(' + Real2Str(A^.Re, Width, Places) + ' + ' +
  142.                        Real2Str(A^.Im, Width, Places) + 'i)'
  143.     else
  144.       Cmp2Str := '(' + Real2Str(A^.Re, Width, Places) + ' - ' +
  145.                        Real2Str(Abs(A^.Im), Width, Places) + 'i)';
  146.     end;
  147.  
  148. function CmpP2StrD(A : Complex; Width, Places : byte) : string;
  149. {Converts a complex in polar form to a string with the angle in degrees.}
  150.   begin
  151.     CmpP2StrD := '('+Real2Str(A^.Re,Width,Places)+' at '+
  152.       Real2Str((A^.Im*F180OverPi),Width,Places)+#248')';
  153.     end;
  154.  
  155. function CmpP2Str(A : Complex; Width, Places : byte) : string;
  156. {Converts a complex in polar form to a string with the angle in radians.}
  157.   begin
  158.     CmpP2Str := '('+Real2Str(A^.Re,Width,Places)+' at '+
  159.       Real2Str((A^.Im),Width,Places)+' rad)';
  160.     end;
  161.  
  162. procedure IncPtr(var StackPtr : byte);
  163.   begin
  164.     StackPtr := (StackPtr + 1) mod StackTop;
  165.     end;
  166.  
  167. function CAbsF(A  : Complex) : ComplexElement;
  168.   begin
  169.     CmpError := 0;
  170.     CAbsF := sqrt(A^.Re * A^.Re + A^.Im * A^.Im);
  171.     end;
  172.  
  173. procedure CAbs(A  : Complex; var Result : ComplexElement);
  174.   begin
  175.     Result := CAbsF(A);
  176.     end;
  177.  
  178. function CConjF(A : Complex) : Complex;
  179.   begin
  180.     CmpError := 0;
  181.     Conj[SpConj]^.Re := A^.Re;
  182.     Conj[SpConj]^.Im := -A^.Im;
  183.     CConjF := Conj[SpConj];
  184.     IncPtr(SpConj);
  185.     end;
  186.  
  187. procedure CConj(A : Complex; Result : Complex);
  188.   begin
  189.     Result^ := CConjF(A)^;
  190.     end;
  191.  
  192. function CAddF(A, B : Complex) : Complex;
  193.   begin
  194.     CmpError := 0;
  195.     Sum[SpSum]^.Re := A^.Re + B^.Re;
  196.     Sum[SpSum]^.Im := A^.Im + B^.Im;
  197.     CAddF := Sum[SpSum];
  198.     IncPtr(SpSum);
  199.     end;
  200.  
  201. procedure CAdd(A, B : Complex; Result : Complex);
  202.   begin
  203.     Result^ := CAddF(A, B)^;
  204.     end;
  205.  
  206. function RxCF(A : ComplexElement; B : Complex) : Complex;
  207.   begin
  208.     CmpError := 0;
  209.     RProd[SpRProd]^.Re := A * B^.Re;
  210.     RPRod[SpRProd]^.Im := A * B^.Im;
  211.     RxCF := RProd[SpRProd];
  212.     IncPtr(SpRProd);
  213.     end;
  214.  
  215. procedure RxC(A : ComplexElement; B : Complex; Result : Complex);
  216.   begin
  217.     Result^ := RxCF(A, B)^;
  218.     end;
  219.  
  220. function CSubF(A, B : Complex) : Complex;
  221.   begin
  222.     CmpError := 0;
  223.     Diff[SpDiff]^ := CAddF(A, RxCF(-1.0,B))^;
  224.     CSubF := Diff[SpDiff];
  225.     IncPtr(SpDiff);
  226.     end;
  227.  
  228. procedure CSub(A, B : Complex; Result : Complex);
  229.   begin
  230.     Result^ := CSubF(A, B)^;
  231.     end;
  232.  
  233. function CMulF(A, B : Complex) : Complex;
  234.   begin
  235.     CmpError := 0;
  236.     Prod[SpProd]^.Re := A^.Re * B^.Re - A^.Im * B^.Im;
  237.     Prod[SpProd]^.Im := A^.Im * B^.Re + A^.Re * B^.Im;
  238.     CMulF := Prod[SpProd];
  239.     IncPtr(SpProd);
  240.     end;
  241.  
  242. procedure CMul(A, B : Complex; Result : Complex);
  243.   begin
  244.     Result^ := CMulF(A, B)^;
  245.     end;
  246.  
  247. function CDivF(A, B : Complex) : Complex;
  248.   var
  249.     Denom : ComplexElement;
  250.   begin
  251.     CmpError := 0;
  252.     Denom := B^.Re * B^.Re + B^.Im * B^.Im;
  253.     if Denom = 0.0 then begin
  254.       CmpError := -1;
  255.       Quot[SpQuot]^.Re := 0.0;
  256.       Quot[SpQuot]^.Im := 0.0;
  257.       CDivF := Quot[SpQuot];
  258.       IncPtr(SpQuot);
  259.       exit;
  260.       end;
  261.     Quot[SpQuot]^ := RxCF(1.0 / Denom, CMulF(A, CConjF(B)))^;
  262.     CDivF := Quot[SpQuot];
  263.     IncPtr(SpQuot);
  264.     end;
  265.  
  266. procedure CDiv(A, B : Complex; Result : Complex);
  267.   begin
  268.     Result^ := CDivF(A, B)^;
  269.     end;
  270.  
  271. procedure C2P(A : Complex; Result : Complex);
  272. {Transforms a complex in cartesian form into polar form}
  273. {The magnitude will be stored in Result^.Re and the angle in Result^.Im}
  274.   begin
  275.     CmpError := 0;
  276.     Result^.Re := CAbsF(A);
  277.     if abs(A^.Re) > abs(A^.Im) then begin
  278.       {Use the tangent formulation}
  279.       Result^.Im := ArcTan(abs(A^.Im) / abs(A^.Re));
  280.       end
  281.     else begin
  282.       {Use the cotangent formulation}
  283.       Result^.Im := PiOver2 - ArcTan(abs(A^.Re) / abs(A^.Im));
  284.       end;
  285.     if A^.Im > 0 then begin {first or second quadrant}
  286.       if A^.Re < 0.0 then {Second quadrant}
  287.         Result^.Im := Pi - abs(Result^.Im);
  288.       end
  289.     else begin {third or fourth quadrant}
  290.       if A^.Re > 0.0 then {Fourth quadrant}
  291.         Result^.Im := TwoPi - abs(Result^.Im)
  292.       else begin {Third Quadrant}
  293.         Result^.Im := Pi + abs(Result^.Im);
  294.         end;
  295.       end;
  296.     if Result^.Im = TwoPi then Result^.Im := 0.0;
  297.     end;
  298.  
  299. function C2PF(A : Complex) : Complex;
  300.   begin
  301.     C2P(A, Polar[SpPolar]);
  302.     C2PF := Polar[SpPolar];
  303.     IncPtr(SpPolar);
  304.     end;
  305.  
  306. procedure P2C(A : Complex; Result : Complex);
  307. {Transforms a complex in polar form into cartesian form}
  308. {The magnitude is stored in A^.Re and the angle in A^.Im}
  309.   begin
  310.     CmpError := 0;
  311.     Result^.Re := A^.Re * cos(A^.Im);
  312.     Result^.Im := A^.Re * sin(A^.Im);
  313.     end;
  314.  
  315. function P2CF(A : Complex) : Complex;
  316.   begin
  317.     P2C(A, Cartsn[SpCartsn]);
  318.     P2CF := Cartsn[SpCartsn];
  319.     IncPtr(SpCartsn);
  320.     end;
  321.  
  322. function CpPwrRF(A : Complex; R : ComplexElement) : Complex;
  323. {Raises a complex (in polar form) to a real power, returning the result
  324.  also in polar form}
  325.   begin
  326.     CmpError := 0;
  327.     if A^.Re <= 0.0 then begin
  328.       CmpError := -2;
  329.       Power[SpPower]^.Re := 0.0;
  330.       Power[SpPower]^.Im := 0.0;
  331.       exit;
  332.       end;
  333.     Power[SpPower]^.Re := exp(R * ln(A^.Re));
  334.     Power[SpPower]^.Im := R * A^.Im;
  335.     CpPwrRF := Power[SpPower];
  336.     IncPtr(SpPower);
  337.     end;
  338.  
  339. procedure CpPwrR(A : Complex; R : ComplexElement; Result : Complex);
  340.   begin
  341.     Result^ := CpPwrRF(A, R)^;
  342.     end;
  343.  
  344. var
  345.   T1  : byte;
  346.  
  347. begin
  348.   for T1 := 0 to StackTop do begin
  349.     New(Conj[T1]); New(Sum[T1]); New(Diff[T1]);
  350.     New(RProd[T1]); New(Prod[T1]); New(Quot[T1]);
  351.     New(Polar[T1]); New(Cartsn[T1]); New(Power[T1]);
  352.     end;
  353.   end.
  354.