program FindBulge4;  {Includes bitmap utility
        Looks for bulges in the distance curve vs force,
        which indicate an unfolding transition.To use this program,
        the data columns must include CycleCount, Y_force, A_dist-Y,
        B_dist-Y, and Status.  Status reporting must be modified in the Host
        ltposix.c to read:
          data->quicktime.recent[recentIndex][kChannelTweezStatus] =
		(float) (data->instrument.runMode     * 1000
		+ (data->quicktime.switches /2)       * 10
		+ data->instrument.pausing            * 4
               	+ data->instrument.stiffTestHighA     * 2
		+ data->instrument.movingUpward       * 1  ); }
uses crt,graph;
const FXsize = 320;  {number of force gradations in display}
      maxForce = 32;
      minForce = 0;
      scale = 3;  {controls display height of bulge}
      recent = 4;  {smoothing length for rip detection}
      nRips = 10; {number of rips to remember}
var  ch:char;
     peakFound,valleyFound,showingMenu,atTheEnd,readyToGo:boolean;
     i,j,k,L,n,f,u,column,index,xMax, yMax, xMid, yMid,CycleCntCol,YforceCol,
     AdistYCol,BdistYCol,StatusCol,IntStatus,oldIntStatus,rowCount,switches,
     pullNum,hue,code: integer;
     peakCnt,noiseCnt: longInt;
     fileName,dataLine,dataWord,titleLine,message:string;
     inFile, outFile, bFile:text;
     APsdY,BPsdY,Yforce,AdistY,BdistY,CycleCnt,Status,avgF,dFmax,dFmin,
     unfForce,refForce,HiLim,LoLim,mean,delta,cross,x,w,noise,drift: real;
     pullRip: array[0..nRips+1, 1..2] of real;  {1=dF, 2=yForce}
     relaxRip: array[0..nRips+1, 1..2] of real;  {1=dF, 2=yForce}
     Position: array[0..FXsize,1..2] of real;
     avSum: array[0..FXsize,1..2] of integer; {number of points to average}
     pastF: array[1..recent] of real;

{$I Bitmap.pas}

procedure openInputFile;
begin
clrscr;
writeln;
writeln('What text file to open?  ".TXT" added automatically');
readln(filename);
assign(infile,filename+'.txt');
reset(infile);
readln(infile,titleLine);
end;

Procedure detectColumns; {finds needed columns}
begin
repeat readln(infile,dataLine) until length(dataLine) > 40;{skip further intro}
i:=0; column:=1; L:=length(dataLine); dataWord:=''; status:=0;
repeat
      if dataLine[i]=chr(9) then  begin inc(column); dataWord:=''; end;
      repeat inc(i); dataWord:=dataWord + dataLine[i];
      until (i=L) or (dataLine[i+1]=chr(9));
      if dataWord='CycleCount' then CycleCntCol := column;
      if dataWord='Y_force' then YforceCol := column;
      if dataWord='A_dist-Y' then AdistYCol := column;
      if dataWord='B_dist-Y' then BdistYCol := column;
      if dataWord='Status' then StatusCol := column;
      inc(i);
until (i>L) or keypressed;    {done parsing row}
writeln;
if  CycleCntCol = 0 then writeln('CycleCount column missing, dataWord=',dataWord,'  L=',L);
if  YforceCol = 0 then writeln('Y_force column missing');
if  AdistYCol = 0 then writeln('A_dist-Y column missing');
if  BdistYCol = 0 then writeln('B_dist-Y column missing');
if  StatusCol = 0 then writeln('Status column missing');
if  (CycleCntCol * YforceCol * AdistYCol *  BdistYCol * StatusCol = 0)
    then begin
    writeln('program aborted.....press any key to exit');
    repeat until keypressed;
    halt;
    end;
end;

procedure openOutputFile;
begin
assign(outFile,filename+'.fit');
rewrite(outFile);
writeln(outFile,titleLine);
writeln(outFile,'CycCnt':10,chr(9),'pull',chr(9),'unfF',chr(9),'refF',chr(9),
   'loL',chr(9),'loF',chr(9),'hiL',chr(9),'hiF',chr(9),'sw',char(9),'noise',
   chr(9),'drift',chr(9),'code');
end;

procedure drawAxes;
begin
setColor(white);
{setColor(switches div 2 + 8);}
line(0,yMid,xMax,yMid);
for i:= 1 to 32 do line(i*20,yMid-5,i*20,yMid+5);
for i:= 1 to 6 do line(i*100,yMid-10,i*100,yMid+10);
for i:= 1 to 3 do line(i*200,yMid-15,i*200,yMid+15);
outTextXY(180,yMid+25,'10pN');
outTextXY(380,yMid+25,'20pN');
outTextXY(580,yMid+25,'30pN');
for i:= 4 to 64 do putPixel(i*10,yMid-scale*10,white);
for i:= 4 to 64 do putPixel(i*10,yMid-scale*20,white);
outTextXY(0,yMid-scale*10-3,'10nm');
outTextXY(0,yMid-scale*20-3,'20nm');
end;

procedure showMenu;
begin
setColor(white);
Str(pullNum,message);
outTextXY(10,yMid+ 150,'pull number = '+message);
Str(peakCnt,message);
outTextXY(10,yMid+ 160,'cycleCount = '+message);
if atTheEnd then
   begin
   outTextXY(10,yMid+ 170,'lastLine= '+ dataLine);
   outTextXY(10,yMid+ 180,'END OF FILE.......push q');
   end
else begin
     Str(switches,message);
     outTextXY(10,yMid+ 170,'switches = '+message);
     Str(rowCount,message);
     outTextXY(10,yMid+ 180,'points recorded = '+message);
     message:='0';
     if noiseCnt >0 then Str(noise/noiseCnt*1000:0:0,message);
     outTextXY(10,yMid+ 190,'noise = '+message);
     outTextXY(xMid-100,yMid+60,'Ins - get next pull');
     outTextXY(xMid-100,yMid+80,'f,u = switch fold, unfold forces');
     outTextXY(xMid-100,yMid+100,'F,U - set new Fold, Unfold forces');
     outTextXY(xMid-100,yMid+120,'arrows - move brackets');
     outTextXY(xMid-100,yMid+140,'b - step bracketed fit');
     outTextXY(xMid-100,yMid+160,'spacebar - missing plateau');
     outTextXY(xMid-100,yMid+180,'z/c- problems left/right');
     end;
outTextXY(xMid-100,yMid+ 200,'r - reset input file');
outTextXY(xMid-100,yMid+ 220,'q - quit');
end;

procedure getPullRips;   {positive rips ordered by size}
var dF:real; i,j:integer;
begin
dF := yForce - avgF;
for i:= nRips downto 1 do for j:=1 to 2 do
    if dF > pullRip[i,1] then
       begin
       pullRip[i+1,j] := pullRip[i,j]; {push down old dF}
       if dF < pullRip[i-1,1] then
          begin
          pullRip[i,1]:= dF;  {fill in new dF}
          pullRip[i,2]:= pastF[recent]; {unfold force, we hope}
          end;
       end;
unfForce:= pullRip[1,2];
HiLim:= unfForce;
noise:=noise + dF*dF;
inc(noiseCnt);
end;

procedure getRelaxRips;   {negative rips ordered by size}
var dF:real; i,j:integer;
begin
dF := yForce - avgF;
for i:= nRips downto 1 do for j:=1 to 2 do
    if dF < relaxRip[i,1] then
       begin
       relaxRip[i+1,j] := relaxRip[i,j]; {push down old dF}
       if dF > relaxRip[i-1,1] then
          begin
          relaxRip[i,1]:= dF;  {fill in new dF}
          relaxRip[i,2]:= pastF[recent]; {unfold force, we hope}
          end;
       end;
refForce:= relaxRip[1,2];
loLim:= refForce;
noise:=noise + dF*dF;
inc(noiseCnt);
end;

procedure plotMolecule;
begin
moveTo(0,yMid);
hue := (switches div 2) + 2;
if (hue=0) or (hue=16) then hue:=8;
if (code=0) then hue :=1;
if code=2 then hue :=7;
if code=3 then hue :=6;
setColor(hue);
for index:=0 to FXsize do {plot difference between pull & relax positions}
    if (position[index,1] <> 0 ) and (position[index,2] <> 0) then
       lineTo(index*2,yMid-round(scale*(position[index,2]-position[index,1])));
setColor(yellow);
circle(2*round(-unfForce*10),yMid,5);
line  (2*round(-hiLim *10),yMid,2*round(-hiLim *10),yMid-50);
for i:=1 to nRips do circle(2*round(-pullRip[i,2]*10),yMid,2);
setColor(12);
for i:=1 to nRips do circle(2*round(-relaxRip[i,2]*10),yMid,2);
circle(2*round(-refForce *10),yMid,5);
line  (2*round(-loLim *10),yMid,2*round(-loLim *10),yMid-50);
for i:= 4 to 64 do putPixel(i*10,yMid,white);
end;

procedure getTheDrift;
begin
drift:=0.0;
for index:=40 to 80 do {for  4 pN to 8 pN}
    if (position[index,1] <> 0 ) and (position[index,2] <> 0) then
       drift:=drift+ position[index,2]-position[index,1];
end;

procedure fitRegion;   {linear fit to bulge between loLim and hiLim}
var i,low,high : integer;
begin
mean:=0.0; n:=0; cross:=0.0;
low:=round(-loLim*10); high:=round(-hiLim*10);
if high > low then
begin
for i:= low to high do
    if (position[i,1] <> 0 ) and (position[i,2] <> 0) then
       begin inc(n); mean:=mean+position[i,2]-position[i,1] ; end;
if n>0 then mean := mean / n;
setColor(8); {gray}
line  (2*round(-loLim *10),yMid-round(scale*(mean-delta)),
       2*round(-hiLim *10),yMid-round(scale*(mean+delta)));
Str(mean-delta:1:2,message);
outTextXY(20*round(-loLim)-40,yMid-60,message);
Str(mean+delta:1:2,message);
outTextXY(20*round(-hiLim)-5,yMid-60,message);
w := (high - low) / 2;
for i:= low to high do
    begin
    x := i - (high+low)/2;
    if (position[i,1] <> 0 ) and (position[i,2] <> 0) and (w<>0) then
       cross := cross+(position[i,2]-position[i,1]-mean-(x/w)*delta) * x;
    end;
if w<>0 then cross:=cross / w;
end; {if high > low}
Str(mean:1:2,message);
outTextXY(10*round(-loLim-hiLim)-10,yMid+80,'mean='+message);
Str(cross:1:2,message);
outTextXY(10*round(-loLim-hiLim)-10,yMid+120,'cross='+message);
Str(n,message);
outTextXY(10*round(-loLim-hiLim)-10,yMid+160,'n='+message);
Str(delta:1:3,message);
outTextXY(10*round(-loLim-hiLim)-10,yMid+200,'delta='+message);
end;

procedure checkAdjustments;
begin {check that user has pressed certain keys before printing to file}
readyToGo := true;
if HiLim = unfForce then readyToGo := false; {user not moved brackets}
if loLim = refForce then readyToGo := false;
if code <> 1 then readyToGo := true;
if ord(ch) <> 82 then readyToGo := false;  {Ins - key}
if pullNum = 0 then readyToGo :=true;
if ch ='q' then readyToGo := true;   {last word wins here}
end;

procedure report;
begin

Str(f,message);
outTextXY(20,20,'f='+message);
Str(relaxRip[f,2]:1:2,message);
outTextXY(60,20,'force='+message);
Str(relaxRip[f,1]:1:3,message);
outTextXY(160,20,'dF='+message);
setColor(cyan);
Str(u,message);
outTextXY(20,40,'u='+message);
Str(pullRip[u,2]:1:1,message);
outTextXY(60,40,'force='+message);
end;

procedure makeBitmap;
begin
assign (outFile,filename+'.bmp');
rewrite (outFile);
writeBitmap;
close(outFile);
end;


begin  {main}
openInputFile;
openOutputFile;
detectColumns;
i:=detect;
initgraph(i,j,'C:\PASCAL\BGI');
xMax:=getMaxX; yMax:= getMaxY; xMid:=xMax div 2; yMid:=yMax div 2;
repeat  {for many stretch-relax pairs}
    for index:=0 to FXsize do for j:= 1 to 2 do
        begin avSum[index,j]:=0; position[index,j]:=0.0; end;
    for i:=1 to nRips+1 do for j:= 1 to 2 do
        begin pullRip[i,j]:=0.0; relaxRip[i,j]:=0.0; end;
    pullRip[0,1]:=999; relaxRip[0,1]:=-999;
    rowCount:=0;
    peakFound:=false; valleyFound:=false;
    clearDevice;
    drawAxes;
    moveTo(0,yMid);
    dFmax:=0; dFmin:=0;
    delta:=0;
    noise:=0.0; noiseCnt:=0;
    code:= 1; {default indicates normal unfold / refold}
    f:=1; u:=1;
    repeat  {for all data rows within a stretch-relax pair}
       repeat readln(infile,dataLine)
       until (length(dataLine)>20){skips "Skipped"} or EOF(infile) or keyPressed;
       atTheEnd := eof(infile);
       i:=0; column:=1; L:=length(dataLine); dataWord:=''; status:=0;
          repeat     {read data columns}
          if dataLine[i]=chr(9) then  begin inc(column); dataWord:=''; end;
             repeat
             inc(i);
             dataWord:=dataWord + dataLine[i];
             until (i=L) or (dataLine[i+1]=chr(9));
          if column = CycleCntCol then val(dataword,CycleCnt,j);
          if column = YforceCol then val(dataword,Yforce,j);
          if column = AdistYCol then val(dataword,AdistY,j);
          if column = BdistYCol then val(dataword,BdistY,j);
          if column = statusCol then val(dataword,Status,j);
          inc(i);
          until (i>L) or keypressed;    {done parsing row}
       inc(rowCount);
       IntStatus := (round(Status)) mod 1000; {convert real into integer}
       switches := IntStatus div 10;
       if ((IntStatus-oldIntStatus)=-1) and (Yforce<-20) then peakFound:=true;
       if ((IntStatus-oldIntStatus)=1) and (Yforce<-1) then valleyFound:=true;
          {force checks rule out oscillations caused by broken molecule}
       oldIntStatus :=IntStatus;
       for k:= recent downto 2 do pastF[k]:=pastF[k-1];
       pastF[1] := yForce;
       avgF:=0; for k:=1 to recent do avgF:=avgF+pastF[k]; avgF:=avgF/recent;
       if yForce<-5 then if odd(IntStatus) then getPullRips else getRelaxRips;
       if peakFound then peakCnt := round(cycleCnt);
       peakFound:=false;
       index:=round(-Yforce*10);
       if index>FXsize then index:=FXsize;  if index<0 then index:=0;
       if odd(IntStatus) then j:=1 else j:=2; {pull column vs relax}
       position[index,j]:= position[index,j] + (AdistY + BdistY) / 2;
       inc(avSum[index,j]);
    until valleyFound or eof(infile);
    for index:=0 to FXsize do for j:=1 to 2 do
       if avSum[index,j]>0 then {average the positions in each force bin}
       position[index,j] :=position[index,j] / avSum[index,j];
    plotMolecule;
    showMenu;

    repeat
    ch:=readkey;
    case ch of
    'r': begin reset(infile); pullNum := -1; atTheEnd:=false; end;
    'b': begin clearDevice; plotMolecule; fitRegion;
                            if n>0 then delta:=delta+cross/n;  end;
    ' ': begin code:=0; plotMolecule; end;
    'z': begin if code<>1 then code:=1 else code:=2;plotMolecule;end;
    'c': begin if code<>1 then code:=1 else code:=3;plotMolecule;end;
    'u':begin inc(u); if u>nRips then u:=1; unfForce:=pullRip[u,2];
          hiLim:= unfForce; clearDevice; plotMolecule; report; showMenu;end;
    'U':begin unfForce:=hiLim;clearDevice; plotMolecule; report; showMenu;end;
    'f':begin inc(f); if f>nRips then f:=1; refForce:=relaxRip[f,2];
          loLim:= refForce; clearDevice; plotMolecule; report; showMenu;end;
    'F':begin refForce:=loLim; clearDevice; plotMolecule; report;showMenu;end;
     'B': makeBitmap;
 chr(0): begin ch:=readkey;
         case ord(ch) of
{L Arrow}   75: begin clearDevice; loLim:=loLim-0.1; plotMolecule; ShowMenu; end;
{R Arrow}   77: begin clearDevice; hiLim:=hiLim+0.1; plotMolecule; ShowMenu; end;
{up arrow}  72: begin clearDevice; hiLim:=hiLim-0.2; plotMolecule; ShowMenu; end;
{dn arrow}  80: begin clearDevice; loLim:=loLim+0.1; plotMolecule; ShowMenu; end;
{insert}    82: begin end; {insert key will not self repeat}
         end;{inner case}
         end;{second readkey}
    end; {outer case}
    checkAdjustments;
    until readyToGo;
    if (pullNum>0) and (not atTheEnd) then
       begin
       repeat
         clearDevice; plotMolecule;
         fitRegion;
         if n>0 then delta:=delta+cross/n;
       until abs(cross) < 0.0001;
       noise := 1000 * noise / noiseCnt;
       getTheDrift;
       writeln(outFile,peakCnt:10,chr(9),pullNum,chr(9),-unfForce:1:1,chr(9),
        -refForce:1:1, chr(9), (mean-delta):1:2, chr(9), -loLim:1:1, chr(9),
        (mean+delta):1:2, chr(9), -hiLim:1:1, chr(9), switches,chr(9),
        noise:0:0, chr(9), drift:0:0, chr(9),code);
       end;
    inc(pullNum);
until ch = 'q';
writeln(outFile,'END');
close(infile); close(outFile);
end.