function [mat,var,skew,kurt,m3,m2,m1,m0]=tmoments(t,f);

% remove background concentration
f=f-min(f(:)); 

%find first four moments
for j=1:size(f,1);
m0(j)=0;
for i=1:size(f,2)-1;
    m0(j)=m0(j)+mean(mean(f(j,i:i+1))).*(t(i+1)-t(i));
end;

m1(j)=0;
for i=1:size(f,2)-1;
    m1(j)=m1(j)+(mean(t(i:i+1))).*mean(f(j,i:i+1)).*(t(i+1)-t(i));
end;

m2(j)=0;
for i=1:size(f,2)-1;
    m2(j)=m2(j)+(mean(t(i:i+1))).^2.*mean(f(j,i:i+1)).*(t(i+1)-t(i));
end;

m3(j)=0;
for i=1:size(f,2)-1;
    m3(j)=m3(j)+(mean(t(i:i+1))).^3.*mean(f(j,i:i+1)).*(t(i+1)-t(i));
end;

m4(j)=0;
for i=1:size(f,2)-1;
    m4(j)=m4(j)+(mean(t(i:i+1))).^4.*mean(f(j,i:i+1)).*(t(i+1)-t(i));
end;
end

% calculate semi-invariants
mu1=m1./m0;
mu2=m2./m0;
mu3=m3./m0;
mu4=m4./m0;

% estimate properties of interest
mat=mu1; % mean arrival time
var=sqrt(m2./m0 - (m1./m0).^2);  % plotting STANDARD DEVIATION HERE
skew=mu3./mu2.^(3/2); % skewness (tailing)
kurt=mu4./mu2.^2; % kurtosis (peakedness)