home *** CD-ROM | disk | FTP | other *** search
/ Shareware Overload / ShartewareOverload.cdr / games / leastq.zip / LEASTQ.PAS < prev    next >
Pascal/Delphi Source File  |  1990-02-17  |  42KB  |  1,197 lines
  1. Program LEASTQ;   { This program is for TURBO PASCAL 5.0 }
  2. {Copywrite (c) July 1989
  3.      John D. Fox
  4.      309 N. W. 24th St.
  5.      Gainesville, Florida 32607 }
  6. {$N+,E+}                     { simulate 80x87 if not available }
  7. {$D+,L+,S+,V+}            { turn these on when changing code }
  8. uses Crt,Dos,Graph,Printer,Drivers;
  9. type float = double;
  10. const bombthreshold = 1.0E-16;
  11. const epsilon = 1.0E-10;
  12. const copywt = 'Copywrite (c) 07/24/89, John D. Fox';
  13. const
  14.       max = 11; maxpts = 1024;     { change these if you like }
  15.       size: array[1..19] of float = (0.01,0.02,0.05,0.1,0.2,0.5,1.0,2.0,5.0,
  16.             10.0,20.0,50.0,100.0,200.0,500.0,1000.0,2000.0,5000.0,10000.0);
  17. { This program uses some (shudder . . gasp . . gurgle . .) LABELS!!! }
  18. label datinput,start,xlimits,terms;
  19. type  matrx = array[0..max] of array[0..max] of float;
  20.       species =
  21.   (powers,sines,cosines,fourier,oddsines,cheby,legendre,bessels,expon,gauss);
  22.       vector = array[1..maxpts] of float;
  23.       coeff = array[0..max] of float;
  24.       cpoint = ^coeff;
  25.       textline = string[80];
  26.       phrase = string[14];
  27. var   H,HINV,TEMP           : matrx;
  28.       x,y,yc,err,sigsq,ysave: vector;
  29.       a,T,P,JJ,errc,pct     : coeff;
  30.       funct                 : species;
  31.       fudged                : array[1..maxpts] of boolean;
  32.       fpoint                : array[1..maxpts] of cpoint;
  33.       dsource,fname,inpt,
  34.       pswitch               : phrase;
  35.       filename              : array[1..20] of phrase;
  36.       data                  : text;
  37.       grDriver,grMode,
  38.       ErrorCode             : integer;    {initialize grapics mode}
  39.       dd                    : string[4];  {for goof proof data input}
  40.       xd,yd                 : word;       {number of decimal points}
  41.       caption               : textline;
  42.       chisq,                    {what we are trying to mininmize}
  43.       xmax,ymax,xmin,ymin,      {limits of plot}
  44.       xtick,                    {spacing of horiz tics}
  45.       yleast,                   {smallest data point}
  46.       asave1,asave2,asave0,     {restore a[i] for Gaussian fudge}
  47.       confid,                   {confidence level of fit }
  48.       len,cy,lmin,lmax,cent : float;
  49.       ix1,ix2,iy1,iy2,                   { coordinates on plot }
  50.       nticks,                            { for frame }
  51.       code,DOF,yscale       : integer;
  52.       bombed,logplt,plotted,
  53.       filefound             : boolean;
  54.       i,                                  {row}
  55.       j,                                  {column}
  56.       k,                                  {data point index}
  57.       Npts,Dim,N,nfiles     : word;
  58.       ch                    : char;
  59.  
  60. Function YN: boolean;
  61. begin
  62. repeat ch := ReadKey; ch := Upcase(ch) until ch in['Y','N'];
  63. if ch = 'Y' then YN := True else YN := False;
  64. end;
  65.  
  66. Procedure ABORT(Msg : string);
  67. begin
  68.    Writeln(msg, ':',GraphErrorMsg(ErrorCode));
  69.    Halt(1);
  70. end;
  71.  
  72. Procedure CWRITE(caption:textline);        { Highlight capital letters}
  73. var ch : char; k : integer;
  74. begin
  75. for k := 1 to length(caption) do begin
  76.    ch  := caption[k];
  77.    if ch in ['A'..'Z'] then begin
  78.       Textcolor(14); write(ch); Textcolor(7); end
  79.    else write(ch);
  80.    end;
  81. end;
  82.  
  83. Procedure CAPITALIZE(var upword : phrase);
  84. begin
  85. for j := 1 to length(upword) do upword[j] := Upcase(upword[j]);
  86. end;
  87.  
  88. Procedure DATAFILES;                    { find the available data files }
  89. const datafl = '*.DAT';                 {    and select one             }
  90. var
  91.   Filerecord : SearchRec;
  92.   i,j        : word;
  93. begin
  94. TextColor(7);
  95. TextBackground(0);
  96. ClrScr;
  97. writeln('Files on this directory with extension ".DAT" are:'); writeln;
  98. j := 0;
  99. FindFirst(datafl,AnyFile, Filerecord);
  100. if (DosError = 0) then begin
  101.    Inc(j); filename[j] := Filerecord.name end;
  102. while (DosError = 0) and (j < 20) do begin
  103.   FindNext(Filerecord);
  104.   if DosError = 0 then begin
  105.      Inc(j); filename[j] := Filerecord.name end;
  106.   end;
  107. nfiles := j;
  108. if nfiles > 0 then
  109.   for j := 1 to nfiles do writeln
  110.   ('                ',j:2,'    ',filename[j]);
  111. writeln('                ',nfiles+1:2,' . .  None of the above . . . .');
  112. writeln;
  113. repeat
  114.   write('Enter number of data file to use: ');
  115.   readln(dd);
  116.   Val(dd,i,code);
  117.   until (code = 0) and (i in [1..nfiles + 1]);
  118. if i = nfiles + 1 then dsource := 'NONAME' else dsource := filename[i];
  119. end;  { of DATAFILES }
  120.  
  121. Procedure SORT;                   { sort data points in order of ascending x }
  122. var xtemp,ytemp,errtemp : float;
  123. begin
  124. for k:= 2 to Npts do begin
  125.    xtemp := x[k]; ytemp := y[k]; errtemp := err[k];
  126.    j := k;
  127.    while (x[j-1] > xtemp) and (j > 1) do begin       { point out of place ?   }
  128.       x[j] := x[j-1];                                { then move the          }
  129.       y[j] := y[j-1];                                { other points up        }
  130.       err[j] := err[j-1];                            { like a deck of cards   }
  131.       j := j-1;                                      { to make a hole for it  }
  132.       end;
  133.    x[j] := xtemp; y[j] := ytemp; err[j] := errtemp;  { and insert point }
  134.    end;
  135. end;
  136.  
  137. Procedure LIMITS;                   { Finds values of xmax,ymax,xmin,ymin }
  138. var xmost,xleast,ymost,span : float;
  139.     ispan : word;
  140. begin
  141.               (****** find x limits ******)
  142. xleast := x[1]; xmost := x[Npts];          { data has been sorted on x }
  143. span := xmost - xleast;
  144. while span > 100.0 + epsilon do span := span/10.0;
  145. while span < (10.0 - epsilon) do span := 10.0*span;   { 10 < span < 100 }
  146. ispan := Trunc(span);
  147. if ispan mod 5 = 0 then begin                      {perfect fit to 5 tics?}
  148.     xmin := xleast; xmax := xmost end
  149. else begin                            { nope, then fit with standard sizes }
  150.    k := 0;
  151.    repeat inc(k) until size[k] > (xmost-xleast)/20.0 - epsilon;
  152.    xmin := 0;
  153.    if (xleast + epsilon) < 0 then
  154.         repeat xmin := xmin - size[k] until xmin < xleast - epsilon
  155.     else begin
  156.        repeat xmin := xmin + size[k] until xmin > xleast + epsilon;
  157.        xmin := xmin - size[k];
  158.        end;
  159.     xmax := xmin;
  160.     repeat xmax := xmax + size[k] until xmax > xmost - epsilon
  161.     end;
  162.                (******  find y limits *****)
  163.  
  164. ymost := -100000.0;
  165. yleast := 100000.0;          { y s not sorted, so find yleast,ymost }
  166. for k := 1 to Npts do begin
  167.    ysave[k] := y[k];
  168.    if y[k] > ymost then ymost := y[k];
  169.    if y[k] < yleast then yleast := y[k];
  170.    end;
  171. k := 0;
  172. repeat inc(k) until size[k] > (ymost - yleast)/10.0 - epsilon;
  173. ymin := 0;
  174. if (yleast + epsilon) < 0 then
  175.         repeat ymin := ymin - size[k] until ymin < yleast - epsilon
  176.    else begin
  177.    repeat ymin := ymin + size[k] until ymin > yleast + epsilon;
  178.    ymin := ymin - size[k]
  179.    end;
  180. ymax := ymin;
  181. repeat ymax := ymax + size[k] until ymax > (ymost + epsilon);
  182. end; { of LIMITS }
  183.  
  184. Procedure READFILE;  {Read in the selected file; must be on the same disk}
  185. begin
  186. filefound := True;
  187. Assign(data,dsource);
  188. {$I-} Reset(data); {$I+}
  189. if IOresult = 0 then begin
  190.    Npts := 0;
  191.    k := 0;
  192.    while not EOF(data) and (k < maxpts) do begin
  193.       Inc(k); readln(data,x[k],y[k],err[k]);
  194.       if err[k] = 0 then k := k-1;   {no points with zero error allowed}
  195.       sigsq[k] := err[k]*err[k];
  196.       end;
  197.    Npts := k;
  198.    SORT; LIMITS;
  199.    Close(data);
  200.    end else filefound := False;
  201. end;
  202.  
  203. Procedure WRITEFILE;               { Save data entered from keyboard }
  204. var DFILE : phrase;
  205. begin
  206. repeat
  207.    write('Name of File (no extension): '); readln(dsource);
  208. until length(dsource) in [1..8];
  209. CAPITALIZE(dsource);
  210. DFILE := dsource + '.DAT';
  211. Assign(data,DFILE);
  212. Rewrite(data);
  213. writeln('Writing data to disk as ',DFILE);
  214. for k := 1 to Npts do writeln(data,x[k]:12:xd,y[k]:12:yd,err[k]:12:yd);
  215. Close(data);
  216. end;
  217.  
  218. Procedure ERASE(nn:integer);                  { Blank out previous text}
  219. var xx,yy :integer;
  220. begin
  221. xx := GetX; yy := GetY;
  222. SetFillStyle(1,1);
  223. Bar(xx,yy,xx+nn,yy+8);
  224. MoveTo(xx,yy);
  225. end;
  226.  
  227. Procedure BLINK;                     { Graphics equiv of blinking cursor }
  228. var x,y : integer;
  229. begin
  230. x :=  GetX; y := GetY;
  231. SetFillStyle(1,1);
  232. repeat
  233.    OutText('_');
  234.    Delay(150);
  235.    MoveTo(x,y);
  236.    Bar(x,y,x+10,y+10);
  237.    Delay(150);
  238.    until KeyPressed;
  239. end;
  240.  
  241. Procedure REALIN(var xx:float);         { Graphics equiv of Readln }
  242. begin
  243. inpt := '';
  244. repeat
  245.     BLINK;
  246.     ch := ReadKey;
  247.     if ch in ['+'..'9'] then begin
  248.        inpt := inpt + ch;
  249.        OutText (ch);
  250.        end;
  251.     if ch = #8 then begin
  252.        MoveTo(Getx-8,Gety);
  253.        Delete(inpt,length(inpt),1);
  254.        ERASE(8)
  255.        end;
  256.     until ch = #13;
  257. if inpt <> '' then Val(inpt,xx,code);
  258. end;
  259.  
  260. Procedure WORDIN(var kk:word);          { Input word }
  261. begin
  262. inpt := '';
  263. repeat
  264.     BLINK;
  265.     ch := ReadKey;
  266.     if ch in ['-','0'..'9'] then begin
  267.        inpt := inpt + ch;
  268.        OutText (ch)
  269.        end;
  270.     if ch = #8 then begin
  271.        MoveTo(GetX-8,Gety);
  272.        Delete(inpt,length(inpt),1);
  273.        ERASE(8)
  274.        end;
  275.     until ch = #13;
  276. if inpt <> '' then Val(inpt,kk,code);
  277. end;
  278.  
  279. Procedure OUTWR(x:float);          {graphics equivalent of Write(x)}
  280. var msg : string[5];
  281. begin
  282. if x > 100.0 then Str(x:5:1,msg)              { e.g. 182.0 }
  283.    else if x > 10.0 then Str(x:5:2,msg)       { e.g.  18.2 }
  284.       else if x > 1.0 then Str(x:5:3,msg)     { e.g.  1.82 }
  285.          else Str(x:5:4,msg);                 { e.g. 0.182 }
  286. OutText(msg);
  287. end;
  288.  
  289. Procedure XTICS;           { calculate horizontal spacing for 4 or 5 tics }
  290. var xscale : float;
  291.     iscale : integer;
  292. begin                       { calculate iscale between 10 and 100 }
  293. xscale := xmax - xmin;
  294. while xscale > 100.001 do xscale := xscale/10.0;
  295. while xscale < 9.9999 do xscale := xscale*10.0;
  296. iscale := Trunc(xscale +0.001);
  297. if (iscale mod 10 = 0) and (iscale > 30) then nticks := iscale div 10 else
  298. if iscale mod 5 = 0 then nticks := 5
  299.    else if iscale mod 4 = 0 then nticks := 4
  300.       else if iscale mod 4 < iscale mod 5 then nticks := 4
  301.          else nticks := 5;
  302. xtick := 500.0*(iscale div nticks)/xscale;
  303. end;
  304.  
  305.  
  306. Procedure DECIMALS;     {Calculate number of decimals for output}
  307. var xxd,yyd : float;
  308. begin
  309. if abs(xmax) > abs(xmin) then xxd := abs(xmax) else xxd := abs(xmin);
  310. if xxd > 10000.0 then xd := 0
  311.    else if xxd > 100.0 then xd := 1
  312.       else if xxd > 10.0 then xd := 2
  313.          else xd := 3;
  314. if abs(ymax) > abs(ymin) then yyd := abs(ymax) else yyd := abs(ymin);
  315. if yyd > 10000.0 then yd := 0
  316.    else if yyd > 100.0 then yd := 1
  317.       else if yyd > 10.0 then yd := 2
  318.          else if yyd > 0.1 then yd := 3
  319.             else yd := 4;
  320. end;
  321.  
  322. Procedure LOGTICS;                          { make semi-log graph paper}
  323. var ii : integer;
  324.     logsize : float;
  325. begin
  326. for k := 1 to 19 do begin
  327.    logsize := ln(size[k]);
  328.    if (logsize > ymin) and (logsize < ymax) then begin
  329.       ii := 10 + Round(yscale*(ymax - logsize)/(ymax - ymin));
  330.       Line(100,ii,110,ii); Line(590,ii,600,ii);           { draw tics }
  331.       MoveTo(70,ii); OUTWR(size[k]);                      { label tics }
  332.       end;
  333.    end;
  334. MoveTo(70,10); OUTWR(exp(ymax));
  335. MoveTo(70,yscale + 10); OUTWR(exp(ymin));
  336. SetTextStyle(0,1,0);
  337. MoveTo(25,yscale div 2); OutText('L O G A R I T M I C   S C A L E');
  338. SetTextStyle(0,0,0);
  339. end;
  340.  
  341. Procedure DRAWFRAME;
  342. var vtick,itick : integer;
  343.     xlabel,ylabel : float;
  344. begin
  345. XTICS; DECIMALS;
  346. SetGraphMode(GetGraphMode);
  347. SetBkColor(1); SetColor(7);
  348. SetTextJustify(CenterText,CenterText);
  349. Rectangle(100,10,600,(yscale+10));                    {   draw frame     }
  350. for j := 0 to nticks do begin
  351.    itick := 100 + Round(j*xtick);                 {  draw horizontal tics }
  352.    Line(itick,yscale,itick,yscale+10);
  353.    Line(itick,10,itick,20);
  354.    xlabel := xmin +j*(xmax-xmin)*xtick/500.0;
  355.    Moveto(itick,yscale+18); OUTWR(xlabel);        { label horizontal tics }
  356.    end;
  357. if logplt then LOGTICS else begin
  358.    if ymax/(ymin+0.001) > 2.6 then ymin := 0.0;
  359.    for i := 0 to 4 do begin                       { draw vertical tics    }
  360.       vtick := 10 + ((4-i)*yscale div 4);
  361.       Line(100,vtick,110,vtick);
  362.       Line(590,vtick,600,vtick);
  363.       ylabel := ymin + i*(ymax-ymin)/4.0;
  364.       MoveTo(70,vtick); OUTWR(ylabel);            { label vertical tics   }
  365.       end;
  366.    end;
  367. MoveTo(350,5); SetColor(12); OutText(dsource);    {   lable picture }
  368. SetTextJustify(LeftText,TopText);
  369. end;  { of DRAWFRAME }
  370.  
  371. Procedure PLOTPOINT(pcolor:integer);         { data point is a little square }
  372. begin                                        {    with error flag  }
  373. SetColor(pcolor);
  374. ix1 := Round(500.0*(x[k] - xmin)/(xmax-xmin)) + 100; ix2 := ix1;
  375. iy1 := 10 + Round(yscale*(ymax - y[k] - err[k])/(ymax - ymin));
  376. iy2 := 10 + Round(yscale*(ymax - y[k] + err[k])/(ymax - ymin));
  377. Line(ix1,iy1,ix2,iy2);                                  {  Draw error flag }
  378. ix1 := ix1 - 1; iy1 := 1 + (iy1 + iy2) div 2;   { fatten flag in the middle }
  379. for i := 0 to 2 do Putpixel(ix1,iy1-i,pcolor);  {       to make point       }
  380. ix1 := ix1 + 2;
  381. for i := 0 to 2 do Putpixel(ix1,iy1-i,pcolor);
  382. end;
  383.  
  384. Procedure PLOTPOINTS;
  385. var pcolor : integer;
  386. begin
  387. for k := 1 to Npts do begin
  388.    if fudged[k] then pcolor := 10 else pcolor := 12;
  389.    PLOTPOINT(pcolor);
  390.    end;
  391. end;
  392.  
  393. Procedure NEATWR(xx:float);    {write floating point number in 7 spaces}
  394. begin
  395. if abs(xx) >= 1000.0 then write(xx:7:1)                { e. g. -1234.5 }
  396.   else if abs(xx) >= 100 then write(xx:7:2)            { e. g. -123.45 }
  397.      else if abs(xx) >= 10.0 then write(xx:7:3)        { e. g. -12.345 }
  398.        else if abs(xx) >= 1.0 then write(xx:7:4)       { e. g. -1.2345 }
  399.            else write(xx:7:4);                         { e. g. -0.1234 }
  400. end;
  401.  
  402. Procedure READIN(var xx:float);               { goof proof data input }
  403. var  x,y : integer;            {input terminated by <CR> or arrow keys}
  404. begin
  405. x := WhereX; y := WhereY;
  406. repeat
  407.   code := 0;
  408.   inpt := '';
  409.   repeat
  410.     ch := Readkey;
  411.     if (ch = #0) and Keypressed then ch := Readkey;
  412.     if ch in ['-'..'9'] then inpt := inpt + ch;
  413.     if ch = #8 then Delete(inpt,Length(inpt),1);
  414.     gotoXY(x,y); write(inpt)
  415.     until ch in [#13,#27,#72,#75,#77,#80];
  416.   if inpt <> '' then begin
  417.     Val(inpt,xx,code);
  418.     if code <> 0 then begin
  419.       gotoXY(x,y); write ('       '); write(#7); gotoXY(x,y) end;
  420.     end;
  421.   until code = 0;
  422. gotoXY(x,y); NEATWR(xx);
  423. end; { of READIN }
  424.  
  425. Procedure REPAIR;               { view data on screen and correct points }
  426. label Start; { <===== yet another loathesome LABEL! }
  427. var xnew,ynew,errnew : float;
  428. begin
  429. Setcolor(12);
  430. Start:
  431. repeat
  432.    MoveTo(150,yscale+26); ERASE(300);
  433.    OutText('Point to fix: ');
  434.    WORDIN(k);
  435.    if (k < 1) or (k > Npts) then begin   { not a possible point }
  436.       OutText(' No such point');
  437.       write(#7);
  438.       Delay(400);
  439.       end;
  440.    until k in [1..Npts];
  441. PLOTPOINT(15);
  442. MoveTo(300,yscale + 26); OutText('This it? ');
  443. repeat PLOTPOINT(1); Delay(150); PLOTPOINT(15); delay(150) until Keypressed;
  444. if not YN then begin PLOTPOINT(12); goto Start end;    { wrong point }
  445.       { point accepted for alteration; now adjust its values}
  446. PLOTPOINT(15);
  447. if plotted then fudged[k] := true;   { henceforth, will be plotted in green }
  448. xnew := x[k]; ynew := y[k]; errnew := err[k];
  449. MoveTo(150,yscale + 38); ERASE(300); OutText('x[k] now =');
  450. OUTWR(x[k]); OutText('  set to: ');
  451. REALIN(xnew);
  452. MoveTo(150,yscale + 38); ERASE(300); OutText('y[k] now =');
  453. if logplt then OUTWR(exp(y[k])) else OUTWR(y[k]);
  454. OutText('  set to: ');
  455. REALIN(ynew);
  456. MoveTo(150,yscale + 38); ERASE(300);
  457. OutText('err[k] now ='); OUTWR(err[k]); OutText(' set to: ');
  458. REALIN(errnew);
  459. PLOTPOINT(1);     {erase old point by plotting it in the background color }
  460. x[k] := xnew; ysave[k] := ynew; err[k] := errnew;
  461. if logplt then y[k] := ln(ynew) else y[k] := ynew;
  462. PLOTPOINT(10);
  463. SetColor(13);
  464. MoveTo(450,yscale + 38); OutText('More?'); BLINK;
  465. if YN then begin
  466.    MoveTo(150,yscale+38); ERASE(380); goto Start end;
  467. RestoreCrtMode;
  468. end;  { of REPAIR }
  469.  
  470. Procedure KEYBOARD;                           { read data from keyboard }
  471. var steps,count       : boolean;
  472.     x0,step,numbin    : float;
  473.     i,j,ix,iy,kk,top  : word;
  474. begin
  475. ClrScr;
  476. writeln('                        OK, Input from Keyboard');
  477. steps := False; count := False;
  478. write('Uniform steps in x ? ');
  479. if YN then begin
  480.    writeln('Yes');
  481.    x0 := 0; step := 1.0;                      { default values }
  482.    steps := True;
  483.    write('     First x = '); READIN(x0);
  484.    write('  step = '); READIN(step);
  485.    writeln('   Program will supply the x coodinates');
  486.    writeln
  487.    end
  488. else writeln('No');
  489. write('Is y a count? ');
  490. if YN then begin
  491.    write(' Yes ');
  492.    writeln(' Program will supply the errors');
  493.    count := True
  494.    end;
  495. k := 1;
  496. gotoXY(1,8);
  497. writeln('   Use cursor keys to correct mistakes, <Esc> to exit');
  498. Textbackground(4);
  499. gotoXY(4,9); write(' Point #      x[k]        y[k]       err[k]      ');
  500. Window(4,10,52,24);
  501. for i := 1 to 16 do begin gotoxY(1,i); ClrEol end;
  502. k := 1; Npts := 1; top := 1;
  503. if not steps then j := 1 else
  504.    begin x[k] := x0; gotoXY(12,1); NEATWR(x[k]); j := 2 end;
  505. repeat
  506.   ix := 12 + 12*(j-1); iy := k-top+1;
  507.   if k < Npts then begin
  508.      case j of
  509.      1 : numbin := x[k];
  510.      2 : numbin := y[k];
  511.      3 : numbin := err[k];
  512.      end;
  513.     gotoXY(ix,iy); Textbackground(0); NEATWR(numbin);
  514.     gotoXY(ix,iy); Textbackground(4);
  515.     end
  516.   else numbin := 0;
  517.   gotoXY(ix,iy);
  518.   READIN(numbin);
  519.   case j of
  520.     1: begin
  521.      x[k] := numbin;
  522.      case ch of
  523.         #13,#77 : inc(j);
  524.         #72 : if k > 1 then dec(k);
  525.         #80 : if k < Npts then inc(k);
  526.         #75 : WRITE(#7);
  527.         end;
  528.      end;
  529.     2 : begin
  530.         y[k] := numbin;
  531.         if count and (ch = #13) then begin
  532.            if (numbin > 1.0 - epsilon) then begin
  533.                err[k] := sqrt(numbin); gotoXY(36,iy); NEATWR(err[k]); end
  534.            else begin
  535.               ch := #20;
  536.               write(#7);
  537.               gotoXY(18,iy); write('*** NEG COUNT NOT ALLOWED ***');
  538.               Delay(1000); gotoXY(18,iy); ClrEOL;
  539.               end;
  540.            end;
  541.         case ch of
  542.         #13,#77 : if not count then inc(j) else begin
  543.            inc(k); if not steps then dec(j) end;
  544.         #80 : begin
  545.                 if not count then err[k] := err[k-1];
  546.                 gotoXY(36,iy); NEATWR(err[k]);
  547.                 inc(k);
  548.                 end;
  549.         #72 : if k > 1 then dec(k) else WRITE(#7);
  550.         #75 : if not steps then dec(j) else WRITE(#7);
  551.         end;
  552.       end;
  553.     3 : begin
  554.         if numbin > epsilon then err[k] := numbin else begin
  555.            err[k] := err[k-1]; gotoXY(ix,iy); NEATWR(err[k]); end;
  556.         case ch of
  557.          #13,#77 : begin if steps then j := 2 else j := 1; inc(k) end;
  558.          #80     : if k < Npts then inc(k);
  559.          #72     : if k > 1 then dec(k);
  560.          #75     : dec(j);
  561.          end;
  562.        end;
  563.     end;  {case j}
  564.   if k > Npts then begin
  565.     Npts := k;
  566.     top := 1; if k > 15 then top := k - 14;
  567.     for kk := top to k-1 do begin
  568.       iy := kk - top + 1;
  569.       gotoXY(5,iy); write(kk);
  570.       gotoXY(12,iy); NEATWR(x[kk]);
  571.       gotoXY(24,iy); NEATWR(y[kk]);
  572.       gotoXY(36,iy); NEATWR(err[kk]);
  573.       end;
  574.     gotoXY(1,iy+1); ClrEOL;
  575.     if steps then begin
  576.       x[k] := x0 + (k-1)*step; gotoXY(12,iy+1); NEATWR(x[k]); end;
  577.     err[k] := err[k-1];
  578.     end;
  579.    until ch = #27;
  580. npts := npts-1;
  581.  
  582. (*repeat                                   { start entering data points }
  583.    write('Point #',k:2,'   x = ');
  584.    if steps then begin
  585.       x[k] := x0 + (k-1)*step;
  586.       gotoXY(16,whereY); write(x[k]:7:3);
  587.       end
  588.    else READIN(x[k]);
  589.    gotoXY(28,whereY); ClrEol; write('y = ');
  590.    repeat
  591.       READIN(y[k]);
  592.       if count and (y[k] < 0) then begin
  593.          write('***NEG COUNT NOT ALLOWED***');
  594.          write(#7);
  595.          Delay(1000);
  596.          gotoXY(32,whereY); ClrEol;
  597.          end;
  598.       until not count or (y[k] > 0);
  599.    gotoXY(44,whereY);
  600.    write('error = ');
  601.    if count then begin
  602.       err[k] := sqrt(y[k]);
  603.       gotoXY(52,whereY);
  604.       writeln(err[k]:8:4);
  605.       end
  606.    else begin err[k] := err[k-1]; READIN(err[k]); writeln end;
  607.    write('More? (<CR>=Yes)');
  608.    ch := ReadKey; ch := Upcase(ch);
  609.    gotoXY(1,whereY); ClrEol;
  610.    if (ch = #13) or (ch = 'Y') and (err[k] > 0) then inc(k);
  611.    until  (k > maxpts) or (Upcase(ch) = 'N') or (Upcase(ch) = 'Q');
  612. Npts := k;*)
  613.  
  614. Window(1,1,80,25);
  615. Textbackground(0);
  616.              { plot data and correct mistakes }
  617. if Npts > 2 then begin
  618.    SORT; LIMITS; DRAWFRAME; PLOTPOINTS;
  619.    MoveTo(250,yscale+ 26); OutText('Is this data OK?'); BLINK;
  620.    if not YN then REPAIR;
  621.    RestoreCrtMode; Textcolor(7);
  622.    for k := 1 to Npts do sigsq[k] := err[k]*err[k];
  623.    writeln; write(' Save input data to disk? ');
  624.    if YN then WRITEFILE;
  625.    filefound := True;
  626.    end;
  627. end;  { of KEYBOARD }
  628.  
  629. Procedure LEGEND(xx:float; nn:word);         {  Legendre polynomials }
  630. var  m : word;                               {  by forward recursion }
  631. begin
  632. P[0] := 1;
  633. P[1] := xx;
  634. P[2] := (3.0*xx*xx - 1.0)/2.0;
  635. if nn > 2 then
  636.    for m := 3 to nn do P[m] := ((2*m-1)*xx*P[m-1] - (m-1)*P[m-2])/m;
  637. end;
  638.  
  639. Procedure CHEB(xx: float; nn:word);          { Chebychev polynomials }
  640. var  m : word;                               { by forward recursion }
  641. begin
  642. T[0] := 1;
  643. T[1] := xx;
  644. T[2] := 2.0*xx*xx - 1;
  645. if nn > 2 then for  m := 3 to nn do T[m] := 2.0*xx*T[m-1] - T[m-2];
  646. end;
  647.  
  648. Procedure BESSY(xx:float; nn:word);          { Bessel functions }
  649. var  n,nmax : word;                          { by backward recursion }
  650.      J : array[0..50] of float;
  651.      sum : float;
  652. begin
  653. for n := 0 to nn do JJ[n] := 0;
  654. if xx > 0.001 then begin
  655.    nmax := Trunc(9.0*sqrt(xx+1.0));
  656.    J[nmax] := 1.0;
  657.    J[nmax+1] := 0.0;
  658.    sum := 0;
  659.    for n := nmax downto 1 do begin
  660.       J[n-1] := (2.0*n/xx)*J[n] - J[n+1];
  661.       if  n mod 2 = 0 then sum := sum + 2.0*J[n];
  662.       end;
  663.    sum := sum + J[0];
  664.    for n := 0 to nn do JJ[n] := J[n]/sum;
  665.    end
  666. else JJ[0] := 1.0;
  667. end;
  668.  
  669. Procedure RANGE;   { set range in x over which functions will operate }
  670. begin
  671. Textcolor(7);
  672. writeln;
  673. if funct in [powers,cheby,legendre,expon] then begin
  674.       lmin := xmin; lmax := xmax;       {Default values}
  675.       write('xmin,',xmin:8:2,' equivalent to: ');
  676.       if funct in [cheby,legendre] then write('(must be >= -1.0 )');
  677.       Textcolor(14); READIN(lmin); Textcolor(7); writeln;
  678.       write('xmax,',xmax:8:2,' equivalent to: ');
  679.       if funct in [cheby,legendre] then write('(must be <= +1.0 )');
  680.       Textcolor(14); READIN(lmax); Textcolor(7); writeln;
  681.       len := (xmax - xmin)/(lmax - lmin);
  682.    end;
  683. if funct in [sines..oddsines] then begin
  684.    lmin := 0; cy := 1.0;                             {default values}
  685.    write('x range represents how many cycles? ');
  686.    Textcolor(14); READIN(cy); Textcolor(7); writeln;
  687.    len := (xmax-xmin)/(2*Pi*cy);
  688.    end;
  689. if funct = bessels then begin          { special treatment for Bessels}
  690.    lmin := 0; lmax := 1.0;  {default}
  691.    repeat
  692.       write('xmin, ',xmin:8:2,' equivalent to: (must be >= 0)');
  693.       Textcolor(14); READIN(lmin); Textcolor(7); writeln;
  694.       write('xmax, ',xmax:8:2,' equivalent to: ');
  695.       Textcolor(14); READIN(lmax); Textcolor(7); writeln;
  696.       if lmax > 25.0 then writeln
  697.          ('Program can"t calculate Bessel funct for arguments this big');
  698.       until (lmax > lmin) and (lmax < 25.0);
  699.    len := (xmax-xmin)/(lmax-lmin);
  700.    end;
  701. if funct = gauss then begin
  702.    lmin := 0; len := 1.0 end;  {start at 1st data point}
  703. end; { of RANGE }
  704.  
  705. Procedure DIMENSION;                   { Input number of terms in fit }
  706. begin
  707. if logplt then begin
  708.   if funct = expon then Dim := 2;
  709.   if funct = gauss then Dim := 3
  710.   end
  711. else repeat
  712.    writeln; write(' Number of terms: ');
  713.    TextColor(14);
  714.    Readln(Dim);
  715.    if (Dim < 2) or (Dim > 12) or (Dim >= Npts) then write(#7);
  716.    if Dim < 2 then writeln('***** MUST HAVE AT LEAST TWO TERMS *****');
  717.    if Dim > 12 then writeln('***** NO MORE THAN 12 TERMS, PLEASE *****');
  718.    if Dim >= Npts then writeln('*****TOO MANY TERMS, NOT ENOUGH DATA *****');
  719.    TextColor(7);
  720.    until (Dim < Npts) and (Dim in [2..max+1]);
  721. N := Dim-1; DOF := Npts - Dim;
  722. end;
  723.  
  724. Function FUNC(x:float; nn:word):float;  { Calculate f[nn](x) according to }
  725. var f,xx : float; m: word;              {    species of function          }
  726. begin
  727. f := 1.0;
  728. xx := lmin + (x - xmin)/len;
  729. case funct of
  730.     powers,expon,gauss : for m := 1 to nn do f := f*xx;
  731.     fourier : begin m := (nn+1) div 2;
  732.               if nn mod 2 = 0 then f := cos(m*xx) else f := sin(m*xx); end;
  733.     sines   :  if nn > 0 then f := sin(nn*xx);
  734.     cosines :  f := cos(nn*xx);
  735.     oddsines :  begin  m := 2*nn+1; f := sin(m*xx); end;
  736.     cheby   :  begin CHEB(xx,nn); f := T[nn]; end;
  737.     legendre:  begin LEGEND(xx,nn); f := P[nn]; end;
  738.     bessels :  begin BESSY(xx,nn); f := JJ[nn]; end;
  739.   end;
  740. FUNC := f;
  741. end;
  742.  
  743. Procedure CONVERT;   { Convert data to log of data for expon and gauss fit }
  744. begin
  745. for k := 1 to Npts do begin
  746.    ysave[k] := y[k];
  747.    err[k] := err[k]/y[k];
  748.    sigsq[k] := err[k]*err[k];
  749.    y[k] := ln(y[k]);
  750.    end;
  751. ymax := ln(ymax);
  752. ymin := ln(yleast);     { subroutine LIMITS may have altered ymin }
  753. end;
  754.  
  755. Procedure RESTORE;      { Undo CONVERT }
  756. begin
  757. for k := 1 to Npts do begin
  758.    y[k] := ysave[k];
  759.    err[k] := err[k]*y[k];
  760.    sigsq[k] := err[k]*err[k];
  761.    end;
  762. ymax := exp(ymax);
  763. ymin := exp(ymin);
  764. end;
  765.  
  766.  
  767. Procedure GETFUN;                      { input species of function to fit}
  768. const  name : array[1..11] of phrase = ('Powerseries','Sines','Cosines',
  769.            'Fourier','Oddsines','Bessels','Tchebychev','Legendre',
  770.            'Exponential','Gaussian','View . . .');
  771. var ii : word;
  772.     capt : array[0..11] of textline;
  773.  
  774. Procedure VIEWPOINTS;      { choose function species from graphics display}
  775. begin
  776. DRAWFRAME; PLOTPOINTS;
  777. SetColor(2);
  778. for i := 1 to 10 do begin
  779.    MoveTo(5,21+i*10); ERASE(8*length(name[i]));
  780.    OutTextXY(5,20+i*10,name[i]);
  781.    Setcolor(10); OutTextXY(5,20+10*i,name[i,1]);
  782.    Setcolor(2);
  783.    end;
  784. MoveTo(5,yscale div 2 + 30); OutText('your choice?');
  785. SetColor(10); BLINK;
  786. ch := ReadKey; ch := Upcase(ch);
  787. RestoreCrtMode;
  788. end;
  789.  
  790. begin
  791. capt[1] := ' y = a0 + a1*x/len + a2*(x/len)^2 + . . . ';
  792. capt[2] := ' y = a0 + a1*sin(x/len) + a2*sin(2*x/len) + . . . ';
  793. capt[3] := ' y = a0 + a1*cos(x/len) + a2*cos(2*x/len) + . . . ';
  794. capt[4] := ' y = a0 + a1*sin(x/len) + a2*cos(x/len) + . . . ';
  795. capt[5] := ' y = a0*sin(x/len) + a1*sin(3*x/len) + . . .';
  796. capt[6] := ' y = a0*J0(x/len) + a1*J1(x/len) + . . . ';
  797. capt[7] := ' y = a0*T0(x/len) + a1*T1(x/len) + . . . ';
  798. capt[8] := ' y = a0*P0(x/len) + a1*P1(x/len) + . . . .';
  799. capt[9] := ' y = a0*exp(x/a1)  ';
  800. capt[10] := ' y = a0*exp[-(x-a1)^2/2*a2^2]  ';
  801. capt[11] := 'a plot of the data';
  802. repeat
  803.    ClrScr; logplt := False;
  804.    writeln('Source file: ',dsource); writeln;
  805.    writeln('                       Catalog of Functions'); writeln;
  806.    for i := 1 to 11 do begin
  807.       write('        '); CWRITE(name[i]);
  808.       gotoXY(22,whereY); writeln(capt[i]);
  809.       end;
  810.    writeln; Textcolor(14); write(' your choice: ');
  811.    ch := ReadKey; ch := Upcase(ch);
  812.    if ch = 'V' then VIEWPOINTS;
  813.    until ch in ['P','S','C','F','O','T','L','B','E','G'];
  814. case ch of
  815.    'P' : begin funct := powers;  ii := 1 end;
  816.    'S' : begin funct := sines; ii := 2 end;
  817.    'C' : begin funct := cosines; ii := 3 end;
  818.    'F' : begin funct := fourier; ii := 4 end;
  819.    'O' : begin funct := oddsines; ii := 5 end;
  820.    'B' : begin funct := bessels; ii := 6 end;
  821.    'T' : begin funct := cheby; ii := 7 end;
  822.    'L' : begin funct := legendre; ii := 8 end;
  823.    'E' : begin funct := expon; ii := 9; logplt := true; end;
  824.    'G' : begin funct := gauss; logplt := true; ii := 10 end;
  825.    end;
  826. fname := name[ii]; caption := capt[ii];
  827. if logplt and (yleast < 0.001) then begin
  828.    writeln('There are negative or zero values of y in your data.');
  829.    writeln('Can"t use this species of function with your data');
  830.    write(#7);
  831.    Delay(2000);
  832.    GETFUN
  833.    end;
  834. Textcolor(14); writeln(' ',fname); writeln;
  835. end;  { of GETFUN }
  836.  
  837. Procedure PRECALC;        { to save time, calculate f[n](x) at all x[k] }
  838. var xx : float;
  839. begin;                    { put precalculated functions on the heap}
  840.                                      { to save space in data segment}
  841. for n := 1 to Npts do GetMem(fpoint[n],sizeof(coeff));
  842. writeln('Wait . . precalculating functions');
  843. if funct = bessels then begin                   { save time }
  844.    for k := 1 to Npts do begin                  { do them 10 at a time}
  845.       xx := lmin + (x[k]-xmin)/len;
  846.       BESSY(xx,max);
  847.       for n := 0 to max do fpoint[k]^[n] := JJ[n];
  848.       end;
  849.    end
  850. else for k := 1 to Npts do
  851.       for n := 0 to max do fpoint[k]^[n] := FUNC(x[k],n);
  852. end;
  853.  
  854. Function YCALC(x:float):float;    {  with recursive functions save time }
  855. var xx,y: float;                  {  by calculating a bunch of them at once }
  856. begin
  857. y := 0;
  858. xx := lmin + (x-xmin)/len;
  859. case funct of
  860.    cheby : begin CHEB(xx,N); for i := 0 to N do y := y + a[i]*T[i]; end;
  861.    legendre : begin LEGEND(xx,N); for i := 0 to N do y := y + a[i]*P[i]; end;
  862.    bessels : begin BESSY(xx,N); for i := 0 to N do y := y + a[i]*JJ[i]; end;
  863.    else for i := 0 to N do y := y + a[i]*FUNC(x,i);
  864.    end;
  865. YCALC := y;
  866. end;
  867.  
  868. Procedure CHISQUARE;
  869. begin
  870. for k := 1 to Npts do yc[k] := YCALC(x[k]);
  871. chisq := 0;
  872. for k := 1 to Npts do chisq := chisq + sqr(y[k] - yc[k])/sigsq[k]
  873. end;
  874.  
  875. Procedure INVERSE;                   { by Gauss-Jordan elimination }
  876. var   fact: float;
  877. begin
  878. bombed := False;
  879. for i := 0 to N do for j := 0 to N do temp[i,j] := 0;
  880. for i := 0 to N do temp[i,i] := 1.0;
  881. for i := 0 to N do begin
  882.    if abs(h[i,i]) > bombthreshold then begin
  883.       if i > 0 then begin
  884.          for j := 0 to i-1 do begin  {zero terms below the diagonal}
  885.             fact := h[j,i]/h[i,i];
  886.             for k := 0 to N do h[j,k] := h[j,k] - fact*h[i,k];
  887.             for k := 0 to N do temp[j,k] := temp[j,k] -fact*temp[i,k];
  888.             end;
  889.          end;
  890.       if i < N then begin
  891.          for j := i+1 to N do begin {zero terms above the diagonal}
  892.             fact := h[j,i]/h[i,i];
  893.             for k := 0 to N do h[j,k] := h[j,k] -fact*h[i,k];
  894.             for k := 0 to N do temp[j,k] := temp[j,k] - fact*temp[i,k];
  895.             end;
  896.          end;
  897.       end
  898.       else bombed := True;
  899.    end;
  900. if bombed then writeln('MATRIX INVERSION BOMBED') else
  901.    for i := 0 to N do
  902.         for j := 0 to N do hinv[i,j] := temp[i,j]/h[i,i];
  903. end;  { of INVERSE }
  904.  
  905. Procedure MATRIX;                 { this procedure finds the a[i] }
  906. var  wy  : vector;
  907. begin
  908. writeln('Building matrix');
  909. for i := 0 to N do begin          {  build matrix on  }
  910.    for j := i to N do begin       {  diagonal and above  }
  911.       h[i,j] := 0;
  912.       for k := 1 to Npts do
  913.           h[i,j] := h[i,j] + fpoint[k]^[i]*fpoint[k]^[j]/sigsq[k];
  914.       h[j,i] := h[i,j];    { to save time, since matrix is symmetric }
  915.       end;
  916.    end;
  917. writeln('Inverting matrix');
  918. INVERSE;
  919. if not bombed then begin
  920.    writeln('Calculating coefficients');
  921.    for k := 1 to Npts do wy[k] := y[k]/sigsq[k];
  922.    for i := 0 to N do begin
  923.       a[i] := 0;
  924.       for j := 0 to N do
  925.          for k := 1 to Npts do
  926.             a[i] := a[i] + wy[k]*fpoint[k]^[j]*hinv[j,i];
  927.       end;
  928.    end; {if not bombed}
  929. end; { Procedure MATRIX }
  930.  
  931. Procedure PLOTCURVE;
  932. var  x2,y2,delta : float;
  933. begin
  934. Setcolor(10);
  935. ix2 := 100; x2 := xmin;
  936. delta := (xmax - xmin)/125.0;
  937. y2 := YCALC(xmin); iy2 := 10 + Round(yscale*(ymax - y2)/(ymax - ymin));
  938. for k := 1 to 125 do begin;
  939.     ix1 := ix2; iy1 := iy2;
  940.     ix2 := 4*k + 100;
  941.     x2 := x2 + delta; y2 := YCALC(x2);
  942.     iy2 := 10 + Round(yscale*(ymax - y2)/(ymax - ymin));
  943.     Line(ix1,iy1,ix2,iy2);
  944.     end;
  945. Str(Dim,dd); dd := ', ' + dd;
  946. MoveTo(250,yscale + 26); OutText(fname); OutText(dd); OutText(' terms');
  947. end;
  948.  
  949. Procedure QUERY;
  950. var caption : textline;
  951. begin
  952. plotted := true;
  953. if logplt then OutTextXY(250,yscale+35,'Press any key to continue ')
  954. else begin
  955.    MoveTo(120,yscale+38); OutText('Chisq = '); OUTWR(chisq);
  956.    Str(DOF,dd);
  957.    caption := ' for ' + dd + ' degrees of freedom.   Accept fit? ';
  958.    OutText(caption); BLINK
  959.    end;
  960. ch := ReadKey;
  961. RestoreCrtMode;
  962. Textcolor(7); ClrScr;
  963. end;
  964.  
  965. Procedure FUDGE;
  966. begin
  967. if funct = expon then begin
  968.    a[0] := ln(a[0]); a[1] := 1/a[1]; CONVERT end;
  969. if funct = gauss then begin
  970.    a[0] := asave0; a[1] := asave1; a[2] := asave2; CONVERT end;
  971. {could have saved the plot image and restored here; this way uses less memory}
  972. DRAWFRAME; PLOTPOINTS; PLOTCURVE; REPAIR;
  973. if logplt then RESTORE;
  974. end;
  975.  
  976. Function CL(Chisq:float; DOF:integer):float;      { find confidence level }
  977. var pchi,h,twoexph,gammah,x,y,step,sum : float;
  978. begin
  979. if chisq/DOF > 12.0 then begin  CL := 0.0; Exit end;
  980. if DOF/(chisq +0.0001) > 12.0 then begin  CL := 1.0; Exit end;
  981. h := DOF/2.0;            { h for half }
  982. if (chisq < 50.0) and (h <= 20.0) then begin
  983.    twoexph := exp(h*ln(2.0));      { = 2^h }
  984.    gammah := 2.507*exp((h-0.5)*ln(h))*(1.0+0.0833/h)/exp(h);   { Gamma(h) by }
  985.    step := Chisq/100.0;                                 { Stirling's approx }
  986.    x := 0; sum := 0;
  987.    for i := 1 to 100 do begin                    { Simpson's rule integration }
  988.       x := x + step;
  989.       pchi := exp((h-1.0)*ln(x))/exp(x/2.0);
  990.       sum := sum + 2.0*step*pchi;
  991.       if (i mod 2) > 0 then sum := sum + 2.0*step*pchi;
  992.       end;
  993.    CL := 1.0 - (sum - step*pchi)/(3.0*twoexph*gammah);
  994.    end
  995. else begin           { DOF > 40, assume p(chisq) normally distributed }
  996.    x := sqrt(2.0*chisq) - sqrt(4.0*h - 1.0);
  997.    y := abs(x);
  998.    {power series approx to integrated error funct calculated by LEASTQ }
  999.    {this expression apparently does not (quite) overflow the 8087 stack }
  1000.    sum := 5000.16 + y*(3978.21 + y*(73.0309 -y*(824.04 -y*(135.194 +
  1001.               y*(75.357 - y*(28.247 -2.73*y))))));
  1002.    if x < 0 then CL := sum/10000.0 else CL := 1.0 - sum/10000.0;
  1003.    if x > 3.5 then CL := 0.0;
  1004.    if x < -3.5 then CL := 1.0;
  1005.    end;
  1006. end;  { of CL }
  1007.  
  1008. Procedure COMPARE;
  1009. var
  1010.     column,k1,k2,k3 : word;
  1011. begin
  1012. if logplt then begin     {convert yc[k] back to arithmetic and recalc chisq}
  1013.    chisq := 0;
  1014.    for k := 1 to Npts do begin
  1015.      yc[k] := exp(yc[k]);
  1016.      chisq := chisq + Sqr(y[k]-yc[k])/sigsq[k];
  1017.      end;
  1018.      DECIMALS;
  1019.    end;
  1020. writeln;
  1021. writeln
  1022. ('   ',dsource,': COMPARISON OF DATA WITH CALCULATED VALUES  (* = fudged) ');
  1023. for k := 1 to 3 do write('  x      ydata   ycalc    '); writeln;
  1024. column := Npts div 3;
  1025. Window(1,Dim+6,80,22);
  1026. gotoXY(1,1);
  1027. if Npts mod 3 > 0 then Inc(column);           { write in snake format }
  1028. for k := 1 to column do begin
  1029.    k1 := k; k2 := column + k; k3 := 2*column + k;
  1030.    write(x[k1]:6:xd,y[k1]:8:yd,yc[k1]:8:yd);
  1031.    if fudged[k1] then write('*   ') else write ('    ');
  1032.    write(x[k2]:6:xd,y[k2]:8:yd,yc[k2]:8:yd);
  1033.    if fudged[k2] then write('*   ') else write ('    ');
  1034.    if k3 <= Npts then begin
  1035.       write(x[k3]:6:xd,y[k3]:8:yd,yc[k3]:8:yd);
  1036.       if fudged[k3] then write('*')
  1037.       end;
  1038.    writeln;
  1039.    if k mod (16-Dim) = 0 then begin
  1040.        write('. . . . hit <CR> for more . . . . ');
  1041.        readln; gotoXY(1,WhereY-1) end;
  1042.    end;
  1043. Window(1,1,80,25);
  1044. end;
  1045.  
  1046. Procedure GAUSSCONV;
  1047. var del0,del1,del2 : float;
  1048.     b : coeff;
  1049. begin
  1050. b[0] := a[0] - a[1]*a[1]/(4*a[2]);           { convert to parabola }
  1051. b[1] := -a[1]/(2*a[2]);                      { centered on b[1]    }
  1052. b[2] := -a[2];
  1053. del0 := errc[0]; del1 := errc[1]; del2 := errc[2];
  1054. errc[0] := sqrt(del0*del0 +sqr(del1*b[1]) + sqr(del2*b[1]*b[1]));
  1055. errc[1] := sqrt(sqr(del1/(2*a[2]))+ sqr(del2*a[1]/(2*a[2]*a[2])));
  1056. errc[2] := del2/(4*sqrt(b[2]*b[2]*b[2]));
  1057. asave0 := a[0]; asave1 := a[1]; asave2 := a[2];    { save for FUDGE}
  1058. a[0] := exp(b[0]); errc[0] := exp(errc[0]);     { to fit standard }
  1059. a[1] := xmin + b[1];                            { definition      }
  1060. a[2] := 1.0/sqrt(2.0*b[2]);                     { of a Gaussian   }
  1061. end;
  1062.  
  1063. Procedure OUTWRITE;  { display coefficents, a[i] and calculate CL }
  1064. begin
  1065. ClrScr; Textcolor(7);
  1066. for i := 0 to N do errc[i] := sqrt(hinv[i,i]);
  1067. if funct = gauss then GAUSSCONV;
  1068. if funct = expon then begin
  1069.    a[0] := exp(a[0]); a[1] := 1.0/a[1];
  1070.    errc[0] := a[0]*errc[0];
  1071.    errc[1] := a[1]*a[1]*errc[1]
  1072.    end;
  1073. write('a[i] for ',caption);
  1074. if funct = gauss then writeln( ' center (a1) = ',a[1]:7:3)
  1075.    else writeln(' len =',len:7:4);
  1076. writeln;
  1077. for i := 0 to N do begin
  1078.    pct[i] := 100*errc[i]/abs(a[i]+0.0001*errc[i]);
  1079.    writeln
  1080. ('          a',i,' = ',a[i]:13,'   +/- ',errc[i]:10,'  (',pct[i]:5:1,' % )');
  1081.    end;
  1082. COMPARE;
  1083. gotoXY(1,22);
  1084. write('Chisq = ',chisq:7:3,' for ',DOF,' degrees of freedom. ');
  1085. write('Calculating CL, wait . . .');
  1086. confid := CL(chisq,DOF);
  1087. gotoXY(42,whereY);  writeln('     Confidence level = ',confid:6:4);
  1088. end; { of OUTWRITE }
  1089.  
  1090. Procedure PRINTOUT;
  1091. var column,k1,k2,k3 : word;
  1092.     year,month,day,dayofweek,hour,minute,second,sec100 : word;
  1093. begin
  1094. Write(Lst,'      Data file: ',dsource,'  fitted at ');
  1095. GetDate(year,month,day,dayofweek);
  1096. getTime(hour,minute,second,sec100);
  1097. writeln(Lst,hour:2,':',minute:2,'  on ', month:2,'/',day:2,'/',(year-1900):2);
  1098. writeln(Lst,'                       function ',fname);
  1099. writeln(Lst);
  1100. write(Lst,'a[i] for ',caption);
  1101. if funct = gauss then writeln(Lst,' center (a1) = ',a[1]:7:3)
  1102.    else writeln(Lst,' len =',len:7:3);
  1103. writeln(Lst);
  1104. for i := 0 to N do writeln(Lst,
  1105. '          a',i,' = ',a[i]:13,'   +/- ',errc[i]:10,'  (',pct[i]:5:1,' % )');
  1106. writeln(Lst);
  1107. writeln(Lst,'                 COMPARISON OF DATA WITH CALCULATED VALUES ');
  1108. for k := 1 to 3 do write(Lst,'  x      ydata   ycalc    '); writeln(Lst);
  1109. column := Npts div 3;
  1110. if Npts mod 3 > 0 then Inc(column);           { write in snake format }
  1111. for k := 1 to column do begin
  1112.    k1 := k; k2 := column + k; k3 := 2*column + k;
  1113.    write(Lst,x[k1]:6:xd,y[k1]:8:yd,yc[k1]:8:yd);
  1114.    if fudged[k1] then write(Lst,'*   ') else write(Lst,'    ');
  1115.    write(Lst,x[k2]:6:xd,y[k2]:8:yd,yc[k2]:8:yd);
  1116.    if fudged[k2] then write(Lst,'*   ') else write(Lst,'    ');
  1117.    if k3 <= Npts then write(Lst,x[k3]:6:xd,y[k3]:8:yd,yc[k3]:8:yd);
  1118.    if fudged[k3] then writeln(Lst,'*') else writeln(Lst);
  1119.    end;
  1120. writeln(Lst,'                                   points marked "*" are fudged');
  1121. writeln(Lst);
  1122. write(Lst,'Chisq = ',chisq:7:3,' for ',DOF,' degrees of freedom. ');
  1123. writeln(Lst,'     Confidence level = ',confid:6:4);
  1124. end;
  1125.  
  1126. begin                    { MAIN PROGRAM }
  1127. if copywt = '' then;
  1128. if length(paramstr(1)) > 2 then begin
  1129.    dsource := paramstr(1); pswitch := paramstr(2) end
  1130.       else begin dsource := ''; pswitch := paramstr(1)  end;
  1131.  
  1132. (* makes EXE file with built in BGI driver for EGA/VGA *)
  1133. if RegisterBGIdriver(@EGAVGADriverProc) < 0 then ABORT('EGA/VGA');
  1134. (* you can compile without this, and without 'drivers' in the uses
  1135.    statement. You will then need the EGAVGA.BGI file in the same
  1136.    directory for LEASTQ.EXE to run *)
  1137.  
  1138. GrDriver := Detect;
  1139. if (pswitch = '/e') or (pswitch = '/E') then begin
  1140.    GrDriver := EGA; GrMode := 1 end;          { To force EGA }
  1141. InitGraph(GrDriver, GrMode, '');
  1142. ErrorCode := GraphResult;                     { Check for errors }
  1143. if ErrorCode <> grOK then ABORT('No Graphics');
  1144. yscale := GetMaxY - 50;  {Vertical scale factor = 300 for EGA, 430 for VGA}
  1145. RestoreCrtMode;
  1146. datinput:                                 { new Data option starts here }
  1147. plotted := false; logplt := False;
  1148. for k := 1 to 255 do fudged[k] := false;
  1149. if dsource = '' then DATAFILES else CAPITALIZE(dsource);
  1150. if dsource = 'NONAME' then KEYBOARD else READFILE;
  1151. if not filefound then begin
  1152.    writeln('******* file not found *******');
  1153.    write(#7);
  1154.    Delay(1000);
  1155.    dsource := '';
  1156.    goto datinput
  1157.    end;
  1158. if Npts < 3 then begin
  1159.    write('Not enough valid points');
  1160.    write(#7);
  1161.    Delay(1000);
  1162.    goto datinput
  1163.    end;
  1164. start:                                    { new Species option starts here }
  1165. logplt := False;
  1166. GETFUN;
  1167. xlimits:                                  { new Range option starts here }
  1168. if logplt then CONVERT;
  1169. RANGE; PRECALC;
  1170. terms:                                    { choose number of terms in fit }
  1171. DIMENSION;
  1172. MATRIX;                                   { calculate coefficients }
  1173. if bombed then
  1174.    if Dim = 2 then goto start      { else no way out of program }
  1175.       else goto terms;
  1176. CHISQUARE;
  1177. DRAWFRAME; PLOTPOINTS; PLOTCURVE; QUERY;  { display curve fit }
  1178. if logplt then RESTORE
  1179.   else if UpCase(ch) = 'N' then goto terms; { user dissatisfied, try again}
  1180. OUTWRITE;                                 { record coefficients }
  1181. CWRITE
  1182. ('new Data  new function Species  new Range in x   Fudge data  Quit program ');
  1183. writeln; Textcolor(14); write('your choice:'); Textcolor(7);
  1184. repeat
  1185.    ch := Readkey; ch := Upcase(ch);
  1186.    until ch in ['D','S','R','F','Q'];
  1187. case ch of
  1188.    'D' : begin dsource := ''; goto datinput end;
  1189.    'S' : goto start;
  1190.    'R' : goto xlimits;
  1191.    'F' : begin FUDGE; goto start end;
  1192.    end;
  1193. CloseGraph;
  1194. Write(' Do you want a printed record of your last fit?');
  1195. if YN then PRINTOUT;
  1196. end.
  1197.