{$N+}
PROGRAM dataplot (inputfile, output);

USES
  Graph, Crt, Dos;

TYPE
  expdt = array[0..2049] OF longint;
  fitdt = array[0..2049] of real;
  fitpar = array[1..8] OF real;

VAR
  x, y       : expdt;
  fity,fitp1,fitp2   : fitdt;
  aa, deltaa : fitpar;

  inf : text;
  ch  : char;

  i, j, nx, ny, ncn, nco, limit     : longint;
  chini, chfin, chrange, ichi       : longint;
  GraphDriver, GraphMode, Errorcode : integer;
  maxx, maxy, maxcolor              : word;
  test1, test2, test3, derr, scale  : real;
  del, xx, chitot, chipt            : real;
  area1, area2, sum, sum1, sum2     : real;
  fname, s1, s2, s3, s4, s5, s6     : string;
  s7, s8, s9                        : string;
  running, stop, stopp, twocurve    : boolean;
  ctot, chctot                      : longint;
  a20,a30,delx,dely                 : real;
{---------------------------------------------------------------------------}

PROCEDURE startgraf;
BEGIN
  Directvideo := false;
  Graphdriver := detect;
  Initgraph(Graphdriver, Graphmode, '');
  maxx := getmaxx;
  maxy := getmaxy;
END;

{---------------------------------------------------------------------------}

FUNCTION xl(ch : integer) : integer;
BEGIN
  xl := 40 + trunc((maxx - 41) * ch / limit);
END;

FUNCTION yl(dat : longint) : longint;
BEGIN
  yl := maxy - 30 - trunc(scale * dat);
END;

FUNCTION xl1(ch : longint) : longint;
BEGIN
  IF chrange < 210 THEN
    xl1 := 40 + 2 * (ch - chini)
  ELSE
    xl1 := 40 + trunc(420.0 * (ch - chini) / chrange);
END;

{---------------------------------------------------------------------------}
{-- Plot expermiental data --}
{---------------------------------------------------------------------------}

PROCEDURE Curser;
BEGIN
  IF ncn > limit THEN
    ncn := 1;
  IF ncn < 1 THEN
    ncn := limit;

  IF y[nco] > 1 THEN
  BEGIN
    setcolor(0);
    nx := xl(nco);
    ny := yl(y[nco] - 1);
    line(nx, maxy - 31, nx, ny);
    setcolor(11);
    rectangle(40, 10, maxx, maxy - 30);
    FOR i := nco - 4 TO nco + 4 DO
    BEGIN
      IF (i > 0) AND (i <= limit) THEN
      BEGIN
        nx := xl(i);
        ny := yl(y[i]);
        putpixel(nx, ny, 15);
      END;
    END;
  END;

  IF y[ncn] > 1 THEN
  BEGIN
    setcolor(10);
    nx := xl(ncn);
    ny := yl(y[ncn] - 1);
    line(nx, maxy - 31, nx, ny);
  END;

  nco := ncn;

  setfillstyle(solidfill, 0);
  setcolor(10);
  str(ncn, s1);
  str(y[ncn], s2);
  ny := maxy - 8;
  bar(330, ny,374, maxy);
  outtextxy(330, ny, s1);
  bar(437, ny, 600, maxy);
  outtextxy(437, ny, s2);
 END;

{---------------------------------------------------------------------------}

PROCEDURE Dataplott;   {plots the experimental data}
BEGIN
  cleardevice;
  setcolor(11);
  rectangle(40, 10, maxx, maxy - 30);

  settextjustify(lefttext, toptext);
  outtextxy(450,0,'Filename : ' + fname);

  settextjustify(centertext, toptext);
  j := trunc(limit / 10.24);
  FOR i := 0 TO 10 DO
  BEGIN
    nx := xl(i * j);
    line(nx, maxy - 25, nx, maxy - 30);
    str(i * j, s1);
    outtextxy(nx, maxy - 20, s1);
  END;

  settextjustify(righttext, centertext);
  FOR i := 0 TO 4 DO
  BEGIN
    j := trunc(i * 100 / scale);
    ny := yl(j);
    line(35, ny, 40, ny);
    str(j, s1);
    outtextxy(32, ny, s1);
  END;

  ny := maxy - 8;
  setcolor(10);
  settextjustify(lefttext, toptext);
  outtextxy(205, ny, 'channel number:');
  outtextxy(375, ny, 'counts:');

  FOR i := 1 TO limit DO            {Plot the data}
  BEGIN
    nx := xl(i);
    ny := yl(y[i]);
    putpixel(nx, ny, 15);
  END;

  setfillstyle(solidfill,7);
  bar(450, 20, 630, 100);
  setcolor(11);
  rectangle(450, 20, 630, 100);
  setcolor(0);
  outtextxy(465, 28, '  decrease scale');
  outtextxy(465, 38, '  increase scale');
  outtextxy(465, 48, '  move left (fast)');
  outtextxy(465, 58, '  move left (slow)');
  outtextxy(465, 68, '  move right (slow)');
  outtextxy(465, 78, '  move right (fast)');
  outtextxy(465, 88, '    exit');
  setcolor(4);
  outtextxy(465, 28, '[');
  outtextxy(465, 38, ']');
  outtextxy(465, 48, 'M');
  outtextxy(465, 58, '<');
  outtextxy(465, 68, '>');
  outtextxy(465, 78, '/');
  outtextxy(465, 88, 'Esc');

  curser;
END;

{---------------------------------------------------------------------------}

PROCEDURE plotdat;
BEGIN
  nco := 1;
  ncn := 1;
  stopp := false;

  startgraf;
  Dataplott;

  REPEAT
  IF keypressed THEN
  BEGIN
    ch := readkey;
    IF Ord(ch) = 27 THEN
      stopp := true
    ELSE
      CASE ch OF
        '[','{' : BEGIN
                    scale := scale / 2.0;
                    Dataplott;
                  END;
        ']','}' : BEGIN
                    scale := scale * 2.0;
                    Dataplott;
                  END;
        '.','>' : BEGIN
                    ncn := ncn + 1;
                    curser;
                  END;
        ',','<' : BEGIN
                    ncn := ncn - 1;
                    curser;
                  END;
        '/','?' : BEGIN
                    ncn := ncn + 10;
                    curser;
                  END;
        'm','M' : BEGIN
                    ncn := ncn - 10;
                    curser;
                  END;
      END;  {CASE}
  END;
  UNTIL stopp;

  closegraph;
END;  {dataplot}

{---------------------------------------------------------------------------}
{-- Curve fit one curve --}
{---------------------------------------------------------------------------}

PROCEDURE manual; {Calculates total chisq for current parameters}
BEGIN
  chitot := 0.0;

  FOR ichi := chini TO chfin DO
  BEGIN
    xx := sqr((ichi - aa[1]) / aa[2]);
    IF y[ichi] = 0 THEN
      derr := 1.0
    ELSE
      derr := sqrt(y[ichi]);
    sum := aa[4];
    IF ichi > (aa[1] - 2.0 * aa[2]) THEN
      sum := (aa[4] - aa[5]) * (aa[1] + 2.0 * aa[2] - ichi) / 4.0 /
             aa[2] + aa[5];
    IF ichi > (aa[1] + 2.0 * aa[2]) THEN
      sum := aa[5];
    IF xx < 10.0 THEN
      sum := sum + (aa[3] - aa[4]) * exp(-xx);
    chitot := chitot + sqr((sum - y[ichi]) / derr);
  END;

  chipt := chitot / chrange;
  area1 := (aa[3] - aa[4]) * aa[2] * sqrt(3.1416);
END;

{---------------------------------------------------------------------------}

PROCEDURE automatic;
BEGIN
  setcolor(11);
  outtextxy(165, 80, 'doing grid search...');

  deltaa[1] := 1.0;
  FOR i := 2 TO 5 DO
    deltaa[i] :=  aa[i] / 100.0;
  deltaa[3] := 1.0;
  deltaa[4] := 1.0;
  deltaa[5] := 1.0;

  FOR j := 1 TO 5 DO
  BEGIN
    aa[j] := aa[j] - deltaa[j];
    manual;
    test1 := chitot;
    aa[j] := aa[j] + deltaa[j];
    manual;
    test2 := chitot;
    IF test2 > test1 THEN
    BEGIN
      deltaa[j] := -deltaa[j];
      test2 := test1;
      test1 := chitot;
    END;
    test3 := test2;

    WHILE test3 <= test2 DO
    BEGIN
      aa[j] := aa[j] + deltaa[j];
      manual;
      test3 := chitot;
      IF test3 < test2 THEN
      BEGIN
        test1 := test2;
        test2 := test3
      END;
    END;

    aa[j] := aa[j] - deltaa[j] *
             ((test3 - test2) / (test3 - 2 * test2 + test1) + 1.0 / 2.0);
  END;
END; {automatic curvefit}

{---------------------------------------------------------------------------}

PROCEDURE autoerror(jjj:integer);
BEGIN
  setcolor(11);
  str(jjj:4,s1);
  outtextxy(165, 80, 'doing grid search without');
  outtextxy(190,90,'parameter '+s1);

  deltaa[1] := 1.0;
  FOR i := 2 TO 5 DO
    deltaa[i] :=  aa[i] / 100.0;
  deltaa[3] := 1.0;
  deltaa[4] := 1.0;
  deltaa[5] := 1.0;

  FOR j := 1 TO 5 DO
  begin
  if j<>jjj then
  BEGIN
    aa[j] := aa[j] - deltaa[j];
    manual;
    test1 := chitot;
    aa[j] := aa[j] + deltaa[j];
    manual;
    test2 := chitot;
    IF test2 > test1 THEN
    BEGIN
      deltaa[j] := -deltaa[j];
      test2 := test1;
      test1 := chitot;
    END;
    test3 := test2;

    WHILE test3 <= test2 DO
    BEGIN
      aa[j] := aa[j] + deltaa[j];
      manual;
      test3 := chitot;
      IF test3 < test2 THEN
      BEGIN
        test1 := test2;
        test2 := test3
      END;
    END;

    aa[j] := aa[j] - deltaa[j] *
             ((test3 - test2) / (test3 - 2 * test2 + test1) + 1.0 / 2.0);
  END;
 end;
END; {autoerror curvefit}

PROCEDURE areaerror;
BEGIN
  setcolor(11);
  outtextxy(165, 80, 'doing grid search without');
  outtextxy(190,90,'changing the area ');

  deltaa[1] := 1.0;
  FOR i := 2 TO 5 DO
    deltaa[i] :=  aa[i] / 100.0;
  deltaa[3] := 1.0;
  deltaa[4] := 1.0;
  deltaa[5] := 1.0;

  FOR j := 1 TO 5 DO
  begin
  if j<>3 then
  BEGIN
    aa[j] := aa[j] - deltaa[j];
    aa[3]:=area1/aa[2]/sqrt(3.1416)+aa[4];
    manual;
    test1 := chitot;
    aa[j] := aa[j] + deltaa[j];
    aa[3]:=area1/aa[2]/sqrt(3.1416)+aa[4];
    manual;
    test2 := chitot;
    IF test2 > test1 THEN
    BEGIN
      deltaa[j] := -deltaa[j];
      test2 := test1;
      test1 := chitot;
    END;
    test3 := test2;

    WHILE test3 <= test2 DO
    BEGIN
      aa[j] := aa[j] + deltaa[j];
      aa[3]:=area1/aa[2]/sqrt(3.1416)+aa[4];
      manual;
      test3 := chitot;
      IF test3 < test2 THEN
      BEGIN
        test1 := test2;
        test2 := test3
      END;
    END;

    aa[j] := aa[j] - deltaa[j] *
             ((test3 - test2) / (test3 - 2 * test2 + test1) + 1.0 / 2.0);
  END;
 end;
END; {autoerror curvefit}

{---------------------------------------------------------------------------}


PROCEDURE borders;
BEGIN
  setcolor(11);
  rectangle(40, 0, 460, maxy - 30);

  settextjustify(centertext, toptext);
  nx := xl1(chini);
  line(nx, maxy - 25, nx, maxy - 30);
  str(chini, s1);
  outtextxy(nx, maxy - 20, s1);

  nx := xl1(chfin);
  line(nx, maxy - 25, nx, maxy - 30);
  str(chfin, s1);
  outtextxy(nx, maxy - 20, s1);

  FOR i := 0 TO 20 DO
  BEGIN
    j := i * 50;
    IF (j >= chini + 25) AND (j <= chfin - 25) THEN
    BEGIN
      nx := xl1(j);
      line(nx, maxy - 25, nx, maxy - 30);
      str(j, s1);
      outtextxy(nx, maxy - 20, s1);
    END;
  END;

  setcolor(10);
  rectangle(470, 0, maxx, 90);
  settextjustify(lefttext,toptext);
  outtextxy(475, 8, 'peak   = ');
  outtextxy(475, 18, 'sigma  = ');
  outtextxy(475, 28, 'height = ');
  outtextxy(475, 48, 'b1     = ');
  outtextxy(475, 58, 'b2     = ');
  outtextxy(475, 78, 'del    = ');

  setcolor(14);
  rectangle(470, 110, maxx, 180);
  outtextxy(475, 118, 'area   = ');
  outtextxy(475, 138, 'total chisq =');
  outtextxy(475, 158, 'chi/data point =');

  setfillstyle(solidfill, 7);
  bar(470, 250, maxx, 449);
  setcolor(11);
  rectangle(470, 250, maxx, 449);
  setcolor(0);
  outtextxy(475, 258, '  decrease scale');
  outtextxy(475, 268, '  increase scale');
  outtextxy(475, 288, '  increase peak');
  outtextxy(475, 298, '  decrease peak');
  outtextxy(475, 308, '  increase sigma');
  outtextxy(475, 318, '  decrease sigma');
  outtextxy(475, 328, '  increase height');
  outtextxy(475, 338, '  decrease height');
  outtextxy(475, 348, '  increase b1');
  outtextxy(475, 358, '  decrease b1');
  outtextxy(475, 368, '  increase b2');
  outtextxy(475, 378, '  decrease b2');
  outtextxy(475, 398, '  change del');
  outtextxy(475, 408, '  redraw');
  outtextxy(475, 428, '       automatic fit');
  outtextxy(475, 438, '    exit');
  setcolor(4);
  outtextxy(475, 258, '[');
  outtextxy(475, 268, ']');
  outtextxy(475, 288, 'A');
  outtextxy(475, 298, 'Z');
  outtextxy(475, 308, 'S');
  outtextxy(475, 318, 'X');
  outtextxy(475, 328, 'D');
  outtextxy(475, 338, 'C');
  outtextxy(475, 348, 'F');
  outtextxy(475, 358, 'V');
  outtextxy(475, 368, 'G');
  outtextxy(475, 378, 'B');
  outtextxy(475, 398, 'L');
  outtextxy(475, 408, 'R');
  outtextxy(475, 428, 'Return');
  outtextxy(475, 438, 'Esc');
END;

{---------------------------------------------------------------------------}

PROCEDURE plotfit;   {Plots the final fit and the experimental data}
BEGIN
  setfillstyle(solidfill, 0);
  bar(0, 0, 460, maxy - 30);
  setcolor(11);
  rectangle(40, 0, 460, maxy - 30);

  settextjustify(righttext, centertext);
  FOR i := 0 TO 4 DO
  BEGIN
    j := trunc(i * 100 / scale);
    ny := yl(j);
    line(35, ny, 40, ny);
    str(j, s1);
    outtextxy(32, ny, s1);
  END;

  FOR i := chini TO chfin DO            {Plot the data}
  BEGIN
    nx := xl1(i);
    ny := yl(y[i]);
    putpixel(nx, ny, 15);
  END;

  manual;

  FOR i := chini TO chfin DO            {plot the experimental data}
  BEGIN
    nx := xl1(i);
    xx := sqr((i - aa[1]) / aa[2]);
    sum := aa[4];
    IF i > (aa[1] - 2.0 * aa[2]) THEN
      sum := (aa[4] - aa[5]) * (aa[1] + 2.0 * aa[2] - i) / 4.0 /
             aa[2] + aa[5];
    IF i > (aa[1] + 2.0 * aa[2]) THEN
      sum := aa[5];
    IF xx < 10 THEN
      sum := sum + (aa[3] - aa[4]) * exp(-xx);
    fity[i]:=sum;
    ny := yl(trunc(sum));
    putpixel(nx,ny,14);
  END;

  str(chitot:10:3, s1);
  str(chipt:10:3, s2);
  str(area1:10:3, s3);

  setcolor(14);
  settextjustify(lefttext, toptext);
  bar(555, 118, maxx - 1, 127);
  outtextxy(555, 118, s3);
  bar(555, 148, maxx - 1, 157);
  outtextxy(555, 148, s1);
  bar(555, 168, maxx - 1, 177);
  outtextxy(555, 168, s2);
END;

{---------------------------------------------------------------------------}

PROCEDURE writeout;
BEGIN
  str(aa[1]:10:3, s1);
  str(aa[2]:10:3, s2);
  str(aa[3]:10:3, s3);
  str(aa[4]:10:3, s4);
  str(aa[5]:10:3, s5);
  str(del:10:3, s6);

  setfillstyle(solidfill, 0);
  settextjustify(lefttext, toptext);
  setcolor(10);
  bar(555, 8, maxx - 1, 17);
  outtextxy(555, 8, s1);
  bar(555, 18, maxx - 1, 27);
  outtextxy(555, 18, s2);
  bar(555, 28, maxx - 1, 37);
  outtextxy(555, 28, s3);
  bar(555, 48, maxx - 1, 57);
  outtextxy(555, 48, s4);
  bar(555, 58, maxx - 1, 67);
  outtextxy(555, 58, s5);
  bar(555, 78, maxx - 1, 87);
  outtextxy(555, 78, s6);
END;

{---------------------------------------------------------------------------}

PROCEDURE Manfit;
BEGIN
  write('input the initial channel for the fit? ');
  readln(chini);
  write('input the final channel for the fit? ');
  readln(chfin);
  writeln;
  writeln('input the starting values: peak sigma height b1 b2');
  writeln('     (hit <space> or <return> between each)');
  write('?');
  readln(aa[1], aa[2], aa[3], aa[4], aa[5]);

  chrange := chfin - chini;
  del := 10.0;
  stop := false;

  startgraf;
  borders;
  plotfit;                 {Plot the fit, data, and parameters}
  writeout;

  REPEAT
    IF keypressed THEN
    BEGIN
      ch := readkey;
      IF Ord(ch) = 27 THEN
        stop := true
      ELSE IF Ord(ch) = 13 THEN
      BEGIN
        automatic;
        plotfit;
      END
      ELSE
        CASE ch OF
          'L','l' : BEGIN
                      del := del / 10.0;
                      IF del < 0.001 THEN
                        del := 100.0;
                    END;
          'R','r' : plotfit;
          'a','A' : aa[1] := aa[1] + del;
          'Z','z' : aa[1] := aa[1] - del;
          'S','s' : aa[2] := aa[2] + del;
          'X','x' : aa[2] := aa[2] - del;
          'D','d' : aa[3] := aa[3] + del;
          'C','c' : aa[3] := aa[3] - del;
          'F','f' : aa[4] := aa[4] + del;
          'V','v' : aa[4] := aa[4] - del;
          'G','g' : aa[5] := aa[5] + del;
          'B','b' : aa[5] := aa[5] - del;
          '[','{' : BEGIN
                      scale := scale / 2.0;
                      plotfit;
                    END;
          ']','}' : BEGIN
                      scale := scale * 2.0;
                      plotfit;
                    END;
          '1','!' : begin
                     autoerror(1);
                     plotfit;
                    end;
          '2','@' : begin
                     autoerror(2);
                     plotfit;
                    end;
          '3','#' : begin
                     areaerror;
                     plotfit;
                    end;
        END;  {CASE}
      writeout;
    END;
  UNTIL stop;

  closegraph;
END; {manfit}

{---------------------------------------------------------------------------}
{-- Curve fit two curves  --}
{---------------------------------------------------------------------------}

PROCEDURE manual2; {Calculates total chisq for current parameters}
BEGIN
  chitot := 0.0;

  FOR ichi := chini TO chfin DO
  BEGIN
    xx := sqr((ichi - aa[1]) / aa[2]);
    IF y[ichi] = 0 THEN
      derr := 1.0
    ELSE
      derr := sqrt(y[ichi]);
    sum := aa[7];
    IF ichi > (aa[1] - 2.0 * aa[2]) THEN
      sum := (aa[7] - aa[8]) * (aa[4] + 2.0 * aa[5] - ichi) /
             (aa[4] - aa[1] + 2.0 * aa[2] + 2.0 * aa[5]) + aa[8];
    IF ichi > (aa[4] + 2.0 * aa[5]) THEN
      sum := aa[8];
    IF xx < 10.0 THEN
      sum := sum + (aa[3] - aa[7]) * exp(-xx);
    xx := sqr((ichi - aa[4]) / aa[5]);
    IF xx < 10.0 THEN
      sum := sum + (aa[6] - aa[7]) * exp(-xx);
    chitot := chitot + sqr((sum - y[ichi]) / derr);
  END;

  chipt := chitot / chrange;
  area1 := (aa[3] - aa[7]) * aa[2] * sqrt(3.1416);
  area2 := (aa[6] - aa[7]) * aa[5] * sqrt(3.1416);
END;

{---------------------------------------------------------------------------}

PROCEDURE automatic2;
BEGIN
  setcolor(11);
  outtextxy(165, 80, 'doing grid search...');

  FOR i := 1 TO 8 DO
    deltaa[i] :=  aa[i] / 100.0;
  deltaa[1] := 1.0;
  deltaa[7] := 1.0;
  deltaa[8] := 1.0;
  deltaa[4] := 1.0;

  FOR j := 1 TO 8 DO
  BEGIN
    aa[j] := aa[j] - deltaa[j];
    manual2;
    test1 := chitot;
    aa[j] := aa[j] + deltaa[j];
    manual2;
    test2 := chitot;
    IF test2 > test1 THEN
    BEGIN
      deltaa[j] := -deltaa[j];
      test2 := test1;
      test1 := chitot;
    END;
    test3 := test2;

    WHILE test3 <= test2 DO
    BEGIN
      aa[j] := aa[j] + deltaa[j];
      manual2;
      test3 := chitot;
      IF test3 < test2 THEN
      BEGIN
        test1 := test2;
        test2 := test3
      END;
    END;

    aa[j] := aa[j] - deltaa[j] *
             ((test3 - test2) / (test3 - 2 * test2 + test1) + 1.0 / 2.0);
  END;
END;  {automatic curvefit}

{---------------------------------------------------------------------------}

PROCEDURE autoerror2(jjj:integer);
BEGIN
  setcolor(11);
  str(jjj:4,s1);
  outtextxy(165, 80, 'doing grid search without');
  outtextxy(190,90,'parameter '+s1);

  FOR i := 1 TO 8 DO
    deltaa[i] :=  aa[i] / 100.0;
  deltaa[1] := 1.0;
  deltaa[4] := 1.0;

  FOR j := 1 TO 8 DO
  begin
  if j<>jjj then
  BEGIN
    aa[j] := aa[j] - deltaa[j];
    manual2;
    test1 := chitot;
    aa[j] := aa[j] + deltaa[j];
    manual2;
    test2 := chitot;
    IF test2 > test1 THEN
    BEGIN
      deltaa[j] := -deltaa[j];
      test2 := test1;
      test1 := chitot;
    END;
    test3 := test2;

    WHILE test3 <= test2 DO
    BEGIN
      aa[j] := aa[j] + deltaa[j];
      manual2;
      test3 := chitot;
      IF test3 < test2 THEN
      BEGIN
        test1 := test2;
        test2 := test3
      END;
    END;

    aa[j] := aa[j] - deltaa[j] *
             ((test3 - test2) / (test3 - 2 * test2 + test1) + 1.0 / 2.0);
  END;
 end;
END;  {autoerror2 curvefit}

{---------------------------------------------------------------------------}

PROCEDURE borders2;
BEGIN
  setcolor(11);
  rectangle(40, 0, 460, maxy - 30);

  settextjustify(centertext, toptext);
  nx := xl1(chini);
  line(nx, maxy - 25, nx, maxy - 30);
  str(chini, s1);
  outtextxy(nx, maxy - 20, s1);

  nx := xl1(chfin);
  line(nx, maxy - 25, nx, maxy - 30);
  str(chfin, s1);
  outtextxy(nx, maxy - 20, s1);

  FOR i := 0 TO 20 DO
  BEGIN
    j := i * 50;
    IF (j >= chini + 25) AND (j <= chfin - 25) THEN
    BEGIN
      nx := xl1(j);
      line(nx, maxy - 25, nx, maxy - 30);
      str(j, s1);
      outtextxy(nx, maxy - 20, s1);
    END;
  END;

  setcolor(10);
  rectangle(470, 0, maxx, 130);
  settextjustify(lefttext, toptext);
  outtextxy(475, 8, 'p1    =');
  outtextxy(475, 18, 'sig1  =');
  outtextxy(475, 28, 'h1    =');
  outtextxy(475, 48, 'p2    =');
  outtextxy(475, 58, 'sig2  =');
  outtextxy(475, 68, 'h2    =');
  outtextxy(475, 88, 'b1    =');
  outtextxy(475, 98, 'b2    =');
  outtextxy(475, 118, 'del   =');

  setcolor(14);
  rectangle(470, 140, maxx, 220);
  outtextxy(475, 148, 'area1 =');
  outtextxy(475, 158, 'area2 =');
  outtextxy(475, 178, 'total chisq =');
  outtextxy(475, 198, 'chi/data point =');

  setfillstyle(solidfill, 7);
  bar(470, 230, maxx, 449);
  setcolor(11);
  rectangle(470, 230, maxx, 449);
  setcolor(0);
  outtextxy(475, 238, '  increase p1');
  outtextxy(475, 248, '  decrease p1');
  outtextxy(475, 258, '  increase sig1');
  outtextxy(475, 268, '  decrease sig1');
  outtextxy(475, 278, '  increase h1');
  outtextxy(475, 288, '  decrease h1');
  outtextxy(475, 298, '  increase p2');
  outtextxy(475, 308, '  decrease p2');
  outtextxy(475, 318, '  increase sig2');
  outtextxy(475, 328, '  decrease sig2');
  outtextxy(475, 338, '  increase h2');
  outtextxy(475, 348, '  decrease h2');
  outtextxy(475, 358, '  increase b1');
  outtextxy(475, 368, '  decrease b1');
  outtextxy(475, 378, '  increase b2');
  outtextxy(475, 388, '  decrease b2');
  outtextxy(475, 398, '  change del');
  outtextxy(475, 408, '  redraw');
  outtextxy(475, 418, '  show both curves');
  outtextxy(475, 428, '       automatic fit');
  outtextxy(475, 438, '    exit');
  setcolor(4);
  outtextxy(475, 238, 'A');
  outtextxy(475, 248, 'Z');
  outtextxy(475, 258, 'S');
  outtextxy(475, 268, 'X');
  outtextxy(475, 278, 'D');
  outtextxy(475, 288, 'C');
  outtextxy(475, 298, 'F');
  outtextxy(475, 308, 'V');
  outtextxy(475, 318, 'G');
  outtextxy(475, 328, 'B');
  outtextxy(475, 338, 'H');
  outtextxy(475, 348, 'N');
  outtextxy(475, 358, 'J');
  outtextxy(475, 368, 'M');
  outtextxy(475, 378, 'K');
  outtextxy(475, 388, ',');
  outtextxy(475, 398, 'L');
  outtextxy(475, 408, 'R');
  outtextxy(475, 418, 'T');
  outtextxy(475, 428, 'Return');
  outtextxy(475, 438, 'Esc');
END;

{---------------------------------------------------------------------------}

PROCEDURE plotfit2;   {Plots the final fit and the experimental data}
BEGIN
  setfillstyle(solidfill, 0);
  bar(0, 0, 460, maxy - 30);
  setcolor(11);
  rectangle(40, 0, 460, maxy - 30);

  settextjustify(righttext, centertext);
  FOR i := 0 TO 4 DO
  BEGIN
    j := trunc(i * 100 / scale);
    ny := yl(j);
    line(35, ny, 40, ny);
    str(j, s1);
    outtextxy(32, ny, s1);
  END;

  FOR i := chini TO chfin DO            {Plot the data}
  BEGIN
    nx := xl1(i);
    ny := yl(y[i]);
    putpixel(nx, ny, 15);
  END;

  manual2;

  FOR i := chini TO chfin DO
  BEGIN
    nx := xl1(i);
    xx := sqr((i - aa[1]) / aa[2]);
    sum := aa[7];
    IF i > (aa[1] - 2.0 * aa[2]) THEN
      sum := (aa[7] - aa[8]) * (aa[4] + 2.0 * aa[5] - i) /
             (aa[4] - aa[1] + 2.0 * aa[2] + 2.0 * aa[5]) + aa[8];
    IF i > (aa[4] + 2.0 * aa[5]) THEN
      sum := aa[8];
    sum1 := (aa[3] - aa[7]) * exp(-xx);
    IF xx < 10.0 THEN
      sum := sum + (aa[3] - aa[7]) * exp(-xx);
    xx := sqr((i - aa[4]) / aa[5]);
    sum2 := (aa[6] - aa[7]) * exp(-xx);
    IF xx < 10.0 THEN
      sum := sum + (aa[6] - aa[7]) * exp(-xx);
    fity[i] := sum;
    fitp1[i] := sum1;
    fitp2[i] := sum2;
    ny := yl(trunc(sum));
    putpixel(nx, ny, 14);
    IF (twocurve) THEN
    BEGIN
      ny := yl(trunc(sum1));
      putpixel(nx, ny, 12);
      ny := yl(trunc(sum2));
      putpixel(nx, ny, 9);
    END;
  END;

  str(chitot:10:3, s1);
  str(chipt:10:3, s2);
  str(area1:10:3, s3);
  str(area2:10:3, s4);

  setcolor(14);
  settextjustify(lefttext, toptext);
  bar(555, 148, maxx - 1, 157);
  outtextxy(555, 148, s3);
  bar(555, 158, maxx - 1, 167);
  outtextxy(555, 158, s4);
  bar(555, 188, maxx - 1, 197);
  outtextxy(555, 188, s1);
  bar(555, 208, maxx - 1, 217);
  outtextxy(555, 208, s2);
END;

{---------------------------------------------------------------------------}

PROCEDURE writeout2;
BEGIN
  str(aa[1]:10:3, s1);
  str(aa[2]:10:3, s2);
  str(aa[3]:10:3, s3);
  str(aa[4]:10:3, s4);
  str(aa[5]:10:3, s5);
  str(aa[6]:10:3, s6);
  str(aa[7]:10:3, s7);
  str(aa[8]:10:3, s8);
  str(del:10:3, s9);

  setfillstyle(solidfill, 0);
  settextjustify(lefttext, toptext);
  setcolor(10);
  bar(555, 8, maxx - 1, 17);
  outtextxy(555, 8, s1);
  bar(555, 18, maxx - 1, 27);
  outtextxy(555, 18, s2);
  bar(555, 28, maxx - 1, 37);
  outtextxy(555, 28, s3);
  bar(555, 48, maxx - 1, 57);
  outtextxy(555, 48, s4);
  bar(555, 58, maxx - 1, 67);
  outtextxy(555, 58, s5);
  bar(555, 68, maxx - 1, 77);
  outtextxy(555, 68, s6);
  bar(555, 88, maxx - 1, 97);
  outtextxy(555, 88, s7);
  bar(555, 98, maxx - 1, 107);
  outtextxy(555, 98, s8);
  bar(555, 118, maxx - 1, 127);
  outtextxy(555, 118,s9);
END;

{---------------------------------------------------------------------------}

PROCEDURE Manfit2;
BEGIN
  write('input the initial channel for the fit? ');
  readln(chini);
  write('input the final channel for the fit? ');
  readln(chfin);
  writeln;
  writeln('input the starting values : peak1 sig1 h1 peak2 sig2 h2 b1 b2');
  writeln('     (hit <space> or <return> between each)');
  readln(aa[1], aa[2], aa[3], aa[4], aa[5], aa[6], aa[7], aa[8]);

  chrange := chfin - chini;
  del := 10.0;
  twocurve := false;
  stop := false;

  startgraf;
  borders2;
  plotfit2;                 {Plot the fit, data, and parameters}
  writeout2;

  REPEAT
    IF keypressed THEN
    BEGIN
      ch := readkey;
      IF Ord(ch) = 27 THEN
        stop := true
      ELSE IF Ord(ch) = 13 THEN
      BEGIN
        automatic2;
        plotfit2;
      END
      ELSE
        CASE ch OF
          'L','l' : BEGIN
                      del := del / 10.0;
                      IF del < 0.001 THEN
                        del := 100.0;
                    END;
          'R','r' : plotfit2;
          'a','A' : aa[1] := aa[1] + del;
          'Z','z' : aa[1] := aa[1] - del;
          'S','s' : aa[2] := aa[2] + del;
          'X','x' : aa[2] := aa[2] - del;
          'D','d' : aa[3] := aa[3] + del;
          'C','c' : aa[3] := aa[3] - del;
          'F','f' : aa[4] := aa[4] + del;
          'V','v' : aa[4] := aa[4] - del;
          'G','g' : aa[5] := aa[5] + del;
          'B','b' : aa[5] := aa[5] - del;
          'H','h' : aa[6] := aa[6] + del;
          'N','n' : aa[6] := aa[6] - del;
          'J','j' : aa[7] := aa[7] + del;
          'm','M' : aa[7] := aa[7] - del;
          'K','k' : aa[8] := aa[8] + del;
          ',','<' : aa[8] := aa[8] - del;
          't','T' : BEGIN
                      IF (twocurve) THEN
                        twocurve := false
                      ELSE
                        twocurve := true;
                      plotfit2;
                    END;
          '[','{' : BEGIN
                      scale := scale / 2.0;
                      plotfit2;
                    END;
          ']','}' : BEGIN
                      scale := scale * 2.0;
                      plotfit2;
                    END;
          '1','!' : begin
                     autoerror2(1);
                     plotfit2;
                    end;
          '2','@' : begin
                     autoerror2(2);
                     plotfit2;
                    end;
          '3','#' : begin
                     autoerror2(3);
                     plotfit2;
                    end;
          '4','$' : begin
                     autoerror2(4);
                     plotfit2;
                    end;
          '5','%' : begin
                     autoerror2(5);
                     plotfit2;
                    end;
          '6','^' : begin
                     autoerror2(6);
                     plotfit2;
                    end;
        END;  {CASE}
      writeout2;
    END;
  UNTIL stop;

  closegraph;
END; {manfit}

Procedure errorupdate;
begin
  setcolor(10);
  rectangle(470, 0, maxx, 90);
  settextjustify(lefttext,toptext);
  outtextxy(475, 8, 'peak   = ');
  outtextxy(475, 18, 'sigma  = ');
  outtextxy(475, 28, 'height = ');
  outtextxy(475, 48, 'b1     = ');
  outtextxy(475, 58, 'b2     = ');
  outtextxy(475, 78, 'del    = ');

  setcolor(14);
  rectangle(470, 110, maxx, 180);
  outtextxy(475, 118, 'area   = ');
  outtextxy(475, 138, 'total chisq =');
  outtextxy(475, 158, 'chi/data point =');

end;

Procedure ploterror(ncol:integer);

begin
nx:=260 + trunc((aa[2]-a20)*420/delx);
ny:=trunc((maxy-30)/2)-trunc((aa[3]-a30)*(maxy-30)/dely);
putpixel(nx,ny,ncol);
end;


Procedure errorhs;
var
 a2min, a2max, a3min, a3max : real;
 cc :integer;

begin
 startgraf;
 cleardevice;
 rectangle(40,0,460,maxy-30);
 a20:=aa[2];
 a30:=aa[3];
 a2min:=aa[2]-2;
 a2max:=aa[2]+2;
 a3min:=aa[3]-3.0*sqrt(aa[3]);
 a3max:=aa[3]+3.0*sqrt(aa[3]);
 delx:=4.0;
 dely:=6.0*sqrt(aa[3]);
 ploterror(15);
 errorupdate;
 writeout;
 stop:=false;
 repeat
  if keypressed then
   begin
    cc:=ord(readkey);
     case cc of
      27 : stop:=true;
      72 : begin
            ploterror(0);
            aa[3]:=aa[3]-sqrt(aa[3])*0.01;
            ploterror(15);
           end;
      80 : begin
            ploterror(0);
            aa[3]:=aa[3]+sqrt(aa[3])*0.01;
            ploterror(15);
           end;
      75 : begin
            Ploterror(0);
            aa[2]:=aa[2]-0.01;
            ploterror(15);
           end;
      77 : begin
            Ploterror(0);
            aa[2]:=aa[2]+0.01;
            ploterror(15);
           end;

      end; {case}
    writeout;
    manual;
    str(chitot:10:3, s1);
    str(chipt:10:3, s2);
    str(area1:10:3, s3);

    setcolor(14);
    settextjustify(lefttext, toptext);
    bar(555, 118, maxx - 1, 127);
    outtextxy(555, 118, s3);
    bar(555, 148, maxx - 1, 157);
    outtextxy(555, 148, s1);
    bar(555, 168, maxx - 1, 177);
    outtextxy(555, 168, s2);

   end;
until stop;
closegraph;
end; {errorhs}

{---------------------------------------------------------------------------}

PROCEDURE announcement;
BEGIN
  str(limit:4, s1);
  writeln('current file open : ' + fname + '  (' + s1 + ' channels)');
  writeln('Total number of counts = ',ctot);
  writeln('Total number of counts*channel = ',chctot);
  writeln;
  writeln('press a key to continue: ');
  writeln('   P is for plotting the data');
  writeln('   1 is for single curvefit');
  writeln('   2 is for double curvefit');
  writeln('   E is for error analysis on h and sigma');
  writeln('   O is for opening different file');
  writeln('   W is for writing a spreadsheet file of the fit');
  writeln('   Q is for quit');
  writeln;
  write('? ');
END;

{---------------------------------------------------------------------------}

PROCEDURE title;
BEGIN
  ClrScr;
  writeln;
  writeln;
  writeln('Welcome to the MCA curvefitting program');
  writeln('    written by Byron Curry and Peter Siegel');
  writeln('         modified by TK');
  writeln;
END;

{---------------------------------------------------------------------------}

PROCEDURE readfile;
BEGIN
  write('name of data file? ');
  readln(fname);
  for i := 1 to Length(fname) do
  fname[i] := UpCase(fname[i]);
  assign(inf, fname);
  reset(inf);

  i := -1;
  WHILE NOT EOF(inf) DO
  BEGIN
    i := i + 1;
    x[i]:=i;
  {  read(inf,x[i]);  }
    read(inf, y[i]);
  END;
  limit := i;
  ctot:=0;
  chctot:=0;
  for j:=1 to i do
    ctot:=ctot+y[j];
  for j:=1 to i do
    chctot:=chctot+j*y[j];

  Close(inf);
  writeln;
  scale := 1.0;
  title;
END;

Procedure writefile;
begin
  write('name of data file? ');
  readln(fname);
  assign(inf, fname);
  rewrite(inf);
  for i:=chini to chfin do
  writeln(inf,i,chr(9),y[i],chr(9),fity[i]:10:4,fitp1[i]:10:4,fitp2[i]:10:4);
  close(inf);
end;

{---------------------------------------------------------------------------}

{main program}
BEGIN
  title;
  readfile;
  running := true;
  announcement;

  WHILE (running) DO
  BEGIN
    IF keypressed THEN
    BEGIN
      ch := readkey;
      IF Ord(ch) = 27 THEN
        running := false
      ELSE
        CASE ch OF
          'P','p' : BEGIN
                      plotdat;
                      title;
                      announcement;
                    END;
          '1','!' : BEGIN
                      writeln('Single curve fit');
                      writeln;
                      manfit;
                      title;
                      announcement;
                    END;
          '2','@' : BEGIN
                      writeln('Double curve fit');
                      writeln;
                      manfit2;
                      title;
                      announcement;
                    END;
          'e','E' : begin
                     writeln('Error Analysis for h and sigma');
                     writeln;
                     errorhs;
                     title;
                     announcement;
                    end;
          'O','o' : BEGIN
                      writeln('Open file');
                      writeln;
                      readfile;
                      announcement;
                    END;
          'W','w' : begin
                     writeln('write file');
                     writeln;
                     writefile;
                     announcement;
                    end;
          'Q','q' : running := false;
        END;  {CASE}
    END;
  END;

END.   {main program}

{---------------------------------------------------------------------------}