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.