--- a
+++ b/ApproxLib.def
@@ -0,0 +1,329 @@
+DEFINITION MODULE ApproxLib;
+
+ (*------------------------------------------------------------------------*)
+ (* Verschiedene Routinen um Polynme oder Rationale Funktionen an eine *)
+ (* Funktion oder einen Datensatz anzupassen. *)
+ (* Module for fitting polynomals and splines to data sets *)
+ (*------------------------------------------------------------------------*)
+ (* Implementation : Michael Riedl *)
+ (* Licence : GNU Lesser General Public License (LGPL) *)
+ (*------------------------------------------------------------------------*)
+
+ (* $Id: ApproxLib.def,v 1.4 2018/04/09 20:48:43 mriedl Exp mriedl $ *)
+
+FROM LMathLib IMPORT Funktion1;
+
+PROCEDURE Ausgleichsgerade(VAR A,B : LONGREAL;
+ VAR dA,dB : LONGREAL;
+ VAR Xq,Yq : LONGREAL;
+ VAR R : LONGREAL;
+ dim : CARDINAL;
+ T : LONGREAL;
+ VAR X,Y : ARRAY OF LONGREAL);
+
+ (*----------------------------------------------------------------*)
+ (* Berechnet eine Ausgleichsgerade durch X,Y nach der Methode *)
+ (* der kleinsten Fehlerquadratsumme. *)
+ (* *)
+ (* A,B : f(x) = Ax + B *)
+ (* *)
+ (* R : Korellationsfaktor. *)
+ (* T : Student-Faktor *)
+ (* Xq,Yq : Mittelwerte von X,Y. *)
+ (*----------------------------------------------------------------*)
+
+CONST FitStraightLine = Ausgleichsgerade; (* For english users ... *)
+
+ (*----------------------------------------------------------------*)
+ (* Calculates the coefficients of a least square fitted straigt *)
+ (* through the data provided in X,Y *)
+ (* *)
+ (* A,B : f(x) = Ax + B *)
+ (* *)
+ (* R : Corellationsfactor. *)
+ (* T : Student-Faktor *)
+ (* Xq,Yq : Mean values of X,Y *)
+ (*----------------------------------------------------------------*)
+
+PROCEDURE TschebKoeff( U,O : LONGREAL;
+ VAR TKoeff : ARRAY OF LONGREAL;
+ Grad : CARDINAL;
+ Funkt : Funktion1);
+
+ (*----------------------------------------------------------------*)
+ (* Berechnet die Tschebytscheffkoeffizienten TKoeff einer T-Ent- *)
+ (* wicklung des Grades Grad f"ur die Funktion Funkt im Intervall *)
+ (* [U,O]. *)
+ (*----------------------------------------------------------------*)
+
+PROCEDURE PolyKoeff(VAR PKoeff : ARRAY OF LONGREAL; (* Polynomkoeffizienten *)
+ VAR TKoeff : ARRAY OF LONGREAL; (* Tschebytscheffkoeff. *)
+ U,O : LONGREAL; (* Intervall d. Polynom *)
+ Grad : CARDINAL);
+
+ (*----------------------------------------------------------------*)
+ (* Transformiert die Koeffizienten einer Tschebytscheffentwick- *)
+ (* lung im Intervall [-1,1] in die entsprechenden Polynomkoef- *)
+ (* fizienten f"ur das Intervall [U,O] mit *)
+ (* f(x) = \sum_{i=0}^Grad PKoeff_i x^i *)
+ (*----------------------------------------------------------------*)
+
+PROCEDURE RatKoeff( Funkt : Funktion1;
+ U,O : LONGREAL;
+ VAR Ak,Bk : ARRAY OF LONGREAL;
+ n,m : CARDINAL);
+
+ (*----------------------------------------------------------------*)
+ (* Berechnet f"ur das Intervall [U,O] der Eingabefunktion *)
+ (* Funkt(x) die Koeffizienten der rationalen Funktion *)
+ (* R^{n,m}(x) = {\sum_{k=0}^n A_k x^k *)
+ (* \over *)
+ (* 1.0 + \sum_{k=1}^m B_k x^k} *)
+ (* die die Funkt(x) im quadratischen Mittel bestapproximiert. *)
+ (*----------------------------------------------------------------*)
+
+PROCEDURE FitTscheb(dim : CARDINAL;
+ Grad : CARDINAL;
+ VAR X,Y : ARRAY OF LONGREAL;
+ VAR TKoeff : ARRAY OF LONGREAL;
+ VAR U,O : LONGREAL;
+ VAR StdAbw : LONGREAL;
+ VAR IFehl : CARDINAL);
+
+ (*----------------------------------------------------------------*)
+ (* Berechnet die Entwicklungkoeffizienten TKoeff einer *)
+ (* Tschebytscheffentwicklung nach der Methode der minimalen *)
+ (* Fehlerquadratsumme als Anpassung an die Datenpaare X,Y *)
+ (* *)
+ (* dim : Anzahl der Datenpaare (X,Y) *)
+ (* Grad : Grad der T-Entwicklung *)
+ (* U,O : Intervall der T-Entwicklung (R"uckgabeparameter !) *)
+ (* IFehl : Siehe PROCEDURRE PolyFit. *)
+ (*----------------------------------------------------------------*)
+
+PROCEDURE FitPoly( dim : CARDINAL; (* Anzahl Datenpaare X,Y *)
+ VAR X : ARRAY OF LONGREAL;
+ VAR Y : ARRAY OF LONGREAL;
+ Grad : CARDINAL; (* Grad des Polynoms *)
+ VAR PKoeff : ARRAY OF LONGREAL;(* Koeff. des Polynoms *)
+ VAR StdAbw : LONGREAL;
+ VAR IFehl : CARDINAL);
+
+ (*----------------------------------------------------------------*)
+ (* Fittet ein Polynom des Grades Grad an die Datenpaare (X,Y) *)
+ (* nach der Methode der minimalen Fehlerquadratsumme. *)
+ (* *)
+ (* IFehl : *)
+ (* 1 : Anzahl der Datenpunkte kleiner zwei. *)
+ (* 2 : Anzahl der Datenpunkte zu gro3. *)
+ (* 3 : Grad des Polynoms kleiner eins. *)
+ (* 4 : Grad des Polynoms gr"o\3er als die Anzahl der *)
+ (* Datenpunkte. *)
+ (* 5 : Fehlergleichungssystem zu gro\3. *)
+ (* 6 : Keine L"osung gefunden. *)
+ (*----------------------------------------------------------------*)
+
+PROCEDURE FitRat(VAR X,Y : ARRAY OF LONGREAL; (* Datenpunkte *)
+ dim : CARDINAL; (* Anzahl Datenpunkte *)
+ VAR Ak,Bk : ARRAY OF LONGREAL;
+ n,m : CARDINAL;
+ VAR StdAbw : LONGREAL); (* Standartabweichung *)
+
+ (*----------------------------------------------------------------*)
+ (* Berechent die Koeffizenten einer rationalen Funktion *)
+ (* R^{n,m}(x) = {\sum_{k=0}^n A_k x^k *)
+ (* \over *)
+ (* 1.0 + \sum_{k=1}^m B_k x^k} *)
+ (* die die Datenpunkt X,Y im Sinne der minimierten Fehler- *)
+ (* quadratsumme approximiert. *)
+ (* R_i = Y_i - Q^{m,n}(x) (Residuenvektor) *)
+ (*----------------------------------------------------------------*)
+
+PROCEDURE SplineKoeff(VAR X : ARRAY OF LONGREAL;
+ VAR Y : ARRAY OF LONGREAL;
+ VAR B : ARRAY OF LONGREAL;
+ VAR C : ARRAY OF LONGREAL;
+ VAR D : ARRAY OF LONGREAL;
+ N : CARDINAL);
+
+ (*----------------------------------------------------------------*)
+ (* Berechnet die kubischen Splinkoeffizeienten der in X,Y ueber- *)
+ (* gebenen Wertepaare in den Feldern B,C und D. Die Y[i] werden *)
+ (* als der konstanter Teil des Spline-Polynoms behandelt, d.h. *)
+ (* *)
+ (* $ f(x) = Y_i + B_i x + C_i x^2 + D x^3 $ *)
+ (* *)
+ (* Der Parameter N gibt die Anzahl der Datenpunkte in X,Y an. *)
+ (* *)
+ (* Calculates the coefficients of the cubic spline interpolation *)
+ (* through the N data points provided in the arrays X and Y. *)
+ (* The coefficients are returned in arrays B,C and D. The values *)
+ (* of Y form the constant part of the polinominals of degree 3 *)
+ (* *)
+ (* Ref.: Forsythe, George; Malcolm, Michael; Moler, Cleve; *)
+ (* "Computer Methods for Mathematical Computations", *)
+ (* Prentice-Hall, Upper Saddle River, US (1977) *)
+ (*----------------------------------------------------------------*)
+
+PROCEDURE EvalSpline( x : LONGREAL;
+ VAR X : ARRAY OF LONGREAL;
+ VAR Y : ARRAY OF LONGREAL;
+ VAR B : ARRAY OF LONGREAL;
+ VAR C : ARRAY OF LONGREAL;
+ VAR D : ARRAY OF LONGREAL;
+ N : CARDINAL) : LONGREAL;
+
+ (*----------------------------------------------------------------*)
+ (* Auswerten einer Spline-Interpolation am Punkt x [X_1<=x<=X_N] *)
+ (* Die Koeffizienten B,C und D koennen z.B. mit der Prozedur *)
+ (* SplineKoeff berechnet werden. Die Y[i] werden als der *)
+ (* konstanter Teil des Spline-Polynoms behandelt. *)
+ (* *)
+ (* $ f(x) = Y_i + B_i x + C_i x^2 + D x^3 $ *)
+ (* *)
+ (* Evaluate the cubic spline interpolation at point x where *)
+ (* *)
+ (* X[i] <= x <= X[i+1] , i = 0..N-1 *)
+ (* *)
+ (* The spline coefficients B,C and D should have been calculated *)
+ (* by procedure "SplineKoeff". *)
+ (* *)
+ (* Ref.: See SplineKoeff *)
+ (*----------------------------------------------------------------*)
+
+PROCEDURE EvalSplineDerivative( x : LONGREAL;
+ VAR X : ARRAY OF LONGREAL;
+ VAR Y : ARRAY OF LONGREAL;
+ VAR B : ARRAY OF LONGREAL;
+ VAR C : ARRAY OF LONGREAL;
+ VAR D : ARRAY OF LONGREAL;
+ N : CARDINAL;
+ VAR y : LONGREAL;
+ VAR dy : LONGREAL;
+ VAR d2y : LONGREAL;
+ VAR iFehl : INTEGER);
+
+ (*----------------------------------------------------------------*)
+ (* Auswerten einer Spline-Interpolation am Punkt x [X_1<=x<=X_N] *)
+ (* Es wird in y der Interpolationswert zurueckgegeben, in dy die *)
+ (* erste Ableitung am Punkt x und in d2y die zweite Ableitung. *)
+ (* iFehl = 1 wenn x ausserhalb des Wertebereichs liegen sollte. *)
+ (* *)
+ (* Calculates the spline-interpolation as in function EvalSpline *)
+ (* and returns the value in y. In dy the first and in d2y the *)
+ (* second derivative at point x are returend as well. *)
+ (* iFehl = 1 if x is outside of the interval [X_1 ... X_N] *)
+ (*----------------------------------------------------------------*)
+
+PROCEDURE Spline(VAR Xa,Ya : ARRAY OF LONGREAL; (* Eingabevektoren *)
+ N : CARDINAL; (* Anzahl der Eingabewerte *)
+ X1,Xm : LONGREAL;
+ VAR Xn,Yn : ARRAY OF LONGREAL; (* Ausgabevektor *)
+ M : CARDINAL); (* Anzahl Splinepunkte *)
+
+ (*----------------------------------------------------------------*)
+ (* Berechnet aus den Werten in Xa,Ya kubisch splineinterpolierte *)
+ (* Werte in Xn,Yn im Bereich X1,XM. *)
+ (* *)
+ (* Xa,Ya : Vorgegebene Datenpunte / original data points *)
+ (* N : Anzahl der Datenpunkte / number of data points *)
+ (* X1,Xm : Bereich, in dem die berechneten Spline-Punkte *)
+ (* liegen sollen / range for the new values *)
+ (* Xn,Yn : Splineinterpolierte Wertepaare / new data points *)
+ (* M : Anzahl der erzeugten Datenpunkte in den Feldern *)
+ (* Xn,Yn / number of values to be calculated in Xn,Yn *)
+ (* *)
+ (* Calculates M cubic spline interpolated values in the range *)
+ (* X1 .. Xm. These interpolated values fit to the N original data *)
+ (* pairs provided in Xa,Ya. This routine is a driver routine for *)
+ (* SplineKoeff and EvalSpline. *)
+ (*----------------------------------------------------------------*)
+
+CONST MaxR = 32; (* Max degree of nominator polynominal A[0..MaxL-1] *)
+ MaxL = 32; (* Max degree of denominator polynominal B[0..MaxL-1] *)
+
+PROCEDURE RatTsch1( l,r : CARDINAL;
+ ap,ep : LONGREAL;
+ eps : LONGREAL;
+ Fct : Funktion1;
+ W : Funktion1;
+ VAR A,B : ARRAY OF LONGREAL;
+ VAR etabar : LONGREAL;
+ VAR iFehl : INTEGER);
+
+ (* ---------------------------------------------------------------*)
+ (* RatTsch1 determines the coefficients A[0:l] and B[0:r] of the *)
+ (* best fit rational function *)
+ (* *)
+ (* ra(x) := (A[0]+A[1]*x + ... + A[l]*x**l) / *)
+ (* (B[0]+B[1]*x + ....+ B[r]*x**r) *)
+ (* *)
+ (* for the function Fct and the weight function W, which *)
+ (* minimizes the maximum modulus of the weighted error function *)
+ (* *)
+ (* (Fct(x) - ra(x)) * W(x) *)
+ (* *)
+ (* on the interval [ap,ep]. Input parameters are the integer *)
+ (* degree l of the numerator of ra, the integer degree r of the *)
+ (* denominator of ra, the real parameters ap and ep specifying *)
+ (* the left and right endpoint of the approximation interval, *)
+ (* respectively, the real procedure Fct giving the function to be *)
+ (* approximated and the weight function W, also of type real *)
+ (* procedure. Fct must be continuous, W continuous and positive *)
+ (* on [ap,ep]. The real input parameter eps should be 10^(-t+1) *)
+ (* if t is the number of decimal digits in the mantissa of the *)
+ (* arithmetic used. Output parameters are the real arrays *)
+ (* *)
+ (* A[0:l], B[0:r] (with B[0]= 1), *)
+ (* *)
+ (* the real etabar with the property *)
+ (* *)
+ (* abs(etabar) = max(abs(Fct(x) - ra(x))*W(x) : ap <= x <= ep). *)
+ (* *)
+ (* The integer iFehl specifies, whether the approximation was *)
+ (* successful (iFehl=0) or not (iFehl=1,2,3,4,5): if iFehl=1, *)
+ (* then it was impossible to obtain a sufficiently accurate ra(x) *)
+ (* within 10 iterations. The values iFehl=2,3 or 4, indicate the *)
+ (* presence of a pole of the approximant ra(x) within [ap,ep]. *)
+ (* If iFehl=5 then the rational interpolation failed *)
+ (* *)
+ (* ==> l,r : Degree of rational approximation *)
+ (* See also definition of constans MaxR and MaxL *)
+ (* ==> ap,ep : Anfangspunkt,Endpunkt (start and end point) *)
+ (* ==> eps : see text *)
+ (* ==> Fct : Function to be approximated *)
+ (* ==> W : Weight function, if unsure us w(x) = 1 *)
+ (* <== A,B : Computes coefficients of the rational approx. *)
+ (* <== etabar : see text *)
+ (* <== iFehl : Error indicator *)
+ (* 0 : Approximation was successful *)
+ (* 1 : The number of iterations admitted by the *)
+ (* limit ito dit not suffice to complete the *)
+ (* Remes-Algorithm. Increasing the number *)
+ (* ito might help. *)
+ (* 2 : A rational function with a pole has *)
+ (* occured - deteted by hitting a point in *)
+ (* the search procedure where error function *)
+ (* > 10**6 or rational function > 10**20 *)
+ (* 3 : A pole of a rational function had been *)
+ (* detected - denominator changes sign *)
+ (* within the interval *)
+ (* 4 : Case 1 and 2 or 3 are combined *)
+ (* 5 : Rational interpolation fails *)
+ (* *)
+ (* HINT: If iFehl = -1 on ENTRY the nuber of iterations and the *)
+ (* return code of iFehl are printed to the console at the *)
+ (* end of the routine. *)
+ (* *)
+ (* This procedure is a Modula-2 translation of the Algol60 *)
+ (* procedure rattsch1 - a Algol60 version is available as well *)
+ (* *)
+ (* Ref: Werner, H; Stoer, J.; Bommas, W.; "Rational Chebyshev *)
+ (* Approximation", Handbook Series Methods of Approximation, *)
+ (* Numerische Mathematik, Springer, Berlin (DE) (1967) *)
+ (* *)
+ (* For L2 approximations the procedure ApproxLib.RatKoeff can *)
+ (* be used. *)
+ (* ---------------------------------------------------------------*)
+
+END ApproxLib.