PROGRAM d14r11(input,output);
(* driver for routine MEDFIT *)
CONST
   npt=100;
   spread=0.1;
   ndata=npt;
TYPE
   glndata = ARRAY [1..ndata] OF real;
   glsarray = glndata;
VAR
   glinext,glinextp : integer;
   glma : ARRAY [1..55] OF real;
   gliset : integer;
   glgset : real;
   a,aa,abdev,abdevt,b,chi2,q,siga,sigb : real;
   i,idum,mwt : integer;
   sig,x,y,arr : glndata;

(*$I MODFILE.PAS *)
(*$I RAN3.PAS *)

(*$I GASDEV.PAS *)

(*$I GAMMLN.PAS *)

(*$I GSER.PAS *)

(*$I GCF.PAS *)

(*$I GAMMQ.PAS *)

(*$I SORT.PAS *)

(*$I FIT.PAS *)

(*$I ROFUNC.PAS *)

(*$I MEDFIT.PAS *)

BEGIN
   gliset := 0;
   idum := -1984;
   FOR i := 1 to npt DO BEGIN
      x[i] := 0.1*i;
      y[i] := -2.0*x[i]+1.0+spread*gasdev(idum);
      sig[i] := spread
   END;
   mwt := 1;
   writeln;
   fit(x,y,npt,sig,mwt,a,b,siga,sigb,chi2,q);
   writeln('According to routine FIT the result is:');
   writeln('a = ':7,a:8:4,'  uncertainty: ':13,siga:8:4);
   writeln('b = ':7,b:8:4,'  uncertainty: ':13,sigb:8:4);
   writeln('chi-squared: ':16,chi2:8:4,' for ',npt:4,' points');
   writeln('goodness-of-fit: ':20,q:8:4);
   writeln;
   writeln('According to routine MEDFIT the result is:');
   medfit(x,y,npt,a,b,abdev);
   writeln('a = ':7,a:8:4);
   writeln('b = ':7,b:8:4);
   writeln(' ':3,'absolute deviation (per data point): ',abdev:8:4);
   writeln(' ':3,'note: gaussian spread is',spread:8:4,')')
END.