{$N+}
Program mcaenergycal;

uses
 crt, dos, graph;

const
 np = 5;

type
 integerarraynp = array[1..np] of integer;
 Realarraynp = array[1..np] of extended;
 realarraynpbynp = array[1..np,1..np] of extended;
 data   =  array[1..10] of extended;

var
 ip,n,i,j,k,ns,nn,j1,i1,nunk,jj, ny, ni : integer;
 ch, energy, d, offset,dele,delch,x,dele2 : extended;
 indx,nsi : integerarraynp;
 col,estrd,fit,delfit,chcal,errcal : realarraynp;
 estrds,chcals,errcals,yest :realarraynp;
 dela,a,y,check,asave: realarraynpbynp;
 running : boolean;
 graphdriver, graphmode, error :integer;
 maxx,maxy,yp   : integer;
 cho   : char;
 chunk, errunk : data;
 inf, outf      :text;
 fname,s1,s2,s3,s4,s5          :string;
 sumx,sumy,sumxy,sumx2,sumx_sq,slope,intercept,oset   :real;
 xscale,yscale,chmax   :real;

function pow(var x:extended;var j:integer):extended;
begin
pow:=exp(j*ln(x));
end;

Procedure lubksb(var a: realarraynpbynp; n:integer;
                 var indx: integerarraynp; var b: realarraynp);
var
 j,ip,ii,i :integer;
 sum : real;

begin
 ii:=0;
 for i:=1 to n do
  begin
   ip:=indx[i];
   sum:=b[ip];
   b[ip]:=b[i];
   if ii<>0 then
    for j:=ii to i-1 do
     sum := sum-a[i,j]*b[j]
    else if sum <> 0.0 then
     ii :=i;
    b[i]:=sum;
   end;
  for i:=n downto 1 do
   begin
    sum:=b[i];
    for j:= i+1 to n do
     sum := sum-a[i,j]*b[j];
    b[i]:=sum/a[i,i];
   end;
end;

Procedure ludcmp(var a: realarraynpbynp; n:integer;
                 var indx : integerarraynp; var d:extended);
const
 tiny = 1.0e-20;

var
 k,j,imax,i :integer;
 sum,dum,big :real;
 vv :^realarraynp;

begin
 new(vv);
 d:=1.0;
 for i:=1 to n do begin
   big :=0.0;
   for j:=1 to n do
    if abs(a[i,j])>big then big :=abs(a[i,j]);
    if big=0.0 then begin
      writeln('pause in ludcmp - singular matrix');
      readln
     end;
    vv^[i]:=1.0/big;
 end;
  for j:=1 to n do begin
    for i:=1 to j-1 do begin
      sum :=a[i,j];
      for k:=1 to i-1 do
       sum :=sum-a[i,k]*a[k,j];
      a[i,j]:=sum
     end;
    big:=0.0;
    for i:=j to n do begin
      sum:=a[i,j];
       for k:=1 to j-1 do
        sum := sum - a[i,k]*a[k,j];
        a[i,j]:=sum;
    dum := vv^[i]*abs(sum);
    if dum >= big then begin
      big :=dum;
      imax:=i;
     end;
   end;
  if j<> imax then begin
   for k:=1 to n do begin
    dum := a[imax,k];
    a[imax,k] := a[j,k];
    a[j,k]:=dum;
   end;
  d:=-d;
  vv^[imax]:=vv^[j];
 end;
indx[j] := imax;
if a[j,j]=0.0 then a[j,j] := tiny;
if j<> n then begin
 dum:=1.0/a[j,j];
 for i:=j+1 to n do
  a[i,j]:=a[i,j]*dum;
 end;
 end;
 dispose(vv);
 end;

Procedure fitter;
begin
writeln;
writeln;
writeln(' number       channel number        energy  ');
for i:= 1 to nn do
 writeln('   ',i,'           ',chcals[i]:8:3,'           ',estrds[i]:8:3);
writeln;
writeln(' How many standards from the above list do you ');
writeln(' want to use?  ');
readln(ns);
writeln(' Which ones?  Enter the numbers in order');
for i:=1 to nn do
 nsi[i]:=0;
for i:=1 to ns do
 read(nsi[i]);
for i:=1 to ns do
 begin
  chcal[i]:=chcals[nsi[i]];
  estrd[i]:=estrds[nsi[i]];
 end;
for i:=1 to 3 do
 for j:=1 to 3 do
   a[i,j]:=0.0;

for i:= 1 to 3 do
 begin
   for j:=i to 3 do
    begin
     a[i,j]:=0.0;
      for k:=1 to ns do
       begin
       i1:=i+j-2;
       a[i,j]:=a[i,j] + pow(chcal[k],i1);
       end;
        asave[i,j]:=a[i,j];
        a[j,i]:=a[i,j];
        asave[j,i]:=a[j,i];
     end;
    end;
 n:=3;
 ludcmp(a,n,indx,d);
 for j:=1 to n do
  begin
   for i:=1 to n do col[i] :=0.0;
   col[j]:=1.0;
   lubksb(a,n,indx,col);
   for i:=1 to n do y[i,j] := col[i];
  end;

for i:=1 to 3 do
begin
 yest[i]:=0.0;
  for j:=1 to ns do
   begin
    i1:=i-1;
    yest[i]:=yest[i] + estrd[j]*pow(chcal[j],i1);
   end;
end;

for i:=1 to ns do
 fit[i]:=0.0;

for i:=1 to 3 do
 for k:=1 to 3 do
  fit[i]:=fit[i]+y[i,k]*yest[k];

for i:=1 to 3 do
 writeln(' Coefficient ', i,' is ',fit[i]);

graphdriver:=detect;
initgraph(graphdriver,graphmode,'');
maxx:=getmaxx;
maxy:=getmaxy;
yp:=340;
chmax:=chcals[1];
for i:=2 to nn do
 if (chcals[i]<chmax) then chmax:=chcals[i];
line(5,0,5,300);
line(5,300,maxx,300);
xscale:=0.001*(maxx-5);
yscale:=300.0/1500.0;
yp:=yp+10;
outtextxy(20,yp,'Standards:   Channel Number      Energy          Difference ');
ip:=1;
for jj:=1 to nn do
begin
  ch:=chcals[jj];
  energy:=0.0;
  dele:=0.0;
  dele2:=0.0;
  for i:=1 to 3 do
   begin
    i1:=i-1;
    energy := energy + fit[i]*pow(ch,i1);
   end;
  str(chcals[jj]:15:2,s1);
  str(energy:15:2,s2);
  str((estrds[jj]-energy):15:2,s3);
  if jj=nsi[ip] then
   begin
    setcolor(13);
    outtextxy(70,trunc(yp+jj*10),s1+'   '+s2+'    '+s3);
    circle(trunc(5+chcals[jj]*xscale),trunc(300-energy*yscale),3);
    ip:=ip+1;
   end
    else
    begin
      setcolor(10);
      outtextxy(70,trunc(yp+jj*10),s1+'   '+s2+'    '+s3);
      circle(trunc(5+chcals[jj]*xscale),trunc(300-energy*yscale),3);
    end;
 end;
yp:=yp+jj*10+20;

setcolor(14);
outtextxy(20,yp,'Unknowns:   Channel Number         Energy');
settextstyle(3,0,1);

for jj:=1 to nunk do
begin
  ch:=chunk[jj];
  energy:=0.0;
  dele:=0.0;
  dele2:=0.0;
  for i:=1 to 3 do
   begin
    i1:=i-1;
    energy := energy + fit[i]*pow(ch,i1);
   end;
  dele:=sqrt(dele);
  str(chunk[jj]:15:2,s1);
  str(energy:15:2,s2);
  outtextxy(60,yp+jj*15,s1+'   '+s2);
  circle(trunc(5+chunk[jj]*xscale),trunc(300-energy*yscale),3);
end;

 for ni:=1 to 1000 do
  begin
   energy:=0.0;
    for i:=1 to 3 do
     begin
      ch:=ni;
      i1:=i-1;
      energy := energy + fit[i]*pow(ch,i1);
     end;
   putpixel(trunc(5+(ni+offset)*xscale),trunc(300-energy*yscale),13);
  end;

 setcolor(15);
 settextstyle(0,0,0);
 outtextxy(10,maxy-10,' f: for a quadratic fit;  l: for linear fit; q: quit');
end;

Procedure linfit;
begin
writeln;
writeln(' number       channel number        energy  ');
for i:= 1 to nn do
 writeln('   ',i,'           ',chcals[i]:8:3,'           ',estrds[i]:8:3);
writeln;
writeln(' How many standards from the above list do you ');
writeln(' want to use?  ');
readln(ns);
writeln(' Which ones?  Enter the numbers in order');
for i:=1 to nn do
 nsi[i]:=0;
for i:=1 to ns do
 read(nsi[i]);
for i:=1 to ns do
 begin
  chcal[i]:=chcals[nsi[i]];
  estrd[i]:=estrds[nsi[i]];
  errcal[i]:=errcals[nsi[i]];
 end;
sumx:=0;
sumy:=0;
sumx2:=0;
sumxy:=0;
for i:=1 to ns do
 begin
  sumx:=sumx+chcal[i];
  sumy:=sumy+estrd[i];
  sumx2:=sumx2+chcal[i]*chcal[i];
  sumxy:=sumxy+chcal[i]*estrd[i];
 end;
sumx_sq:=sumx*sumx;
slope:=(ns*sumxy-sumx*sumy)/(ns*sumx2-sumx_sq);
intercept:=(sumy*sumx2-sumxy*sumx)/(ns*sumx2-sumx_sq);
oset:=-intercept/slope;
graphdriver:=detect;
initgraph(graphdriver,graphmode,'');
maxx:=getmaxx;
maxy:=getmaxy;
chmax:=chcals[1];
for i:=2 to nn do
 if (chcals[i]>chmax) then chmax:=chcals[i];
yp:=340;
line(5,0,5,300);
line(5,300,maxx,300);
xscale:=1.0*(maxx-5)/(chmax+10);
yscale:=300.0/1500.0;
str(slope:15:4,s1);
str(oset:15:4,s2);
outtextxy(20,yp,'The slope is ' + s1 + ' the offset is ' + s2);
ip:=1;
yp:=yp+10;
outtextxy(20,yp,'Standards:   Channel Number      Energy          Difference ');
for jj:=1 to nn do
 begin
  energy:=slope*chcals[jj]+intercept;
  str(chcals[jj]:15:2,s1);
  str(energy:15:2,s2);
  str((estrds[jj]-energy):15:2,s3);
  if jj=nsi[ip] then
   begin
    setcolor(13);
    outtextxy(70,trunc(yp+jj*10),s1+'   '+s2+'    '+s3);
    circle(trunc(5+chcals[jj]*xscale),trunc(300-energy*yscale),3);
    ip:=ip+1;
   end
    else
     begin
      setcolor(10);
      outtextxy(70,trunc(yp+jj*10),s1+'   '+s2+'    '+s3);
      circle(trunc(5+chcals[jj]*xscale),trunc(300-energy*yscale),3);
     end;
 end;
yp:=yp+jj*10+20;

setcolor(14);
outtextxy(20,yp,'Unknowns:   Channel Number         Energy');
settextstyle(3,0,1);

for i:=1 to nunk do
 begin
  energy:=slope*chunk[i]+intercept;
  str(chunk[i]:15:2,s1);
  str(energy:15:2,s2);
  outtextxy(60,yp+i*15,s1+'   '+s2);
  circle(trunc(5+chunk[i]*xscale),trunc(300-energy*yscale),3);
 end;
 for i:=0 to 1000 do
  begin
   ny:=trunc(slope*i+intercept);
   putpixel(trunc(5+i*xscale),trunc(300-ny*yscale),13);
  end;
 setcolor(15);
 settextstyle(0,0,0);
 outtextxy(10,maxy-10,' f: for a quadratic fit;  l: for linear fit; q: quit');

end;

Procedure screenin;
begin
writeln(' input the number of calibration energies (limit=5)');
readln(nn);
for i:=1 to nn do
 begin
  writeln(' input the calibration energy and channel No. for energy ',i);
  readln(estrds[i],chcals[i]);
  writeln;
 end;
writeln(' input the number of unknown energies to be determined (limit=10)');
readln(nunk);
for i:=1 to nunk do
 begin
  writeln(' input the channel No. for unknown energy ',i);
  readln(chunk[i]);
 end;
writeln;
writeln('filename to save file:');
readln(fname);
assign(outf,fname);
rewrite(outf);
writeln(outf,nn,'  ',nunk);
for i:=1 to nn do
  writeln(outf,estrds[i],chcals[i]);
for i:=1 to nunk do
  writeln(outf,chunk[i]);
close(outf);
end;

Procedure readfl;
 begin
 writeln;
 writeln('file to be read in:');
 readln(fname);
 assign(inf,fname);
 reset(inf);
readln(inf,nn,nunk);
for i:=1 to nn do
  readln(inf,estrds[i],chcals[i]);
for i:=1 to nunk do
  readln(inf,chunk[i]);
close(inf);
writeln;
writeln('Calibration energys, channel number');
for i:=1 to nn do
 writeln(estrds[i]:15:4,'  ',chcals[i]:15:4);
writeln;
writeln('Unknown channel numbers');
for i:=1 to nunk do
writeln(chunk[i]:15:4);
writeln;
 end;

begin
writeln(' Welcome to the MCA calibration program');
writeln;
writeln('      This program determines the coefficients for a polynomial');
writeln(' expansion of the E(C) relationship: the energy E as a function');
writeln(' of channel number C');
writeln;
running:=true;
writeln(' press s to input data from the screen');
writeln(' press r to read data from a file');
while (running) do
 begin
  if keypressed then
   begin
   cho:=readkey;
    case cho of
     's' : begin
          screenin;
          running:=false;
         end;
     'r' : begin
          readfl;
          running:=false;
         end;
     end;{case}
    end;
  end;
running:=true;
writeln;
writeln(' press f to do a least-square quadratic fit ');
writeln(' press l to do a least-square linear fit ');
writeln(' press q to exit the program');
while (running) do
 begin
  if keypressed then
   begin
    cho:=readkey;
     case cho of
      'Q', 'q' : begin
                  closegraph;
                  running :=false;
                 end;
      'F', 'f' : begin
                 closegraph;
                 fitter;
                 end;
      'L','l' :  begin
                 closegraph;
                 linfit;
                 end;
     end;{case}
    end;
   end;
end.