function [msac,eyetrial] = iscan_detect_saccades(xrawo,yrawo,micthresh,...
                                        sangnum,RFx,RFy,plotit) 
%** Note: revision 1/4/2012, by Jude Mitchell for Sina Salehi data 
%**** [msac,eyetrial] = iscan_detect_saccade(xrawo,yrawo,micthresh,sangnum,RFx,RFy)
%****** main idea (detection algorithm):  
%******   A) filter out any points of missing data in recording
%******        (as indicated by steps in position > 5 deg in 5ms)
%******   B) search every point surrounded by a 400ms window
%******      in eye position by doing parametric curve fitting:
%******        a) fit a 6th order spline over that window (7 params with mean)
%******        Then fit a saccade model, spline plus single sigmoid step
%******        b) the saccade model is 2nd order spline + 
%******                               sigmoid (width + amp parameter)
%******                               1st derivative of sigmoid (gauss)
%******                               2nd derivative of sigmoid (dog)
%******        c) if percentage improvement greater than 30% for (b) over (a)
%******           then label it a saccade and substitute the discontinuity
%******           model fit in place of the spline fit
%******        d) compute the velocity and acceleration based on the 
%******           traces fit by the models, if a saccade does not pass
%******           a 10 deg/s velocity threshold and a 1000 deg/s^2 
%******           acceleration threshold then do not count it
%**********************************************************************
%*****input: 
%******  xrawo - x position over time (at 1000Hz), in visual degrees
%******  yrawo - y position over time (at 1000Hz), in visaul degrees
%******  micthresh - reject saccades of amplitude smaller than this
%******  sangnum - build histogram of saccade sangnum directions
%******  RFx - x coordinate of neuron RF, set to 1 if not relevent
%******  RFy - y coordinate of neuron RF, set to 0 if not relevent
%******  plotit - if 1, then plot out the raw trace, fit, and saccades
%*****output:
%****** msac - list of microsaccade infomation, each row being
%******        [time of peak vel, start time, end time, size, peak velocity, 
%******         direction relative to RF, direction relative to real space]
%**************
%******eyetrial.xrawa = raw x trace filtered by FIR -3dB at 50Hz
%******eyetrial.yrawa = raw y trace filtered by gauss 
%******eyetrial.vrawa = velocity computed off xrawa and yrawa
%************* other data fit via splines to raw data **********
%******eyetrial.xfit = spline fit x position trace 
%******eyetrial.yfit = spline fit y position trace
%******eyetrial.dxfit = derivative of fit x
%******eyetrial.dyfit = derivative of fit y
%******eyetrial.vfit = velocity of the spline fits
%******eyetrial.acfit = acceleration of the split fits
%******eyetrial.xfita = same as xfit, but rotated to RF coordinates;
%******eyetrial.yfita,dxfita,dyfita
%********************************
%******** other internal parameters:
%***** model_thresh: percentage improvement of sigmoid model over
%*****       spline model necessary to declare an actual saccade event
%***** vel_thresh: min velocity threshold for saccade, 10 deg/s
%***** acc_thresh: min acceleration, 1000 deg/s, marks saccade start+end
%***** TM: downsampling from 1000Hz to lower, sample every 4ms
%*********************************************************  

  %**************************
  vel_thresh = 10;    % must have velocity > 10 deg/s (FIR filtered signal)
  acc_thresh = 1000;    % must have acceleration > 1000 deg/s^2 
  model_thresh = 30;   % at least at this % improvement to accept sigmoid
                       % step model over alternative spline model
  %***************************
  NM = size(xrawo,2);
  TM = 4;   % down sample traces since i-scan max sampling is 240 Hz
  downsample = TM:TM:NM;
  %*********** go through whole series fitting splines
  TT = 400;          % window for fitting raw data 
  TT2 = floor(TT/2);
  T = floor(TT/TM);  % look at plus and minus 200 ms from current point        
  T2 = floor(T/2);  % half-width of smoothing 
  ET = 1:T;         % evaluate error in this rage to decide if a saccade
  
  %************* filter out any huge fast steps due to signal failures
  xrawa = zeros(1,NM);
  yrawa = zeros(1,NM);
  oldx = xrawo(1);
  oldy = yrawo(1);
  xrawa(1) = oldx;
  yrawa(1) = oldy;
  for jj = 2:NM
      %******************
      if ( abs(xrawo(jj)-oldx) > 2) | (abs(yrawo(jj)-oldy) > 2)
          xrawa(jj) = oldx;
          yrawa(jj) = oldy;
      else
          xrawa(jj) = xrawo(jj);
          yrawa(jj) = yrawo(jj);
          oldx = xrawo(jj);
          oldy = yrawo(jj);
      end
  end
  
  %************ identify the mean (discretization of data)
  xd = xrawa(2:NM) - xrawa(1:(NM-1));
  yd = yrawa(2:NM) - yrawa(1:(NM-1));
  xd(NM) = 0;
  yd(NM) = 0;
  %********** mean non-zero step size from ms to ms ********
  yy = find( (xd>0.02) & (xd<0.4) );  
  minxd = mean( xd(yy) );
  yy = find( (yd>0.02) & (yd<0.4) );
  minyd = mean( yd(yy) );
  %**********************************************************
  
  %***************** smooth the raw eye traces **************
  firfilt = fir1(31,0.14);  % FIR filter with -3dB at 50 Hz  
  xreal = smoothfilt(xrawa,firfilt)';
  yreal = smoothfilt(yrawa,firfilt)';
  %**************** compute the velocity from smooth trace ****
  dxreal = zeros(1,NM);
  dxreal(1:(NM-1)) = 1000 * (xreal(2:NM) - xreal(1:(NM-1)));
  dyreal = zeros(1,NM);
  dyreal(1:(NM-1)) = 1000 * (yreal(2:NM) - yreal(1:(NM-1)));
  vreal = sqrt( (dxreal .^ 2) + (dyreal .^ 2) );
  %********** smooth velocity and compute the acceleration ****
  x2real = smoothfilt(dxreal,firfilt)';
  y2real = smoothfilt(dyreal,firfilt)';
  dx2real = zeros(1,NM);
  dx2real(1:(NM-1)) = 1000 * (x2real(2:NM) - x2real(1:(NM-1)));
  dy2real = zeros(1,NM);
  dy2real(1:(NM-1)) = 1000 * (y2real(2:NM) - y2real(1:(NM-1)));
  areal = sqrt( (dx2real .^ 2) + (dy2real .^ 2) );
  %**************** down sample for curve fitting ***********
  M = size(downsample,2);
  %*********** add noise of discretization error ... why, suppose
  %*********** a value is constant for several time points, with
  %*********** what accuracy do you really know that value is constant,
  %*********** no better than your discretization.  Adding continuous
  %*********** noise here insures the bar is raised to account for this.
  xrawdat = xreal(downsample);
  yrawdat = yreal(downsample);
  vrawa = vreal(downsample);
  arawa = areal(downsample);
  
  %************* flag all 50 ms intervals where vel > 10deg/s and
  %************* where the acc > 1000 deg/s^2 ******************
  vsearch = (T2+1):(M-T2);
  
  %**************** now fit spline + sigmoid step models ***
  %*********************************************************
  xa = ones(1,T);  % mean of fit
  xb = 1:T;
  xb = (xb - mean(xb)) * (1/T);  %1st order, linear
  xc = xb .^ 2;                  %2nd order, quadratic
  xc = (xc - mean(xc));
  xc = xc / sum( xc .^ 2 );
  xd = xc .* xb;
  xd = (xd - mean(xd));
  xd = xd / sum( xd .^ 2 );
  xd2 = xd .* xb;
  xd2 = (xd2 - mean(xd2));
  xd2 = xd2 / sum(xd2 .^ 2);
  xd3 = xd2 .* xb;
  xd3 = (xd3 - mean(xd3));
  xd3 = xd3 / sum(xd3 .^ 2);
  xd4 = xd3 .* xb;
  xd4 = (xd4 - mean(xd4));
  xd4 = xd4 / sum(xd4 .^ 2);
  XX = [xa ; xb ; xc ; xd ; xd2 ; xd3 ; xd4 ];  %mean,linear,2nd,3rd,4th,5th order terms
  Auto = pinv( XX * XX' );
  
  %*********** fit bursts and sigmoid steps of different width
  sigmas = (exp(0:0.5:4.0)-0.9);  % 5 coarse search of sigma parameter in sigmoid
  sigmas = sigmas / TM;  % resize for tighter time scale
  SN = size(sigmas,2);
  ZZ = cell(1,SN);
  AutoZZ = cell(1,SN);
  xe = cell(1,SN);
  xf = cell(1,SN);
  xg = cell(1,SN);
  for jj = 1:SN
   xe{jj} = 1:T;
   xe{jj} = xe{jj} - mean(xe{jj});
   xe{jj} = xe{jj} / sigmas(jj);
   xe{jj} = 1 + exp( -xe{jj} );
   xe{jj} = 1 ./ xe{jj};  % a sigmoid function
   xf{jj} = xe{jj} .* (1 - xe{jj});  % derivative is sigmoid resembles gauss
   xg{jj} = xf{jj} .* (1 - (2*xe{jj}));  % 2nd derivative resembles diff of gauss
   xe{jj} = xe{jj} - mean(xe{jj});
   xf{jj} = xf{jj} - mean(xf{jj});
   xg{jj} = xg{jj} - mean(xg{jj});   
  end
  %******************************************
  for kk = 1:SN
   ZZ{1,kk} = [xa ; xb ; xc ; xe{kk} ; xf{kk} ; xg{kk} ];   
   AutoZZ{1,kk} = pinv( ZZ{1,kk} * ZZ{1,kk}' );
  end
  
  %************** a second fit is done with higher time precision
  sigmas2 = (exp( 0:0.15:4 ) - 0.9);  % 21 log spaced values
  SN2 = size(sigmas2,2);
  ZZ2 = cell(1,SN2);
  AutoZZ2 = cell(1,SN2);
  %***********************
  xa2 = ones(1,TT);
  xb2 = 1:TT;
  xb2 = (xb2 - mean(xb2)) * (1/TT);
  xc2 = xb2 .^ 2;
  xc2 = (xc2 - mean(xc2));
  xc2 = xc2 / sum( xc2 .^ 2 );
  %********************************
  xe = cell(1,SN2);
  xf = cell(1,SN2);
  for jj = 1:SN2
   xe{jj} = 1:TT;
   xe{jj} = xe{jj} - mean(xe{jj});
   xe{jj} = xe{jj} / sigmas2(jj);
   xe{jj} = 1 + exp( -xe{jj} );
   xe{jj} = 1 ./ xe{jj};  % a sigmoid with sigma 10ms
   xf{jj} = xe{jj} .* (1 - xe{jj});
   xg{jj} = xf{jj} .* (1 - (2*xe{jj}));  % 2nd derivative resembles diff of gauss 
   xe{jj} = xe{jj} - mean(xe{jj});
   xf{jj} = xf{jj} - mean(xf{jj});
   xg{jj} = xg{jj} - mean(xg{jj});
  end
  %*********************
  for kk = 1:SN2
      ZZ2{1,kk} = [xa2 ; xb2 ; xc2 ; xe{kk} ; xf{kk} ; xg{kk} ]; 
      AutoZZ2{1,kk} = pinv( ZZ2{1,kk} * ZZ2{1,kk}' );
  end 
  
  %******** compute a smoothed version of eye trace with splines
  xsmo = zeros(1,M);
  ysmo = zeros(1,M);
  %********* fill in start by spline smoothing 
  xda = xrawdat(1:T);
  yda = yrawdat(1:T);
  Wx = Auto * (XX * xda');
  Wy = Auto * (XX * yda');
  xh = Wx' * XX;
  yh = Wy' * XX;
  xsmo(1:(T2+2)) = xh(1:(T2+2));
  ysmo(1:(T2+2)) = yh(1:(T2+2));
  %********* fill in end by spline smoothing
  xda = xrawdat((M-T+1):M);
  yda = yrawdat((M-T+1):M);
  Wx = Auto * (XX * xda');
  Wy = Auto * (XX * yda');
  xh = Wx' * XX;
  yh = Wy' * XX;
  xsmo((M-T2-1):M) = xh((T2-1):T);
  ysmo((M-T2-1):M) = yh((T2-1):T);
  %************************************
  
  %********* now search for model fits of candidate saccade points
  stepval = zeros(1,M);
  for t = vsearch  
    %*********** for every point compute its smoothed spline  
    xda = xrawdat((t-T2):(t+T2-1));
    yda = yrawdat((t-T2):(t+T2-1));
    Wx = Auto * (XX * xda');
    Wy = Auto * (XX * yda');
    xh = Wx' * XX;
    yh = Wy' * XX;
    xsmo(t) = xh(T2);
    ysmo(t) = yh(T2);
         
    %************* find best model with 2nd order spline, sigmoid and pulse 
    minerr = 10000;
    %**************** fit discontinuity model then ***************
    for kk = 1:SN
       ZWx = AutoZZ{1,kk} * (ZZ{1,kk} * xda');
       ZWy = AutoZZ{1,kk} * (ZZ{1,kk} * yda');
       zxh = ZWx' * ZZ{1,kk};
       zyh = ZWy' * ZZ{1,kk};
       %**** when computing error, realize zeros a limited by discretization limits
       val = rms_error(zxh,xda,zyh,yda,minxd,minyd);
       if (val<minerr)
           bzxh = zxh; 
           bzyh = zyh;
           minerr = val;
           zerrkk(t) = sigmas(kk);
       end
    end
     
    %******* compute the errors of two best models taking into
    %******* account limitations due to data discretization
    serr = discretized_error(xh(ET),xda(ET),yh(ET),yda(ET),minxd,minyd);
    szerr = discretized_error(bzxh(ET),xda(ET),bzyh(ET),yda(ET),minxd,minyd);
    %************** error ratio: percentage improvement
    stepval(t) = 200*(serr-szerr)/(serr+szerr);
  end
  
  %***** go back, find large improvements of step model *******
  %***** and mark the points of peak improvement as saccades **
  imsac = [];   % locate any peaks above threshold
  ii = 1;
  bump = floor(30/TM);  % search +/- this ms ahead
  bumpup = floor(100/TM); % min distance between events
  while (ii < M)
     if (stepval(ii) > model_thresh) 
        xii = max( (ii-bump), 1);
        xoo = min( (ii+bump), M);
        yy = find( stepval(xii:xoo) == max(stepval(xii:xoo)) );  % find peak sigmoid fit
        bxii = xii+yy(1);
        imsac = [imsac bxii];  % mark saccade location
        ii = bxii + bumpup;
     else
        ii = ii + 1;
     end
  end
  
  %********** now return to saccade points for a final fit of higher
  %********** temporal precision, fit the sigmoid to them and then
  %********** insert that fit over the spline fit at the saccade 
  %********** this is to provide best smooth estimate of eye position
  xfit = spline(downsample,xsmo,1:NM);
  yfit = spline(downsample,ysmo,1:NM);
  pweight = 1:TT;
  pweight = pweight - ((TT+1)/2);
  pweight = pweight / 50;
  pweight = - (pweight .^ 2);
  pweight = exp(pweight);
  %*************** now go back, fit the step as accurately as possible
  dipwin = 15;  %allow fit to shift in time slightly to better center sigmoid
  for ii = 1:size(imsac,2)
      xii = imsac(ii)*TM;
      %************** relocate the peak to best fit point ****
      xiia = max(TT2,(xii-dipwin));
      xiib = min((NM-TT2),(xii+dipwin));
      minerr = 100000;
      besttt = xii;
      cowo = [];
      %**************************************
      for tt = xiia:xiib
       ta = (tt-TT2);
       tb = (tt+TT2-1);
       %********* if hitting a boundary shift back
       htt = tt;
       while (ta<1)
           ta = ta + 1;
           tb = tb + 1;
           htt = htt + 1;
       end
       while (tb>NM)
           ta = ta - 1;
           tb = tb - 1;
           htt = htt - 1;
       end
       %**************************************
       cobo = [];
       xda = xreal(ta:tb); 
       yda = yreal(ta:tb);
       for kk = 1:SN2   
         ZWx = AutoZZ2{1,kk} * (ZZ2{1,kk} * xda');
         ZWy = AutoZZ2{1,kk} * (ZZ2{1,kk} * yda');
         zxh = ZWx' * ZZ2{1,kk};
         zyh = ZWy' * ZZ2{1,kk};
         val = rms_error(zxh,xda,zyh,yda,minxd,minyd);
         if (val<minerr)
           minerr = val;
           bestkk = kk;
           bestjj = jj;
           bzxh = zxh;
           bzyh = zyh;
           besttt = htt;
         end
         cobo = [cobo val];
       end
       cowo = [cowo ; cobo];
      end
      %****** modify time of saccade based on best fit ******
      imsac(ii) = besttt;
      %**********************************
      ta = (besttt-TT2);
      tb = (besttt+TT2-1);
      xda = xreal(ta:tb);
      yda = yreal(ta:tb);
      %**********************************
      xfit(ta:tb) = ((1-pweight) .* xfit(ta:tb)) + (pweight .* bzxh);
      yfit(ta:tb) = ((1-pweight) .* yfit(ta:tb)) + (pweight .* bzyh);
      %**********************************
      
      %********** plot individual fits : debugging ****
      if (0)  %set to 1 if you wish to plot raw fits to each event
       figure(30);
       subplot(2,2,1); hold off;
       plot(xda,'r:'); hold on;
       plot(yda,'b:');
       plot(bzxh,'r-');
       plot(bzyh,'b-');
       subplot(2,2,2); hold off;
       imagesc(cowo');
       subplot(2,2,3); hold off;
       plot(ta:tb,xrawa(ta:tb),'r:'); hold on;
       plot(ta:tb,yrawa(ta:tb),'b:');
       plot((ta:tb),xfit(ta:tb),'r-');
       plot((ta:tb),yfit(ta:tb),'b-');
       axis tight;
       [minxd,minyd]
       input('check');
      end
  end
  
  %*************** compute velocity and acceleration off spline fits
  dxfit = zeros(1,NM);
  dyfit = zeros(1,NM);
  dxfit(1:(NM-1)) = 1000 * (xfit(2:NM) - xfit(1:(NM-1)));
  dyfit(1:(NM-1)) = 1000 * (yfit(2:NM) - yfit(1:(NM-1)));
  vfit = sqrt( (dxfit .^ 2) + (dyfit .^ 2) );
  dx2fit = zeros(1,NM);
  dy2fit = zeros(1,NM);
  dx2fit(1:(NM-1)) = 1000 * (dxfit(2:NM) - dxfit(1:(NM-1)));
  dy2fit(1:(NM-1)) = 1000 * (dyfit(2:NM) - dyfit(1:(NM-1)));
  acfit = sqrt( (dx2fit .^ 2) + (dy2fit .^ 2) );
  %*************** tranform signal into RF coordinates ********
  angol = atan2(RFy,RFx);
  xfita = (cos(angol) * xfit) + (sin(angol) * yfit);
  yfita = (-sin(angol) * xfit) + (cos(angol) * yfit);
  dxfita = (cos(angol) * dxfit) + (sin(angol) * dyfit);
  dyfita = (-sin(angol) * dxfit) + (cos(angol) * dyfit);
  
  %******** take the identified saccades and get their metrics **
  msac = [];
  %**************************
  for ii = 1:size(imsac,2)
      xii = imsac(ii); %time of saccade fit  
      %********* find the peak velocity *****
      nearbo = 30;
      aa = (xii-nearbo);
      bb = (xii+nearbo);
      velpeak = max(max(vfit(aa:bb)));
      y = find( vfit(aa:bb) == velpeak );
      xii = (aa-1+y(1));
      %********* find the peak acceleration **
      XXY = aa:bb;
      y = find( acfit(XXY) == max( acfit(XXY) ) );
      acpeak = acfit( (aa+y(1)) );   % peak acceleration within 75ms
      
      %****** reinforce original contrainsts on velocity and acceleration
      %******* slow movements near noise events are not acceptable either
      if (acpeak<acc_thresh) | (velpeak<vel_thresh)
          disp(sprintf('Saccade removed below threshold %6.2f %7.1f',velpeak,acpeak));
          continue;
      end
      
      acthresh = min(acc_thresh,(acpeak/3));
      %**************** search for two consecutive point below that
      astart = max( (xii-75), 2);
      while (astart<=xii)
          if ( (acfit(astart)>acthresh) & (acfit(astart+1)>acthresh) )
                 break;
          end
          astart = astart + 1;
      end
      %*********************************************
      aend = min( (xii+75), (NM-2));
      while (aend>=xii)
          if ( (acfit(aend)>acthresh) & (acfit(aend-1)>acthresh) )
                  break;
          end
          aend = aend - 1;
      end

      %********* given start and end, compute saccade amplitude
      xgo = xfit(astart);
      xfi = xfit(aend);
      ygo = yfit(astart);
      yfi = yfit(aend);
      %********** find the max deviations in position ******
      xdev = range(xfit(astart:aend));
      xdev = xdev * sign((xfi-xgo));
      ydev = range(yfit(astart:aend));
      ydev = ydev * sign((yfi-ygo));
      distv = sqrt( xdev^2 + ydev^2 );
      %******************************
      dist = sqrt( (xfi-xgo)^2 + (yfi-ygo)^2 );   
      %**************************
      if (dist<micthresh)    
          continue;
      end
      sacdur = aend-astart;
      sacsize = dist;
      sacsizev = distv;
      
      %********** saccade direction, space relative**********
      xx = xdev; 
      yy = ydev; 
      if (xx < 0)
          ango = atan( (yy/xx) );
      else
          ango = pi - atan( (-yy/xx) );
      end
      ango = ango + (pi/sangnum);
      if (ango<0)
          ango = ango + (2*pi);
      end
      if (ango>=(2*pi))
          ango = ango - (2*pi);
      end
      angar2 = 1 + floor( ango / (2*pi/sangnum) );
      
      %********** saccade direction, relative to RF location
      xdeva = (cos(angol) * xx) + (sin(angol) * yy);
      ydeva = (-sin(angol) * xx) + (cos(angol) * yy);
      
      %***********************************
      xx = xdeva;
      yy = ydeva; 
      if (xx < 0)
           ango = atan( (yy/xx) );
      else
           ango = pi - atan( (-yy/xx) );
      end
      ango = ango + (pi/sangnum);
      if (ango<0)
           ango = ango + (2*pi);
      end
      if (ango>(2*pi))
           ango = ango - (2*pi);
      end
      angar = 1 + floor( ango / (2*pi/sangnum) );
      
      %***** store summary of saccade metrics ************
     
      % revision 1/4/12 by Jude Mitchell to deal with common
      %    problem appearing in Sina Salehi's data:
      %    it is flagging the same saccade twice
      %*** hack to fix this:
      %    check here the current saccade is not overlapping
      %    with the start and end of the previous one, and
      %    if so, then do not include it
      overlap = 0;
      if (~isempty(msac))
         %****** if a previous saccade, check how much it
         %****** overlaps the current one in time
         last = size(msac,1);
         prev_astart = msac(last,2);
         prev_aend = msac(last,3);
         oversum = 0;
         for zi = astart:aend
             if ( (zi >= prev_astart) & (zi <= prev_aend) )
                 oversum = oversum + 1;
             end
         end
         overlap = (oversum/(aend-astart));
         %********* if overlap less than 25%, then include it
         if ( overlap < 0.25)   
               msac = [msac ; [xii,astart,aend,sacdur,sacsize,...
                         velpeak,angar,angar2,sacsizev]];  
         end
      else
         msac = [msac ; [xii,astart,aend,sacdur,sacsize,...
                         velpeak,angar,angar2,sacsizev]];
      end
      
      if (0) % for debugging of saccade end and starts
        XX = (xii-100):(xii+100);
        figure(40);
        subplot(3,1,1); hold off;
        plot(XX,xrawa(XX),'r:'); hold on;
        plot(XX,yrawa(XX),'b:');
        plot(XX,xreal(XX),'r.'); hold on;
        plot(XX,yreal(XX),'b.');
        plot(XX,xfit(XX),'r-');
        plot(XX,yfit(XX),'b-');
        plot(astart,xfit(astart),'ro');
        plot(astart,yfit(astart),'bo');
        plot(aend,xfit(aend),'rs');
        plot(aend,yfit(aend),'bs');
        axis tight;
        subplot(3,1,2); hold off;
        plot(XX,vfit(XX),'k-'); hold on;
        plot(XX,vreal(XX),'k:');
        axis tight;
        V = axis;
        plot([V(1),V(2)],[vel_thresh,vel_thresh],'k-');
        plot([astart,astart],[V(3),V(4)],'b-');
        plot([aend,aend],[V(3),V(4)],'b-');
        subplot(3,1,3); hold off;
        plot(XX,acfit(XX),'k-'); hold on;
        plot(XX,areal(XX),'k:');
        axis tight;
        V = axis;
        plot([V(1),V(2)],[acc_thresh,acc_thresh],'k-');
        plot([astart,astart],[V(3),V(4)],'b-');
        plot([aend,aend],[V(3),V(4)],'b-');
        title(sprintf('%6.3f %6.3f range %6.3f overlap(%f)',velpeak,dist,distv,overlap));
        input('check');
      end    
  end % loop of all sac
  
  %*****************************************************
  if (0)   % plot prettier depiction of results 
     figure(10);
     subplot('position',[0.1 0.65 0.4 0.25]); hold off;
     goa = 1; %700;
     hoa = NM; %3100;
     Xt = goa:5:hoa;
     %******************************
     maxa = max(max([xrawa(Xt),yrawa(Xt)]));
     mina = min(min([xrawa(Xt),yrawa(Xt)]));
     hob = mina;
     gob = maxa;
     if ((gob-hob) < 1.5)
         dv = (1.5-(gob-hob))*0.5;
         gob = gob + dv;
         hob = hob - dv;
     end
     
     gcol = [0,0.5,0];
     bcol = [0,0.65,0.65];
     grycol = [0.6,0.6,0.6];
     
     H = plot([goa,hoa],[0,0],'k-'); hold on;
     set(H,'Linewidth',1.5);
     set(H,'Color',grycol);
     H = plot([1450,1450],[hob,gob],'k-');
     set(H,'Linewidth',1.5);
     set(H,'Color',grycol)
     H = plot([2450,2450],[hob,gob],'k-');
     set(H,'Linewidth',1.5);
     set(H,'Color',grycol);
     
     for jj = 1:size(msac,1)
        H = rectangle('Position',[msac(jj,2) hob (msac(jj,3)-msac(jj,2)) (gob-hob)]);
        set(H,'FaceColor',[1,0.6,0.6]); hold on;
        set(H,'EdgeColor',[1,0,0]); hold on;
     end
     
     H = plot(Xt,xrawa(Xt),'g-'); hold on;
     set(H,'Linewidth',1.0);
     set(H,'Color',gcol);
     H = plot(Xt,yrawa(Xt),'b-'); hold on;
     set(H,'Linewidth',1.0);
     set(H,'Color',bcol);
     H = plot(Xt,xfit(Xt),'k-');
     set(H,'Linewidth',1.0);
     H = plot(Xt,yfit(Xt),'k-');
     set(H,'Linewidth',1.0);
     axis([goa hoa hob gob]);
     
     xlabel('Time (ms)');
     ylabel('Visual Degrees');
     
     input('Check');
     close(10);
     
  end
  
  if (plotit == 1)
     X = 1:M;
     Xt = TM * X;
     XN = 1:NM;
  
     figure(10);
     subplot(4,1,1); hold off;
     H = plot(XN,xrawa(XN),'r:'); hold on;
     set(H,'Linewidth',1.5);
     H = plot(XN,yrawa(XN),'b:');
     set(H,'Linewidth',1.5);
     plot(xfit(XN),'r-'); hold on;
     plot(yfit(XN),'b-');
     axis tight;
     V = axis;
     V(3) = V(3) - 0.1;
     V(4) = V(4) + 0.1;
     axis([0 NM V(3) V(4)]);
     xlabel('Time (ms)');
     ylabel('Eyepos (degs)');
     title('Raw Traces and Spline Fits');
     
     subplot(4,1,2); hold off;
     plot(Xt,stepval(X),'k-'); hold on;
     plot([0,NM],[model_thresh,model_thresh],'b-');
     axis tight;
     V = axis;
     axis([0 NM -10 V(4)]);
     xlabel('Time (ms)');
     ylabel('Simoid Step % Improvement');
     title('Improvement by step discontinuity');
     
     if (1)
     subplot(4,1,3); hold off;
     plot(vfit(XN),'k-'); hold on;
     plot([1,NM],[vel_thresh,vel_thresh],'k-');
     %plot(Xt,vrawa(X),'g-');
     axis tight;
     V = axis;
     for jj = 1:size(msac,1)
        H = plot([msac(jj,2),msac(jj,2)],[0,V(4)],'b-');
        set(H,'Linewidth',1);
        H = plot([msac(jj,3),msac(jj,3)],[0,V(4)],'b-');
        set(H,'Linewidth',1);
     end
     V = axis;
     axis([0 NM V(3) V(4)]);
     xlabel('Time (ms)');
     ylabel('Velocity (deg/s)');
     title('Spline Fit Velocity and Threshold');
     
     subplot(4,1,4); hold off;
     plot(acfit(XN),'k-'); hold on;
     plot([1,NM],[acc_thresh,acc_thresh],'k-');
     % plot(areal(XN),'g-');
     axis tight;
     V = axis;
     for jj = 1:size(msac,1)
        H = plot([msac(jj,2),msac(jj,2)],[0,V(4)],'b-');
        set(H,'Linewidth',1);
        H = plot([msac(jj,3),msac(jj,3)],[0,V(4)],'b-');
        set(H,'Linewidth',1);
     end
     axis tight;
     V = axis;
     axis([0 NM V(3) V(4)]);
     xlabel('Time (ms)');
     ylabel('Accel (deg^2/s)');
     title('Spline Fit Acceleration and Threshold');
     
     end
     
     input('check');
  end

  eyetrial.minxd = minxd;  % discretization of x signal
  eyetrial.minyd = minyd;  % discretization of y signal
  eyetrial.xrawa = xrawa;  % raw signal, data failures removed
  eyetrial.yrawa = yrawa;
  eyetrial.xreal = xreal;  % smoothed version of signal
  eyetrial.yreal = yreal;  % smoothed version of signal
  eyetrial.vrawa = vreal;
  eyetrial.arawa = areal;
  %************************** other fits via the spline fit **********
  eyetrial.xfit = xfit;
  eyetrial.yfit = yfit;
  eyetrial.dxfit = dxfit;
  eyetrial.dyfit = dyfit;
  eyetrial.vfit = vfit;
  eyetrial.acfit = acfit;
  eyetrial.xfita = xfita;
  eyetrial.yfita = yfita;
  eyetrial.dxfita = dxfita;
  eyetrial.dyfita = dyfita;
  
return;


function output = smooth(input, window, kernel_type)
% Smoothing function: 
% output = smooth(input, window, type)
% "Type" should be set to 'gauss' for a gaussian kernel or to 'boxcar'
% for a simple moving window average.  "Window" is the total kernel width.
% Input array must be one-dimensional.
input_dims = ndims(input);
input_size = size(input);
if input_dims > 2 | min(input_size) > 1,
   disp('Input array is too large.');
   return
end

if input_size(2) > input_size(1),
   input = input';
   toggle_dims = 1;
else
   toggle_dims = 0;
end

if window/2 ~= round(window/2),
   window = window + 1;
end
halfwin = window/2;

input_length = length(input);

if kernel_type(1:5) == 'gauss',
   x = -halfwin:1:halfwin;
   kernel = exp(-x.^2/(window/2)^2);
else
   kernel = ones(window, 1);
end
kernel = kernel/sum(kernel);

padded(halfwin+1:input_length+halfwin) = input;
padded(1:halfwin) = ones(halfwin, 1)*input(1);
padded(length(padded)+1:length(padded)+halfwin) = ones(halfwin, 1)*input(input_length);

output = conv(padded, kernel);
output = output(window:input_length+window-1);

if toggle_dims == 1,
   output = output';
end


function output = smoothfilt(input, windowtemp)
% Smoothing function: 
% output = smooth(input, window)
% window is a filter template for smoothing
% Input array must be one-dimensional.
input_dims = ndims(input);
input_size = size(input);
if input_dims > 2 | min(input_size) > 1,
   disp('Input array is too large.');
   return
end

if input_size(2) > input_size(1),
   input = input';
   toggle_dims = 1;
else
   toggle_dims = 0;
end

window = size(windowtemp,2);
halfwin = floor(window/2);


input_length = length(input);
kernel = windowtemp;

padded(halfwin+1:input_length+halfwin) = input;
padded(1:halfwin) = ones(halfwin, 1)*input(1);
padded(length(padded)+1:length(padded)+halfwin) = ones(halfwin, 1)*input(input_length);

output = conv(padded, kernel);
output = output(window:input_length+window-1);

if toggle_dims == 1,
   output = output';
end

return;

%*********************************************************************
function val = discretized_error(zxh,xda,zyh,yda,minxd,minyd);
%********* suppose you observe an error of zero in which the true
%********* underlying value is in the middle of a discretization bin,size d,
%********* and is thus consistent, compare this to the case where the
%********* true value lies on the discretization boundary such that the
%********* value oscillates from + to - d/2, so an error of d/2 and 0 may
%********* be no different from each other.  On average, any error less
%********* than d/2 should be treated equally, as an error of d/4 
  xd = abs( zxh - xda );
  yy = find( xd < (minxd/2) );
  xd(yy) = (minxd/4);
  
  yd = abs( zyh - yda );
  yy = find( yd < (minyd/2) );
  yd(yy) = (minyd/4);
  
  val = mean( xd .^ 2) + mean( yd .^ 2);

return;

%*********************************************************************
function val = rms_error(zxh,xda,zyh,yda,minxd,minyd);
  
  xd = (zxh - xda);
  yd = (zyh - yda);
  
  val = sum( xd .^ 2) + sum( yd .^ 2);

return;