function [tau,R,f,S]=aint(t,u,w,dt,K2,varargin)
%
% for further information see
% http://ldvproc.nambis.de/programs
%
% Inputs: (call the function e.g. like aint(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 in the
%    present version
% 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 aint(t,u,w,dt,K2,'meanrem',
%          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
%
% 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;
  p = inputParser;
  addParameter(p,'meanrem',defaultmeanrem);
  parse(p,varargin{:});
  meanrem=p.Results.meanrem;

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

    Je=ceil(t(N)/dt)+K2+1;
    Je=2^ceil(log(Je)/log(2));
    
    ue=zeros([1,Je])+0j;
    we=zeros([1,Je]);
    for i=1:N-1
      for j=ceil(t(i)/dt)+1:ceil(t(i+1)/dt)
        we(j)=1.0;
        ue(j)=u(i);
      end
    end
    if meanrem
      me=sum(we.*ue)/sum(we);
      ue=(ue-me).*we;
      ne=zeros([1,Je]);
      for i=1:N-1
        for j=ceil(t(i)/dt)+1:ceil(t(i+1)/dt)
          ne(ceil(t(i)/dt))=ne(ceil(t(i)/dt))+1;
        end
      end
    end
    U1=fft(ue);
    R1=ifft(conj(U1).*U1);
    U0=fft(we);
    R0=real(ifft(conj(U0).*U0));
    if meanrem
      Un=fft(ne);
      Rn=real(ifft(conj(Un).*Un));
    end

    dn=(N-1)/(t(N)-t(1));
    c=exp(-dn*dt)/(1-exp(-dn*dt))^2;
    tau=circshift((K1:K2),K1)*dt;
    R=zeros([1,K])+0j;
    for k=K1:K2
      ksm=mod(k-1,length(R1))+1;
      ksc=mod(k,length(R1))+1;
      ksp=mod(k+1,length(R1))+1;
      kd=mod(k,length(R))+1;
      if (k~=0) && (R0(ksc)>0) && (R0(ksm)>0) && (R0(ksp)>0)
        R(kd)=(2*c+1)*R1(ksc)/R0(ksc)-c*(R1(ksm)/R0(ksm)+R1(ksp)/R0(ksp));
      end
    end
    if R0(1)>0
      R(1)=R1(1)/R0(1);
    end
    if meanrem
      R=R+(R(1)*Rn(1)+2*sum(real(R(2:K2+1)).*Rn(2:K2+1)))/(sum(ne)^2-Rn(1)-2*sum(Rn(2:K2+1)));
    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 and the data vector must have the same length.")
  end
end
