diff --git a/Project.toml b/Project.toml
index 5cdf77f..9711847 100644
--- a/Project.toml
+++ b/Project.toml
@@ -5,13 +5,14 @@ version = "1.0.0-DEV"
 
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 QuadGK = "1fd47b50-473d-5c70-9696-f719f8f3bcdc"
 Roots = "f2b01f46-fcfa-551c-844a-d8ac1e96c665"
 
 [compat]
-Roots = "2"
-QuadGK = "2"
 Documenter = "1"
+QuadGK = "2"
+Roots = "2"
 julia = "1.8"
 
 [extras]
diff --git a/src/Background.jl b/src/Background.jl
index 8398e62..9acdb60 100644
--- a/src/Background.jl
+++ b/src/Background.jl
@@ -24,11 +24,6 @@
 ##################################################################################
 
 
-# First define some basic constants (the only one we will need)
-const G_NEWTON::Float64 = 4.517103049894965e-48 # in Mpc^3 / Msun / s^(2)
-const C_LIGHT::Float64 = 9.715611890180198e-15 # Mpc / s 
-const MPC_TO_KM::Float64 = 3.085677581491367e19
-const KM_TO_MPC::Float64 = 1.0/MPC_TO_KM
 
 export BkgCosmology, FLRW, FlatFLRW, planck18_bkg, edsPlanck18_bkg, convert_cosmo
 export hubble_H0, hubble_E, hubble_E2
@@ -248,7 +243,7 @@ z_eq_Λm(bkg_cosmo::BkgCosmology = planck18_bkg) = bkg_cosmo.z_eq_Λm
 a_eq_mr(bkg_cosmo::BkgCosmology = planck18_bkg) = z_to_a(bkg_cosmo.z_eq_mr)
 a_eq_Λm(bkg_cosmo::BkgCosmology = planck18_bkg) = z_to_a(bkg_cosmo.z_eq_Λm)
 
-δt_s(a0::Real, a1::Real, bkg_cosmo::BkgCosmology = planck18_bkg; kws...) = QuadGK.quadgk(a -> 1 / hubble_E(a_to_z(a), bkg_cosmo) / a, a0, a1, rtol=1e-3; kws...)[1] / (hubble_H0(bkg_cosmo) * km_to_Mpc)
+δt_s(a0::Real, a1::Real, bkg_cosmo::BkgCosmology = planck18_bkg; kws...) = QuadGK.quadgk(a -> 1 / hubble_E(a_to_z(a), bkg_cosmo) / a, a0, a1, rtol=1e-3; kws...)[1] / (hubble_H0(bkg_cosmo) * KM_TO_MPC)
 
 """ age of the universe (s)"""
 universe_age(z::Float64 = 0.0, bkg_cosmo::BkgCosmology = planck18_bkg; kws...) = δt_s(0, z_to_a(z), bkg_cosmo; kws...)
diff --git a/src/CosmoTools.jl b/src/CosmoTools.jl
index b168093..6da8f48 100644
--- a/src/CosmoTools.jl
+++ b/src/CosmoTools.jl
@@ -18,8 +18,9 @@
 
 module CosmoTools
 
-import QuadGK, Roots
+import QuadGK, Roots, HypergeometricFunctions
 
+include("Units.jl")
 include("Background.jl")
 include("TransferFunction.jl")
 include("PowerSpectrum.jl")
diff --git a/src/Halos.jl b/src/Halos.jl
index 41cd8b3..cebf608 100644
--- a/src/Halos.jl
+++ b/src/Halos.jl
@@ -18,9 +18,9 @@
 
 
 export Halo, nfwProfile, αβγProfile, HaloProfile, coreProfile, plummerProfile
-export halo_from_ρs_and_rs, halo_from_mΔ_and_cΔ
+export halo_from_mΔ_and_cΔ
 export mΔ_from_ρs_and_rs, mΔ, rΔ_from_ρs_and_rs, rΔ, ρ_halo, μ_halo, m_halo
-export velocity_dispersion, gravitational_potential, escape_velocity
+export velocity_dispersion, gravitational_potential, escape_velocity, orbital_frequency
 
 abstract type HaloProfile{T<:Real} end
 
@@ -69,8 +69,8 @@ end
 
 ## Definition of relationships between concentration, mass, scale density and scale radius
 function cΔ_from_ρs(ρs::Real, hp::HaloProfile = nfwProfile, Δ::Real = 200, ρ_ref::Real = planck18_bkg.ρ_c0)
-    g(c::Real) = c^3 / μ(c, hp) - 3 * ρs / Δ / ρ_ref
-    Roots.find_zero(g, (1e-10, 1e+10), Bisection()) 
+    g(c::Real) = c^3 / μ_halo(c, hp) - 3 * ρs / Δ / ρ_ref
+    Roots.find_zero(g, (1e-10, 1e+10), Roots.Bisection()) 
 end
 
 mΔ_from_ρs_and_rs(ρs::Real, rs::Real, hp::HaloProfile = nfwProfile, Δ::Real = 200, ρ_ref::Real = planck18_bkg.ρ_c0) = 4 * pi * ρs * rs^3 * μ_halo(cΔ_from_ρs(ρs, hp, Δ, ρ_ref), hp)
@@ -89,7 +89,7 @@ Base.iterate(iter::Halo) = (iter, nothing)
 Base.iterate(::Halo, state::Nothing) = nothing
 
 
-halo_from_ρs_and_rs(hp::HaloProfile, ρs::Real, rs::Real) = Halo(hp, promote(ρs, rs)...)
+Halo(hp::HaloProfile, ρs::Real, rs::Real) = Halo(hp, promote(ρs, rs)...)
 
 function halo_from_mΔ_and_cΔ(hp::HaloProfile, mΔ::Real, cΔ::Real;  Δ::Real = 200, ρ_ref::Real = planck18_bkg.ρ_c0)
     ρs = ρs_from_cΔ(cΔ, hp, Δ, ρ_ref)
@@ -132,6 +132,8 @@ circular_velocity(r::Real, h::Halo) = sqrt(G_NEWTON * m(r, h) / r)
 """ circular velocity in (km / s)"""
 circular_velocity_kms(r::Real, h::Halo) = sqrt(G_NEWTON * m(r, h) / r) * MPC_TO_KM
 
+""" orbital frequency in (1 / s) """
+orbital_frequency(r::Real, h::Halo) = circular_velocity(r, h) / (2*π*r)
 
 
 
diff --git a/src/Units.jl b/src/Units.jl
new file mode 100644
index 0000000..38092a5
--- /dev/null
+++ b/src/Units.jl
@@ -0,0 +1,36 @@
+##################################################################################
+# This file is part of CosmoTools.jl
+#
+# Copyright (c) 2024, Gaétan Facchinetti
+#
+# CosmoTools.jl is free software: you can redistribute it and/or modify it 
+# under the terms of the GNU General Public License as published by 
+# the Free Software Foundation, either version 3 of the License, or any 
+# later version. CosmoTools.jl is distributed in the hope that it will be useful, 
+# but WITHOUT ANY WARRANTY; without even the implied warranty of 
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. 
+# See the GNU General Public License for more details.
+#
+# You should have received a copy of the GNU 
+# General Public License along with 21cmCAST. 
+# If not, see <https://www.gnu.org/licenses/>.
+##################################################################################
+#
+# Contains functions related to the transfer functions
+#
+# author: Gaetan Facchinetti
+# email: gaetanfacc@gmail.com
+#
+##################################################################################
+
+
+export G_NEWTON, C_LIGHT, MPC_TO_KM, KM_TO_MPC, YR_TO_S, MYR_TO_S
+
+# First define some basic constants (the only one we will need)
+const G_NEWTON::Float64 = 4.517103049894965e-48 # in Mpc^3 / Msun / s^(2)
+const C_LIGHT::Float64 = 9.715611890180198e-15 # Mpc / s 
+const MPC_TO_KM::Float64 = 3.085677581491367e19
+const KM_TO_MPC::Float64 = 1.0/MPC_TO_KM
+const YR_TO_S::Float64 = 31557600 
+const MYR_TO_S::Float64 = 31557600e+6
+