diff options
Diffstat (limited to 'plugins/Libs/CplxMath.pas')
| -rw-r--r-- | plugins/Libs/CplxMath.pas | 278 |
1 files changed, 278 insertions, 0 deletions
diff --git a/plugins/Libs/CplxMath.pas b/plugins/Libs/CplxMath.pas new file mode 100644 index 0000000000..7cd180af9e --- /dev/null +++ b/plugins/Libs/CplxMath.pas @@ -0,0 +1,278 @@ +unit CplxMath;
+{* This unit contains functins for working with complex numbers. To use with
+ KOL library and its kolmath.pas unit instead of standard math.pas, define
+ synmbol KOL in project options, or uncomment its definition below. }
+
+interface
+
+//{$DEFINE KOL}
+
+{$IFNDEF KOL}
+ {$IFDEF KOL_MCK}
+ {$DEFINE KOL}
+ {$ENDIF}
+{$ENDIF}
+
+uses {$IFDEF KOL} kolmath, kol {$ELSE} math, sysutils {$ENDIF};
+
+type
+ {$IFDEF CPLX_EXTENDED}
+ Double = Extended;
+ {$ENDIF}
+
+ Complex = record Re, Im: double end;
+ {* }
+
+ function CfromReIm( Re, Im: Double ): Complex;
+ {* Re + i * Im }
+
+ function Cadd( const X, Y: Complex ): Complex;
+ {* X + Y }
+
+ function Cneg( const X: Complex ): Complex;
+ {* -X }
+
+ function Csub( const X, Y: Complex ): Complex;
+ {* X - Y }
+
+ function Cmul( const X, Y: Complex ): Complex;
+ {* X * Y }
+
+ function CmulD( const X: Complex; D: Double ): Complex;
+ {* X * D }
+
+ function CmulI( const X: Complex ): Complex;
+ {* i * X }
+
+ function Cdiv( const X, Y: Complex ): Complex;
+ {* X / Y }
+
+ function Cmod( const X: Complex ): Double;
+ {* Q( X.Re^2 + X.Im^2 ) }
+
+ function Carg( const X: Complex ): Double;
+ {* arctg( X.Im / X.Re ) }
+
+ function CfromModArg( R, Arg: Double ): Complex;
+ {* R * ( cos Arg + i * sin Arg ) }
+
+ function Cpow( const X: Complex; Pow: Double ): Complex;
+ {* X ^ Pow }
+
+ function Cpower( const X, Pow: Complex ): Complex;
+ {* X ^ Pow }
+
+ function CIntPower( const X: Complex; Pow: Integer ): Complex;
+ {* X ^ Pow}
+
+ function Csqrt( const X: Complex ): Complex;
+ {* Q( X ) }
+
+ function Cexp( const X: Complex ): Complex;
+ {* exp( X ) }
+
+ function Cln( const X: Complex ): Complex;
+ {* ln( X ) }
+
+ function Ccos( const X: Complex ): Complex;
+ {* cos( X ) }
+
+ function Csin( const X: Complex ): Complex;
+ {* sin( X ) }
+
+ function C2Str( const X: Complex ): String;
+ {* }
+
+ function C2StrEx( const X: Complex ): String;
+ {* experimental }
+
+implementation
+
+ function CfromReIm( Re, Im: Double ): Complex;
+ begin
+ Result.Re := Re;
+ Result.Im := Im;
+ end;
+
+ function Cadd( const X, Y: Complex ): Complex;
+ begin
+ Result.Re := X.Re + Y.Re;
+ Result.Im := X.Im + Y.Im;
+ end;
+
+ function Cneg( const X: Complex ): Complex;
+ begin
+ Result.Re := -X.Re;
+ Result.Im := -X.Im;
+ end;
+
+ function Csub( const X, Y: Complex ): Complex;
+ begin
+ Result := Cadd( X, Cneg( Y ) );
+ end;
+
+ function Cmul( const X, Y: Complex ): Complex;
+ begin
+ Result.Re := X.Re * Y.Re - X.Im * Y.Im;
+ Result.Im := X.Re * Y.Im + X.Im * Y.Re;
+ end;
+
+ function CmulD( const X: Complex; D: Double ): Complex;
+ begin
+ Result.Re := X.Re * D;
+ Result.Im := X.Im * D;
+ end;
+
+ function CmulI( const X: Complex ): Complex;
+ begin
+ Result.Re := -X.Im;
+ Result.Im := X.Re;
+ end;
+
+ function Cdiv( const X, Y: Complex ): Complex;
+ var Z: Double;
+ begin
+ Z := 1.0 / ( Y.Re * Y.Re + Y.Im * Y.Im );
+ Result.Re := (X.Re * Y.Re + X.Im * Y.Im ) * Z;
+ Result.Im := (X.Im * Y.Re - X.Re * Y.Im ) * Z;
+ end;
+
+ function Cmod( const X: Complex ): Double;
+ begin
+ Result := sqrt( X.Re * X.Re + X.Im * X.Im );
+ end;
+
+ function Carg( const X: Complex ): Double;
+ begin
+ Result := ArcTan2( X.Im, X.Re );
+ end;
+
+ function CfromModArg( R, Arg: Double ): Complex;
+ begin
+ Result.Re := R * cos( Arg );
+ Result.Im := R * sin( Arg );
+ end;
+
+ function Cpow( const X: Complex; Pow: Double ): Complex;
+ var R, A: Double;
+ begin
+ R := power( Cmod( X ), Pow );
+ A := Pow * Carg( X );
+ Result := CfromModArg( R, A );
+ end;
+
+ function Cpower( const X, Pow: Complex ): Complex;
+ begin
+ Result := Cexp( Cmul( X, Cln( Pow ) ) );
+ end;
+
+ function CIntPower( const X: Complex; Pow: Integer ): Complex;
+ begin
+ if (Pow < 0) or (Pow > 100) then Result := Cpow( X, Pow )
+ else if Pow = 0 then
+ begin
+ Result.Re := 1;
+ Result.Im := 0;
+ end
+ else
+ begin
+ Result := X;
+ while Pow > 1 do
+ begin
+ Result := Cmul( Result, X );
+ dec( Pow );
+ end;
+ end;
+ end;
+
+ function Csqrt( const X: Complex ): Complex;
+ begin
+ Result := Cpow( X, 0.5 );
+ end;
+
+ function Cexp( const X: Complex ): Complex;
+ var Z: Double;
+ begin
+ Z := exp( X.Re );
+ Result.Re := Z * cos( X.Im );
+ Result.Im := Z * sin( X.Im );
+ end;
+
+ function Cln( const X: Complex ): Complex;
+ begin
+ Result := CfromModArg( ln( Cmod( X ) ), Carg( X ) );
+ end;
+
+ function Ccos( const X: Complex ): Complex;
+ begin
+ Result := CmulI( X );
+ Result := CmulD( Cadd( Cexp( Result ), Cexp( Cneg( Result ) ) ),
+ 0.5 );
+ end;
+
+ function Csin( const X: Complex ): Complex;
+ begin
+ Result := CmulI( X );
+ Result := CmulD( Csub( Cexp(Result), Cexp( Cneg(Result) ) ),
+ 0.5 );
+ end;
+
+ {$IFDEF KOL}
+ function Abs( X: Double ): Double;
+ begin
+ Result := EAbs( X );
+ end;
+ {$ENDIF}
+
+ {$IFNDEF KOL}
+ function Double2Str( D: Double ): String;
+ begin
+ Result := DoubleToStr( D );
+ end;
+ {$ENDIF}
+
+ function C2Str( const X: Complex ): String;
+ begin
+ if Abs( X.Im ) < 1e-307 then
+ begin
+ Result := Double2Str( X.Re );
+ end
+ else
+ begin
+ Result := '';
+ if Abs( X.Re ) > 1e-307 then
+ begin
+ Result := Double2Str( X.Re );
+ if X.Im > 0.0 then
+ Result := Result + ' + ';
+ end;
+ if X.Im < 0.0 then
+ Result := Result + '- i * ' + Double2Str( -X.Im )
+ else
+ Result := Result + 'i * ' + Double2Str( X.Im );
+ end;
+ end;
+
+ function C2StrEx( const X: Complex ): String;
+ begin
+ if Abs( X.Im ) < 1e-307 then
+ begin
+ Result := Double2StrEx( X.Re );
+ end
+ else
+ begin
+ Result := '';
+ if Abs( X.Re ) > 1e-307 then
+ begin
+ Result := Double2StrEx( X.Re );
+ if X.Im > 0.0 then
+ Result := Result + ' + ';
+ end;
+ if X.Im < 0.0 then
+ Result := Result + '- i * ' + Double2StrEx( -X.Im )
+ else
+ Result := Result + 'i * ' + Double2StrEx( X.Im );
+ end;
+ end;
+
+end.
|