 end
    
    dSdi(i,j,:,:) = Z*t - reshape(M(i)*dMdi(j,:,:) + dMdi(i,:,:)*M(j),n,D);
    dSdi(j,i,:,:) = dSdi(i,j,:,:);
    dSdX(j,i,:,:) = dSdX(i,j,:,:);
  end % j
  
  S(i,i) = S(i,i) + 1e-06;    % add small diagonal jitter for numerical reasons
end % i

dSdX(:,:,D+1,:) = -dSdX(:,:,D+2,:);
dSdX(:,:,D+1,:) = -dSdX(:,:,D+2,:);

% 4) centralize moments
S = S - M*M';
%S(diag(S)<0,diag(S)<0) = 1e-6;

% 5) Vectorize derivatives
dMds=reshape(dMds,[E D*D]);
dSdm=reshape(dSdm,[E*E D]); dSds=reshape(dSds,[E*E D*D]