%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function J = calcJ(Wf,hT,K);
global Rx dsdshT
% calculate the error term for diagonalization
J = 0;
for k = 1:K
Eomega = fast_ip(fast_ip(Wf,Rx(:,:,:,k)),hermite3(Wf));
Eomega(dsdshT) = 0; % this uses the optimal lambdaS
% note that trace(A*A') == A(:)' * A(:)
J = J + Eomega(:)' * Eomega(:);
end
J = real(J)/K/hT;