function [] = testeh()
E = @(z) 0.23.*(1+z).^3+0.77
rhoc = @(z) E(z).*(2.775.*10.^(11))
omegam = @(z) 0.23.*(1+z).^3./(E(z))
x = @(z) omegam(z) - 1
deltac = @(z) 177.65+82.*x(z)-39.*x(z).^2
rvir = @(M,z) (3.*M./(12.566.*rhoc(z).*deltac(z))).^(1/3)
c = @(M,z) (9./(1+z)).*(M./((1.23).*10.^(12))).^(-0.13)
rhos = @(M,z) (c(M,z)).^3.*(M./(12.566.*(rvir(M,z)).^(3))).*(1./(log(1+c(M,z))-(c(M,z))./(1+c(M,z))))
rho = @(M,z,r) rhos(M,z)./((c(M,z).*r./(rvir(M,z))).*(1+(c(M,z).*r./(rvir(M,z)))).^2)
rho1 = @(r) rho(10.^14,0,r)
rho2 = @(r) rho(10.^14,1,r)
rho3 = @(r) rho(10.^15,0,r)
rho4 = @(r) rho(10.^15,1,r)
hold on
ezplot(rho1, [10.^-2 10.^2])
ezplot(rho2, [10.^-2 10.^2])
ezplot(rho3, [10.^-2 10.^2])
ezplot(rho4, [10.^-2 10.^2])
set(gca, 'XScale', 'log')
set(gca, 'YScale', 'log')
ylim([10.^5 10.^17])
