function [hk,dh_dx,dh_du,Rk,sn,T] = diag_TCV_Psib(xk,gk,vk,~,RAPTORmodel,varargin)
%#codegen

% defaults
sn   = 10e-3; % noise standard deviation

% input processing
if numel(varargin) >= 1 && ~isempty(varargin{1})
  sn = varargin{1};
end

%% Psi %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% diagnostic parameters
scal = RAPTORmodel.psi.scal; % typical scale

% get psi
[psi,dpsi_dx] = eval_psi(xk,gk,vk,RAPTORmodel,false);
yk = psi(end); dyk_dx = dpsi_dx(end,:);

% scaling matrix
T = 1/scal*eye(numel(yk));
%iT = scal*eye(numel(yk));

% h (scaled output value)
hk = T*yk; % scaled measured output

% dh_dx
dh_dx = T*dyk_dx;

% dh_du
dh_du = T*zeros(numel(hk),RAPTORmodel.dims.nu);

% noise covariance
Rk = sn.^2*T*eye(numel(hk))*T';
return
