Programmer 7500/MAX_PROGRAMMERS.iso/PASCAL/PNL006.ZIP/SHCMPLX1.ZIP/SHCMPLX.PAS

Pascal/Delphi Source File | 1991-03-08 | 10.0 KB | 357 lines

{$N+,E+}
Unit ShCmplx;
{
                                ShCmplx

       A Complex Arithmetic Unit for Turbo Pascal 5.0 and above.

                                   by

                   W. G. Madison and Associates, Ltd.
                          8425 Greenbelt Road
                              P.O. Box 898
                          Greenbelt, MD 20770
                             [phone]
                             CIS 73240,342

                  Copyright 1991 Madison & Associates
                          All Rights Reserved

        This unit and the associated .DOC and TEST*.* files may
        be freely copied and distributed, provided only that no
        fee is charged for the package beyond a nominal copying
        charge, and provided that the package is distributed IN
        UNALTERED FORM. The sole exception to  this  latter re-
        striction is that bona-fide clubs and  user  groups may
        append text material to the documentation file, provid-
        ed that any material appended is clearly  identified as
        to its source, its beginning, and its end.
}

interface

type
  ComplexElement  = extended;
  ComplexBaseType = record
                      Re,
                      Im  : ComplexElement;
                      end;
  Complex         = ^ComplexBaseType;

function CmplxError  : integer;

function Cmp2Str(A : Complex; Width, Places : byte) : string;
{Converts a complex value to a string of the form (Re + Im i)}

function CmpP2Str(A : Complex; Width, Places : byte) : string;
{Converts a complex in polar form to a string with the angle in radians.}

function CmpP2StrD(A : Complex; Width, Places : byte) : string;
{Converts a complex in polar form to a string with the angle in degrees.}

procedure CAbs(A  : Complex; var Result : ComplexElement);
function CAbsF(A  : Complex)  : ComplexElement;
{Returns the absolute value of a complex number}

procedure CConj(A : Complex; Result : Complex);
function CConjF(A : Complex) : Complex;
{Returns the complex conjugate of a complex number}

procedure CAdd(A, B : Complex; Result : Complex);
function CAddF(A, B : Complex) : Complex;
{Returns the sum of two complex numbers A + B}

procedure RxC(A : ComplexElement; B : Complex; Result : Complex);
function RxCF(A : ComplexElement; B : Complex) : Complex;
{Returns the complex product of a real and a complex}

procedure CSub(A, B : Complex; Result : Complex);
function CSubF(A, B : Complex) : Complex;
{Returns the difference between two complex numbers A - B}

procedure CMul(A, B : Complex; Result : Complex);
function CMulF(A, B : Complex) : Complex;
{Returns the product of two complex numbers A * B}

procedure CDiv(A, B : Complex; Result : Complex);
function CDivF(A, B : Complex) : Complex;
{Returns the quotient of two complex numbers A / B}

procedure C2P(A : Complex; Result : Complex);
function C2PF(A : Complex) : Complex;
{Transforms a complex in cartesian form into polar form}
{The magnitude will be stored in Result^.Re and the angle in Result^.Im}

procedure P2C(A : Complex; Result : Complex);
function P2CF(A : Complex) : Complex;
{Transforms a complex in polar form into cartesian form}
{The magnitude is stored in A^.Re and the angle in A^.Im}

procedure CpPwrR(A : Complex; R : ComplexElement; Result : Complex);
function CpPwrRF(A : Complex; R : ComplexElement) : Complex;
{Raises a complex (in polar form) to a real power}

{===========================================================}

implementation

const
  MaxStack  = 5;
  StackTop  = MaxStack - 1;
  Pi        = 3.14159265358979;
  PiOver2   = Pi / 2.0;
  TwoPi     = Pi * 2.0;
  F180OverPi= 180.0 / Pi;
  SpConj  : byte = 0;
  SpSum   : byte = 0;
  SpDiff  : byte = 0;
  SpRProd : byte = 0;
  SpProd  : byte = 0;
  SpQuot  : byte = 0;
  SpPolar : byte = 0;
  SpCartsn: byte = 0;
  SpPower : byte = 0;

var
  Conj,
  Sum,
  Diff,
  RProd,
  Prod,
  Quot,
  Polar,
  Cartsn,
  Power    : array[0..StackTop] of Complex;
  CmpError : integer;

function CmplxError  : integer;
  begin
    CmplxError := CmpError;
    CmpError := 0;
    end;

function Real2Str(A : ComplexElement; Width, Places : byte) : string;
  var
    T1  : string;
  begin
    Str(A:Width:Places, T1);
    Real2Str := T1;
    end;

function Cmp2Str(A : Complex; Width, Places : byte) : string;
  begin
    if A^.Im >= 0.0 then
      Cmp2Str := '(' + Real2Str(A^.Re, Width, Places) + ' + ' +
                       Real2Str(A^.Im, Width, Places) + 'i)'
    else
      Cmp2Str := '(' + Real2Str(A^.Re, Width, Places) + ' - ' +
                       Real2Str(Abs(A^.Im), Width, Places) + 'i)';
    end;

function CmpP2StrD(A : Complex; Width, Places : byte) : string;
{Converts a complex in polar form to a string with the angle in degrees.}
  begin
    CmpP2StrD := '('+Real2Str(A^.Re,Width,Places)+' at '+
      Real2Str((A^.Im*F180OverPi),Width,Places)+#248')';
    end;

function CmpP2Str(A : Complex; Width, Places : byte) : string;
{Converts a complex in polar form to a string with the angle in radians.}
  begin
    CmpP2Str := '('+Real2Str(A^.Re,Width,Places)+' at '+
      Real2Str((A^.Im),Width,Places)+' rad)';
    end;

procedure IncPtr(var StackPtr : byte);
  begin
    StackPtr := (StackPtr + 1) mod StackTop;
    end;

function CAbsF(A  : Complex) : ComplexElement;
  begin
    CmpError := 0;
    CAbsF := sqrt(A^.Re * A^.Re + A^.Im * A^.Im);
    end;

procedure CAbs(A  : Complex; var Result : ComplexElement);
  begin
    Result := CAbsF(A);
    end;

function CConjF(A : Complex) : Complex;
  begin
    CmpError := 0;
    Conj[SpConj]^.Re := A^.Re;
    Conj[SpConj]^.Im := -A^.Im;
    CConjF := Conj[SpConj];
    IncPtr(SpConj);
    end;

procedure CConj(A : Complex; Result : Complex);
  begin
    Result^ := CConjF(A)^;
    end;

function CAddF(A, B : Complex) : Complex;
  begin
    CmpError := 0;
    Sum[SpSum]^.Re := A^.Re + B^.Re;
    Sum[SpSum]^.Im := A^.Im + B^.Im;
    CAddF := Sum[SpSum];
    IncPtr(SpSum);
    end;

procedure CAdd(A, B : Complex; Result : Complex);
  begin
    Result^ := CAddF(A, B)^;
    end;

function RxCF(A : ComplexElement; B : Complex) : Complex;
  begin
    CmpError := 0;
    RProd[SpRProd]^.Re := A * B^.Re;
    RPRod[SpRProd]^.Im := A * B^.Im;
    RxCF := RProd[SpRProd];
    IncPtr(SpRProd);
    end;

procedure RxC(A : ComplexElement; B : Complex; Result : Complex);
  begin
    Result^ := RxCF(A, B)^;
    end;

function CSubF(A, B : Complex) : Complex;
  begin
    CmpError := 0;
    Diff[SpDiff]^ := CAddF(A, RxCF(-1.0,B))^;
    CSubF := Diff[SpDiff];
    IncPtr(SpDiff);
    end;

procedure CSub(A, B : Complex; Result : Complex);
  begin
    Result^ := CSubF(A, B)^;
    end;

function CMulF(A, B : Complex) : Complex;
  begin
    CmpError := 0;
    Prod[SpProd]^.Re := A^.Re * B^.Re - A^.Im * B^.Im;
    Prod[SpProd]^.Im := A^.Im * B^.Re + A^.Re * B^.Im;
    CMulF := Prod[SpProd];
    IncPtr(SpProd);
    end;

procedure CMul(A, B : Complex; Result : Complex);
  begin
    Result^ := CMulF(A, B)^;
    end;

function CDivF(A, B : Complex) : Complex;
  var
    Denom : ComplexElement;
  begin
    CmpError := 0;
    Denom := B^.Re * B^.Re + B^.Im * B^.Im;
    if Denom = 0.0 then begin
      CmpError := -1;
      Quot[SpQuot]^.Re := 0.0;
      Quot[SpQuot]^.Im := 0.0;
      CDivF := Quot[SpQuot];
      IncPtr(SpQuot);
      exit;
      end;
    Quot[SpQuot]^ := RxCF(1.0 / Denom, CMulF(A, CConjF(B)))^;
    CDivF := Quot[SpQuot];
    IncPtr(SpQuot);
    end;

procedure CDiv(A, B : Complex; Result : Complex);
  begin
    Result^ := CDivF(A, B)^;
    end;

procedure C2P(A : Complex; Result : Complex);
{Transforms a complex in cartesian form into polar form}
{The magnitude will be stored in Result^.Re and the angle in Result^.Im}
  begin
    CmpError := 0;
    Result^.Re := CAbsF(A);
    if abs(A^.Re) > abs(A^.Im) then begin
      {Use the tangent formulation}
      Result^.Im := ArcTan(abs(A^.Im) / abs(A^.Re));
      end
    else begin
      {Use the cotangent formulation}
      Result^.Im := PiOver2 - ArcTan(abs(A^.Re) / abs(A^.Im));
      end;
    if A^.Im > 0 then begin {first or second quadrant}
      if A^.Re < 0.0 then {Second quadrant}
        Result^.Im := Pi - abs(Result^.Im);
      end
    else begin {third or fourth quadrant}
      if A^.Re > 0.0 then {Fourth quadrant}
        Result^.Im := TwoPi - abs(Result^.Im)
      else begin {Third Quadrant}
        Result^.Im := Pi + abs(Result^.Im);
        end;
      end;
    if Result^.Im = TwoPi then Result^.Im := 0.0;
    end;

function C2PF(A : Complex) : Complex;
  begin
    C2P(A, Polar[SpPolar]);
    C2PF := Polar[SpPolar];
    IncPtr(SpPolar);
    end;

procedure P2C(A : Complex; Result : Complex);
{Transforms a complex in polar form into cartesian form}
{The magnitude is stored in A^.Re and the angle in A^.Im}
  begin
    CmpError := 0;
    Result^.Re := A^.Re * cos(A^.Im);
    Result^.Im := A^.Re * sin(A^.Im);
    end;

function P2CF(A : Complex) : Complex;
  begin
    P2C(A, Cartsn[SpCartsn]);
    P2CF := Cartsn[SpCartsn];
    IncPtr(SpCartsn);
    end;

function CpPwrRF(A : Complex; R : ComplexElement) : Complex;
{Raises a complex (in polar form) to a real power, returning the result
 also in polar form}
  begin
    CmpError := 0;
    if A^.Re <= 0.0 then begin
      CmpError := -2;
      Power[SpPower]^.Re := 0.0;
      Power[SpPower]^.Im := 0.0;
      exit;
      end;
    Power[SpPower]^.Re := exp(R * ln(A^.Re));
    Power[SpPower]^.Im := R * A^.Im;
    CpPwrRF := Power[SpPower];
    IncPtr(SpPower);
    end;

procedure CpPwrR(A : Complex; R : ComplexElement; Result : Complex);
  begin
    Result^ := CpPwrRF(A, R)^;
    end;

var
  T1  : byte;

begin
  for T1 := 0 to StackTop do begin
    New(Conj[T1]); New(Sum[T1]); New(Diff[T1]);
    New(RProd[T1]); New(Prod[T1]); New(Quot[T1]);
    New(Polar[T1]); New(Cartsn[T1]); New(Power[T1]);
    end;
  end.