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.
