%computes the skewness of a field (up to three dimensions, skewness
%calculated over the first dimension)- NORMALIZED BY TH CUBE OF THE STD

function [out] = skew(in_field);

dim=length(find(size(in_field) > 1));

if dim==1;
    out = mean((in_field-mean(in_field)).^3)/((std(in_field))^3);
end

%the re-ordering convention here follows that used in the matlab mean
%command
if dim == 2;
    [NN,N]=size(in_field);
    for n=1:N;
        out(1,n)=mean((in_field(:,n)-mean(in_field(:,n))).^3)/((std(in_field(:,n)))^3);
    end
end

%the re-ordering convention here is to compress the time dimensions
if dim == 3
    [NT, NN,N]=size(in_field);
    for n=1:N;
        for nn=1:NN
            out(nn,n)=mean((in_field(:,nn,n)-mean(in_field(:,nn,n))).^3)/((std(in_field(:,nn,n)))^3);
        end
    end
end
