function rho = fcn_rho(z)
global X
[Nx,~] = size(X);
[Nz,~] = size(z);
rho = 2+repmat(sin(X),[1,Nz])./ repmat(exp(z'),[Nx,1]);


