function [Tp] = pert(T);
% Calculates the perturbation about the mean of the first direction.

N=size(T);
T=reshape(T,[N(1) prod(N)/N(1)]);

q=T;q(isnan(T))=0;
qm=0*q+1;qm(find(q==0))=0;
n=sum(qm,1); n(find(n==0))=1;
tx=sum(q,1)./n;
n=sum(qm,2); n(find(n==0))=1;
ty=sum(q,2)./n;
[Ty,Tx]=ndgrid(ty,tx);

Tp=reshape( T-Tx.*qm ,N);