IMPLEMENTATION MODULE LngCmplxMath;
(*------------------------------------------------------------------------*)
(* Baisisrechenoperationen f端r den Record-Typ Deklera.LONGCMPLX *)
(* und abgeleitete Funktionen *)
(* Mathematical functions for the type LONGCMPLX *)
(*------------------------------------------------------------------------*)
(* Letzte Bearbeitung: *)
(* *)
(* 12.07.93, MRi: Erstellen der ersten Version *)
(* 18.07.15, MRi: Umbenennung in LngCmplxMath *)
(* 04.08.15, MRi: Anpassung abs (Eleminierung der Vertauschung) *)
(* 24.09.17, MRi: Umbenennung von LONGCOMPLX in LONGCMPLX *)
(* Zusammenlegen von CMathLib mit LngCmplxMath *)
(*------------------------------------------------------------------------*)
(* Implementation : Michael Riedl *)
(* Licence : GNU Lesser General Public License (LGPL) *)
(*------------------------------------------------------------------------*)
(* $Id: LngCmplxMath.mod,v 1.2 2017/09/24 10:14:24 mriedl Exp mriedl $ *)
FROM Errors IMPORT FatalError;
FROM Deklera IMPORT LONGCMPLX;
IMPORT LongMath;
FROM LMathLib IMPORT sinh,cosh;
PROCEDURE CMPLX(real,imag : LONGREAL) : LONGCMPLX;
VAR z : LONGCMPLX;
BEGIN
WITH z DO re:=real; im:=imag; END; RETURN z;
END CMPLX;
PROCEDURE CAdd(x,y : LONGCMPLX) : LONGCMPLX;
VAR z : LONGCMPLX;
BEGIN
WITH z DO re:=x.re + y.re; im:=x.im + y.im; END; RETURN z;
END CAdd;
PROCEDURE CSub(x,y : LONGCMPLX) : LONGCMPLX;
VAR z : LONGCMPLX;
BEGIN
WITH z DO re:=x.re - y.re; im:=x.im - y.im; END; RETURN z;
END CSub;
PROCEDURE CMul(x,y : LONGCMPLX) : LONGCMPLX;
VAR z : LONGCMPLX;
BEGIN
WITH z DO re:=x.re*y.re - x.im*y.im; im:=x.re*y.im + x.im*y.re; END;
RETURN z;
END CMul;
PROCEDURE CDiv(x,y : LONGCMPLX) : LONGCMPLX;
VAR Zw : LONGREAL;
z : LONGCMPLX;
BEGIN
IF (y.re = 0.0) AND (y.im = 0.0) THEN
FatalError('Division durch Null in LngCmplxMath.CDiv.');
(* Eventuell besser auf NAN stellen *)
END;
IF (ABS(y.re) > ABS(y.im)) THEN
Zw:=y.im/y.re; y.re:=Zw*y.im + y.re;
z.re:=(x.re + Zw*x.im) / y.re;
z.im:=(x.im - Zw*x.re) / y.re;
ELSE
Zw:=y.re/y.im; y.im:=Zw*y.re + y.im;
z.re:=(Zw*x.re + x.im) / y.im;
z.im:=(Zw*x.im - x.re) / y.im;
END;
RETURN z;
END CDiv;
PROCEDURE CRez(x : LONGCMPLX) : LONGCMPLX;
VAR z : LONGCMPLX;
Zw : LONGREAL;
BEGIN
IF (x.re = 0.0) AND (x.im = 0.0) THEN
FatalError('Division durch Null in LngCmplxMath.CRez.');
END;
IF (ABS(x.re) > ABS(x.im)) THEN
Zw:=x.im / x.re; x.re:=Zw*x.im + x.re;
z.re:= 1.0 / x.re;
z.im:= - Zw / x.re;
ELSE
Zw:=x.re / x.im; x.im:=Zw*x.re + x.im;
z.re:= Zw / x.im;
z.im:= - 1.0 / x.im;
END;
RETURN z;
END CRez;
PROCEDURE abs(x : LONGCMPLX) : LONGREAL;
VAR Zw : LONGREAL;
BEGIN
x.re:=ABS(x.re); x.im:=ABS(x.im);
IF (x.im > x.re) THEN
Zw:=(x.re / x.im);
RETURN x.im*LongMath.sqrt(1.0 + Zw*Zw);
ELSIF (x.im = 0.0) THEN
RETURN x.re;
ELSE
Zw:=(x.im / x.re);
RETURN x.re*LongMath.sqrt(1.0 + Zw*Zw);
END;
END abs;
PROCEDURE conj(x : LONGCMPLX) : LONGCMPLX;
BEGIN
x.im:= - x.im; RETURN x;
END conj;
PROCEDURE sqrt(x : LONGCMPLX) : LONGCMPLX;
VAR Zw : LONGREAL;
BEGIN
Zw:=LongMath.sqrt(0.5*(ABS(x.re) + abs(x)));
IF (x.im # 0.0) THEN x.im:=x.im/(2.0*Zw); END;
IF (x.re >= 0.0) THEN
x.re:=Zw;
ELSIF (x.im >= 0.0) THEN
x.re:=x.im; x.im:=Zw;
ELSE
x.re:= - x.im; x.im:= - Zw;
END;
RETURN x;
END sqrt;
PROCEDURE exp(x : LONGCMPLX) : LONGCMPLX;
VAR interim : LONGREAL;
z : LONGCMPLX;
BEGIN
interim:=LongMath.exp(x.re);
z.re:=interim*LongMath.cos(x.im);
z.im:=interim*LongMath.sin(x.im);
RETURN z;
END exp;
PROCEDURE ln(x : LONGCMPLX) : LONGCMPLX;
VAR z : LONGCMPLX;
BEGIN
IF (x.re = 0.0) AND (x.im = 0.0) THEN
FatalError('Argument Null (LngCmplxMath.ln) !');
END;
z.re:=0.5*LongMath.ln(x.re*x.re + x.im*x.im);
z.im:=LongMath.arctan(x.im / x.re);
RETURN z;
END ln;
PROCEDURE sin(x : LONGCMPLX) : LONGCMPLX;
VAR z : LONGCMPLX;
BEGIN
z.re:=LongMath.sin(x.re)*cosh(x.im);
z.im:=LongMath.cos(x.re)*sinh(x.im);
RETURN z;
END sin;
PROCEDURE cos(x : LONGCMPLX) : LONGCMPLX;
VAR z : LONGCMPLX;
BEGIN
z.re := LongMath.cos(x.re)*cosh(x.im);
z.im := - LongMath.sin(x.re)*sinh(x.im);
RETURN z;
END cos;
PROCEDURE tan(x : LONGCMPLX) : LONGCMPLX;
VAR z : LONGCMPLX;
yr,yi,Nenner : LONGREAL;
BEGIN
yr:=2.0*x.re; yi:=2.0*x.im;
Nenner:=1.0 / (LongMath.cos(yr) + cosh(yi));
z.re:=LongMath.sin(yr)*Nenner;
z.im:=sinh(yi)*Nenner;
RETURN z;
END tan;
PROCEDURE arctan(x : LONGCMPLX) : LONGCMPLX;
VAR z : LONGCMPLX;
v,w : LONGCMPLX;
BEGIN
v.re := 1.0 - x.im; v.im := x.re;
w.re := 1.0 + x.im; w.im := - x.re;
IF (w.re = 0.0) AND (w.im = 0.0) THEN
FatalError('Unzulaessiges Argument in LngCmplxMath.arctan.');
END;
z:=CDiv(v,w);
IF (z.re = 0.0) AND (z.im = 0.0) THEN
FatalError('Unzulaessiges Argument in LngCmplxMath.arctan.');
END; (* v:=CLn(z) *)
v.re:=0.5*LongMath.ln(z.re*z.re + z.im*z.im);
v.im:=LongMath.arctan(z.im / z.re);
z.re:=0.5*v.im; z.im:=-0.5*v.re;
RETURN z;
END arctan;
PROCEDURE scalarMult(scalar: LONGREAL; z: LONGCMPLX) : LONGCMPLX;
BEGIN
z.re:=scalar*z.re; z.im:=scalar*z.im; RETURN z;
END scalarMult;
END LngCmplxMath.