{ BEGINENVIRON}
PROGRAM L2main;

USES Gamma;

CONST
   ndatap = 50;
   map =   20;
TYPE
   RealArrayNDATA = ARRAY [1..ndatap] OF real;
   RealArrayNDAT2 = ARRAy [1..2,1..ndatap] of REAL;
   RealArrayMA = ARRAY [1..map] OF real;
   IntegerArrayMFIT = ARRAY [1..map] OF integer;
   RealArrayMAbyMA = ARRAY [1..map,1..map] OF real;
   RealArrayNPbyMP = ARRAY [1..map,1..1] OF real;
   IntegerArrayNP = ARRAY [1..map] OF integer;

PROCEDURE gaussj(VAR a: RealArrayMAbyMA;
                     n: integer;
                 VAR b: RealArrayNPbyMP;
                     m: integer);
VAR
   big,dum,pivinv: real;
   i,icol,irow,j,k,l,ll: integer;
   indxc,indxr,ipiv: ^IntegerArrayNP;
BEGIN
   new(indxc);
   new(indxr);
   new(ipiv);
   FOR j := 1 TO n DO ipiv^[j] := 0;
   FOR i := 1 TO n DO BEGIN
      big := 0.0;
      FOR j := 1 TO n DO
         IF ipiv^[j] <> 1 THEN
            FOR k := 1 TO n DO
               IF ipiv^[k] = 0 THEN
                  IF abs(a[j,k]) >= big THEN BEGIN
                     big := abs(a[j,k]);
                     irow := j;
                     icol := k
                  END
               ELSE IF ipiv^[k] > 1 THEN BEGIN
                  writeln('pause 1 in GAUSSJ - singular matrix');
                  readln
               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
         writeln('pause 2 in GAUSSJ - singular matrix');
         readln
      END;
      pivinv := 1.0/a[icol,icol];
      a[icol,icol] := 1.0;
      FOR l := 1 TO n DO
         a[icol,l] := a[icol,l]*pivinv;
      FOR l := 1 TO m DO
         b[icol,l] := b[icol,l]*pivinv;
      FOR ll := 1 TO n DO
         IF ll <> icol THEN BEGIN
            dum := a[ll,icol];
            a[ll,icol] := 0.0;
            FOR l := 1 TO n DO
               a[ll,l] := a[ll,l]-a[icol,l]*dum;
            FOR l := 1 TO m DO
               b[ll,l] := b[ll,l]-b[icol,l]*dum
         END
   END;
   FOR l := n DOWNTO 1 DO
      IF indxr^[l] <> indxc^[l] THEN
         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;
   dispose(ipiv);
   dispose(indxr);
   dispose(indxc)
END;
{END GAUSSJ}
PROCEDURE covsrt(VAR covar: RealArrayMAbyMA;
                        ma: integer;
                 VAR lista: IntegerArrayMFIT;
                      mfit: integer);
VAR
   j,i: integer;
   swap: real;
BEGIN
   FOR j := 1 TO ma-1 DO
      FOR i := j+1 TO ma DO covar[i,j] := 0.0;
   FOR i := 1 TO mfit-1 DO BEGIN
      FOR j := i+1 TO mfit DO
         IF lista[j] > lista[i] THEN
            covar[lista[j],lista[i]] := covar[i,j]
         ELSE
            covar[lista[i],lista[j]] := covar[i,j]
   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
      covar[lista[j],lista[j]] := covar[1,j];
   FOR j := 2 TO ma DO
      FOR i := 1 TO j-1 DO
         covar[i,j] := covar[j,i]
END;
{END COVSRT}

PROCEDURE funcs(x1,x2:Real;VAR af:RealArrayMA;mf:integer);
BEGIN
  If (mf<1) OR (mf>4) THEN  BEGIN
     Writeln(' Pocet funkci musi byt 1-4 PAUSE in funcs');readln;end;
  af[1] := 1;
  af[2] := x1;
  af[3] := x2;
  af[4] := x1 * x2;
END;

PROCEDURE l2fit(VAR x:RealArrayNDAT2; y,sig: RealArrayNDATA;
                     ndata: integer;
                     VAR a: RealArrayMA;
                        ma: integer;
                 VAR lista: IntegerArrayMFIT;
                      mfit: integer;
                 VAR covar: RealArrayMAbyMA;
                 VAR chisq: real);
VAR
   k,kk,j,ihit,i: integer;
   ym,wt,sum,sig2i: real;
   beta: ^RealArrayNPbyMP;
   afunc: ^RealArrayMA;
BEGIN
   new(beta);
   new(afunc);
   kk := mfit+1;
{   writeln('zacatek l2fit');}
   FOR j := 1 TO ma DO BEGIN
      ihit := 0;
      FOR k := 1 TO mfit DO
         IF lista[k] = j THEN ihit := ihit+1;
      IF ihit = 0 THEN BEGIN
         lista[kk] := j;
         kk := kk+1
      END
      ELSE IF ihit > 1 THEN BEGIN
         writeln('pause in routine LFIT');
         writeln('improper permutation in LISTA');
         readln
      END
   END;
   IF kk <> ma+1 THEN BEGIN
      writeln('pause in routine LFIT');
      writeln('improper permutation in LISTA');
      readln
   END;
   FOR j := 1 TO mfit DO BEGIN
      FOR k := 1 TO mfit DO covar[j,k] := 0.0;
      beta^[j,1] := 0.0
   END;
   FOR i := 1 TO ndata DO BEGIN
      funcs(x[1,i],x[2,i],afunc^,ma);
      ym := y[i];
      IF mfit < ma THEN
         FOR j := mfit+1 TO ma DO
            ym := ym-a[lista[j]]*afunc^[lista[j]];
      sig2i := 1.0/sqr(sig[i]);
      FOR j := 1 TO mfit DO BEGIN
         wt := afunc^[lista[j]]*sig2i;
        FOR k := 1 TO j DO
            covar[j,k] := covar[j,k]+wt*afunc^[lista[k]];
         beta^[j,1] := beta^[j,1]+ym*wt
      END
   END;
   IF mfit > 1 THEN BEGIN
      FOR j := 2 TO mfit DO
         FOR k := 1 TO j-1 DO covar[k,j] := covar[j,k]
   END;
   gaussj(covar,mfit,beta^,1);
   FOR j := 1 TO mfit DO
      a[lista[j]] := beta^[j,1];
   chisq := 0.0;
   FOR i := 1 TO ndata DO BEGIN
      funcs(x[1,i],x[2,i],afunc^,ma);
      sum := 0.0;
      FOR j := 1 TO ma DO
         sum := sum+a[j]*afunc^[j];
      chisq := chisq+sqr((y[i]-sum)/sig[i])
   END;
   covsrt(covar,ma,lista,mfit);
   dispose(afunc);
   dispose(beta)
END;
{END LFIT}

VAR
   xd,yd,zd,sigz,chisq,q,nu: Real;
   nbod,nx,ny,i: Integer;
   f1,f2 : text;
   xpole: RealArrayNDAT2;
   zpole,sigpole: RealArrayNDATA;
   a,sa: RealArrayMA;
   ma: integer;
   lista: IntegerArrayMFIT;
   mfit: integer;
   covar: RealArrayMAbyMA;

BEGIN
   ASSIGN(f1,'l2fit.dat');
   REset(f1);
{   ASSIGN(f2,'zkus2.vys');
   REwrite(f2);}
   readln(f1,nbod);
   for i:=1 to nbod do begin
       read(f1,xpole[1,i],xpole[2,i],zpole[i]);
       sigpole[i] := 1; END;
   ma := 4;
   writeln(' zadej mfit  4 (3 - bez x*y) >');
   readln(mfit);
   {mfit := 4;}
   for i := 1 to 4 do lista[i] := i;
   l2fit (xpole,zpole,sigpole,nbod,a,ma,lista,mfit,covar,chisq);
   writeln('   a        b        c      d');
   writeln(a[1]:8:3,a[2]:8:3,a[3]:8:3,a[4]:8:3);
   FOR I := 1 TO 4 DO sa[i] := sqrt(covar[i,i]); 
   writeln('  sa       sb       sc     sd');
   writeln(sa[1]:8:3,sa[2]:8:3,sa[3]:8:3,sa[4]:8:3);
   nu := nbod - mfit;
   q := gammq(nu/2,chisq/2);
   writeln('chisq =',chisq,'  quality',q);
   writeln('i             covar[i,j]');
   for i := 1 to 4 do writeln(i:4,covar[i,1]:8:3,
              covar[i,2]:8:3,covar[i,3]:8:3,covar[i,4]:8:3);
END.