home *** CD-ROM | disk | FTP | other *** search
/ DP Tool Club 15 / CD_ASCQ_15_070894.iso / vrac / sk210f.zip / SHCMPLX.PAS < prev    next >
Pascal/Delphi Source File  |  1994-05-15  |  13KB  |  468 lines
  1. {$I SHDEFINE.INC}
  2.  
  3. {$I SHUNITSW.INC}
  4.  
  5. {$D-,L-}
  6. unit ShCmplx;
  7. {
  8.                                 ShCmplx
  9.  
  10.                        A Complex Arithmetic Unit
  11.  
  12.                                    by
  13.  
  14.                               Bill Madison
  15.  
  16.                    W. G. Madison and Associates, Ltd.
  17.                           13819 Shavano Downs
  18.                             P.O. Box 780956
  19.                        San Antonio, TX 78278-0956
  20.                              (512)492-2777
  21.                              CIS 73240,342
  22.                 Internet bill.madison@lchance.sat.tx.us
  23.  
  24.                 Copyright 1990, '94 Madison & Associates
  25.                           All Rights Reserved
  26.  
  27.         This file may  be used and distributed  only in accord-
  28.         ance with the provisions described on the title page of
  29.                   the accompanying documentation file
  30.                               SKYHAWK.DOC
  31. }
  32.  
  33. interface
  34.  
  35. uses
  36. TPCRT,
  37. TPDOS,
  38. TPSTRING,
  39.   TpInline,
  40.   ShErrMsg;
  41.  
  42. const
  43.   Copyr = 'Copyright 1990, 1994 by W.G. Madison';
  44.  
  45. type
  46. {$IFNDEF Gen87}
  47.   extended        = real;
  48.   Float           = real;
  49. {$ELSE}
  50.   Float           = extended;
  51. {$ENDIF}
  52.  
  53.   ComplexElement  = Float;
  54.   ComplexBaseType = record
  55.                       Re,
  56.                       Im  : ComplexElement;
  57.                       end;
  58.   Complex         = ^ComplexBaseType;
  59.  
  60.   PFreeRec        = ^TFreeRec;
  61.   TFreeRec        = record
  62.                       Next  : PFreeRec;
  63.                       Size  : pointer;
  64.                       end;
  65.  
  66. const
  67.   cxHaltErr     : boolean        = true;    {Stop program execution on
  68.                                               non-I/O errors}
  69.   cxOK          = 0;
  70.   cxDivByZero   = 200;
  71.   cxBadMagnitude= 201;
  72.  
  73. procedure CmplxInit;
  74. {Initializes the complex arithmetic package.}
  75.  
  76. procedure CmplxDeInit;
  77. {De-initializes the complex arithmetic package, releasing all the ring
  78.  buffer heap space.}
  79.  
  80. function CmplxError  : integer;
  81.  
  82. function Cmp2Str(A : Complex; Width, Places : byte) : string;
  83. {Converts a complex value to a string of the form (Re + Im i)}
  84.  
  85. function CmpP2Str(A : Complex; Width, Places : byte) : string;
  86. {Converts a complex in polar form to a string with the angle in radians.}
  87.  
  88. function CmpP2StrD(A : Complex; Width, Places : byte) : string;
  89. {Converts a complex in polar form to a string with the angle in degrees.}
  90.  
  91. procedure CAbs(A  : Complex; var Result : ComplexElement);
  92. function CAbsF(A  : Complex)  : ComplexElement;
  93. {Returns the absolute value of a complex number}
  94.  
  95. procedure CConj(A : Complex; Result : Complex);
  96. function CConjF(A : Complex) : Complex;
  97. {Returns the complex conjugate of a complex number}
  98.  
  99. procedure CAdd(A, B : Complex; Result : Complex);
  100. function CAddF(A, B : Complex) : Complex;
  101. {Returns the sum of two complex numbers A + B}
  102.  
  103. procedure RxC(A : ComplexElement; B : Complex; Result : Complex);
  104. function RxCF(A : ComplexElement; B : Complex) : Complex;
  105. {Returns the complex product of a real and a complex}
  106.  
  107. procedure CSub(A, B : Complex; Result : Complex);
  108. function CSubF(A, B : Complex) : Complex;
  109. {Returns the difference between two complex numbers A - B}
  110.  
  111. procedure CMul(A, B : Complex; Result : Complex);
  112. function CMulF(A, B : Complex) : Complex;
  113. {Returns the product of two complex numbers A * B}
  114.  
  115. procedure CDiv(A, B : Complex; Result : Complex);
  116. function CDivF(A, B : Complex) : Complex;
  117. {Returns the quotient of two complex numbers A / B}
  118.  
  119. procedure C2P(A : Complex; Result : Complex);
  120. function C2PF(A : Complex) : Complex;
  121. {Transforms a complex in cartesian form into polar form}
  122. {The magnitude will be stored in Result^.Re and the angle in Result^.Im}
  123.  
  124. procedure P2C(A : Complex; Result : Complex);
  125. function P2CF(A : Complex) : Complex;
  126. {Transforms a complex in polar form into cartesian form}
  127. {The magnitude is stored in A^.Re and the angle in A^.Im}
  128.  
  129. procedure CpPwrR(A : Complex; R : ComplexElement; Result : Complex);
  130. function CpPwrRF(A : Complex; R : ComplexElement) : Complex;
  131. {Raises a complex (in polar form) to a real power}
  132.  
  133. {===========================================================}
  134.  
  135. implementation
  136.  
  137. VAR OUTPUT : TEXT;
  138.  
  139. const
  140.   MaxStack            = 5;
  141.   StackTop            = MaxStack - 1;
  142.   Pi                  = 3.14159265358979;
  143.   PiOver2             = Pi / 2.0;
  144.   TwoPi               = Pi * 2.0;
  145.   F180OverPi          = 180.0 / Pi;
  146.   SpConj    : byte    = 0;
  147.   SpSum     : byte    = 0;
  148.   SpDiff    : byte    = 0;
  149.   SpRProd   : byte    = 0;
  150.   SpProd    : byte    = 0;
  151.   SpQuot    : byte    = 0;
  152.   SpPolar   : byte    = 0;
  153.   SpCartsn  : byte    = 0;
  154.   SpPower   : byte    = 0;
  155.   CmpInited : boolean = false;
  156.  
  157. type
  158.   StackArray  = array[0..StackTop] of Complex;
  159.  
  160. var
  161.   Conj,
  162.   Sum,
  163.   Diff,
  164.   RProd,
  165.   Prod,
  166.   Quot,
  167.   Polar,
  168.   Cartsn,
  169.   Power    : StackArray;
  170.   CmpError : integer;
  171.  
  172. procedure ChkInit;
  173.   begin
  174.     if not CmpInited then
  175.       RunErrorMsg(400, 'Unit ShCmplx not initialized.');
  176.     end;
  177.  
  178. procedure CmplxInit;
  179.   var
  180.     T1  : byte;
  181.   begin {CmplxInit}
  182.     for T1 := 0 to StackTop do begin
  183.       New(Cartsn[T1]);
  184.       New(Conj[T1]);
  185.       New(Diff[T1]);
  186.       New(Polar[T1]);
  187.       New(Power[T1]);
  188.       New(Prod[T1]);
  189.       New(Quot[T1]);
  190.       New(RProd[T1]);
  191.       New(Sum[T1]);
  192.       end;
  193.     CmpInited := true;
  194.     end; {CmplxInit}
  195.  
  196. procedure CmplxDeInit;
  197. {De-initializes the complex arithmetic package, releasing all the ring
  198.  buffer heap space.}
  199.   var
  200.     T1  : byte;
  201.   begin {CmplxDeInit}
  202.     for T1 := StackTop downto 0 do begin
  203.       {Flush the heap space}
  204.       if    Sum[T1] <> nil then
  205.         Dispose(Sum[T1]);
  206.       if  RProd[T1] <> nil then
  207.         Dispose(RProd[T1]);
  208.       if   Quot[T1] <> nil then
  209.         Dispose(Quot[T1]);
  210.       if   Prod[T1] <> nil then
  211.         Dispose(Prod[T1]);
  212.       if  Power[T1] <> nil then
  213.         Dispose(Power[T1]);
  214.       if  Polar[T1] <> nil then
  215.         Dispose(Polar[T1]);
  216.       if   Diff[T1] <> nil then
  217.         Dispose(Diff[T1]);
  218.       if   Conj[T1] <> nil then
  219.         Dispose(Conj[T1]);
  220.       if Cartsn[T1] <> nil then
  221.         Dispose(Cartsn[T1]);
  222.       {Reset the heap pointers}
  223.       Conj[T1] := nil;  Sum[T1] := nil;     Diff[T1] := nil;
  224.       RProd[T1] := nil; Prod[T1] := nil;    Quot[T1] := nil;
  225.       Polar[T1] := nil; Cartsn[T1] := nil;  Power[T1] := nil;
  226.       end;
  227.     {Reset the ring buffer pointers}
  228.     SpConj := 0;  SpSum := 0;     SpDiff := 0;
  229.     SpRProd := 0; SpProd := 0;    SpQuot := 0;
  230.     SpPolar := 0; SpCartsn := 0;  SpPower := 0;
  231.     CmpInited := false;
  232.     end; {CmplxDeInit}
  233.  
  234. function CmplxError  : integer;
  235.   begin
  236.     CmplxError := CmpError;
  237.     CmpError := cxOK;
  238.     end;
  239.  
  240. function Real2Str(A : ComplexElement; Width, Places : byte) : string;
  241.   var
  242.     T1  : string;
  243.   begin
  244.     Str(A:Width:Places, T1);
  245.     Real2Str := T1;
  246.     end;
  247.  
  248. function Cmp2Str(A : Complex; Width, Places : byte) : string;
  249.   begin
  250.     if A^.Im >= 0.0 then
  251.       Cmp2Str := '(' + Real2Str(A^.Re, Width, Places) + ' + ' +
  252.                        Real2Str(A^.Im, Width, Places) + 'i)'
  253.     else
  254.       Cmp2Str := '(' + Real2Str(A^.Re, Width, Places) + ' - ' +
  255.                        Real2Str(Abs(A^.Im), Width, Places) + 'i)';
  256.     end;
  257.  
  258. function CmpP2StrD(A : Complex; Width, Places : byte) : string;
  259. {Converts a complex in polar form to a string with the angle in degrees.}
  260.   begin
  261.     CmpP2StrD := '('+Real2Str(A^.Re,Width,Places)+' at '+
  262.       Real2Str((A^.Im*F180OverPi),Width,Places)+#248')';
  263.     end;
  264.  
  265. function CmpP2Str(A : Complex; Width, Places : byte) : string;
  266. {Converts a complex in polar form to a string with the angle in radians.}
  267.   begin
  268.     CmpP2Str := '('+Real2Str(A^.Re,Width,Places)+' at '+
  269.       Real2Str((A^.Im),Width,Places)+' rad)';
  270.     end;
  271.  
  272. procedure IncPtr(var StackName : StackArray; var StackPtr : byte);
  273.   begin
  274.     StackPtr := (StackPtr + 1) mod StackTop;
  275.     if StackName[StackPtr] <> nil then begin
  276.       Dispose(StackName[StackPtr]);
  277.       end;
  278.     New(StackName[StackPtr]);
  279.     end;
  280.  
  281. function CAbsF(A  : Complex) : ComplexElement;
  282.   begin
  283.     CmpError := cxOK;
  284.     CAbsF := sqrt(A^.Re * A^.Re + A^.Im * A^.Im);
  285.     end;
  286.  
  287. procedure CAbs(A  : Complex; var Result : ComplexElement);
  288.   begin
  289.     Result := CAbsF(A);
  290.     end;
  291.  
  292. function CConjF(A : Complex) : Complex;
  293.   begin
  294.     CmpError := cxOK;
  295.     ChkInit;
  296.     Conj[SpConj]^.Re := A^.Re;
  297.     Conj[SpConj]^.Im := -A^.Im;
  298.     CConjF := Conj[SpConj];
  299.     IncPtr(Conj, SpConj);
  300.     end;
  301.  
  302. procedure CConj(A : Complex; Result : Complex);
  303.   begin
  304.     Result^ := CConjF(A)^;
  305.     end;
  306.  
  307. function CAddF(A, B : Complex) : Complex;
  308.   begin
  309.     CmpError := cxOK;
  310.     ChkInit;
  311.     ChkInit;
  312.     Sum[SpSum]^.Re := A^.Re + B^.Re;
  313.     Sum[SpSum]^.Im := A^.Im + B^.Im;
  314.     CAddF := Sum[SpSum];
  315.     IncPtr(Sum, SpSum);
  316.     end;
  317.  
  318. procedure CAdd(A, B : Complex; Result : Complex);
  319.   begin
  320.     Result^ := CAddF(A, B)^;
  321.     end;
  322.  
  323. function RxCF(A : ComplexElement; B : Complex) : Complex;
  324.   begin
  325.     CmpError := cxOK;
  326.     ChkInit;
  327.     RProd[SpRProd]^.Re := A * B^.Re;
  328.     RPRod[SpRProd]^.Im := A * B^.Im;
  329.     RxCF := RProd[SpRProd];
  330.     IncPtr(RProd, SpRProd);
  331.     end;
  332.  
  333. procedure RxC(A : ComplexElement; B : Complex; Result : Complex);
  334.   begin
  335.     Result^ := RxCF(A, B)^;
  336.     end;
  337.  
  338. function CSubF(A, B : Complex) : Complex;
  339.   begin
  340.     CmpError := cxOK;
  341.     ChkInit;
  342.     Diff[SpDiff]^ := CAddF(A, RxCF(-1.0,B))^;
  343.     CSubF := Diff[SpDiff];
  344.     IncPtr(Diff, SpDiff);
  345.     end;
  346.  
  347. procedure CSub(A, B : Complex; Result : Complex);
  348.   begin
  349.     Result^ := CSubF(A, B)^;
  350.     end;
  351.  
  352. function CMulF(A, B : Complex) : Complex;
  353.   begin
  354.     CmpError := cxOK;
  355.     ChkInit;
  356.     Prod[SpProd]^.Re := A^.Re * B^.Re - A^.Im * B^.Im;
  357.     Prod[SpProd]^.Im := A^.Im * B^.Re + A^.Re * B^.Im;
  358.     CMulF := Prod[SpProd];
  359.     IncPtr(Prod, SpProd);
  360.     end;
  361.  
  362. procedure CMul(A, B : Complex; Result : Complex);
  363.   begin
  364.     Result^ := CMulF(A, B)^;
  365.     end;
  366.  
  367. function CDivF(A, B : Complex) : Complex;
  368.   var
  369.     Denom : ComplexElement;
  370.   begin
  371.     CmpError := cxOK;
  372.     Denom := B^.Re * B^.Re + B^.Im * B^.Im;
  373.     if Denom = 0.0 then begin
  374.       CmpError := cxDivByZero;
  375.       ChkInit;
  376.       Quot[SpQuot]^.Re := 0.0;
  377.       Quot[SpQuot]^.Im := 0.0;
  378.       CDivF := Quot[SpQuot];
  379.       IncPtr(Quot, SpQuot);
  380.       exit;
  381.       end;
  382.     Quot[SpQuot]^ := RxCF(1.0 / Denom, CMulF(A, CConjF(B)))^;
  383.     CDivF := Quot[SpQuot];
  384.     IncPtr(Quot, SpQuot);
  385.     end;
  386.  
  387. procedure CDiv(A, B : Complex; Result : Complex);
  388.   begin
  389.     Result^ := CDivF(A, B)^;
  390.     end;
  391.  
  392. procedure C2P(A : Complex; Result : Complex);
  393. {Transforms a complex in cartesian form into polar form}
  394. {The magnitude will be stored in Result^.Re and the angle in Result^.Im}
  395.   begin
  396.     CmpError := cxOK;
  397.     Result^.Re := CAbsF(A);
  398.     if abs(A^.Re) > abs(A^.Im) then begin
  399.       {Use the tangent formulation}
  400.       Result^.Im := ArcTan(abs(A^.Im) / abs(A^.Re));
  401.       end
  402.     else begin
  403.       {Use the cotangent formulation}
  404.       Result^.Im := PiOver2 - ArcTan(abs(A^.Re) / abs(A^.Im));
  405.       end;
  406.     if A^.Im > 0 then begin {first or second quadrant}
  407.       if A^.Re < 0.0 then {Second quadrant}
  408.         Result^.Im := Pi - abs(Result^.Im);
  409.       end
  410.     else begin {third or fourth quadrant}
  411.       if A^.Re > 0.0 then {Fourth quadrant}
  412.         Result^.Im := TwoPi - abs(Result^.Im)
  413.       else begin {Third Quadrant}
  414.         Result^.Im := Pi + abs(Result^.Im);
  415.         end;
  416.       end;
  417.     if Result^.Im = TwoPi then Result^.Im := 0.0;
  418.     end;
  419.  
  420. function C2PF(A : Complex) : Complex;
  421.   begin
  422.     ChkInit;
  423.     C2P(A, Polar[SpPolar]);
  424.     C2PF := Polar[SpPolar];
  425.     IncPtr(Polar, SpPolar);
  426.     end;
  427.  
  428. procedure P2C(A : Complex; Result : Complex);
  429. {Transforms a complex in polar form into cartesian form}
  430. {The magnitude is stored in A^.Re and the angle in A^.Im}
  431.   begin
  432.     CmpError := cxOK;
  433.     Result^.Re := A^.Re * cos(A^.Im);
  434.     Result^.Im := A^.Re * sin(A^.Im);
  435.     end;
  436.  
  437. function P2CF(A : Complex) : Complex;
  438.   begin
  439.     ChkInit;
  440.     P2C(A, Cartsn[SpCartsn]);
  441.     P2CF := Cartsn[SpCartsn];
  442.     IncPtr(Cartsn, SpCartsn);
  443.     end;
  444.  
  445. function CpPwrRF(A : Complex; R : ComplexElement) : Complex;
  446. {Raises a complex (in polar form) to a real power, returning the result
  447.  also in polar form}
  448.   begin
  449.     CmpError := cxOK;
  450.     if A^.Re <= 0.0 then begin
  451.       CmpError := cxBadMagnitude;
  452.       ChkInit;
  453.       Power[SpPower]^.Re := 0.0;
  454.       Power[SpPower]^.Im := 0.0;
  455.       exit;
  456.       end;
  457.     Power[SpPower]^.Re := exp(R * ln(A^.Re));
  458.     Power[SpPower]^.Im := R * A^.Im;
  459.     CpPwrRF := Power[SpPower];
  460.     IncPtr(Power, SpPower);
  461.     end;
  462.  
  463. procedure CpPwrR(A : Complex; R : ComplexElement; Result : Complex);
  464.   begin
  465.     Result^ := CpPwrRF(A, R)^;
  466.     end;
  467.   end.
  468.