
function fretG40P
%@copyright Taekjip Ha, Michael Schlierf 1999-2010 new features are added




% warning off MATLAB:divideByZero
% pack;
close all;clear all;

folder='/Users/Michael/Documents/Urbana/DataPostdoc/G40P/';
%folder= input('enter the location of fret trace file \n','s');
cd(folder);

%fname=input('enter the name of the trace file without the suffix \n','s');
fname='hel1';
acceptorcheckstart=60;  %% statistics check survival of acceptor after xx seconds; adjust accoridng to data set
% define time unit
%timeunit=input('enter the integration time used in the experiment in
%seconds \n');
timeunit=0.03;
intervals=[];
% minimum intensity
minTotal=40;
% read data


fid=fopen([fname '.traces'],'r');
len=fread(fid,1,'int32');
Ntraces=fread(fid,1,'int16')
mal=[num2str(Ntraces/2) ' molecules ' num2str(len) ' data points each ' folder ' ' fname];
disp(mal);
raw=fread(fid,Ntraces*len,'int16');
%raw=abs(raw); %
fclose(fid);
% convert into donor and acceptor traces
index=(1:Ntraces*len);
Data=zeros(Ntraces,len);
donor=zeros(Ntraces/2,len);
acceptor=zeros(Ntraces/2,len);
Data(index)=raw(index);



for i=1:(Ntraces/2),
    donor(i,:)=Data(i*2-1,:);
    acceptor(i,:)=Data(i*2,:);
    dbackground=mean(donor(i,(end-40):end-10));   %standard 0 - MS set here to end-40:end-10
    
    %  30 data points at the end of the measurement. This is not necessarily the background but a precorrection 
    
    abackground=mean(acceptor(i,(end-40):end-10));   %standard 0 - MS set here to end-40:end-10 Important - this should not contain the acceptor test with direct acceptor excitation
 
   % donor(i,:)=donor(i,:)-dbackground;
   % acceptor(i,:)=acceptor(i,:)-abackground;
end

count=0;
for i=1:(Ntraces/2),   % MS statistics acceptor survival at the end
    
    if mean(acceptor(i,ceil(acceptorcheckstart/timeunit):floor((acceptorcheckstart+0.3)/timeunit))) > 40
        count=count+1;
    end
    if i==(Ntraces/2)
        disp(['number of traces with acceptor survival: ' num2str(count) ' => ' num2str(count/(Ntraces/2)*100) ' %']);
        disp(['acceptor test interval= ' num2str(acceptorcheckstart) ' s +0.3 s' ]);
    end
end

donor_transp=donor';
acceptor_transp=acceptor';
donor_linear=donor_transp(:);
acceptor_linear=acceptor_transp(:);
validfretE=[];
med_d=median(donor_linear);
med_a=median(acceptor_linear);

AcceptorTreshold=med_a*0.0;%=75.;
DonorTreshold=med_d*0.0;%=75.;

donor_linear=donor_linear-DonorTreshold;
acceptor_linear=acceptor_linear-AcceptorTreshold;
total_linear=donor_linear+acceptor_linear;
median(donor_linear);
median(acceptor_linear);
fretE_linear=1./(1+donor_linear./acceptor_linear);
%median(fretE_linear);

maxlist=max([max(donor,[],2) max(acceptor,[],2)],[],2);
dummy=hist(maxlist);
dummy=dummy/sum(dummy);
i=0;
total=0;
while total<0.8
    i=i+1;
    total=sum(dummy(1:i));
end
SignalPlotRangeY=max(maxlist)*i/10+min(maxlist);

validfretE=fretE_linear(find(fretE_linear>-2 &fretE_linear<2&(donor_linear>2*DonorTreshold | acceptor_linear>2*AcceptorTreshold) &total_linear>200&total_linear<1000));
outname=['validfretE_' fname '.dat'];
save(outname,'validfretE','-ascii');
starttime=zeros(Ntraces/2,1);
endtime=zeros(Ntraces/2,1);
moleculenumber=zeros(Ntraces/2,1);
time=(0:(len-1))*timeunit;

%calculate, plot and save average traces
dAvg=sum(donor,1)/Ntraces*2;
aAvg=sum(acceptor,1)/Ntraces*2;

elevelAvg=aAvg./(aAvg+dAvg);
figure;
hdl1=gcf;
subplot(2,1,1);
plot(time,dAvg,'g',time,aAvg,'r');
title('Average donor and acceptor signal');
zoom on;
grid on;
subplot(2,1,2);
plot(time,elevelAvg,'b');
zoom on;
grid on;
avgOutput=[time' dAvg' aAvg' ]; avgFileName=[fname '_avg.dat'];

save(avgFileName,'avgOutput','-ascii');
% calculate E level from the first 10 points and plot histograms of E level
% and total intensity. Also save the same info.

j=0;

elevel=zeros(Ntraces/2,1);
total=elevel;

outname=[fname '_elevel_2p.dat'];

for i=1:(Ntraces/2),
    
    tempD=sum(donor(i,(3:3)),2);
    tempA=sum(acceptor(i,(3:3)),2);
    total(i)=(tempA+tempD)/1;
    elevel(i)=tempA/(tempA+tempD);
    if total(i) < minTotal
        elevel(i)=NaN;
        %total(i)=NaN;
    end
    
end

save(outname,'elevel','-ascii');
FRETbin=(-0.2:0.05:1.2);

figure;
hdl2=gcf;

subplot(2,2,1);
FRETbin=(-0.2:0.025:1.3);
hist(elevel,FRETbin);
temp=axis;
temp(1)=-0.2;
temp(2)=1.2;
axis(temp);
title('FRET histogram');

subplot(2,2,2);
hist(total,30);
title('Total intensity histogram');

subplot (2,2,3);
plot(elevel, total,'bd');
temp=axis;
temp(1)=-0.2;
temp(2)=1.2;
axis(temp);
title('intensity vs. FRET');


k=0;
result=zeros(Ntraces/2,3);
starttime=zeros(Ntraces/2,1);
endtime=zeros(Ntraces/2,1);
moleculenumber=zeros(Ntraces/2,1);

rk=0;

r_starttime=zeros(Ntraces/2,1);
r_endtime=zeros(Ntraces/2,1);
stallduration=zeros(Ntraces/2,1);
recoverytime=zeros(Ntraces/2,1);
sumE=[];

figure;
hdl3=gcf;

for i=1:(Ntraces/2),
    
    elevel(i)=-1;
    
end

i=0;

while i<((Ntraces/2) - 1),
    i = i + 1;
    backd(i) = 0.;% min(donor(i,:));
    backa(i) = 0.;%min(acceptor(i,:));
    donor(i,:) = donor(i,:) - backd(i);
    acceptor(i,:) = acceptor(i,:) - backa(i);
end

i=0;
iInc=+1;


while i < (Ntraces/2),
    
    i=i+iInc;
    
    if iInc==50  %jump only once
        iInc=1;
    end
    
    if (total(i) < 10000) && ( i <= (Ntraces/2))%&& (acceptor(i,ceil(acceptorcheckstart/timeunit)) >50) % set criteria for showing molecules & (acceptor(i,:)) MS
        %  total(i) < minTotal &
        if iInc~=1
            iInc=+1;
        end
                
        % Trace window
        figure(hdl3);
        ax(1)=subplot(2,1,1);   %MS
        hold off;
        
        plot(time,donor(i,:),'g-');hold on;%g-o,'MarkerFaceColor','g','MarkerSize',4);hold on;   % 4 is good marker size %MS
        plot(time,acceptor(i,:),'r-');%r-o,'MarkerFaceColor','r','MarkerSize',4);
        
        temp=axis;       
        temp(3)=-50;
        temp(4)=SignalPlotRangeY+500;
        hold on;
        plot(time,(donor(i,:)+acceptor(i,:))+700,'color',[0.7, 0.7, 0.7]);
        grid on;
        axis(temp);
        posBackSlash=findstr(folder,'/');
        firstline=[folder((posBackSlash(end)+1):end) ' Minimum Donor Signal:' num2str(min(backd(i)))  ' Minimum Acceptor Signal:' num2str(min(backa(i)))];
        secondline=[fname ' Molecule=' num2str(i)];
        mal=char(firstline);
        title(mal);
        
        zoom on;
        hold off;
        ax(2)=subplot(2,1,2);   %MS
        fretE=zeros(len,1);
        for m=1:len,
            if (acceptor(i,m)+donor(i,m))==0
                fretE(m)=NaN; % This is to avoid undefined fretE interfering with future analysis
                %elseif (acceptor(i,m)<2*AcceptorTreshold) & (donor(i,m)<2*DonorTreshold)
                % fretE(m)=NaN;
            elseif acceptor(i,m)+donor(i,m) <minTotal %(acceptor(i,m)+donor(i,m) < 200 &  acceptor(i,m)+donor(i,m) > 1000)% minTotal
                fretE(m)=NaN;
            else
                fretE(m)=(acceptor(i,m)+0)./(acceptor(i,m)+donor(i,m)+0);
            end
        end
        linkaxes(ax,'x'); %MS
                
        plot(time,fretE,'k-');%k-o, 'MarkerFaceColor',[0 0 0],'MarkerSize',4)%,[0 0 0]); %colorhange from 0.6 .6 1  %MS
        temp(3)=-0.1;
        temp(4)=1.2;
        axis(temp);
        mal=char(secondline);
        title(mal);
        grid on;
        zoom on;
        mal=['molecule #' num2str(i) '  acceptor avg ' num2str(floor(mean(acceptor(i,ceil(acceptorcheckstart/timeunit):floor((acceptorcheckstart+0.5)/timeunit)))/40))];
        ans=input(mal,'s');
        
        if isempty(ans)
            ans='null'; %just to avoid error messages
        end
        if ans=='z'
            disp(['avg int acceptor=' num2str(mean(acceptor(i,:)))])
        end
        
        if ans=='s'
            if exist('traces','dir')==7
            else
                mkdir('traces');
            end
            
            j=j+1;
            mN(j)=i;
            
            outname=['traces/' fname '-' 'tr' num2str(i) '.dat'];
            % output=[time' donor(i,:)' acceptor(i,:)'];
            output=[time' fretE acceptor(i,:)' donor(i,:)'];
            save(outname,'output','-ascii');
        end
        if ans=='p' %subtract background
            figure(hdl3);
            [X,Y]=ginput(2);
            backgroundrange=(floor(X(1)/timeunit):floor(X(2)/timeunit));
            backavgD=mean(donor(i,backgroundrange));
            backavgA=mean(acceptor(i,backgroundrange));
            donor(i,:)=donor(i,:)-backavgD;
            acceptor(i,:)=acceptor(i,:)-backavgA;
            iInc=0;
            disp(['D_bg=' num2str(backavgD) ' A_bg=' num2str(backavgA)])
        end
        
        if ans=='b' %backward
            iInc=-1;
        end
        
        if ans(1)=='g'
            i=str2num(ans(2:end))-1;
        end
        
        if ans=='j' %jump to trace
            tracenumber=input('Jump to trace number: [if empty just one further]');
            if isempty(tracenumber)
                iInc=1;
            else
                iInc=tracenumber-i;
            end
        end
        
        
        if ans=='a'
            tempD=sum(donor(i,(3:3)),2);
            tempA=sum(acceptor(i,(3:3)),2);
            
            elevel(i)=tempA/(tempA+tempD);           
        end
            
        
        
        %%%MS functions follow
            
        if ans=='c'  %%show a certain x axis section (just for viaulization)
            clear X Y mousebutton;
            display('Please click onthe point where unwinding starts.');
            [X,Y, mousebutton]=ginput;  % at the moment I only use X
            
            %if isempty(sectionstart)
            %     sectionstart=0;  %MS
            %end
            sectionstart=X-2;
            
            figure(hdl3);
            ax(1)=subplot(2,1,1);
            plot(time,donor(i,:),'g',time,acceptor(i,:),'r');
            maxcps=1.2*max(donor(i,:));
            axis([sectionstart (sectionstart+5) -50 maxcps]);
            grid on;
            zoom on;
            
            %First figure, second subplot
            ax(2)=subplot(2,1,2);
            plot(time,fretE,'k-');
            axis([sectionstart (sectionstart+5) -0.1 1.2]);
            grid on;
            zoom on;
            linkaxes(ax,'x'); %MS
            input('press any key to proceed','s');
        end % MS axis section viewer
        
                
        
        if ans=='u'    %MS version of the clicktimes program 
            Times=[];
            Mousebuttons=[];
            clear X Y mousebutton;
            if exist('times','dir')==7 && exist('traces','dir')==7
            else
                mkdir('times');
                mkdir('traces');
            end
            divisions=input('how many divisions do you want the window split into?');
            if isempty(divisions)
                divisions=5;  %MS
            end
            for j=1:divisions  %MS
                if isempty(Mousebuttons)    % This is only good if you just
                    %want to measure average unwindingtime. not valid for stalls!!
                elseif Mousebuttons(end)==3
                    break
                elseif length(Times')==4
                    break
                end
                figure(hdl3);
                ax(1)=subplot(2,1,1);
                plot(time,donor(i,:),'g',time,acceptor(i,:),'r');
                maxcps=1.2*max(donor(i,:)); %MS
                axis([(j-1)*time(end)/divisions j*time(end)/divisions -50 maxcps]);
                grid on;
                zoom on;
                
                %First figure, second subplot
                ax(2)=subplot(2,1,2);
                plot(time,fretE,'bx-');
                axis([(j-1)*time(end)/divisions j*time(end)/divisions -0.2 1.2]);
                grid on;
                zoom on;
                linkaxes(ax,'x'); %MS
                [X,Y, mousebutton]=ginput;  % at the moment I only use X
                Times=[Times; X];
                Mousebuttons=[Mousebuttons; mousebutton];
            end %for loop for j
            
            %ginput takes data in columns (transpose Times)
            Times=Times';
            Mousebuttons=Mousebuttons';
            
            clickindex=round(Times/timeunit);
            
            %timeDD = [timeUnwinding' timeWait' timeDelay'];
            s=length(Times)-1;
            if (s==1)   % timeunwinding is unwinding from high fret to low fret.
                timeUnwinding=Times(2)-Times(1);
                timeWait=Times(2)-Times(1);
                timeDelay=Times(1);
            end
            
            if (s==2)
                timeUnwinding=Times(2)-Times(1);
                timeWait=Times(3)-Times(1);
                timeDelay=Times(1);
            end
            
            if (s==3)  % MS special for ATPgS correlate times with mousebutton for single stalls only. 
                if(Mousebuttons(4)==1) & (Mousebuttons(1)==1)  % the first and last time with left click gives total unwinding time
                    timeUnwinding=Times(4)-Times(1);
                else
                    display('Please click the left mouse buttons for total unwinding time')
                end
                if(Mousebuttons(2)==3) %timestallstart gives time difference from unwinding start to stall start 
                    timestallstart=Times(2)-Times(1);
                    if(Mousebuttons(3)==3)
                        timestalltotal=Times(3)-Times(2); %timestalltotal gives the total time of first stall. 
                        averageFRETplateau=mean(fretE(ceil(Times(2)/timeunit):floor(Times(3)/timeunit)));  %calculates average fret value during stall
                    end
                    
                else
                    display('Please click the right mouse buttons for stall start and end')
                end
                timeWait=Times(4)-Times(1);  % total unwinding time. without test.
                timeDelay=Times(1);
            end
            
            
            lentrace=length(time);
            
            idealized=ones(1,lentrace);
            %offset=(times(1)/inttime)
            offset=0;
            
            counter= 2;
            
            %this idealization assumes first click is for a high state
            for k=1:s
                v=[clickindex(k)+1:clickindex(k+1)+1];
                %plus one comes bc input file starts at 0
                idealized(1,v)=rem(counter,3)*ones(1,length(v));
                counter=counter+1;
            end
            
            idealized=idealized/2.;
            
            %redisplays unzoomed trace
            
            %First figure, first subplot
            figure(hdl3);
            ax(1)=subplot(2,1,1);
            plot(time,donor(i,:),'g',time,acceptor(i,:),'r');
            axis([time(1) time(end) -50 maxcps]);
            grid on;
            zoom on;
            
            %First figure, second subplot
            ax(2)=subplot(2,1,2);
            plot(time,fretE,'b', time, idealized, 'm');
            axis([time(1) time(end) -0.2 1.2]);
            grid on;
            zoom on;
            linkaxes(ax,'x');
            happy=input('is this good, type "n" if you want it skipped?','s');
            
            if isempty(happy)
                happy='y';
            end
            
            if happy~='n'
                if (s==3)
                    fname3=['times/' fname '-' 'tr' num2str(i) 'UnW.dat'];
                    save(fname3, 'timeUnwinding', '-ascii');
                    fname4=['times/' fname '-' 'tr' num2str(i) 'starttoleaving.dat'];
                    save(fname4, 'timeWait', '-ascii');
                    fname5=['times/' fname '-' 'tr' num2str(i) 'Delay.dat'];
                    save(fname5, 'timeDelay', '-ascii');
                    fname6=['times/' fname '-' 'tr' num2str(i) 'stallstart.dat'];
                    save(fname6, 'timestallstart', '-ascii');
                    fname7=['times/' fname '-' 'tr' num2str(i) 'stallduration.dat'];
                    save(fname7, 'timestalltotal', '-ascii');
                    fname8=['times/' fname '-' 'tr' num2str(i) 'stallFRETavg.dat'];
                    save(fname8, 'averageFRETplateau', '-ascii');
                else
                    fname3=['times/' fname '-' 'tr' num2str(i) 'UnW.dat'];
                    save(fname3, 'timeUnwinding', '-ascii');
                    fname4=['times/' fname '-' 'tr' num2str(i) 'starttoleaving.dat'];
                    save(fname4, 'timeWait', '-ascii');
                    fname5=['times/' fname '-' 'tr' num2str(i) 'Delay.dat'];
                    save(fname5, 'timeDelay', '-ascii');
                end
                %save the analyzed trace as ascii file
                savename=['traces/' fname '-' 'tr' num2str(i) '.dat'];
                output=[time' fretE acceptor(i,:)' donor(i,:)'];
                save(savename,'output','-ascii');
            end
        end %MS version of the clicktimes program 
        
        
        if ans=='l'   %slippage detection -  looks at the FRET trace
            %click left-left for slippage without unwinding and
            %left right for slippage with unwinding
            Xpos=[];
            Mousebuttonsslip=[];
            Levels=[];
            clear X Y mousebutton;
            mkdir('slippage');
            divisions=input('how many divisions do you want the window split into?');
            if isempty(divisions)
                divisions=4;  %MS
            end
            for j=1:divisions  %MS
                if isempty(Mousebuttonsslip)
                elseif Mousebuttonsslip(end)==3
                    break
                end
                
                
                figure(hdl3);
                ax(1)=subplot(2,1,1);
                plot(time,donor(i,:),'g',time,acceptor(i,:),'r');
                maxcps=1.2*max(donor(i,:)); %MS
                axis([(j-1)*time(end)/divisions j*time(end)/divisions -50 maxcps]);
                grid on;
                zoom on;
                
                %First figure, second subplot
                ax(2)=subplot(2,1,2);
                % fretE2=smooth(fretE,7,'moving',2);
                windowSize=5; 
                fretE2=filter(ones(1,windowSize)/windowSize,1,fretE);
                plot(time,fretE2,'b');
                axis([(j-1)*time(end)/divisions j*time(end)/divisions -0.2 1.2]);
                grid on;
                zoom on;
                linkaxes(ax,'x'); %MS
                [X,Y, mousebutton]=ginput;  % reads X  and Y position of mouse, leftclick start , right click stop. first Y click is important for level to cut
                Xpos=[Xpos; X];
                Levels=[Levels; Y];
                Mousebuttonsslip=[Mousebuttonsslip; mousebutton];
            end %for loop for j
            sizeXpos=length(Xpos)/2;
            Mousebuttonsslip=Mousebuttonsslip';
            Levels=Levels';
            Xpos=Xpos';
            unsuccounter=0;
            for k=1:sizeXpos
                if (Mousebuttonsslip(2*k-1)==1) & (Mousebuttonsslip(2*k)==3)  %% attempts before sucessfull unwinding -- save name is jsut different, nothing else in the detection
                    slippagestarttime=[];
                    slippagelifetime=[];
                    fp = round(Xpos(2*k-1)/timeunit);
                    lp = round(Xpos(2*k)/timeunit);
                    tmptime = time(fp:lp);
                    tmpFRET = fretE2(fp:lp);
                    subsectionlength=length(tmptime');
                    slip=0;
                    for l=2:subsectionlength
                        if (tmpFRET(l)<Levels(2*k-1)) & (slip==0)
                            slippagestarttime=[slippagestarttime;tmptime(l)];
                            slipstart=tmptime(l);
                            slip=1;
                        elseif (slip==1) & tmpFRET(l)>Levels(2*k-1)
                            sliptime=tmptime(l)-slipstart;
                            slippagelifetime=[slippagelifetime;sliptime];
                            slip=0;
                        else
                            
                        end
                    end
                    slippagelifetime=[slippagelifetime;111]; %%write 111 more in slippage lifetime. important for saving, ignored later for plotting
                    %slippagestarttime(end)=[]; %% make slippagestart one shorter, since the last event was the unwinding.
                    fnameslip=['slippage/' fname '-' num2str(i) '-' 'slip_succesful.dat']; %first colume slipstart, second column slipduration
                    output=[slippagestarttime slippagelifetime];
                    save(fnameslip, 'output', '-ascii');
                    
                elseif (Mousebuttonsslip(2*k-1)==1) & (Mousebuttonsslip(2*k)==1)   % no successful unwinding
                    slippagestarttime=[];
                    slippagelifetime=[];
                    unsuccounter=unsuccounter+1;
                    fp = round(Xpos(2*k-1)/timeunit);
                    lp = round(Xpos(2*k)/timeunit);
                    tmptime = time(fp:lp);
                    tmpFRET = fretE2(fp:lp);
                    subsectionlength=length(tmptime');
                    slip=0;
                    for l=2:subsectionlength
                        if (tmpFRET(l)<Levels(2*k-1)) & (slip==0)
                            slippagestarttime=[slippagestarttime;tmptime(l)];
                            slipstart=tmptime(l);
                            slip=1;
                        elseif (slip==1) & tmpFRET(l)>Levels(2*k-1)
                            sliptime=tmptime(l)-slipstart;
                            slippagelifetime=[slippagelifetime;sliptime];
                            slip=0;
                        else
                        end
                    end
                    if length(slippagestarttime)==length(slippagelifetime)
                        fnameslipun=['slippage/' fname '-' num2str(i) '-slip_unsuccesful-' num2str(unsuccounter) '.dat']; %first colume slipstart, second column slipduration
                        output=[slippagestarttime slippagelifetime];
                        save(fnameslipun, 'output', '-ascii');
                    else
                        display('There were more starting times than end times detected, please check again!')
                    end
                end
            end
            
        end   %ans==l end slippage detection 
        
        
        if ans=='o'    %MS mutiple stalls in unwinding trace
            Times=[];
            averageFRETplateau=[];
            lifetimestalls=[];
            Mousebuttons=[];
            clear X Y mousebutton;
            if exist('multiplestalls','dir')==7
            else
                mkdir('multiplestalls');
            end
            divisions=input('how many divisions do you want the window split into?');
            if isempty(divisions)
                divisions=5;  %MS
            end
            for j=1:divisions  %MS
              
                figure(hdl3);
                ax(1)=subplot(2,1,1);
                plot(time,donor(i,:),'g',time,acceptor(i,:),'r');
                maxcps=1.2*max(donor(i,:)); %MS
                axis([(j-1)*time(end)/divisions j*time(end)/divisions -50 maxcps]);
                grid on;
                zoom on;
                
                %First figure, second subplot
                ax(2)=subplot(2,1,2);
                plot(time,fretE,'bx-');
                axis([(j-1)*time(end)/divisions j*time(end)/divisions -0.2 1.2]);
                grid on;
                zoom on;
                linkaxes(ax,'x'); %MS
                [X,Y, mousebutton]=ginput;  % at the moment I only use X and calcualte the average FRET value, left and right mousebuttons make different stuff
                Times=[Times; X];
                Mousebuttons=[Mousebuttons; mousebutton];
            end %for loop for j
            
            %ginput takes data in columns (transpose Times)
            Times=Times';
            Mousebuttons=Mousebuttons';
            
            clickindex=round(Times/timeunit);
            
            %timeDD = [timeUnwinding' timeWait' timeDelay'];
            s=length(Times);         %% Attention new s definition
            if (mod(s,2))
                display('Please reanalyze, you clicked odd number of times')
            else
                half=s/2;
                for nn=1:half
                    if(Mousebuttons(2*nn-1)==1) && (Mousebuttons(2*nn)==3)   % check left right click
                        avg=mean(fretE(ceil(Times(2*nn-1)/timeunit):floor(Times(2*nn)/timeunit))); % calculate average
                        averageFRETplateau=[averageFRETplateau; avg]; % add average FRET value
                        life=Times(2*nn)-Times(2*nn-1);   %% start end stll lifetime
                        lifetimestalls=[lifetimestalls; life];  %% add to list of stlll lifetimes
                    else
                        display('Please reanalyze: stall start left click, stall end right click')
                    end
                end
            end
            
            
            lentrace=length(time);
            
            idealized=ones(1,lentrace);
            %offset=(times(1)/inttime)
            offset=0;
            
            counter= 2;
            
            %this idealization assumes first click is for a high state
            for k=1:(s-1)
                v=[clickindex(k)+1:clickindex(k+1)+1];
                %plus one comes bc input file starts at 0
                idealized(1,v)=rem(counter,3)*ones(1,length(v));
                counter=counter+1;
            end
            
            idealized=idealized/2.;
            
            %redisplays unzoomed trace
            
            %First figure, first subplot
            figure(hdl3);
            ax(1)=subplot(2,1,1);
            plot(time,donor(i,:),'g',time,acceptor(i,:),'r');
            axis([time(1) time(end) -50 maxcps]);
            grid on;
            zoom on;
            
            %First figure, second subplot
            ax(2)=subplot(2,1,2);
            plot(time,fretE,'b', time, idealized, 'm');
            axis([time(1) time(end) -0.2 1.2]);
            grid on;
            zoom on;
            linkaxes(ax,'x');
            happy=input('is this good, type "n" if you want it skipped?','s');
            
            if isempty(happy)
                happy='y';
            end
            
            if happy~='n'
                fname3=['multiplestalls/' fname '-' 'tr' num2str(i) 'stallFRETavgs.dat'];
                save(fname3, 'averageFRETplateau', '-ascii');
                averageFRETplateau
                fname4=['multiplestalls/' fname '-' 'tr' num2str(i) 'stalldurations.dat'];
                save(fname4, 'lifetimestalls', '-ascii');
                %save the analyzed trace as ascii file
                savename=['multiplestalls/' fname '-' 'tr' num2str(i) '.dat'];
                output=[time' fretE acceptor(i,:)' donor(i,:)'];
                save(savename,'output','-ascii');
            end
        end %end mutliple stalls
        
     
        
    end
    
end



