home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / TURBOPAS / PAS_ENG.ZIP / DIFFUS.PAS < prev    next >
Encoding:
Pascal/Delphi Source File  |  1985-07-18  |  3.5 KB  |  141 lines
  1. { die ln-funktion bei MT+ ist fuer Werte kleiner 1E-5 aeussertst langsam,
  2.   so dass man das Gefuehl bekommt, das dass Programm sich aufhaengt. In-
  3.   folgedessen konnte ich dieses Programm nicht verifizieren. -Juergen }
  4.  
  5. program diffus; { --> 302 }
  6.  
  7. { Pascal program to perform a linear least-squares fit }
  8. { for the diffusion of Zn and Cu }
  9. { with Gauss-Jordan routine }
  10. { Sperate modules needed:
  11.    GAUSSJ,
  12.    PLOT }
  13.  
  14.  
  15. const maxr = 20;  { data prints }
  16.  maxc = 4;  { polynomial terms }
  17.  r = 1.987; { gas constant }
  18.  
  19. type ary = array[1..maxr] of real;
  20.  arys = array[1..maxc] of real;
  21.  ary2 = array[1..maxr,1..maxc] of real;
  22.  ary2s = array[1..maxc,1..maxc] of real;
  23.  
  24. var
  25.  x,y,y_calc : ary;
  26.  t,d,resid : ary;
  27.  coef,sig : arys;
  28.  nrow,ncol : integer;
  29.  correl_coef,srs : real;
  30.  
  31.  
  32.  
  33. procedure get_data(var x,y,t,d: ary;
  34.      var nrow: integer);
  35. { get values for nrow and arrays t,d }
  36.  
  37. var i : integer;
  38. begin
  39.   nrow:=7;
  40.   t[1]:=600.0; d[1]:=1.4E-12;
  41.   t[2]:=650.0; d[2]:=5.5E-12;
  42.   t[3]:=700.0; d[3]:=1.8E-11;
  43.   t[4]:=750.0; d[4]:=6.1E-11;
  44.   t[5]:=800.0; d[5]:=1.6E-10;
  45.   t[6]:=850.0; d[6]:=4.4E-10;
  46.   t[7]:=900.0; d[7]:=1.2E-9;
  47.   for i:=1 to nrow do
  48.     begin
  49.       x[i]:=1.0/(t[i]+273.0);
  50.       y[i]:=ln(d[i])
  51.     end
  52. end;  { procedure get data }
  53.  
  54.  
  55. procedure write_data;
  56. { print out the answers }
  57. var i : integer;
  58. begin
  59.   writeln;
  60.   writeln;
  61.   writeln('  I      TC        D    DCALC');
  62.   for i:=1 to nrow do
  63.     writeln(i:3,t[i]:8:0,d[i],' ',y_calc[i]);
  64.   writeln; writeln(' Coefficients ');
  65.   writeln(coef[1],' constant term');
  66.   for i:=2 to ncol do
  67.     writeln(coef[i]);  { other terms }
  68.   writeln;
  69.   writeln('D0=',(exp(coef[1])):7:2,' cm sq/sec.');
  70.   writeln('Q =',(-r*coef[2]/1000.0):8:2,'kcal/mole');
  71.   writeln;writeln('SRS= ',srs:7:3)
  72. end;  { write_data }
  73.  
  74. {procedure square(x: ary2;
  75.    y: ary;
  76.       var a: ary2s;
  77.       var g: arys;
  78.   nrow,ncol: integer);}
  79. { matrix multiplication routine }
  80. { a= transpose x times x }
  81. { g= y times x }
  82.  
  83. {$I SQUARE.LIB }
  84.  
  85. {external procedure gaussj(var b: ary2s;
  86.          y: arys;
  87.      var coef: arys;
  88.      ncol:  integer;
  89.      var error: boolean);
  90. }
  91. {$I GAUSSJ.LIB }
  92.  
  93. procedure linfit(x,  { independant variable }
  94.    y: ary; { dependent variable }
  95.    var y_calc: ary; { calculated dep. variable }
  96.    var resid:  ary; { array of residuals }
  97.    var coef:   arys; { coefficients }
  98.    var sig:    arys; { error on coefficients }
  99.    nrow:       integer; { length of array }
  100.    var ncol:   integer); { number of terms }
  101.  
  102. { least squares fit to nrow sets of x and y pairs of points }
  103. { Seperate procedures needed:
  104.  SQUARE -> form square coefficient matrix
  105.  GAUSSJ -> Gauss-Jordan elimination }
  106.  
  107. var xmatr  : ary2;  { data matrix }
  108.  a  : ary2s; { coefficient matrix }
  109.  g  : arys;  { constant vector }
  110.  error  : boolean;
  111.  i,j,nm  : integer;
  112.  see,a1  : real;
  113.  
  114. begin  { procedure linfit }
  115.   ncol:=2; { number of terms }
  116.   for i:=1 to nrow do
  117.     begin  { setup matrix }
  118.       xmatr[i,1]:=1.0; { first column }
  119.       xmatr[i,2]:=x[i] { second column }
  120.     end;
  121.   square(xmatr,y,a,g,nrow,ncol);
  122.   gaussj(a,g,coef,ncol,error);
  123.   srs:=0.0;
  124.   a1:=exp(coef[1]);
  125.   for i:=1 to nrow do
  126.     begin
  127.       y_calc[i]:=a1*exp(coef[2]*x[i]);
  128.       if y[i]<>0.0 then resid[i]:=y_calc[i]/y[i]-1.0
  129.       else resid[i]:=y[i]/y_calc[i]-1.0;
  130.       srs:=srs+sqr(resid[i])
  131.     end
  132. end; { linfit }
  133.  
  134.  
  135. begin  { main program }
  136.   clrscr;
  137.   get_data(x,y,t,d,nrow);
  138.   linfit(x,y,y_calc,resid,coef,sig,nrow,ncol);
  139.   write_data
  140. end.
  141.