%% Function - Curvefit to recovery
%------------------------------------------------------------------
%function [Curvefit, goodness, x] = myfitY(y1,makeplot,si,unitsY,fig)
% [Curvefit, goodness] = myfity2(y1,makeplot,si,yUnit)
% Curve Fit
% 
% Input
% y1 - curve to be fitted must be exponential rising curve 
% makeplot - enter plot option PLOT (y), DONT PLOT (n)
% si - sample interval in seconds (si) 
% unitsY - units of y-axis
% fig - figure number
%
% Output
% Curvefit - curve fit to data
% goodness - goodness of fit least square error
% fitx - x data vector with time in seconds


    Y = y1(:);%curve, make column vector
    %Y data index starts at 1

    % Translate Y to origin, eliminating constant terms
    minY = min(Y);
    Y = Y - minY;

    % Create time vector
    x = 0:(length(Y)-1); %start x-axis from zero
    x = x*si; %convert to seconds
    x = x(:); %make column vector
    %x = x+(5*FrameTreat); %shift the curve to match time stamp
    %myc = round(min(Y));
    %myalpha = round(max(Y)-myc);

    myalpha = max(Y);
    %temp = 0.63*myalpha+myc; %find value 63%
    temp = 0.63*myalpha; %find value 63%
    temp = abs(Y-temp); %subtract value 66% of max from curve
    [~, tempx] = min(temp); %get index of min difference = index of 63%myalpha+myc
    mybeta = 1/(x(tempx));
    %r = 5*FrameTreat; %create right shift known constant (seconds)

    % s = fitoptions('Mehtod','NonlinearLeastSquares','Robust','LAR',...
    %                 'Lower',[myalpha-.01,0,myc-.01],...
    %                 'Upper',[myalpha+.01,inf,myc+.01],...
    %                 'Startpoint',[myalpha,mybeta,myc]); %set the fitting parameters
    s = fitoptions('Mehtod','NonlinearLeastSquares','Robust','LAR',...
                    'Lower',[myalpha-1,0],...
                    'Upper',[myalpha+1,inf],...
                    'Startpoint',[myalpha,mybeta]); %set the fitting parameters


    % f = fittype( @(alpha,beta,c,x) alpha*(1 - exp(-beta *(x-r))) +c,...
    %             'coefficients',{'alpha','beta','c'},'options',s ); %set the curve fit function
    % f = fittype( @(alpha,beta,x) alpha*(exp(-beta *(x))),...
    %             'coefficients',{'alpha','beta'},'options',s ); %set the curve fit function

    f = fittype( @(alpha,beta,x) alpha*(1 - exp(-beta *(x))),...
                'coefficients',{'alpha','beta'},'options',s ); %set the curve fit function

    [Curvefit, goodness] = fit(x,Y,f); %fit data to an exponential rise curve    
    
    
    % Plot the original waveform, fit, and fit stats in figure
    % make plot of curve fit if requested
    if strncmpi(makeplot,'y',1) %case insensitive comparison of first character
        %figure;

        figure(fig); plot(Curvefit,'r--'),hold on; %create a plot of the curve fit and original data
        plot(x, Y,'b.','markersize',12); %plot data points
        legend('Fit Curve', 'Data','Location','SouthEast'); %place legend
        xlabel('Time(seconds)','FontSize',20); %x-axis label
        ylabel(['Amplitude (', unitsY, ')'], 'FontSize',20); %y-axis label
        Txticks = get(gca,'XTick'); %x-axis tick array
        Pxticks = round(length(get(gca,'XTick'))/4); %1/4 length of x-axis tick array
        xticks = Txticks(Pxticks); %1/4 length value of x-axis ticks
        tempytick = abs(x-xticks); %x-axis data - x-axis 1/4 length tick value
        [~, tempytick] = min(tempytick); %find position of min difference

        % make strings of equations using tex syntax
        strn1 = '$y={\alpha}\left(1-\exp({-\beta}{x})\right)$';
        strn2 = sprintf('$y={%0.4f}\\left(1-\\exp({-%0.g}{x})\\right)$',Curvefit.alpha, Curvefit.beta);
        strn3 = sprintf('$R^2={%0.4f}$',abs(goodness.rsquare));
        strn4 = sprintf('$\\tau=\\frac{1}{\\beta}={%0.4g}\\left( seconds \\right)$',(1/Curvefit.beta));

        % % For plotting w.r.t. origin    
            text('Interpreter','latex','string',...
            {strn1;strn2;strn3;strn4},...
            'Position',[xticks,Y(tempytick)],...
            'VerticalAlignment','Cap','FontSize',16)

    end

end

%------------------------------------------------------------------
%% Function - Axon binary read
%------------------------------------------------------------------
%function [d,si]=abfload1(fn,varargin)
% 139 & 321
% ** function [d,si]=abfload(fn,varargin)
% loads and returns data in the Axon abf format.
% Data may have been acquired in the following modes:
% (1) event-driven variable-length
% (2) event-driven fixed-length 
% (3) gap-free
% Information about scaling, the time base and the number of channels and 
% episodes is extracted from the header of the abf file (see also abfinfo.m).
% All optional input parameters listed below (= all except the file name) 
% must be specified as parameter/value pairs, e.g. as in 
%          d=abfload('d:\data01.abf','start',100,'stop','e');
%
%
%Credits: Harald Hentschk, University of Tuebingen 
%


%
%                    >>> INPUT VARIABLES >>>
%
% NAME        TYPE, DEFAULT      DESCRIPTION
% fn          char array         abf data file name
% start       scalar, 0          only gap-free-data: start of cutout to be read (unit: sec)
% stop        scalar or char,    only gap-free-data: end of cutout to be read (unit: sec). 
%             'e'                 May be set to 'e' (end of file).
% sweeps      1d-array or char,  only episodic data: sweep numbers to be read. By default, 
%             'a'                 all sweeps will be read ('a')
% channels    cell array         names of channels to be read, like {'IN 0','IN 8'};
%              or char, 'a'       ** make sure spelling is 100% correct (including blanks) **
%                                 if set to 'a', all channels will be read
% chunk       scalar, 0.05       only gap-free-data: the elementary chunk size (Megabytes) 
%                                 to be used for the 'discontinuous' mode of reading data 
%                                 (when less channels shall be read than exist in file)
% machineF    char array,        the 'machineformat' input parameter of the
%              'ieee-le'          matlab fopen function. 'ieee-le' is the correct 
%                                 option for windows; depending on the
%                                 platform the data were recorded/shall be read
%                                 by abfload 'ieee-be' is the alternative.
% 
%                    <<< OUTPUT VARIABLES <<<
%
% NAME        TYPE            DESCRIPTION
% d                           the data read, the format depending on the recording mode
%             1. GAP-FREE:
%             2d array        2d array of size 
%                              '<data pts>' by '<number of chans>'
%                              Examples of access:
%                              d(:,2)       data from channel 2 at full length
%                              d(1:100,:)   first 100 data points from all channels
%             2. EPISODIC FIXED-LENGTH:
%             3d array        3d array of size 
%                              '<data pts per sweep>' by '<number of chans>' by '<number of sweeps>'
%                              Examples of access:
%                              d(:,2,:)            is a matrix which contains all events (at full length) 
%                                                  of channel 2 in its columns
%                              d(1:200,:,[1 11])   contains first 200 data points of events #1 
%                                                  and #11 of all channels
%             3. EPISODIC VARIABLE-LENGTH:
%             cell array      cell array whose elements correspond to single sweeps. Each element is
%                              a (regular) array of size
%                              '<data pts per sweep>' by '<number of chans>'
%                              Examples of access:
%                              d{1}                 a 2d-array which contains sweep #1 (all of it, all channels)
%                              d{2}(1:100,2)        a 1d-array containing first 100 data points of of channel 2 in sweep #1
%
% si          scalar          the sampling interval in usec
%
% Date of this version: Nov 23, 2006

% future improvements:
% - handle expansion of header in transition to file version 1.65 better

% --- defaults   
% gap-free
start=0.0;
stop='e';
% episodic
sweeps='a';
% general
channels='a';
% the size of data chunks (see above) in Mb. 0.05 Mb is an empirical value
% which works well for abf with 6-16 channels and recording durations of 
% 5-30 min
chunk=0.05;
machineF='ieee-le';
verbose=1;

% assign values of optional input parameters, if any were given
pvpmod(varargin);   

%if verbose, disp(['** ' mfilename]); end
d=[]; 
si=[];
if ischar(stop)
  if ~strcmpi(stop,'e')
    error('input parameter ''stop'' must be specified as ''e'' (=end of recording) or as a scalar');
  end
end

% --- obtain vital header parameters and initialize them with -1 
% temporary initializing var
tmp=repmat(-1,1,16);
% for the sake of typing economy set up a cell array (and convert it to a struct below)
% column order is 
%        name, position in header in bytes, type, value)
headPar={
  'fFileVersionNumber',4,'float',-1;  
  'nOperationMode',8,'int16',-1; 
  'lActualAcqLength',10,'int32',-1;
  'nNumPointsIgnored',14,'int16',-1;
  'lActualEpisodes',16,'int32',-1;
  'lFileStartTime',24,'int32',-1;
  'lDataSectionPtr',40,'int32',-1;
  'lSynchArrayPtr',92,'int32',-1;
  'lSynchArraySize',96,'int32',-1; 
  'nDataFormat',100,'int16',-1;            
  'nADCNumChannels', 120, 'int16', -1;
  'fADCSampleInterval',122,'float', -1; 
  'fSynchTimeUnit',130,'float',-1;
  'lNumSamplesPerEpisode',138,'int32',-1;        
  'lPreTriggerSamples',142,'int32',-1;        
  'lEpisodesPerRun',146,'int32',-1;        
  'fADCRange', 244, 'float', -1;
  'lADCResolution', 252, 'int32', -1;
  'nFileStartMillisecs', 366, 'int16', -1;
  'nADCPtoLChannelMap', 378, 'int16', tmp;
  'nADCSamplingSeq', 410, 'int16',  tmp;
  'sADCChannelName',442, 'uchar', repmat(tmp,1,10);
  'fADCProgrammableGain', 730, 'float', tmp;
  'fInstrumentScaleFactor', 922, 'float', tmp;
  'fInstrumentOffset', 986, 'float', tmp;
  'fSignalGain', 1050, 'float', tmp;
  'fSignalOffset', 1114, 'float', tmp;
  'nTelegraphEnable',4512,'int16',tmp;
  'fTelegraphAdditGain',4576,'float',tmp
};

fields={'name','offs','numType','value'};
s=cell2struct(headPar,fields,2);
numOfParams=size(s,1);
clear tmp headPar;

if ~exist(fn,'file'), error(['could not find file ' fn]); end
%if verbose, disp(['opening ' fn '..']); end
[fid,messg]=fopen(fn,'r',machineF); 

if fid == -1,error(messg);end;   
% determine absolute file size
fseek(fid,0,'eof');
fileSz=ftell(fid);
fseek(fid,0,'bof');

% read all vital information in header
% convert names in structure to variables and read value from header
for g=1:numOfParams,
  if fseek(fid, s(g).offs,'bof')~=0, 
    fclose(fid);
    error(['something went wrong locating ' s(g).name]); 
  end;
  sz=length(s(g).value);
  eval(['[' s(g).name ',n]=fread(fid,sz,''' s(g).numType ''');']);
  if n~=sz, 
    fclose(fid);    
    error(['something went wrong reading value(s) for ' s(g).name]); 
  end;
end;

if lActualAcqLength<nADCNumChannels, 
  fclose(fid);
  error('less data points than sampled channels in file'); 
end
% the numerical value of all recorded channels (numbers 0..15)
recChIdx=nADCSamplingSeq(1:nADCNumChannels);
% the corresponding indices into loaded data d
recChInd=1:length(recChIdx);
% the channel names, e.g. 'IN 8'
recChNames=(reshape(char(sADCChannelName),10,16))';
recChNames=recChNames(recChIdx+1,:);

chInd=[];
eflag=0;
if ischar(channels) 
  if strcmp(channels,'a')
    chInd=recChInd;
  else
    fclose(fid);
    error('input parameter ''channels'' must either be a cell array holding channel names or the single character ''a'' (=all channels)');
  end
else
  for i=1:length(channels)
    tmpChInd=strmatch(channels{i},recChNames,'exact');
    if ~isempty(tmpChInd)
      chInd=[chInd tmpChInd];
    else
      % set error flag to 1
      eflag=1;
    end
  end;
end
if eflag
  fclose(fid);
  disp('**** available channels:');
  disp(recChNames);
  disp(' ');
  disp('**** requested channels:');
  disp(strvcat(channels));
  error('at least one of the requested channels does not exist in data file (see above)');
end

% fFileVersionNumber needs a fix - for whatever reason its value is sometimes 
% a little less than what it should be (e.g. 1.6499999xxxx instead of 1.65)
fFileVersionNumber=.001*round(fFileVersionNumber*1000);

% gain of telegraphed instruments, if any
if fFileVersionNumber>=1.65
  addGain=nTelegraphEnable.*fTelegraphAdditGain;
  addGain(addGain==0)=1;
else
  addGain=ones(size(fTelegraphAdditGain));
end

% tell me where the data start
blockSz=512;
switch nDataFormat
  case 0
    dataSz=2;  % bytes/point
    precision='int16';
  case 1
    dataSz=4;  % bytes/point
    precision='float32';
  otherwise
    fclose(fid);
    error('invalid number format');
end;
headOffset=lDataSectionPtr*blockSz+nNumPointsIgnored*dataSz;
% fADCSampleInterval is the TOTAL sampling interval
si=fADCSampleInterval*nADCNumChannels;

if ischar(sweeps) && sweeps=='a'
  nSweeps=lActualEpisodes;
  sweeps=1:lActualEpisodes;
else
  nSweeps=length(sweeps);
end;  

switch nOperationMode
  case 1
    if verbose, 
      disp('data were acquired in event-driven variable-length mode'); 
    end
    warndlg('function abfload has not yet been thorougly tested for data in event-driven variable-length mode - please double-check that the data loaded is correct','Just a second, please');
    if (lSynchArrayPtr<=0 || lSynchArraySize<=0), 
      fclose(fid);
      error('internal variables ''lSynchArraynnn'' are zero or negative'); 
    end
    switch fSynchTimeUnit
      case 0  % time information in synch array section is in terms of ticks
        synchArrTimeBase=1;
      otherwise % time information in synch array section is in terms of usec
        synchArrTimeBase=fSynchTimeUnit;    
    end;  
    % the byte offset at which the SynchArraySection starts
    lSynchArrayPtrByte=blockSz*lSynchArrayPtr;
    % before reading Synch Arr parameters check if file is big enough to hold them
    % 4 bytes/long, 2 values per episode (start and length)
    if lSynchArrayPtrByte+2*4*lSynchArraySize<fileSz, 
      fclose(fid);
      error('file seems not to contain complete Synch Array Section'); 
    end
    if fseek(fid,lSynchArrayPtrByte,'bof')~=0, 
      fclose(fid);
      error('something went wrong positioning file pointer to Synch Array Section'); 
    end
    [synchArr,n]=fread(fid,lSynchArraySize*2,'int32');
    if n~=lSynchArraySize*2,
      fclose(fid);
      error('something went wrong reading synch array section');
    end
    % make synchArr a lSynchArraySize x 2 matrix
    synchArr=permute(reshape(synchArr',2,lSynchArraySize),[2 1]);
    % the length of episodes in sample points
    segLengthInPts=synchArr(:,2)/synchArrTimeBase;
    % the starting ticks of episodes in sample points WITHIN THE DATA FILE
    segStartInPts=cumsum([0 (segLengthInPts(1:end-1))']*dataSz)+headOffset;
    % start time (synchArr(:,1)) has to be divided by nADCNumChannels to get true value
    % go to data portion
    if fseek(fid,headOffset,'bof')~=0, 
      fclose(fid);
      error('something went wrong positioning file pointer (too few data points ?)'); 
    end
    for i=1:nSweeps,
      % if selected sweeps are to be read, seek correct position
      if ~isequal(nSweeps,lActualEpisodes), 
        fseek(fid,segStartInPts(sweeps(i)),'bof'); 
      end;
      [tmpd,n]=fread(fid,segLengthInPts(sweeps(i)),precision);
      if n~=segLengthInPts(sweeps(i)), 
        warning(['something went wrong reading episode ' int2str(sweeps(i)) ': ' segLengthInPts(sweeps(i)) ' points should have been read, ' int2str(n) ' points actually read']); 
      end;
      dataPtsPerChan=n/nADCNumChannels;
      if rem(n,nADCNumChannels)>0, 
        fclose(fid);
        error('number of data points in episode not OK'); 
      end
      % separate channels..
      tmpd=reshape(tmpd,nADCNumChannels,dataPtsPerChan);
      % retain only requested channels
      tmpd=tmpd(chInd,:);
      tmpd=tmpd';
      % if data format is integer, scale appropriately; if it's float, tmpd is fine 
      if ~nDataFormat
        for j=1:length(chInd),
          ch=recChIdx(chInd(j))+1;
          tmpd(:,j)=tmpd(:,j)/(fInstrumentScaleFactor(ch)*fSignalGain(ch)*fADCProgrammableGain(ch)*addGain(ch))...
            *fADCRange/lADCResolution+fInstrumentOffset(ch)-fSignalOffset(ch);
        end;
      end
      % now place in cell array, an element consisting of one sweep with channels in columns
      d{i}=tmpd;
    end;  
  case {2,5}
    if nOperationMode==2
      if verbose
        disp('data were acquired in event-driven fixed-length mode');  
      end
%     else 
%       if verbose
%         disp('data were acquired in waveform fixed-length mode (clampex only)');  
%       end
    end
    % determine first point and number of points to be read 
    startPt=0;
    dataPts=lActualAcqLength;
    dataPtsPerChan=dataPts/nADCNumChannels;
    if rem(dataPts,nADCNumChannels)>0, 
      fclose(fid);
      error('number of data points not OK'); 
    end
    dataPtsPerChanPerSweep=dataPtsPerChan/lActualEpisodes;
    if rem(dataPtsPerChan,lActualEpisodes)>0
      fclose(fid);
      error('number of data points not OK');
    end
    dataPtsPerSweep=dataPtsPerChanPerSweep*nADCNumChannels;
    if fseek(fid,startPt*dataSz+headOffset,'bof')~=0
      fclose(fid);
      error('something went wrong positioning file pointer (too few data points ?)'); 
    end
    d=zeros(dataPtsPerChanPerSweep,length(chInd),nSweeps);
    % the starting ticks of episodes in sample points WITHIN THE DATA FILE
    selectedSegStartInPts=((sweeps-1)*dataPtsPerSweep)*dataSz+headOffset;
    for i=1:nSweeps,
      fseek(fid,selectedSegStartInPts(i),'bof'); 
      [tmpd,n]=fread(fid,dataPtsPerSweep,precision);
      if n~=dataPtsPerSweep, 
        fclose(fid);
        error(['something went wrong reading episode ' int2str(sweeps(i)) ': ' dataPtsPerSweep ' points should have been read, ' int2str(n) ' points actually read']); 
      end;
      dataPtsPerChan=n/nADCNumChannels;
      if rem(n,nADCNumChannels)>0
        fclose(fid);
        error('number of data points in episode not OK'); 
      end;
      % separate channels..
      tmpd=reshape(tmpd,nADCNumChannels,dataPtsPerChan);
      % retain only requested channels
      tmpd=tmpd(chInd,:);
      tmpd=tmpd';
      % if data format is integer, scale appropriately; if it's float, d is fine 
      if ~nDataFormat
        for j=1:length(chInd),
          ch=recChIdx(chInd(j))+1;
          tmpd(:,j)=tmpd(:,j)/(fInstrumentScaleFactor(ch)*fSignalGain(ch)*fADCProgrammableGain(ch)*addGain(ch))...
            *fADCRange/lADCResolution+fInstrumentOffset(ch)-fSignalOffset(ch);
        end;
      end
      % now fill 3d array
      d(:,:,i)=tmpd;
    end;  
    
  case 3
    if verbose, disp('data were acquired in gap-free mode'); end
    % from start, stop, headOffset and fADCSampleInterval calculate first point to be read 
    %  and - unless stop is given as 'e' - number of points
    startPt=floor(1e6*start*(1/fADCSampleInterval));
    % this corrects undesired shifts in the reading frame due to rounding errors in the previous calculation
    startPt=floor(startPt/nADCNumChannels)*nADCNumChannels;
    % if stop is a char array, it can only be 'e' at this point (other values would have 
    % been caught above)
    if ischar(stop),
      dataPtsPerChan=lActualAcqLength/nADCNumChannels-floor(1e6*start/si);
      dataPts=dataPtsPerChan*nADCNumChannels;
    else
      dataPtsPerChan=floor(1e6*(stop-start)*(1/si));
      dataPts=dataPtsPerChan*nADCNumChannels;
      if dataPts<=0 
        fclose(fid);
        error('start is larger than or equal to stop'); 
      end
    end
    if rem(dataPts,nADCNumChannels)>0
      fclose(fid);
      error('number of data points not OK'); 
    end
    if fseek(fid,startPt*dataSz+headOffset,'bof')~=0, 
      fclose(fid);
      error('something went wrong positioning file pointer (too few data points ?)');
    end
    % *** decide on the most efficient way to read data:
    % (i) all (of one or several) channels requested: read, done
    % (ii) one (of several) channels requested: use the 'skip' feature of
    % fread 
    % (iii) more than one but not all (of several) channels requested:
    % 'discontinuous' mode of reading data. Read a reasonable chunk of data
    % (all channels), separate channels, discard non-requested ones (if
    % any), place data in preallocated array, repeat until done. This is
    % faster than reading the data in one big lump, separating channels and
    % discarding the ones not requested
    if length(chInd)==1 && nADCNumChannels>1
      % --- situation (ii)
      % jump to proper reading frame position in file
      if fseek(fid,(chInd-1)*dataSz,'cof')~=0
        fclose(fid);
        error('something went wrong positioning file pointer (too few data points ?)');
      end
      % read, skipping nADCNumChannels-1 data points after each read
      dataPtsPerChan=dataPts/nADCNumChannels;
      [d,n]=fread(fid,dataPtsPerChan,precision,dataSz*(nADCNumChannels-1));
      if n~=dataPtsPerChan,
        fclose(fid);
        error(['something went wrong reading file (' int2str(dataPtsPerChan) ' points should have been read, ' int2str(n) ' points actually read']);
      end
    elseif length(chInd)/nADCNumChannels<1
      % --- situation (iii)
      % prepare chunkwise upload:
      % preallocate d
      d=repmat(nan,dataPtsPerChan,length(chInd));
      % the number of data points corresponding to the maximal chunk size, 
      % rounded off such that from each channel the same number of points is
      % read (do not forget that each data point will by default be made a
      % double of 8 bytes, no matter what the original data format is) 
      chunkPtsPerChan=floor(chunk*2^20/8/nADCNumChannels);
      chunkPts=chunkPtsPerChan*nADCNumChannels;
      % the number of those chunks..
      nChunk=floor(dataPts/chunkPts);
      % ..and the remainder
      restPts=dataPts-nChunk*chunkPts;
      restPtsPerChan=restPts/nADCNumChannels;
      % chunkwise row indices into d
      dix=(1:chunkPtsPerChan:dataPtsPerChan)';
      dix(:,2)=dix(:,1)+chunkPtsPerChan-1;
      dix(end,2)=dataPtsPerChan;
      if verbose && nChunk
        disp(['reading file in ' int2str(nChunk) ' chunks of ~' num2str(chunk) ' Mb']);
      end
      % do it
      for ci=1:size(dix,1)-1
        [tmpd,n]=fread(fid,chunkPts,precision);
        if n~=chunkPts
          fclose(fid);
          error(['something went wrong reading chunk #' int2str(ci) ' (' ...
            int2str(chunkPts) ' points should have been read, ' int2str(n) ' points actually read']); 
        end
        % separate channels..
        tmpd=reshape(tmpd,nADCNumChannels,chunkPtsPerChan);
        d(dix(ci,1):dix(ci,2),:)=tmpd(chInd,:)';
      end
      % collect the rest
      [tmpd,n]=fread(fid,restPts,precision);
      if n~=restPts
        fclose(fid);
        error(['something went wrong reading last chunk (' ...
          int2str(restPts) ' points should have been read, ' int2str(n) ' points actually read']);
      end
      % separate channels..
      tmpd=reshape(tmpd,nADCNumChannels,restPtsPerChan);
      d(dix(end,1):dix(end,2),:)=tmpd(chInd,:)';
    else
      % --- situation (i)
      [d,n]=fread(fid,dataPts,precision);
      if n~=dataPts, 
        fclose(fid);
        error(['something went wrong reading file (' int2str(dataPts) ' points should have been read, ' int2str(n) ' points actually read']); 
      end;
      % separate channels..
      d=reshape(d,nADCNumChannels,dataPtsPerChan);
      d=d';
    end
    % if data format is integer, scale appropriately; if it's float, d is fine 
    if ~nDataFormat
      for j=1:length(chInd),
        ch=recChIdx(chInd(j))+1;
        d(:,j)=d(:,j)/(fInstrumentScaleFactor(ch)*fSignalGain(ch)*fADCProgrammableGain(ch)*addGain(ch))...
          *fADCRange/lADCResolution+fInstrumentOffset(ch)-fSignalOffset(ch);
      end;
    end
  otherwise
    disp('recording mode of data must be event-driven variable-length (1), event-driven fixed-length (2) or gap-free (3) -- returning empty matrix');
    d=[];
    si=[];
end;
  
fclose(fid);

%------------------------------------------------------------------
%% Automatic Population Spike detection
%------------------------------------------------------------------
%function [M1, M2, M3] = markers(waveform)
%Input a typical stimulation + popsipike waveform
%Output M1 = marker indices for 1st marker
%       M2 = marker indices for 2nd marker
%       M3 = marker indices for 3rd marker


%acquire indices and value of max
[y1, p1] =  max(waveform);
%partition waveform starting from the max index to end of file
waveform1 = waveform(p1:length(waveform));

%get indices and value of marker#2 w.r.t. partitioned waveform
[y2, m2] = min(waveform1);
M2 = [waveform(m2 + p1-1), m2 + p1-1]; %marker position on waveform not waveform1

%determine the position of 2nd marker
wavem1m2 = waveform1(1:m2);

%compute the discrete derivative
for k = 1:length(wavem1m2) - 1
    dwavem1m2(k) = wavem1m2(k+1) - wavem1m2(k);
end

%find the edges of the original waveform using the 1st derivative
n = 1;
redge = [];
fedge = [];
dghost = dwavem1m2;
trim = floor(0.2*length(dwavem1m2)); %trim search to avoid false positives
dwavem1m2 = dwavem1m2(trim:length(dwavem1m2)-trim);

for k = 1:length(dwavem1m2) - 1
    if dwavem1m2(k) < 0 && dwavem1m2(k+1) >= 0
        redge(n) = k;
        n = n+1;
    elseif dwavem1m2(k) >= 0 && dwavem1m2(k+1) < 0
        fedge(n) = k;
        n = n+1;
    end
end

%try alternate algorithm if no peaks exist
if isempty(fedge)

    %guess position of marker#1 as halfway point between max and marker#2
    tempM1pos = ceil(m2/2);

    %if the guessed value is the correct value, where dy/dx=0
    %then keep the gussed index
    if waveform1(tempM1pos+1) - waveform1(tempM1pos) <= 0.0005
        M1 = [waveform(tempM1pos + p1-1), tempM1pos + p1-1];

    else %if the guessed position is incorrect then sweep across a
        %region that spans 25% of distance (between max and marker#2)
        %to the left and right of guessed index. It is probable that
        %the dy/dx!=0 so assign marker#2 to the index that has the
        %smallest dy/dx

        dghost(1:tempM1pos-ceil(0.25*m2)) = nan;
        dghost(tempM1pos+ceil(0.25*m2):length(dghost)) = nan;
        [yi, xi] = max(dghost);
        
        M1 = [waveform(xi + p1-1), xi + p1-1];
        
%         tol = .00075;
%         i = tempM1pos - ceil(0.25*m2);
%         count = 0;
%         while waveform(i+1) - waveform1(i) >= tol
%             i = i+1;
%             if i == tempM1pos + ceil(0.25*m2)
%                 i = tempM1pos - ceil(0.25*m2);
%                 tol = tol + 0.00025;
%                 count = count+1;
%             end
%             if count > 500
%                 tempM1pos = i;
%                 break;
%             end
%             tempM1pos = i;
%         end
% 
%         M1 = [waveform(tempM1pos + p1-1), tempM1pos + p1-1];

    end


%try peak detection algorithm to find 1st marker
else
    fedge = trim + fedge;
    redge = trim + redge;    
    if length(fedge) > 1
    [ymax, xmax] = max(fedge);
    fedge = fedge(xmax);
    end
    
    %find the leftmost value of peak. Done to get midpoint of plateu.
    leftedge = fedge;
    boarder = 0;
    while boarder ~= 1
        if wavem1m2(leftedge) == wavem1m2(leftedge-1)
            leftedge = leftedge - 1;
        else
            boarder = 1;
        end
    end

    %find the rightmost value of peak. Done to get midpoint of plateu.
    rightedge = fedge;
    boarder = 0;
    while boarder ~= 1
        if wavem1m2(rightedge) == wavem1m2(rightedge+1)
            rightedge = rightedge + 1;
        else
            boarder = 1;
        end
    end

    %shift the marker acccordingly
    left = floor((leftedge - fedge)/2);
    right = floor((rightedge - fedge)/2);
    fedge = fedge + right + left;

    %assign marker position
    M1 = [waveform(fedge + p1-1), fedge + p1-1];
end

%the position of marker#3 relative to marker#2 and end of sweep is
%the max of partitioned waveform
[y3, m3] = max(waveform1(m2:length(waveform1)));
M3 = [waveform(m3 + p1-1 + m2-1), m3 + p1-1 + m2-1];

%------------------------------------------------------------------
%% Modified questdlg
%------------------------------------------------------------------
%function ButtonName=myquestdlg(Question,Title,Btn1,Btn2,Btn3,Default)
% ln-110: FigPos(1) = 950; %EDITED%get(0,'DefaultFigurePosition');
% ln-114: %EDITEDT%FigPos    = getnicedialoglocation(FigPos, get(0,'DefaultFigureUnits'));
% ln-172: %EDITED%FigPos(3)=max(FigPos(3),MsgTxtXOffset+NumButtons*(BtnWidth+2*DefOffset));
% ln-279: %EDITED%set(QuestFig ,'Position',getnicedialoglocation(FigPos, get(QuestFig,'Units')));
% ln-305-308: 
%EDITED%
% if DefaultValid
%     setdefaultbutton(QuestFig, BtnHandle(DefaultButton));
% end