{ die ln-funktion bei MT+ ist fuer Werte kleiner 1E-5 aeussertst langsam,
so dass man das Gefuehl bekommt, das dass Programm sich aufhaengt. In-
folgedessen konnte ich dieses Programm nicht verifizieren. -Juergen }

program nlin3; { -> 310 }
{ Pascal program to perform a nonlinear least-squares fit for the diffusion
of Zn in CU }

const maxr = 20; { data prints }
maxc = 4; { polynomial terms }
r = 1.987; { gas constant }
type
index = 1..maxr;
ary = array[index] of real;
arys = array[1..maxc] of real;
ary2 = array[1..maxr,1..maxc] of real;

var
x,y,y_calc : ary;
t,d,ex : ary;
coef : arys;
i,n : integer;
nrow,ncol : integer;
done,error : boolean;
correl_coef,srs,
a,b,x2 : real;

external procedure cls;

procedure get_data(var x,y: ary;
var n: integer);
{ get values for n and arrays t,d }

var i : integer;

begin
n:=7;
t[1]:=600.0; d[1]:=1.4E-12;
t[2]:=650.0; d[2]:=5.5E-12;
t[3]:=700.0; d[3]:=1.8E-11;
t[4]:=750.0; d[4]:=6.1E-11;
t[5]:=800.0; d[5]:=1.6E-10;
t[6]:=850.0; d[6]:=4.4E-10;
t[7]:=900.0; d[7]:=1.2E-9;
for i:=1 to n do
begin
x[i]:=1.0/(t[i]+273.0);
y[i]:=d[i]
end
end; { proceddure get data }

procedure write_data;
{ print out the answers }
var i : integer;
begin
writeln;
writeln;
writeln(' I TC D DCALC');
for i:=1 to n do
writeln(i:3,t[i]:8:0,d[i],' ',y_calc[i]);
writeln; writeln(' Coefficients ');
writeln(coef[1],' constant term');
for i:=2 to ncol do
writeln(coef[i]); { other terms }
writeln;
writeln('D0=',a:7:2,' cm sq/sec.');
writeln('Q =',(-r*b/1000.0):8:2,'kcal/mole');
writeln;writeln('SRS= ',srs:8:4)
end; { write_data }

procedure func(b: real;
var fb,dfb: real);
var i : integer;
s1,s2,s3,s4,s5,s6,
ex1,ex2,xi,
x2,yi,y2 : real;
begin
s1:=0.0;
s2:=0.0;
s3:=0.0;
s4:=0.0;
s5:=0.0;
s6:=0.0;
for i:=1 to n do
begin
xi:=x[i];
x2:=xi*xi;
yi:=y[i];
y2:=yi*yi;
ex1:=exp(b*xi);
ex[i]:=ex1;
ex2:=ex1*ex1;
s1:=s1+xi*ex2/y2;
s2:=s2+ex1/yi;
s3:=s3+xi*ex1/yi;
s4:=s4+ex2/y2;
s5:=s5+2.0*x2*ex2/y2;
s6:=s6+x2*ex1/yi
end;
fb:=s1*s2-s3*s4;
dfb:=s2*s5-s1*s3-s4*s6;
a:=s2/s4
end; { func }

procedure newton(var x: real);
const tol = 1.0E-6;
max = 20;
var fx,dfx,dx,x1 : real;
i : integer;

begin { newton }
error:=false;
i:=0;
repeat
i:=i+1;
x1:=x;
func(x,fx,dfx);
if dfx=0.0 then
begin
error:=true;
x:=1.0;
writeln('ERROR: slope zero')
end
else
begin
dx:=fx/dfx;
x:=x1-dx;
end
until
error or
(i>max) or
(abs(dx)<=abs(tol*x));
if i>max then
begin
writeln(chr(7),'ERROR: no convergence in ',max,' loops');
error:=true
end
end; { newton }

procedure nlin(x,y: ary;
var y_calc: ary;
n: integer);
{ fits the diffusion equation through n sets of x and y pairs of points }
var
resid : ary;
matr : ary2;
i : integer;

xi,yi,sum_x,
sum_y,sum_y2,b1,
sum_xy,sum_x2 : real;
begin { nlin }
ncol:=2; { number of terms }
sum_x:=0.0;
sum_y:=0.0;
sum_xy:=0.0;
sum_x2:=0.0;
for i:=1 to n do
begin
xi:=x[i];
yi:=ln(y[i]);
sum_x:=sum_x+xi;
sum_y:=sum_y+yi;
sum_y2:=sum_y2+yi*yi;
sum_xy:=sum_xy+xi*yi;
sum_x2:=sum_x2+xi*xi
end;
b:=(sum_xy-sum_x*sum_y/n)/(sum_x2-sqr(sum_x)/n);
newton(b);
coef[1]:=a;
coef[2]:=b;
srs:=0.0;
for i:=1 to n do
begin
y_calc[i]:=a*ex[i];
if y[i]<>0.0 then
resid[i]:=y_calc[i]/y[i]-1.0
else resid[i]:=y[i]/y_calc[i]-1.0;
srs:=srs+sqr(resid[i])
end
end; { nlin }

begin { main program }
cls;
writeln(' start get_data ');
get_data(x,y,n);
writeln(' start nlin ');
nlin(x,y,y_calc,n);
writeln(' start write_data ');
write_data
end.
