function [tau,R,f,S]=aslot(t,u,w,dt,K2,varargin) 
%
% for further information see
% http://ldvproc.nambis.de/programs
%
% Inputs: (call the function e.g. like adir(t,u,w,dt,K2),
%         arguments in the given order)
%
% t: array of real values, random arrival times (single block
%    or whole LDV data set)
% u: array of real or complex values, velocities (single block
%    or whole LDV data set)
% w: array of real values, individual weights (single block
%    or whole LDV data set), e.g. transit times, not used if
%    atw==True)
% dt: real value, fundamental time step, e.g. temporal resolution
%     of the correlation function
% K2: positive integer, maximum lag time in the estimated correlation
%     function in units of dt, yields a correlation function estimated
%     for lag times from -K2*dt to K2*dt
%
% Options: (call the function e.g. like adir(t,u,w,dt,K2,'locnor',
%          true) or adir(t,u,w,dt,K2,'meanrem',false,'atw',false,
%          'locnor',true,'fuzzy',false,'selfproducts',false), 
%          optional arguments in arbitrary order)
%
% meanrem: boolean, if meanrem is true, the mean is calculated and
%          removed from the data block, this also enables Bessel's
%          correction, default is meanrem=true, set meanrem false
%          to disable mean removal and Bessel's correction
% atw: boolean, set atw true for forward-backward inter-arrival
%      time weighting (for legacy and research reasons only, do
%      not use this option for productive data processing)
% locnor: boolean, set locnor true to enable local normalization,
%         default is locnor=false
% fuzzy: boolean, set fuzzy true to enable fuzzy slotting (for
%        legacy and research reasons only, do not use this option 
%        for productive data processing)
% selfproducts: boolean, set selfproducts true to enable counting
%               also self-products. This may help in the case of
%               significant dead times. However, in the presence
%               of noise, counting self-products will cause a bias.
%
% Outputs: [tau,R,f,S]
%
% tau: array of real values, time lags of the correlation
%      (for lag times from -K2*dt to K2*dt)
% R: array of complex values, correlation values R(tau), with
%    dimension of u^2
% f: array of real values, frequencies of the power density
%    spectrum (between -1/(2*dt) and 1/(2*dt))
% S: array of complex values, power spectral density S(f), with
%    dimension of u^2*dt
%
  defaultmeanrem=true;
  defaultatw=false;
  defaultlocnor=false;
  defaultfuzzy=false;
  defaultselfproducts=false;
  p = inputParser;
  addParameter(p,'meanrem',defaultmeanrem);
  addParameter(p,'atw',defaultatw);
  addParameter(p,'locnor',defaultlocnor);
  addParameter(p,'fuzzy',defaultfuzzy);
  addParameter(p,'selfproducts',defaultselfproducts);
  parse(p,varargin{:});
  meanrem=p.Results.meanrem;
  atw=p.Results.atw;
  locnor=p.Results.locnor;
  fuzzy=p.Results.fuzzy;
  selfproducts=p.Results.selfproducts;

  N=length(t);
  if (length(u)==N) && (length(w)==N)
    K1=-K2;
    K=K2-K1+1;

    if atw
      wf=[t(2:N)-t(1:N-1),0];
      for i=1:length(wf)
        if (wf(i)>5*(t(N)-t(1))/(N-1)) || (wf(i)<0)
          wf(i)=0.0;
        end
      end
      wb=[0,t(2:N)-t(1:N-1)];
      for i=1:length(wb)
        if (wb(i)>5*(t(N)-t(1))/(N-1)) || (wb(i)<0)
          wb(i)=0.0;
        end
      end
    else
      wf=w;
      wb=w;
    end
    if meanrem
      ur=u-sum(wb.*u)/sum(wb);
    else
      ur=u;
    end

    Jr=1000;
    Je=ceil(t(N)/dt)+Jr+K2;
    if locnor
      Je=Je+mod(Je+1,2); %for loc. nor. allways odd
    %else
    %  Je=2^ceil(log(Je)/log(2))
    end
    df=1/(Je*dt);
    fe=(mod((0:Je-1)+fix(Je/2),Je)-fix(Je/2))*df; 
    if atw
      U1=zeros([1,Je])+0j;
      U0=zeros([1,Je])+0j;
      Ub=zeros([1,Je])+0j;
      Wb=zeros([1,Je])+0j;
      if locnor
        Ui=zeros([1,Je])+0j;
        Uj=zeros([1,Je])+0j;
        Qb=zeros([1,Je])+0j;
      end
      for i=1:N
        E=exp(-2j*pi*fe*t(i));
        Uf=wf(i)*ur(i)*E;
        U1=U1+conj(Ub).*Uf;
        Ub=Ub+wb(i)*ur(i)*E;
        Wf=wf(i)*E;
        U0=U0+conj(Wb).*Wf;
        if locnor
          Qf=wf(i)*abs(ur(i))^2*E;
          Ui=Ui+conj(Qb).*Wf;
          Uj=Uj+conj(Wb).*Qf;
          Qb=Qb+wb(i)*abs(ur(i))^2*E;
        end
        Wb=Wb+wb(i)*E;
      end
      if fuzzy
        U1=U1.*sinc(fe*dt).^2;
        U0=U0.*sinc(fe*dt).^2;
        if locnor
          Ui=Ui.*sinc(fe*dt).^2;
          Uj=Uj.*sinc(fe*dt).^2;
        end
      end
      R1=ifft(U1+conj(U1));
      if locnor
        Ri=real(ifft(Ui+conj(Uj)));
        Rj=real(ifft(Uj+conj(Ui)));
      end
      if (~locnor) || meanrem
        R0=real(ifft(U0+conj(U0)));
      end
    else
      U1=zeros([1,Je])+0j;
      U0=zeros([1,Je])+0j;
      if locnor
        U2=zeros([1,Je])+0j;
      end
      for i=1:N
        E=exp(-2j*pi*fe*t(i));
        U1=U1+w(i)*ur(i)*E;
        U0=U0+w(i)*E;
        if locnor
          U2=U2+w(i)*abs(ur(i))^2*E;
        end
      end
      E1=conj(U1).*U1-sum(abs(w.*ur).^2);
      if locnor
        Ei=conj(U2).*U0-sum(abs(w.*ur).^2);
        Ej=conj(U0).*U2-sum(abs(w.*ur).^2);
      end
      if (~locnor) || meanrem
        E0=conj(U0).*U0-sum(abs(w).^2);
      end
      if fuzzy
        E1=E1.*sinc(fe*dt).^2;
        if locnor
          Ei=Ei.*sinc(fe*dt).^2;
          Ej=Ej.*sinc(fe*dt).^2;
        end
        if (~locnor) || meanrem
          E0=E0.*sinc(fe*dt).^2;
        end
      end
      R1=ifft(E1);
      if locnor
        Ri=real(ifft(Ei));
        Rj=real(ifft(Ej));
      end
      if (~locnor) || meanrem
        R0=real(ifft(E0));
      end
    end
    if selfproducts
      R1(1)=R1(1)+sum(wb.*wf.*abs(ur).^2);
      if locnor
        Ri(1)=Ri(1)+sum(wb.*wf.*abs(ur).^2);
        Rj(1)=Rj(1)+sum(wb.*wf.*abs(ur).^2);
      end;
      if (~locnor) || meanrem
        R0(1)=R0(1)+sum(wb.*wf);
      end
    end
    tau=circshift((K1:K2),K1)*dt;
    R=zeros([1,K])+0j;
    if locnor
      s2u=sum(wb.*abs(ur).^2)/sum(wb);
      for k=K1:K2
        ks=mod(k,length(R1))+1;
        kd=mod(k,length(R))+1;
        if Ri(ks)*Rj(ks)>0
          R(kd)=s2u*R1(ks)/sqrt(Ri(ks)*Rj(ks));
        end
      end
    else
      for k=K1:K2
        ks=mod(k,length(R1))+1;
        kd=mod(k,length(R))+1;
        if R0(ks)>0
          R(kd)=R1(ks)/R0(ks);
        end
      end
    end
    if meanrem
      if selfproducts
        R=R+(R(1)*R0(1)+2*sum(real(R1(2:K2+1))))/(sum(wf)*sum(wb)-R0(1)-2*sum(R0(2:K2+1)))
      else
        R=R+(R(1)*R0(1)+2*sum(real(R1(2:K2+1)))+sum(wf.*wb.*abs(ur).^2))/(sum(wf)*sum(wb)-R0(1)-2*sum(R0(2:K2+1))-sum(wf.*wb));
      end
    end
    df=1/(K*dt);
    f=(mod((0:K-1)+fix(K/2),K)-fix(K/2))*df;
    S=dt*fft(R);
  else
    print("The time vector, the data vector and the weights vector must have the same length.")
  end
end
