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 }

ClrScr;

get_data(x,y,n);

nlin(x,y,y_calc,n);

write_data

end.