program cfit4; {164}

{ plot service included }
{ Pascal program to perform a linear least-squares fit }

const max = 20;

type index = 1..max;
ary = array[index] of real;

var x,y,y_calc : ary;
n : integer;
first,done : boolean;
a,b,correl_coef,
sigma_a,sigma_b,
see,seed : real;

external procedure cls;

function random(dummy: integer): real;
{ random number 0-1 }
{ define seed=4.0 as global }

const pi = 3.14159;

var x : real;
i : integer;

begin { RANDOM }
x:=seed+pi;
x:=exp(5.0*ln(x));
seed:=x-trunc(x);
random:=seed
end; { RANDOM }

procedure get_data(var x,y: ary;
var n: integer);
{ get values for n and arrays x,y }
{ y is randomly scattered about a straight line }

const a = 2.0;
b = 5.0;

var i,j : integer;
fudge : real;

begin
write('Fudge? ');
if fudge<0.0 then done:=true
else
begin
repeat
write('How many points? ');
until (n>2) and (n<=max);
if first then first:=false else cls;
for i:=1 to n do
begin
j:=n+1-i;
x[i]:=j;
y[i]:=(a+b*j)*(1.0+(2.0*random(0)-1.0)*fudge)
end { for-loop }
end { if }
end; { procedure get_data }

procedure write_data;
{ print out the answers }
var i,j : integer;

begin
writeln;
writeln(' I X Y YCALC');
for i:=1 to n do
writeln(i:3,x[i]:8:1,y[i]:9:2,y_calc[i]:9:2);
writeln;
writeln('Intercept is ',a:8:3,', sigma is ',sigma_a:8:3);
writeln(' Slope is ',b:8:2,', sigma is ',sigma_b:8:3);
writeln;
writeln('Correlation coeffizient is ',correl_coef:7:4);
for j:=1 to 40 do for i:=1 to 10000 do; cls
end; { write_data }

{\$I C:LINFIT2.LIB }

{\$I PLOT.LIB }

begin { MAIN program }
seed:=4.0;
cls;
first:=true;
done:=false;
repeat
get_data(x,y,n);
if not done then
begin
linfit(x,y,y_calc,a,b,n);
write_data;
plot(x,y,y_calc,n);
end
until done
end.
