class Cosmology(object):
    """
    A class that nicely deals with cosmological parameters.
    
    Most cosmological parameters are merely input and are made available as
    attributes in the class. However, more complicated relations such as 
    the interrelation of omegab, omegac, omegam, omegav for example are dealt
    with in a more robust manner. 
    
    The secondary purpose of this class is to provide simple mappings of the
    parameters to common python cosmology libraries (for now just `cosmolopy`
    and `pycamb`). It has been found by the authors that using more than one
    library becomes confusing in terms of naming all the parameters, so this 
    class helps deal with that.
    
    .. note :: Currently, only several combinations of the density parameters
            are valid:
    
            1. ``omegab`` and ``omegac``
            #. ``omegam`` 
            #. ``omegab_h2`` and ``omegac_h2``
            #. None of them
    
            To this one may add ``omegav`` (dark-energy density) at will. More
            combinations will be added in the future. 
    
    Parameters
    ----------
    
    default : str, {``None``, ``"planck1_base"``}
        Defines a set of default parameters, based on a published set from WMAP
        or Planck. These defaults are applied in a smart way, so as not to 
        override any user-set parameters. 
        
        Current options are
        
        1. ``None``
        #. ``"planck1_base"``: The cosmology of first-year PLANCK mission (with no lensing or WP)
    
    force_flat : bool, default ``False``
        If ``True``, enforces a flat cosmology :math:`(\Omega_m+\Omega_\lambda=1)`.
        This will modify ``omegav`` only, never ``omegam``.
        
    \*\*kwargs : 
        The list of available keyword arguments is as follows:
        
        1. ``sigma_8``: The normalisation. Mass variance in top-hat spheres with :math:`R=8Mpc h^{-1}`
        #. ``n``: The spectral index
        #. ``w``: The dark-energy equation of state
        #. ``cs2_lam``: The constant comoving sound speed of dark energy
        #. ``t_cmb``: Temperature of the CMB
        #. ``y_he``: Helium fraction
        #. ``N_nu``: Number of massless neutrino species
        #. ``N_nu_massive``:Number of massive neutrino species
        #. ``z_reion``: Redshift of reionization
        #. ``tau``: Optical depth at reionization
        #. ``delta_c``: The critical overdensity for collapse
        #. ``h``: The hubble parameter
        #. ``H0``: The hubble constant
        #. ``omegan``: The normalised density of neutrinos
        #. ``omegam``: The normalised density of matter
        #. ``omegav``: The normalised density of dark energy
        #. ``omegab``: The normalised baryon density
        #. ``omegac``: The normalised CDM density
        #. ``omegab_h2``: The normalised baryon density by ``h**2``
        #. ``omegac_h2``: The normalised CDM density by ``h**2``
        
        .. note :: The reason these are implemented as `kwargs` rather than the
                usual arguments, is because the code can't tell *a priori* which
                combination of density parameters the user will input.
                 
    """
    # A dictionary of bounds for each parameter
    # This also forms a list of all parameters possible
    # Note that just because a parameter is within these bounds doesn't mean
    # it will actually work in say CAMB.
    _bounds = {"sigma_8":[0.1, 10],
              "n":[-3, 4],
              "w":[-1.5, 0],
              "cs2_lam":[-1, 2],
              "t_cmb":[0, 10.0],
              "y_he":[0, 1],
              "N_nu":[1, 10],
              "N_nu_massive":[0, 3],
              "z_reion":[2, 1000],
              "tau":[0, 1],
              "omegan":[0, 1],
              "h":[0.05, 5],
              "H0":[5, 500],
              "omegab":[0.0001, 1],
              "omegac":[0, 2],
              "omegav":[0, 2],
              "omegam":[0.0001, 3],
              "omegab_h2":[0.0001, 1],
              "omegac_h2":[0, 2]}

    def __init__(self, default=None, force_flat=False, **kwargs):

        # Map the 'default' cosmology to its dictionary
        self.default = default
        if default == "planck1_base":
            self.__base = dict(planck1_base, **extras)

        # Set some simple parameters
        self.force_flat = force_flat
        self.crit_dens = 27.755e10

        #=======================================================================
        # Check values in kwargs
        #=======================================================================
        for k in kwargs:
            if k not in Cosmology._bounds:
                raise ValueError(k + " is not a valid parameter for Cosmology")

        #=======================================================================
        # First set the "easy" values (no dependence on anything else
        #=======================================================================
        easy_params = ["sigma_8", "n", 'w', 'cs2_lam', 't_cmb', 'y_he', "N_nu",
                       "z_reion", "tau", "omegan", 'delta_c', "N_nu_massive"]
        for p in easy_params:
            if p in kwargs:
                self.__dict__.update({p:kwargs.pop(p)})
            elif default is not None:
                self.__dict__.update({p:self.__base[p]})


        #=======================================================================
        # Now the hard parameters (multi-dependent)
        #=======================================================================
        ################### h/H0 ###############################################
        if "h" in kwargs and "H0" in kwargs:
            if kwargs['h'] != kwargs["H0"] / 100.0:
                print "h and H0 specified inconsistently, using h"



        if "H0" in kwargs:
            self.H0 = kwargs.pop("H0")
            self.h = self.H0 / 100.0

        if "h" in kwargs:
            self.h = kwargs.pop("h")
            self.H0 = 100 * self.h

        if not hasattr(self, "h") and default is not None:
            self.H0 = self.__base["H0"]
            self.h = self.H0 / 100.0

        ################### The omegas #########################################
        if "omegav" in kwargs:
            self.omegav = kwargs.pop("omegav")

        if len(kwargs) == 0: