PROGRAM VYPFIT1;
USES GAMMA;

{ BEGINENVIRON }
CONST
   ndatap =  50;
   nbodobr = 100;
TYPE
   RealArrayNDATA = ARRAY [1..ndatap] OF real;
   RealArrayNOBR = ARRAY [1..nbodobr] OF real;
VAR
   Vstfile, Vystfile : text;
   tex1,tex2:string[12];
   x, y, sig  : RealArrayNDATA;
   mwt, i1, i2 , ndata: integer;
   a,b,siga,sigb,chi2,q,xmin,xmax,dx : real;
   xbod, ybod : real;
 { ENDENVIRON}

PROCEDURE fit(VAR x,y: RealArrayNDATA;
                ndata: integer;
              VAR sig: RealArrayNDATA;
                  mwt: integer;
         VAR a,b,siga: real;
      VAR sigb,chi2,q: real);
VAR
   i: integer;
   wt,t,sy,sxoss,sx,st2,ss,sigdat: real;
BEGIN
   sx := 0.0;
   sy := 0.0;
   st2 := 0.0;
   b := 0.0;
   IF mwt <> 0 THEN BEGIN
      ss := 0.0;
      FOR i := 1 TO ndata DO BEGIN
         wt := 1.0/sqr(sig[i]);
         ss := ss+wt;
         sx := sx+x[i]*wt;
         sy := sy+y[i]*wt
      END
   END
   ELSE BEGIN
      FOR i := 1 TO ndata DO BEGIN
         sx := sx+x[i];
         sy := sy+y[i]
      END;
      ss := ndata
   END;
   sxoss := sx/ss;
   IF mwt <> 0 THEN BEGIN
      FOR i := 1 TO ndata DO BEGIN
         t := (x[i]-sxoss)/sig[i];
         st2 := st2+t*t;
         b := b+t*y[i]/sig[i]
      END
   END
   ELSE BEGIN
      FOR i := 1 TO ndata DO BEGIN
         t := x[i]-sxoss;
         st2 := st2+t*t;
         b := b+t*y[i]
      END
   END;
   b := b/st2;
   a := (sy-sx*b)/ss;
   siga := sqrt((1.0+sx*sx/(ss*st2))/ss);
   sigb := sqrt(1.0/st2);
   chi2 := 0.0;
   IF mwt = 0 THEN BEGIN
      FOR i := 1 TO ndata DO
         chi2 := chi2+sqr(y[i]-a-b*x[i]);
      q := 1.0;
      sigdat := sqrt(chi2/(ndata-2));
{vsuvka Limp}
      Writeln('Sigdat =  ',sigdat);
{konec vsuvky}
      siga := siga*sigdat;
      sigb := sigb*sigdat
   END
   ELSE BEGIN
      FOR i := 1 TO ndata DO
         chi2 := chi2+sqr((y[i]-a-b*x[i])/sig[i]);
      q := gammq(0.5*(ndata-2),0.5*chi2)
   END;
END;
{endfit}
VAR
   sigm,sigm2:real;
BEGIN
   writeln('Zadej nazev vstupni file');
   readln(tex1);
   ASSIGN (Vstfile,tex1);
   writeln('Zadej nazev vystupni file');
   readln(tex2);
   ASSIGN (Vystfile,tex2);
   RESET(Vstfile);
   REWRITE(Vystfile);
   READLN(Vstfile,ndata,mwt);
   i1 := 0;
   WHILE (i1 < ndata) AND (NOT EOF(Vstfile)) DO
   BEGIN
      i1 := i1 + 1;
      IF mwt = 0
         THEN READLN (Vstfile, x[i1],y[i1])
         ELSE READLN (Vstfile, x[i1],y[i1],sig[i1]);
      IF mwt = 0
         THEN WRITELN (i1,x[i1],y[i1])
         ELSE WRITELN (i1,x[i1],y[i1],sig[i1]);
   END;
   IF i1=0 THEN Writeln('nic jsem neprecetl')
           ELSE IF i1 < ndata
                THEN BEGIN
                     Writeln('nacteno',i1,' udaju < ndat =',ndata);
                     ndata := i1
                     END
                 ELSE Writeln('nacteno  ',ndata,' udaju');
                 ;
   xmin := x[1];
   xmax := y[1];
   FOR i1 := 2 to ndata DO
       BEGIN
       IF x[i1] < xmin THEN xmin:=x[i1];
       if x[i1] > xmax THEN xmax:=x[i1];
       END;
   WRITELN ('xmin,xmax ',xmin,xmax);
   fit(x,y,ndata,sig,mwt,a,b,siga,sigb,chi2,q);
   WRITELN('a=',a,' b=',b);
   WRITELN('chi2=',chi2,' q=',q);
   WRITELN('siga=',siga,' sigb=',sigb);
   dx := (xmax - xmin)/(nbodobr - 1);
   xbod := xmin - dx;
   IF (mwt=0) THEN BEGIN
              sigm2 := chi2/(ndata-2);
              sigm  := sqrt(sigm2);
              WRITELN('sigma mereni (pri stejnych abs. chybach) = ',sigm)
              END;
   FOR i1 := 1 TO nbodobr DO
       BEGIN
         xbod := xbod + dx;
         ybod := a + b * xbod;
         IF i1 <= ndata
            THEN WRITELN(Vystfile,xbod,' ',ybod,' ',x[i1],' ',y[i1])
            ELSE WRITELN(Vystfile,xbod,' ',ybod)
       END;
   CLOSE(Vystfile)
END.