{-----------------------------------------------------------------------------
 Unit Name: ComplexNumbers
 Author:    kaetzel
 Purpose:   Contains simple algebraic operations with complex numbers
            of either cartesian or polar type and the conversion between the
            two types.
 Comments:  This unit has been written since the JCL Library and the Complex
            library of Delphi are erroneous and give wrong results when many
            computations are made. This led e.g. to values for a Complex Zero
            with (10^30+i*10^12) which is obviously NOT ZERO!
-----------------------------------------------------------------------------}

unit ComplexNumbers;

interface

uses JclMath,SysUtils,Math;

        {  Two types of Presentation of Complex Numbers can be defined   }

type
        // The Cartesian Type which is z= a + i*b
  TComplexCart=record
    Re:extended;
    Im:extended;
  end;
        // And the Polar or Euler Type which is z= r*exp(i*arg)
  TComplexPolar=record
    r:extended;         //r must per definition be greater or equal to zero.
    arg:extended;       //It is checked in this unit whether the user specifies
  end;

  {r must per definition be greater or equal to zero. It is checked in this unit
   whether the user specifies a correct r, otherwise an error will occur.
   Arguments greater than Pi or smaller -Pi are accepted but will be changed by
   every calculation to the above mentioned range}

        //Error Class for this Unit
  EComplexError=class(Exception);

const
        //These constants are global, they can be used in all units which
        //implement this unit

  ComplexZeroCart:TComplexCart    = (Re:0;Im:0);
  ComplexZeroPolar:TComplexPolar  = (r:0;arg:0);
  ImaginaryUnitCart:TComplexCart  = (Re:0;Im:1);
  ImaginaryUnitPolar:TComplexPolar= (r:1;arg:1.5707963267948966192313216916398);
  RealOne:TComplexCart            = (Re:1;Im:0);
  RealNegOne:TComplexCart         = (Re:-1;Im:0);
  NegImaginaryUnit:TComplexCart   = (Re:0;Im:-1);

{These are all the standard procedures and functions which can be used in all
 units which use this unit. Most of the functions are overloaded so that for
 both types of complex numbers only one function name has to be used.}

     {Type Transformation}

function ComplexCartToPolar (z:TComplexCart):TComplexPolar;
function ComplexPolarToCart (z:TComplexPolar):TComplexCart;
function Conjugate (z:TComplexCart):TComplexCart; overload;
function Conjugate (z:TComplexPolar):TComplexPolar; overload;
function AbsoluteValue (z:TComplexCart):extended;
function RealPart (z:TComplexCart):extended;
function ImaginaryPart (z:TComplexCart) :extended;
function InverseCart (z:TComplexCart) :TComplexCart;

     {Basic Math Operations}

function ComplexAdd (z1,z2:TComplexCart):TComplexCart; overload;
function ComplexAdd (z1,z2:TComplexPolar):TComplexPolar; overload;
function ComplexRealAdd (z1:TComplexCart;addend:extended)
                         :TComplexCart; overload;
function ComplexRealAdd (z1:TComplexPolar;addend:extended)
                         :TComplexPolar; overload;
function ComplexSub (z1,z2:TComplexCart):TComplexCart; overload;
function ComplexSub (z1,z2:TComplexPolar):TComplexPolar; overload;
function ComplexRealSub (z1:TComplexCart;minuend:extended)
                         :TComplexCart; overload;
function ComplexRealSub (z1:TComplexPolar;minuend:extended)
                         :TComplexPolar; overload;
function ComplexMul (z1,z2:TComplexCart):TComplexCart; overload;
function ComplexMul (z1,z2:TComplexPolar):TComplexPolar; overload;
function ComplexRealMul (z1:TComplexCart;factor:extended)
                          :TComplexCart; overload;
function ComplexRealMul (z1:TComplexPolar;factor:extended)
                          :TComplexPolar; overload;
function ComplexDiv (z1,z2:TComplexCart):TComplexCart; overload;
function ComplexDiv (z1,z2:TComplexPolar):TComplexPolar; overload;
function ComplexRealDiv (z1:TComplexCart;divisor:extended)
                         :TComplexCart; overload;
function ComplexRealDiv(z1:TComplexPolar;divisor:extended)
                         :TComplexPolar; overload;

      {Power Function}

function ComplexRealPower (z1:TComplexCart;exponent:extended)
                           :TComplexCart; overload;
function ComplexRealPower (z1:TComplexPolar;exponent:extended)
                           :TComplexPolar; overload;
function ImaginaryIntegerPower (exponent:integer):TComplexCart;
function ComplexCartSqr(z1:TComplexCart):TComplexCart;

      {Display Functions}

function DisplayComplex (z1:TComplexCart):String; overload;
function DisplayComplex (z1:TComplexPolar):String; overload;

      {Mathematical Help functions, no direct use with complex numbers}

function ModReal (x,y:extended):extended;


implementation


{-----------------------------------------------------------------------------
  Procedure: ComplexCartToPolar
  Purpose:   Converts the Cartesian type into the polar type
  Comments:
-----------------------------------------------------------------------------}

function ComplexCartToPolar (z:TComplexCart):TComplexPolar;
begin
        //Absolute Value
  Result.r:=Power((z.Re*z.Re+z.Im*z.Im),0.5);
  if ((z.Re=0)and(z.Im=0)) then
    begin
      Result.arg:=0;
      exit;
    end;
  if ((z.Re=0)and(z.Im>0)) then
    begin
      Result.arg:=PiOn2;
      exit;
    end;
  if ((z.Re=0)and(z.Im<0)) then
    begin
      Result.arg:=-PiOn2;
      exit;
    end;
  if z.Re>0 then
    Result.arg:=arctan(z.Im/z.Re)
  else
    begin
      if z.Im>=0 then
        Result.arg:=arctan(z.Im/z.Re)+Pi
      else
        Result.arg:=arctan(z.Im/z.Re)-Pi;
    end;
end;


{-----------------------------------------------------------------------------
  Procedure: ComplexPolarToCart
  Purpose:   Converts the polar type into the Cartesian one.
  Comments:
-----------------------------------------------------------------------------}

function ComplexPolarToCart (z:TComplexPolar):TComplexCart;
begin
  if (((z.r=0) and (z.arg<>0))or(z.r<0)) then
    Raise EComplexError.Create('Wrong definition in ComplexNumbers.pas: Polar '+
          'Complex Numbers must have a radius>=0 and a argument=0 when r=0!');
  if z.arg=0 then
    begin
      Result.Re:=z.r;
      Result.Im:=0;
      exit;
    end;
        //These conditions are only written to circumvent numerical errors
        //with angular functions
  if (abs(cos(z.arg)))<0.000000000000001 then
    begin
      Result.Re:=0;
      Result.Im:=z.r;
      exit;
    end;
  if (((abs(z.arg/Pi))<1.000000000000001) and ((abs(z.arg/Pi))>0.999999999999999999)) then
    begin
      Result.Re:=-z.r;
      Result.Im:=0;
      exit;
    end;
  Result.Re:=z.r*cos(z.arg);
  Result.Im:=z.r*sin(z.arg);
end;


{-----------------------------------------------------------------------------
  Procedure: Conjugate
  Purpose:   Calculate the conjugate complex number
  Comments:  a+i*b ---> a-i*b
-----------------------------------------------------------------------------}

function Conjugate (z:TComplexCart):TComplexCart; overload;
begin
  Result.Re:=z.Re;
  Result.Im:=-z.Im;
end;


{-----------------------------------------------------------------------------
  Procedure: SimplifyArgument
  Purpose:   Simplifies the argument of a polar complex number so that it will
             be always in the range of [-Pi,+Pi]
  Comments:
-----------------------------------------------------------------------------}

procedure SimplifyArgument (var z:TComplexPolar);
begin
  if (abs((z.arg/Pi)))>1.000000001 then
    begin
      if z.arg>0 then
        begin
          if odd(trunc(z.arg/Pi)) then
            z.arg:=ModReal(z.arg,Pi)-Pi
          else
            z.arg:=ModReal(z.arg,Pi);
        end
      else
        begin
          if odd(trunc(z.arg/(-Pi))) then
            z.arg:=Pi+ModReal(z.arg,Pi)
          else
            z.arg:=ModReal(z.arg,-Pi);
        end;
      if abs(z.arg)<1e-12 then
        z.arg:=0;
    end
  else
    z:=z;
end;


{-----------------------------------------------------------------------------
  Procedure: Conjugate
  Purpose:   Calculate the conjugate complex number
  Comments:  r*exp(i*arg)-->r*exp(i*-arg)
-----------------------------------------------------------------------------}

function Conjugate (z:TComplexPolar):TComplexPolar; overload;
begin
  if (((z.r=0) and (z.arg<>0))or(z.r<0)) then
    Raise EComplexError.Create('Wrong definition in ComplexNumbers.pas: Polar '+
          'Complex Numbers must have a radius>=0 and a argument=0 when r=0!');
  SimplifyArgument(z);
  Result.r:=z.r;
  Result.arg:=-z.arg;
end;


{-----------------------------------------------------------------------------
  Procedure: AbsoluteValue
  Purpose:   Calculate the absolute value of Cartesian complex number
  Comments:
-----------------------------------------------------------------------------}

function AbsoluteValue (z:TComplexCart):extended;
begin
  Result:=Power((z.Re*z.Re+z.Im*z.Im),0.5);
end;


{-----------------------------------------------------------------------------
  Procedure: RealPart
  Purpose:   Calculate the Real Part of a Cartesian complex Number
  Comments:
-----------------------------------------------------------------------------}

function RealPart (z:TComplexCart):extended;
begin
  Result:=z.Re;
end;


{-----------------------------------------------------------------------------
  Procedure: ImaginaryPart
  Purpose:   Calculate the Imaginary Part of a Cartesian Complex Number
  Comments:
-----------------------------------------------------------------------------}

function ImaginaryPart (z:TComplexCart):extended;
begin
  Result:=z.Im;
end;


{-----------------------------------------------------------------------------
  Procedure: InverseCart
  Purpose:   Calculates the inverse complex number for Cartesian type
  Comments:  z2=1/z1
-----------------------------------------------------------------------------}

function InverseCart (z:TComplexCart):TComplexCart;
begin
  Result.Re:=z.Re/(z.Re*z.Re+z.Im*z.Im);
  Result.Im:=-z.Im/(z.Re*z.Re+z.Im*z.Im);
end;


{-----------------------------------------------------------------------------
  Procedure: ComplexAdd
  Purpose:   Calculate a addition of two complex numbers of equal type
  Comments:  Not defined for the polar type, first converted into
             the Cartesian type and reconverted after addition
-----------------------------------------------------------------------------}

function ComplexAdd (z1,z2:TComplexCart):TComplexCart; overload;
begin
  Result.Re:=z1.Re+z2.Re;
  Result.Im:=z1.Im+z2.Im;
end;

function ComplexAdd (z1,z2:TComplexPolar):TComplexPolar; overload;
var
  aid1,aid2,aid3:TComplexCart;
begin
  if (((z1.r=0) and (z1.arg<>0)) or (z1.r<0)) then
    Raise EComplexError.Create('Wrong definition in ComplexNumbers.pas: Polar '+
          'Complex Numbers must have a radius>=0 and a argument=0 when r=0!');
  if (((z2.r=0) and (z2.arg<>0)) or (z2.r<0)) then
    Raise EComplexError.Create('Wrong definition in ComplexNumbers.pas: Polar '+
          'Complex Numbers must have a radius>=0 and a argument=0 when r=0!');
  SimplifyArgument(z1);
  SimplifyArgument(z2);
  aid1:=ComplexPolarToCart(z1);
  aid2:=ComplexPolarToCart(z2);
  aid3:=ComplexAdd(aid1,aid2);
  Result:=ComplexCartToPolar(aid3);
end;


{-----------------------------------------------------------------------------
  Procedure: ComplexRealAdd
  Purpose:   Adds a real value to a complex number
  Comments:
-----------------------------------------------------------------------------}

function ComplexRealAdd (z1:TComplexCart;addend:extended)
                         :TComplexCart; overload;
begin
  Result.Re:=z1.Re+addend;
  Result.Im:=z1.Im;
end;

function ComplexRealAdd (z1:TComplexPolar;addend:extended)
                         :TComplexPolar; overload;
var
  aid:TComplexPolar;
begin
  if (((z1.r=0) and (z1.arg<>0)) or (z1.r<0)) then
    Raise EComplexError.Create('Wrong definition in ComplexNumbers.pas: Polar '+
          'Complex Numbers must have a radius>=0 and a argument=0 when r=0!');
  SimplifyArgument(z1);
  aid.r:=addend;
  aid.arg:=0;
  Result:=ComplexAdd(z1,aid);
end;


{-----------------------------------------------------------------------------
  Procedure: ComplexSub
  Purpose:   Calculate a subtraction of two complex numbers of equal type
  Comments:  Not defined for the polar type, first converted into
             the Cartesian type and reconverted after addition
-----------------------------------------------------------------------------}

function ComplexSub (z1,z2:TComplexCart):TComplexCart; overload;
begin
  Result.Re:=z1.Re-z2.Re;
  Result.Im:=z1.Im-z2.Im;
end;

function ComplexSub (z1,z2:TComplexPolar):TComplexPolar; overload;
var
  aid1,aid2,aid3:TComplexCart;
begin
  if (((z1.r=0) and (z1.arg<>0)) or (z1.r<0)) then
    Raise EComplexError.Create('Wrong definition in ComplexNumbers.pas: Polar '+
          'Complex Numbers must have a radius>=0 and a argument=0 when r=0!');
  if (((z2.r=0) and (z2.arg<>0)) or (z2.r<0)) then
    Raise EComplexError.Create('Wrong definition in ComplexNumbers.pas: Polar '+
          'Complex Numbers must have a radius>=0 and a argument=0 when r=0!');
  SimplifyArgument(z1);
  SimplifyArgument(z2);
  aid1:=ComplexPolarToCart(z1);
  aid2:=ComplexPolarToCart(z2);
  aid3:=ComplexSub(aid1,aid2);
  Result:=ComplexCartToPolar(aid3);
end;


{-----------------------------------------------------------------------------
  Procedure: ComplexRealSub
  Purpose:   Subtracts a real value from a complex number
  Comments:
-----------------------------------------------------------------------------}

function ComplexRealSub (z1:TComplexCart;minuend:extended)
                         :TComplexCart; overload;
begin
  Result.Re:=z1.Re-minuend;
  Result.Im:=z1.Im;
end;

function ComplexRealSub (z1:TComplexPolar;minuend:extended)
                         :TComplexPolar; overload;
var
  aid:TComplexPolar;
begin
  if (((z1.r=0) and (z1.arg<>0)) or (z1.r<0)) then
    Raise EComplexError.Create('Wrong definition in ComplexNumbers.pas: Polar '+
          'Complex Numbers must have a radius>=0 and a argument=0 when r=0!');
  SimplifyArgument(z1);
  aid.r:=minuend;
  aid.arg:=0;
  Result:=ComplexSub(z1,aid);
end;


{-----------------------------------------------------------------------------
  Procedure: ComplexMul
  Purpose:   Declare the multiplication of two complex numbers
  Comments:
-----------------------------------------------------------------------------}

function ComplexMul (z1,z2:TComplexCart):TComplexCart; overload;
begin
  Result.Re:=z1.Re*z2.Re-z1.Im*z2.Im;
  Result.Im:=z1.Re*z2.Im+z1.Im*z2.Re;
end;

function ComplexMul (z1,z2:TComplexPolar):TComplexPolar; overload;
begin
  if (((z1.r=0) and (z1.arg<>0)) or (z1.r<0)) then
    Raise EComplexError.Create('Wrong definition in ComplexNumbers.pas: Polar '+
          'Complex Numbers must have a radius>=0 and a argument=0 when r=0!');
  if (((z2.r=0) and (z2.arg<>0)) or (z2.r<0)) then
    Raise EComplexError.Create('Wrong definition in ComplexNumbers.pas: Polar '+
          'Complex Numbers must have a radius>=0 and a argument=0 when r=0!');
  Result.r:=z1.r*z2.r;
  Result.arg:=z1.arg+z2.arg;
  SimplifyArgument(Result);
end;


{-----------------------------------------------------------------------------
  Procedure: ComplexRealMul
  Purpose:   Multiplies a real value with a complex number
  Comments:
-----------------------------------------------------------------------------}

function ComplexRealMul (z1:TComplexCart;factor:extended)
                         :TComplexCart; overload;
begin
  Result.Re:=z1.Re*factor;
  Result.Im:=z1.Im*factor;
end;

function ComplexRealMul (z1:TComplexPolar;factor:extended)
                         :TComplexPolar; overload;
begin
  if (((z1.r=0) and (z1.arg<>0)) or (z1.r<0)) then
    Raise EComplexError.Create('Wrong definition in ComplexNumbers.pas: Polar '+
          'Complex Numbers must have a radius>=0 and a argument=0 when r=0!');
  if factor=0 then
    begin
      Result.r:=0;
      Result.arg:=0;
      exit;
    end;
  if factor<0 then
    begin
      Result.r:=z1.r*abs(factor);
      SimplifyArgument(z1);
      if z1.arg>0 then
        Result.arg:=z1.arg-Pi
      else
        Result.arg:=z1.arg+Pi;
      exit;
    end;
  Result.r:=z1.r*factor;
  Result.arg:=z1.arg;
end;


{-----------------------------------------------------------------------------
  Procedure: ComplexDiv
  Purpose:   Calculate the division of two complex numbers of the same type
  Comments:
-----------------------------------------------------------------------------}

function ComplexDiv (z1,z2:TComplexCart):TComplexCart; overload;
begin
  Result.Re:=(z1.Re*z2.Re+z1.Im*z2.Im)/(z2.Re*z2.Re+z2.Im*z2.Im);
  Result.Im:=(z2.Re*z1.Im-z1.Re*z2.Im)/(z2.Re*z2.Re+z2.Im*z2.Im);
end;

function ComplexDiv (z1,z2:TComplexPolar):TComplexPolar; overload;
begin
  if (((z1.r=0) and (z1.arg<>0)) or (z1.r<0)) then
    Raise EComplexError.Create('Wrong definition in ComplexNumbers.pas: Polar '+
          'Complex Numbers must have a radius>=0 and a argument=0 when r=0!');;
  if (((z2.r=0) and (z2.arg<>0)) or (z2.r<0)) then
    Raise EComplexError.Create('Wrong definition in ComplexNumbers.pas: Polar '+
          'Complex Numbers must have a radius>=0 and a argument=0 when r=0!');
  SimplifyArgument(z1);
  SimplifyArgument(z2);
  Result.r:=z1.r/z2.r;
  Result.arg:=z1.arg-z2.arg;
  SimplifyArgument(Result);
end;


{-----------------------------------------------------------------------------
  Procedure: ComplexRealDiv
  Purpose:   Calculate the division of complex number by a real number
  Comments:  Divisor=0 will cause an exception
-----------------------------------------------------------------------------}

function ComplexRealDiv (z1:TComplexCart;divisor:extended)
                         :TComplexCart; overload;
begin
  Result.Re:=z1.Re/divisor;
  Result.Im:=z1.Im/divisor;
end;

function ComplexRealDiv(z1:TComplexPolar;divisor:extended)
                        :TComplexPolar; overload;
begin
  if (((z1.r=0) and (z1.arg<>0)) or (z1.r<0)) then
    Raise EComplexError.Create('Wrong definition in ComplexNumbers.pas: Polar '+
          'Complex Numbers must have a radius>=0 and a argument=0 when r=0!');
  if divisor<0 then
    begin
      Result.r:=z1.r/(abs(divisor));
      SimplifyArgument(z1);
      if z1.arg>0 then
        Result.arg:=z1.arg-Pi
      else
        Result.arg:=z1.arg+Pi;
      exit;
      end;
  Result.r:=z1.r/divisor;
  Result.arg:=z1.arg;
end;


{-----------------------------------------------------------------------------
  Procedure: ComplexRealPower
  Purpose:   Calculate the power of a complex number
  Comments:  Exponent must be a real number
-----------------------------------------------------------------------------}

function ComplexRealPower (z1:TComplexCart;exponent:extended)
                           :TComplexCart; overload;
        {The function is not generally defined in the cartesian type,
         therefore the variable is first converted into the polar type
         and reconverted after calculation                             }
var
  aid1,aid2:TComplexPolar;
begin
  aid1:=ComplexCartToPolar(z1);
  aid2:=ComplexRealPower(aid1,exponent);
  Result:=ComplexPolarToCart(aid2);
end;

function ComplexRealPower (z1:TComplexPolar;exponent:extended)
                           :TComplexPolar; overload;
begin
  if z1.r<0 then
    Raise EComplexError.Create('Wrong radius in ComplexNumbers.pas, r must be'+
          ' >=0 !');
  if z1.r=0 then
    begin
      if exponent=0 then
        Result:=ComplexZeroPolar
      else
        if z1.arg<>0 then
          Raise EComplexError.Create('Wrong calculation in ComplexNumbers.pas,'+
                                     ' if r=0 then arg must also be zero!')
        else
          Result:=ComplexZeroPolar;
    end
  else
    begin
      Result.r:=Power(z1.r,exponent);
      Result.arg:=z1.arg*exponent;
    end;
  SimplifyArgument(Result);
end;


{-----------------------------------------------------------------------------
  Procedure: ImaginaryIntegerPower
  Purpose:   Calculates the power of the imaginary unit i
  Comments:  Exponent must be an integer
             Saves calculation time in multisphere scattering
-----------------------------------------------------------------------------}

function ImaginaryIntegerPower (exponent:integer):TComplexCart;
var
  test:integer;
begin
  Test:=exponent mod 4;
  case test of
    0: Result:=RealOne;
    1: Result:=ImaginaryUnitCart;
    2: Result:=RealNegOne;
    3: Result:=NegImaginaryUnit;
  end;
end;


{-----------------------------------------------------------------------------
  Procedure: ComplexCartSqr
  Purpose:   Calculates the square of a Complex Number a+ib
  Comments:  (aČ-bČ)+i(2ab)
-----------------------------------------------------------------------------}

function ComplexCartSqr(z1:TComplexCart):TComplexCart;
begin
  Result.Re:=z1.Re*z1.Re-z1.Im*z1.Im;
  Result.Im:=2*z1.Re*z1.Im;
end;


{-----------------------------------------------------------------------------
  Procedure: DisplayComplex
  Purpose:   Convert a complex number into a string
  Comments:  This can be displayed e.g. in a Edit component with property Text
-----------------------------------------------------------------------------}

function DisplayComplex (z1:TComplexCart):String; overload;
begin
  SetLength(Result,40);
  if z1.Im>=0 then
    Result:=(FloatToStrF(z1.Re,ffGeneral,4,4)+'+'+
             FloatToStrF(z1.Im,ffGeneral,4,4)+'i')
  else
    Result:=(FloatToStrF(z1.Re,ffGeneral,4,4)+'-'+
             FloatToStrF(abs(z1.Im),ffGeneral,4,4)+'i');
end;

function DisplayComplex (z1:TComplexPolar):String; overload;
begin
  if (((z1.r=0) and (z1.arg<>0)) or (z1.r<0)) then
    Raise EComplexError.Create('Wrong definition in ComplexNumbers.pas: Polar '+
          'Complex Numbers must have a radius>=0 and a argument=0 when r=0!');
  SimplifyArgument(z1);
  SetLength(Result,40);
  Result:=(FloatToStrF(z1.r,ffGeneral,4,4)+'*exp(i*'+
           FloatToStrF(z1.arg,ffGeneral,4,4)+')');
end;


{-----------------------------------------------------------------------------
  Procedure: ModReal
  Purpose:   Realizes the mod function (remainder of a division) for real values
  Comments:  Needed in the calculation of the argument for polar complex numbers
-----------------------------------------------------------------------------}

function ModReal (x,y:extended):extended;
begin
  Result:=x-(trunc(x/y))*y;
end;

end.
