--- a
+++ b/LngCmplxMath.mod
@@ -0,0 +1,218 @@
+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.