unit fit;

Interface


const mma=10;
      xvector_size=1000;

type
   x_vector=array[1..xvector_size] of double;
   masktype=array[1..mma] of boolean;
   paramtype=array[1..mma] of double;
   fitfunctype=function(x:double;parameter:paramtype):double;
   fitfuncvector=array[1..mma] of fitfunctype;
   svector=array[0..12] of double;
   matrix=array[0..12] of svector;
   mvector=array[0..90] of double;
   kmatrix=array[0..90] of mvector;
   bmatrix=^kmatrix;
   intvector=array[0..12] of longint;
   intfunctype=function(x:double):double;



var
  fitdfunc:fitfuncvector;
  ffitfunc:fitfunctype;
  axisconvx,axisconvy:fitfunctype;
  linearcombination:boolean;
  mask:masktype;
  message:string;
  fepsilon:double;
  mfit:longint;

Function power(x:double;k:double):double; {Az x k-adik hatvanyat adja vissza x>0}
Function dpower(x:double;k:double):double; {A power k szerinti derivaltja az x helyen}
Function Integral(func:intfunctype;a,b:double):double;
Procedure fitstep(var x,y,sigma:x_vector;ndata:longint;var a:paramtype;var mask:masktype;var chi,alamda:double);
Procedure polifit(var x,y,sigma:x_vector;ndata:longint;s:double;var a:paramtype);
Function dfitfunc(x:double;a:paramtype):double;
procedure order(var x,y:x_vector;ndata:longint);
procedure splineconstant(var x,y:x_vector;m:intvector;s:longint;g:double);
procedure egyenlet(var a:bmatrix;var s:mvector;n:longint);
function sf(xr:double;var x:x_vector;m:intvector):double;
function sdf(xr:double;var x:x_vector;m:intvector):double;
function linear(x:double;a:paramtype):double;
function polinom(x:double;a:paramtype):double;

{*************************************************************}

Implementation


TYPE
   glndata = x_vector;
   glmma = paramtype;
   gllista = ARRAY [1..mma] OF longint;
   gllim=gllista;
   glmmabymma = ARRAY [1..mma,1..mma] OF double;
   glncabynca=glmmabymma;
   glcovar=glncabynca;

Var
   ff                       :fitfunctype;
   funcs:procedure(x:double; var a:glmma;var yfit:double;var dyda:glmma; mma:longint);
   glochisq                 : double;
   glbeta                   : glmma;
   i                        :longint;
   csp                      :mvector;




FUNCTION power(x:double;k:double):double; {Az x k-adik hatvanyat adja vissza}

   var p:double;

   begin
    if abs(x)<1e-10 then power:=0
    else
    if abs(abs(k)-round(abs(k)))<1e-8 then
      begin
        p:=1;
        for i:=1 to round(abs(k)) do p:=p*x;
        if k<0 then p:=1/p;
        power:=p;
      end
    else power:=exp(k*ln(abs(x)));
   end;



FUNCTION dpower(x:double;k:double):double; {A power k szerinti derivaltja az x helyen}

   begin
    if x=0 then dpower:=0 else
    dpower:=ln(x)*power(x,k);
   end;



function dfitfunc(x:double;a:paramtype):double;

 begin
  dfitfunc:=(ffitfunc(x+fepsilon,a)-ffitfunc(x,a))/fepsilon;
 end;


{***************************************************}

Function Integral(func:intfunctype;a,b:double):double;

 var glit:longint;
     s:double;

PROCEDURE trapzd(a,b: double; VAR s: double; n: longint);
{ Programs calling TRAPZD must provide a function
func(x:double):double which is to be integrated. They must
also define the variable
VAR
   glit: longint;
in the main routine. }
VAR
   j: longint;
   x,tnm,sum,del: double;
BEGIN
   IF (n = 1) THEN BEGIN
      s := 0.5*(b-a)*(func(a)+func(b));
      glit := 1
   END
   ELSE BEGIN
      tnm := glit;
      del := (b-a)/tnm;
      x := a+0.5*del;
      sum := 0.0;
      FOR j := 1 to glit DO BEGIN
         sum := sum+func(x);
         x := x+del
      END;
      s := 0.5*(s+(b-a)*sum/tnm);
      glit := 2*glit
   END
END;


PROCEDURE qsimp(a,b: double; VAR s: double);
LABEL 99;
CONST
   eps=1.0e-6;
   jmax=20;
VAR
   j: longint;
   st,ost,os: double;
BEGIN
   ost := -1.0e30;
   os := -1.0e30;
   FOR j := 1 to jmax DO BEGIN
      trapzd(a,b,st,j);
      s := (4.0*st-ost)/3.0;
      IF (abs(s-os) < eps*abs(os)) THEN GOTO 99;
      os := s;
      ost := st
   END;
   writeln ('pause in QSIMP - too many steps'); readln;
99:   END;


Begin
 qsimp(a,b,s);
 Integral:=s;
end;

{***************************************************}

procedure nfuncs(x:double; var a:glmma;var yfit:double;var dyda:glmma; mma:longint);

{A dyda[i] a fittelendo fuggveny a[i] parameter szerinti derivaltja
 Az yfit pegig a fittelendo fuggveny az x helyen}

   var i:longint;


   begin
     yfit:=0;
     for i:=1 to mfit do
       begin
        ff:=fitdfunc[i];
        dyda[i]:=ff(x,a);
        if linearcombination then
          yfit:=yfit+a[i]*dyda[i];
       end;
     if not linearcombination then
	begin
	yfit:=ffitfunc(x,a);
	end;
   end;



PROCEDURE gaussj(VAR a: glncabynca; n,np: longint;
       VAR b: glncabynca; m,mp: longint);
{ Programs using GAUSSJ must define the types
TYPE
   glnpbynp = ARRAY [1..np,1..np] OF double;
   glnpbymp = ARRAY [1..np,1..mp] OF double;
   glnp = ARRAY [1..np] OF longint;
in the main routine. }
VAR
   big,dum,pivinv: double;
   i,icol,irow,j,k,l,ll: longint;
   indxc,indxr,ipiv: gllim;
BEGIN
   FOR j := 1 to n DO BEGIN
      ipiv[j] := 0
   END;
   FOR i := 1 to n DO BEGIN
      big := 0.0;
      FOR j := 1 to n DO BEGIN
         IF (ipiv[j] <> 1) THEN BEGIN
            FOR k := 1 to n DO BEGIN
               IF (ipiv[k] = 0) THEN BEGIN
                  IF (abs(a[j,k]) >= big) THEN BEGIN
                     big := abs(a[j,k]);
                     irow := j;
                     icol := k
                  END
               END ELSE IF (ipiv[k] > 1) THEN BEGIN
                  message:='pause 1 in GAUSSJ - singular matrix';

                  Exit;
               END
            END
         END
      END;
      ipiv[icol] := ipiv[icol]+1;
      IF (irow <> icol) THEN BEGIN
         FOR l := 1 to n DO BEGIN
            dum := a[irow,l];
            a[irow,l] := a[icol,l];
            a[icol,l] := dum
         END;
         FOR l := 1 to m DO BEGIN
            dum := b[irow,l];
            b[irow,l] := b[icol,l];
            b[icol,l] := dum
         END
      END;
      indxr[i] := irow;
      indxc[i] := icol;
      IF (a[icol,icol] = 0.0) THEN BEGIN
         message:='pause 1 in GAUSSJ - singular matrix';

         Exit;
      END;
      pivinv := 1.0/a[icol,icol];
      a[icol,icol] := 1.0;
      FOR l := 1 to n DO BEGIN
         a[icol,l] := a[icol,l]*pivinv
      END;
      FOR l := 1 to m DO BEGIN
         b[icol,l] := b[icol,l]*pivinv
      END;
      FOR ll := 1 to n DO BEGIN
         IF (ll <> icol) THEN BEGIN
            dum := a[ll,icol];
            a[ll,icol] := 0.0;
            FOR l := 1 to n DO BEGIN
               a[ll,l] := a[ll,l]-a[icol,l]*dum
            END;
            FOR l := 1 to m DO BEGIN
               b[ll,l] := b[ll,l]-b[icol,l]*dum
            END
         END
      END
   END;
   FOR l := n DOWNTO 1 DO BEGIN
      IF (indxr[l] <> indxc[l]) THEN BEGIN
         FOR k := 1 to n DO BEGIN
            dum := a[k,indxr[l]];
            a[k,indxr[l]] := a[k,indxc[l]];
            a[k,indxc[l]] := dum
         END
      END
   END
END;



PROCEDURE covsrt(VAR covar: glcovar; ncvm: longint; ma: longint;
       var lista: gllista; mfit: longint);
{* Programs using routine COVSRT must define the types
TYPE
   glcovar = ARRAY [1..ncvm,1..ncvm] OF double;
   gllista = ARRAY [1..mfit] OF longint;
in the calling program. *}
VAR
   j,i: longint;
   swap: double;
BEGIN
   FOR j := 1 to ma-1 DO BEGIN
      FOR i := j+1 to ma DO BEGIN
         covar[i,j] := 0.0
      END
   END;
   FOR i := 1 to mfit-1 DO BEGIN
      FOR j := i+1 to mfit DO BEGIN
         IF (lista[j] > lista[i])  THEN BEGIN
            covar[lista[j],lista[i]] := covar[i,j]
         END ELSE BEGIN
            covar[lista[i],lista[j]] := covar[i,j]
         END
      END
   END;
   swap := covar[1,1];
   FOR j := 1 to ma DO BEGIN
      covar[1,j] := covar[j,j];
      covar[j,j] := 0.0
   END;
   covar[lista[1],lista[1]] := swap;
   FOR j := 2 to mfit DO BEGIN
      covar[lista[j],lista[j]] := covar[1,j]
   END;
   FOR j := 2 to ma DO BEGIN
      FOR i := 1 to j-1 DO BEGIN
         covar[i,j] := covar[j,i]
      END
   END
END;




PROCEDURE mrqcof(var x,y,sig: glndata; ndata: longint;
       VAR a: glmma; mma: longint; var lista: gllista;
       mfit: longint; VAR alpha: glncabynca;
       VAR beta: glmma; nalp: longint; VAR chisq: double);


{ Programs using routine MRQMIN must provide a
PROCEDURE funcs(xx:double; a:glmma; yfit:double; dyda:glmma; mma:longint);
that evaluates the fitting function yfit and its derivatives dyda
with respect to the parameters a at point xx.  Also they
must define the types
TYPE
   glndata = ARRAY [1..ndata] OF double;
   glmma = ARRAY [1..mma] OF double;
   gllista = ARRAY [1..mma] OF longint;
   glnalbynal = ARRAY [1..nalp,1..nalp] OF double;
in the main routine }


VAR
   k,j,i: longint;
   ymod,wt,sig2i,dy: double;
   dyda: glmma;
BEGIN

   FOR j := 1 to mfit DO BEGIN
      FOR k := 1 to j DO BEGIN
         alpha[j,k] := 0.0
      END;
      beta[j] := 0.0
   END;
   chisq := 0.0;
   FOR i := 1 to ndata DO BEGIN
      funcs(x[i],a,ymod,dyda,mma);
      sig2i := 1.0/(sig[i]*sig[i]);
      dy := y[i]-ymod;
      FOR j := 1 to mfit DO BEGIN
         wt := dyda[lista[j]]*sig2i;
         FOR k := 1 to j DO BEGIN
             alpha[j,k] := alpha[j,k]+wt*dyda[lista[k]]
         END;
         beta[j] := beta[j]+dy*wt
      END;
      chisq := chisq+dy*dy*sig2i;
   END;
   FOR j := 2 to mfit DO BEGIN
      FOR k := 1 to j-1 DO BEGIN
         alpha[k,j] := alpha[j,k]
      END
   END;

END;



PROCEDURE mrqmin(var x,y,sig: glndata; ndata: longint;
         VAR a: glmma; mma: longint; var lista: gllista;
         mfit: longint; VAR covar,alpha: glncabynca;
         nca: longint; VAR chisq,alamda: double);

{ Programs using routine MRQMIN must define the types
TYPE
      glndata = ARRAY [1..ndata] OF double;
      glmma = ARRAY [1..mma] OF double;
      gllista = ARRAY [1..mma] OF longint;
      glncabynca = ARRAY [1..nca,1..nca] OF double;
and the variables
VAR
      glochisq: double;  glbeta: glmma;
in the main routine. Also note that this routine calls MRQCOF, which
requires a user-defined procedure FUNCS, described in that routine.	}


VAR
      k,kk,j,ihit: longint;  atry,da: glmma;  oneda: glncabynca;
BEGIN
      IF (alamda < 0.0) THEN BEGIN
        kk := mfit+1;
        FOR j := 1 TO mma DO BEGIN
          ihit := 0;
          FOR k := 1 TO mfit DO BEGIN
            IF (lista[k] = j) THEN ihit := ihit+1
          END;
          IF (ihit = 0) THEN BEGIN
            lista[kk] := j;  kk := kk+1  END
          ELSE IF (ihit > 1) THEN BEGIN
            writeln('pause 1 in routine MRQMIN');
            writeln('Improper permutation in LISTA'); readln
          END  END;
        IF (kk <> (mma+1)) THEN BEGIN
          writeln('pause 2 in routine MRQMIN');
          writeln('Improper permutation in LISTA'); readln
        END;
        alamda := 0.001;
        if linearcombination then alamda:=0.0;
        mrqcof(x,y,sig,ndata,a,mma,lista,mfit,alpha,glbeta,nca,chisq);
        glochisq := chisq;
        FOR j := 1 TO mma DO BEGIN
          atry[j] := a[j]  END; end;
      FOR j := 1 TO mfit DO BEGIN
        FOR k := 1 TO mfit DO BEGIN
          covar[j,k] := alpha[j,k]  END;
        covar[j,j] := alpha[j,j]*(1.0+alamda);  oneda[j,1] := glbeta[j]  END;
      gaussj(covar,mfit,nca,oneda,1,1);
      if message<>'' then halt;
      FOR j := 1 TO mfit DO da[j] := oneda[j,1];
      IF (alamda = 0.0) and (not linearcombination) THEN BEGIN
        covsrt(covar,nca,mma,lista,mfit);
        exit;
      END;
      FOR j := 1 TO mfit DO BEGIN
        atry[lista[j]] := a[lista[j]]+da[j]  END;
      If linearcombination then glochisq:=1.1*chisq Else
       begin
         mrqcof(x,y,sig,ndata,atry,mma,lista,mfit,covar,da,nca,chisq);
       end;
      IF (chisq <= glochisq) THEN BEGIN
        if alamda>1e-29 then alamda := 0.1*alamda;  glochisq := chisq;
        FOR j := 1 TO mfit DO BEGIN
          FOR k := 1 TO mfit DO BEGIN
            alpha[j,k] := covar[j,k]  END;
          glbeta[j] := da[j];  a[lista[j]] := atry[lista[j]]  END  END
      ELSE BEGIN
        alamda := 10.0*alamda;  chisq := glochisq  END;
   END;


Function summ(x:double;a:paramtype):double;

  var i:longint;
      s:double;
      ff:fitfunctype;

  begin
   s:=0;
   for i:=1 to mfit do
     begin
      ff:=fitdfunc[i];
      s:=a[i]*ff(x,a)+s;
     end;
   summ:=s;
  end;

var  alpha,covar:glmmabymma;

Procedure fitstep(var x,y,sigma:x_vector;ndata:longint;var a:paramtype;var mask:masktype;var chi,alamda:double);

  var l          :gllista;
      i,mafit  :longint;

      aa:glmma;

  begin
   mafit:=0;
   message:='';
   funcs:=@nfuncs;
   if linearcombination then
     begin
       ffitfunc:=@summ;
     end;
   for i:=1 to mfit do
    begin
     if mask[i]=true then
       begin
        inc(mafit);
        l[mafit]:=i;
       end;
    end;
   for i:=1 to mma do aa[i]:=a[i];
       mrqmin(x,y,sigma,ndata,aa,mma,l,mafit,covar,alpha,mma,chi,alamda);
      for i:=1 to mma do a[i]:=aa[i];
  end;


procedure pfuncs(x:double; var a:glmma;var yfit:double;var dyda:glmma; mma:longint);

   begin
     dyda[1]:=power(x,a[10]);
     yfit:=a[1]*dyda[1];
     for i:=2 to mfit do
       begin
         dyda[i]:=dyda[i-1]*x;
         yfit:=yfit+a[i]*dyda[i];
       end;
   end;



function polinom(x:double;a:paramtype):double;

  var yfit,p:double;

    begin
     p:=power(x,a[10]);
     yfit:=a[1]*p;
     for i:=2 to mfit do
       begin
        p:=p*x;
        yfit:=yfit+a[i]*p;
       end;
     polinom:=yfit;
    end;


Procedure polifit(var x,y,sigma:x_vector;ndata:longint;s:double;var a:paramtype);

  var l          :gllista;
      i        :longint;
      alpha,covar:glmmabymma;
      alamda,chi:double;

  begin
   a[10]:=s;
   funcs:=@pfuncs;
   ffitfunc:=@polinom;
   linearcombination:=true;
   alamda:=-1;
   for i:=1 to mfit do
        l[i]:=i;
   mrqmin(x,y,sigma,ndata,a,mma,l,mfit,covar,alpha,mma,chi,alamda);
  end;


{****************************************************}

 procedure egyenlet(var a:bmatrix;var s:mvector;n:longint);

 var I,J,K,L:longint;
     x      :double;

 begin
   FOR K:=1 TO N-1 do
     begin
       X:=ABS(a^[K,K]);
       L:=K;
       FOR I:=K+1 TO N do
         begin
           IF X<ABS(a^[I,K]) THEN
             begin
               X:=ABS(a^[I,K]);
               L:=I;
             end;
         end;
       FOR J:=K TO N+1 do
         begin
            X:=a^[K,J];
            a^[K,J]:=a^[L,J];
            a^[L,J]:=X;
         end;
       FOR I:=K+1 TO N do
          FOR J:=K+1 TO N+1 do
             a^[I,J]:=a^[I,J]-a^[I,K]*a^[K,J]/a^[K,K];
     end;
   S[N]:=a^[N,N+1]/a^[N,N];
   FOR I:=N-1 downTO 1 do
     begin
       S[I]:=a^[I,N+1];
       FOR J:=I+1 TO N do
          S[I]:=S[I]-a^[I,J]*S[J];
       S[I]:=S[I]/a^[I,I];
     end;
 end;




 var
    //ii                 :longint;
    Q,Q1,Q2,QQ    :double;
    a                   :bmatrix;
    //ddmin,ddmax         :double;


 procedure splinematrix(var x,y:x_vector;m:intvector;s:longint;g:double) ;

  Var
     I,J,K,L,o    :longint;
     pp:double;

  begin
    FOR I:=1 TO 9*S-4 do FOR J:=1 TO 9*S-3 do a^[I,J]:=0;
    FOR K:=M[0] TO M[1] do
    begin
     pp:=1;
     for o:=1 to 5 do begin a^[1,o]:=a^[1,o]+pp; pp:=pp*x[K];end;
      a^[1,9*S-3]:=a^[1,9*S-3]+y[K];
    end;
    a^[1,6]:=-0.5;
    FOR L:=2 TO S-1 do
      begin
        a^[L,9*(L-2)+6]:=0.5;
        a^[L,9*(L-1)+6]:=-0.5;
        J:=9*(L-1);
        FOR K:=M[L-1] TO M[L] do
          begin
           pp:=1;
           for o:=1 to 5 do begin a^[L,J+o]:=a^[L,J+o]+pp;pp:=pp*x[K];end;
           a^[L,9*S-3]:=a^[L,9*S-3]+y[K];
          end;
        end;
    L:=S;
    J:=9*(L-1);
    FOR K:=M[L-1] TO M[L] do
      begin
        pp:=1;
        for o:=1 to 5 do begin a^[S,J+o]:=a^[S,J+o]+pp;pp:=pp*x[K];end;
        a^[S,9*S-3]:=a^[S,9*S-3]+y[K];
      end;
    a^[S,9*(L-2)+6]:=0.5;
    { ******* S SOR KESZ *******}
    FOR K:=M[0] TO M[1] do
      begin
        pp:=x[K];
        for o:=1 to 5 do begin a^[S+1,o]:=a^[S+1,o]+pp;pp:=pp*x[K];end;
        a^[S+1,9*S-3]:=a^[S+1,9*S-3]+y[K]*x[K];
      end;
    a^[S+1,6]:=-0.5*x[M[1]];
    a^[S+1,7]:=-0.5;
    FOR L:=2 TO S-1 do
      begin
        J:=9*(L-1);
        FOR K:=M[L-1] TO M[L] do
          begin
            pp:=x[K];
            for o:=1 to 5 do begin a^[S+L,J+o]:=a^[S+L,J+o]+pp;pp:=pp*x[K];end;
            a^[S+L,9*S-3]:=a^[S+L,9*S-3]+y[K]*x[K];
          end;
        a^[S+L,9*(L-2)+6]:=0.5*x[M[L-1]];
        a^[S+L,9*(L-2)+7]:=0.5;
        a^[S+L,9*(L-1)+6]:=-0.5*x[M[L]];
        a^[S+L,9*(L-1)+7]:=-0.5;
      end;
    L:=S;
    J:=9*(L-1);
    FOR K:=M[L-1] TO M[L] do
      begin
        pp:=x[K];
        for o:=1 to 5 do begin a^[2*S,J+o]:=a^[2*S,J+o]+pp;pp:=pp*x[K];end;
        a^[2*S,J+6]:=a^[2*S,J+6]+x[K]*y[K];
      end;
    a^[2*S,9*(L-2)+6]:=0.5*x[M[L-1]];
    a^[2*S,9*(L-2)+7]:=0.5;
    {******* 2*S SOR KESZ ********}
    Q2:=x[M[1]];
    Q1:=x[M[0]];
    FOR K:=M[0] TO M[1] do
      begin
        pp:=x[K]*x[K];
        for o:=1 to 5 do begin a^[2*S+1,o]:=a^[2*S+1,o]+pp;pp:=pp*x[K];end;
        a^[2*S+1,3]:=a^[2*S+1,3]+G*(x[M[1]]-x[M[0]]);
        a^[2*S+1,4]:=a^[2*S+1,4]+G*1.5*(Q2*Q2-Q1*Q1);
        a^[2*S+1,5]:=a^[2*S+1,5]+G*2*(Q2*Q2*Q2-Q1*Q1*Q1);
        a^[2*S+1,9*S-3]:=a^[2*S+1,9*S-3]+y[K]*x[K]*x[K];
      end;
    a^[2*S+1,6]:=-0.5*x[M[1]]*x[M[1]];
    a^[2*S+1,7]:=-x[M[1]];
    a^[2*S+1,8]:=-1;
    FOR L:=2 TO S-1 do
      begin
        J:=9*(L-1);
        Q1:=x[M[L-1]];
        Q2:=x[M[L]];
        FOR K:=M[L-1] TO M[L] do
          begin
            pp:=x[K]*x[K];
            for o:=1 to 5 do begin a^[2*S+L,J+o]:=a^[2*S+L,J+o]+pp;pp:=pp*x[K];end;
            a^[2*S+L,J+3]:=a^[2*S+L,J+3]+G*(Q2-Q1);
            a^[2*S+L,J+4]:=a^[2*S+L,J+4]+G*1.5*(Q2*Q2-Q1*Q1);
            a^[2*S+L,J+5]:=a^[2*S+L,J+5]+G*2*(Q2*Q2*Q2-Q1*Q1*Q1);
            a^[2*S+L,9*S-3]:=a^[2*S+L,9*S-3]+y[K]*x[K]*x[K];
           end;
        a^[2*S+L,9*(L-1)+6]:=-0.5*x[M[L]]*x[M[L]];
        a^[2*S+L,9*(L-2)+6]:=0.5*x[M[L-1]]*x[M[L-1]];
        a^[2*S+L,9*(L-1)+7]:=-x[M[L]];
        a^[2*S+L,9*(L-2)+7]:=x[M[L-1]];
        a^[2*S+L,9*(L-1)+8]:=-1;
        a^[2*S+L,9*(L-2)+8]:=1;
      end;
    L:=S;
    J:=9*(L-1);
    Q2:=x[M[L]];
    Q1:=x[M[L-1]];


  end;


 procedure folytatas(var x,y:x_vector;m:intvector;s,L:longint;g:double) ;

  var J,K,o:longint;
      pp:double;

  begin
    J:=0;
    FOR K:=M[L-1] TO M[L] do
      begin
        pp:=x[K]*x[K];
        for o:=1 to 5 do begin a^[3*S,J+o]:=a^[3*S,J+o]+pp;pp:=pp*x[K];end;
        a^[3*S,J+3]:=a^[3*S,J+3]+G*(Q2-Q1);
        a^[3*S,J+4]:=a^[3*S,J+4]+G*(Q2*Q2-Q1*Q1);
        a^[3*S,J+5]:=a^[3*S,J+5]+G*2*(Q2*Q2*Q2-Q1*Q1*Q1);
        a^[3*S,9*S-3]:=a^[3*S,9*S-3]+y[K]*x[K]*x[K];
      end;
    a^[3*S,9*(L-2)+6]:=0.5*x[M[L-1]]*x[M[L-1]];
    a^[3*S,9*(L-2)+7]:=x[M[L-1]];
    a^[3*S,9*(L-2)+8]:=1;

    {******* 3*S SOR KESZ ********}
    Q2:=x[M[1]];
    Q1:=x[M[0]];
    FOR K:=M[0] TO M[1] do
      begin
        pp:=x[K]*x[K]*x[K];
        for o:=1 to 5 do begin a^[3*S+1,o]:=a^[3*S+1,o]+pp;pp:=pp*x[K];end;
        a^[3*S+1,3]:=a^[3*S+1,3]+G*0.5*(Q2*Q2*-Q1*Q1);
        a^[3*S+1,4]:=a^[3*S+1,4]+G*(Q2*Q2*Q2-Q1*Q1*Q1);
        a^[3*S+1,5]:=a^[3*S+1,5]+G*1.5*(Q2*Q2*Q2*Q2-Q1*Q1*Q1*Q1);
        a^[3*S+1,9*S-3]:=a^[3*S+1,9*S-3]+y[K]*x[K]*x[K]*x[K];
      end;
    Q:=x[M[1]];
    a^[3*S+1,6]:=-0.5*Q*Q*Q;
    a^[3*S+1,7]:=-1.5*Q*Q;
    a^[3*S+1,8]:=-3*Q;
    a^[3*S+1,9]:=-3;
    FOR L:=2 TO S-1 do
      begin
        J:=9*(L-1);
        FOR K:=M[L-1] TO M[L] do
          begin
            Q:=x[K];
            Q2:=x[M[L]];
            Q1:=x[M[L-1]];
            pp:=Q*Q*Q;
            for o:=1 to 5 do begin a^[3*S+L,J+o]:=a^[3*S+L,J+o]+pp;pp:=pp*Q;end;
            a^[3*S+L,J+3]:=a^[3*S+L,J+3]+G*0.5*(Q2*Q2-Q1*Q1);
            a^[3*S+L,J+4]:=a^[3*S+L,J+4]+G*(Q2*Q2*Q2*-Q1*Q1*Q1);
            a^[3*S+L,J+5]:=a^[3*S+L,J+5]+G*1.5*(Q2*Q2*Q2*Q2-Q1*Q1*Q1*Q1);
            a^[3*S+L,9*S-3]:=a^[3*S+L,9*S-3]+y[K]*Q*Q*Q;
          end;
        a^[3*S+L,9*(L-2)+6]:=0.5*x[M[L-1]]*x[M[L-1]]*x[M[L-1]];
        a^[3*S+L,9*(L-1)+6]:=-0.5*x[M[L]]*x[M[L]]*x[M[L]];
        a^[3*S+L,9*(L-2)+7]:=1.5*x[M[L-1]]*x[M[L-1]];
        a^[3*S+L,9*(L-1)+7]:=-1.5*x[M[L]]*x[M[L]];
        a^[3*S+L,9*(L-2)+8]:=3*x[M[L-1]];
        a^[3*S+L,9*(L-1)+8]:=-3*x[M[L]];
        a^[3*S+L,9*(L-2)+9]:=3;
        a^[3*S+L,9*(L-1)+9]:=-3;
      end;
    L:=S;
    J:=9*(L-1);
    Q2:=x[M[L]];
    Q1:=x[M[L-1]];
    FOR K:=M[L-1] TO M[L] do
      begin
        Q:=x[K];
        pp:=Q*Q*Q;
        for o:=1 to 5 do begin a^[4*S,J+o]:=a^[4*S,J+o]+pp;pp:=pp*Q;end;
        a^[4*S,J+3]:=a^[4*S,J+3]+G*0.5*(Q2*Q2-Q1*Q1);
        a^[4*S,J+4]:=a^[4*S,J+4]+G*(Q2*Q2*Q2-Q1*Q1*Q1);
        a^[4*S,J+5]:=a^[4*S,J+5]+G*1.5*(Q2*Q2*Q2*Q2+Q1*Q1*Q1*Q1);
        a^[4*S,9*S-3]:=a^[4*S,9*S-3]+y[K]*Q*Q*Q;
      end;
    Q:=x[M[L-1]];
    a^[4*S,9*(L-2)+6]:=0.5*Q*Q*Q;
    a^[4*S,9*(L-2)+7]:=1.5*Q*Q;
    a^[4*S,9*(L-2)+8]:=3*Q;
    a^[4*S,9*(L-2)+9]:=3;
    {******* 4*S SOR KESZ ********}
    Q2:=x[M[1]];
    Q1:=x[M[0]];
    FOR K:=M[0] TO M[1] do
      begin
        Q:=x[K];
        pp:=Q*Q*Q*Q;
        for o:=1 to 5 do begin a^[4*S+1,J+o]:=a^[4*S+1,J+o]+pp;pp:=pp*Q;end;
        a^[4*S+1,3]:=a^[4*S+1,3]+G*1/3*(Q2*Q2*Q2-Q1*Q1*Q1);
        a^[4*S+1,4]:=a^[4*S+1,4]+G*3/4*(Q2*Q2*Q2*Q2-Q1*Q1*Q1*Q1);
        a^[4*S+1,5]:=a^[4*S+1,5]+Q*Q*Q*Q*Q*Q*Q*Q+G*6/5*(Q2*Q2*Q2*Q2*Q2-Q1*Q1*Q1*Q1*Q1);
        a^[4*S+1,9*S-3]:=a^[4*S+1,9*S-3]+y[K]*Q*Q*Q*Q;
      end;
    Q:=x[M[1]];
    a^[4*S+1,6]:=-0.5*Q*Q*Q*Q;
    a^[4*S+1,7]:=-2*Q*Q*Q;
    a^[4*S+1,8]:=-6*Q*Q;
    a^[4*S+1,9]:=-12*Q;
    FOR L:=2 TO S-1 do
      begin
        J:=9*(L-1);
        FOR K:=M[L-1] TO M[L] do
          begin
            Q:=x[K];
            Q2:=x[M[L]];
            Q1:=x[M[L-1]];
            pp:=Q*Q*Q*Q;
            for o:=1 to 5 do begin a^[4*S+L,J+o]:=a^[4*S+L,J+o]+pp;pp:=pp*Q;end;
            a^[4*S+L,J+3]:=a^[4*S+L,J+3]+G*1/3*(Q2*Q2-Q1*Q1);
            a^[4*S+L,J+4]:=a^[4*S+L,J+4]+G*3/4*(Q2*Q2*Q2-Q1*Q1*Q1);
            a^[4*S+L,J+5]:=a^[4*S+L,J+5]+G*6/5*(Q2*Q2*Q2*Q2-Q1*Q1*Q1*Q1);
            a^[4*S+L,9*S-3]:=a^[4*S+L,9*S-3]+y[K]*Q*Q*Q*Q;
          end;
        a^[4*S+L,9*(L-2)+6]:=0.5*x[M[L-1]]*x[M[L-1]]*x[M[L-1]]*x[M[L]];
        a^[4*S+L,9*(L-1)+6]:=-0.5*x[M[L]]*x[M[L]]*x[M[L]]*x[M[L]];
        a^[4*S+L,9*(L-2)+7]:=2*x[M[L-1]]*x[M[L-1]]*x[M[L-1]];
        a^[4*S+L,9*(L-1)+7]:=-2*x[M[L]]*x[M[L]]*x[M[L]];
        a^[4*S+L,9*(L-2)+8]:=6*x[M[L-1]]*x[M[L-1]];
        a^[4*S+L,9*(L-1)+8]:=-6*x[M[L]]*x[M[L]];
        a^[4*S+L,9*(L-2)+9]:=12*x[M[L-1]];
        a^[4*S+L,9*(L-1)+9]:=-12*x[M[L]];
      end;
    L:=S;
    J:=9*(L-1);
    Q2:=x[M[L]];
    Q1:=x[M[L-1]];
    FOR K:=M[L-1] TO M[L] do
      begin
        Q:=x[K];
        pp:=Q*Q*Q*Q;
        for o:=1 to 3 do begin a^[5*S,J+o]:=a^[5*S,J+o]+pp;pp:=pp*Q;end;
        a^[5*S,J+3]:=a^[5*S,J+3]+G*1/3*(Q2*Q2*Q2-Q1*Q1*Q1);
        a^[5*S,J+4]:=a^[5*S,J+4]+Q*Q*Q*Q*Q*Q*Q*G*3/4*(Q2*Q2*Q2*Q2-Q1*Q1*Q1*Q1);
        a^[5*S,J+5]:=a^[5*S,J+5]+Q*Q*Q*Q*Q*Q*Q*Q*G*6/5*(Q2*Q2*Q2*Q2*Q2-Q1*Q1*Q1*Q1*Q1);
        a^[5*S,9*S-3]:=a^[5*S,(9*S-3)]+y[K]*Q*Q*Q*Q;
       end;
    Q:=x[M[L-1]];
    a^[5*S,9*(L-2)+6]:=0.5*Q*Q*Q*Q;
    a^[5*S,9*(L-2)+7]:=2*Q*Q*Q;
    a^[5*S,9*(L-2)+8]:=6*Q*Q;
    a^[5*S,9*(L-2)+9]:=12*Q;
    {******* 5*S SOR KESZ ********}
    FOR L:=1 TO S-1 do
      begin
        QQ:=x[M[L]];
        a^[5*S+L,9*(L-1)+1]:=1;
        a^[5*S+L,9*L+1]:=-1;
        a^[5*S+L,9*(L-1)+2]:=QQ;
        a^[5*S+L,9*L+2]:=-QQ;
        a^[5*S+L,9*(L-1)+3]:=QQ*QQ;
        a^[5*S+L,9*L+3]:=-QQ*QQ;
        a^[5*S+L,9*(L-1)+4]:=QQ*QQ*QQ;
        a^[5*S+L,9*L+4]:=-QQ*QQ*QQ;
        a^[5*S+L,9*(L-1)+5]:=QQ*QQ*QQ*QQ;
        a^[5*S+L,9*L+5]:=-QQ*QQ*QQ*QQ;
      end;
    {****** 6*S-1 SOR KESZ ******}
    FOR L:=1 TO S-1 do
      begin
        QQ:=x[M[L]];
        a^[6*S+L-1,9*(L-1)+2]:=1;
        a^[6*S+L-1,9*L+2]:=-1;
        a^[6*S+L-1,9*(L-1)+3]:=2*QQ;
        a^[6*S+L-1,9*L+3]:=-2*QQ;
        a^[6*S+L-1,9*(L-1)+4]:=3*QQ*QQ;
        a^[6*S+L-1,9*L+4]:=-3*QQ*QQ;
        a^[6*S+L-1,9*(L-1)+5]:=4*QQ*QQ*QQ;
        a^[6*S+L-1,9*L+5]:=-4*QQ*QQ*QQ;
      end;
    {******* 7*S-2 SOR KESZ *******}
    FOR L:=1 TO S-1 do
      begin
        a^[7*S+L-2,9*(L-1)+3]:=1;
        a^[7*S+L-2,9*L+3]:=-1;
        a^[7*S+L-2,9*(L-1)+4]:=3*x[M[L]];
        a^[7*S+L-2,9*L+4]:=-3*x[M[L]];
        a^[7*S+L-2,9*(L-1)+5]:=6*x[M[L]]*x[M[L]];
        a^[7*S+L-2,9*L+5]:=-6*x[M[L]]*x[M[L]];
      end;
    {***** 8*S-3 SOR KESZ ******}
    FOR L:=1 TO S-1  do
      begin
        a^[8*S+L-3,9*(L-1)+4]:=1;
        a^[8*S+L-3,9*L+4]:=-1;
        a^[8*S+L-3,9*(L-1)+5]:=4*x[M[L]];
        a^[8*S+L-3,9*L+5]:=-4*x[M[L]];
      end;
  end;





procedure splineconstant(var x,y:x_vector;m:intvector;s:longint;g:double) ;

 begin
  new(a);
  splinematrix(x,y,m,s,g) ;
  folytatas(x,y,m,s,1,g) ;
  egyenlet(a,csp,9*s-4);
  {egyenlet(a,csp,7*s-3);}
  dispose(a);
 end;


function sf(xr:double;var x:x_vector;m:intvector):double;

   var i:longint;

   begin
    i:=-1;
    if xr<=x[m[0]] then i:=0 else
    repeat
     inc(i);
    until ((x[m[i+1]]>xr)and(x[m[i]]<=xr))or(m[i]>m[i+1]);
    if m[i]>m[i+1] then dec(i);
    sf:=csp[9*i+1]+csp[9*i+2]*xr+csp[9*i+3]*xr*xr+csp[9*i+4]*xr*xr*xr+csp[9*i+5]*xr*xr*xr*xr;
    {sf:=csp[7*i+1]+csp[7*i+2]*xr+csp[7*i+3]*xr*xr+csp[7*i+4]*xr*xr*xr;}
   end;


function sdf(xr:double;var x:x_vector;m:intvector):double;

   var i:longint;

   begin
    i:=-1;
    if xr<=x[m[0]] then i:=0 else
    repeat
     inc(i);
    until ((x[m[i+1]]>=xr)and(x[m[i]]<=xr))or(m[i]>m[i+1]);
    if m[i]>m[i+1] then dec(i);
    sdf:=csp[9*i+2]+csp[9*i+3]*xr*2+csp[9*i+4]*xr*xr*3+csp[9*i+5]*xr*xr*xr*4;
    {sdf:=csp[7*i+2]+csp[7*i+3]*xr*2+csp[7*i+4]*xr*xr*3;}
   end;


procedure order(var x,y:x_vector;ndata:longint);

 var i,j,k,l:longint;
     x0,a:double;

 begin
  j:=ndata;
  for i:=1 to ndata do
    begin
     l:=0;
     a:=-1e+32;
     for k:=1 to j do
      if x[k]>a then
        begin
         a:=x[k];
         l:=k;
        end;
     if l<>0 then
      begin
       x0:=x[l];
       x[l]:=x[j];
       x[j]:=x0;
       x0:=y[l];
       y[l]:=y[j];
       y[j]:=x0;
      end;
     dec(j);
    end;
 end;

{****************************************************}




function linear(x:double;a:paramtype):double;

   begin
     linear:=x;
   end;

Begin
  message:='';
  mfit:=10;
  linearcombination:=true;
  axisconvx:=@linear;
  axisconvy:=@linear;
  for i:=1 to mma do mask[i]:=false;
  fepsilon:=1e-6;
End.