PDF for the Chabrier IMF.
         
        The functional form and coefficients are described in Eq 17
        and Tab 1 of Chabrier (2003):

          m <= 1 Msun: E(log m) = A1*exp(-(log m - log m_c)^2 / 2 sigma^2)
          m  > 1 Msun: E(log m) = A2 * m^-x
         
          A1 = 1.58 : normalization [ log(Msun)^-1 pc^-3]
          m_c = 0.079 [Msun]
          sigma = 0.69
          A2 = 4.43e-2
          x = 1.3
          
        We redefine a = A1, A2 = a * b;

        The normalization is set so that the IMF integrates to 1 over
        the mass range from 0.1 Msun to 100 Msun

        Parameters:
        -----------
        mass: stellar mass (solar masses)
        log_mode[True]: return number per logarithmic mass range, i.e., dN/dlog(M)
         
        Returns:
        --------
        number per (linear or log) mass bin, i.e., dN/dM or dN/dlog(M) where mass unit is solar masses

def pdf(cls, mass, log_mode=True):
        """ PDF for the Chabrier IMF.
         
        The functional form and coefficients are described in Eq 17
        and Tab 1 of Chabrier (2003):

          m <= 1 Msun: E(log m) = A1*exp(-(log m - log m_c)^2 / 2 sigma^2)
          m  > 1 Msun: E(log m) = A2 * m^-x
         
          A1 = 1.58 : normalization [ log(Msun)^-1 pc^-3]
          m_c = 0.079 [Msun]
          sigma = 0.69
          A2 = 4.43e-2
          x = 1.3
          
        We redefine a = A1, A2 = a * b;

        The normalization is set so that the IMF integrates to 1 over
        the mass range from 0.1 Msun to 100 Msun

        Parameters:
        -----------
        mass: stellar mass (solar masses)
        log_mode[True]: return number per logarithmic mass range, i.e., dN/dlog(M)
         
        Returns:
        --------
        number per (linear or log) mass bin, i.e., dN/dM or dN/dlog(M) where mass unit is solar masses
        """
        log_mass = np.log10(mass)

        # Constants from Chabrier 2003
        m_c = 0.079
        sigma = 0.69
        x = 1.3 

        # This value is set to normalize integral from 0.1 to 100 Msun
        a=1.31357499301
        # This value is required so that the two components match at 1 Msun
        b = 0.279087531047 

        dn_dlogm = ((log_mass <= 0) * a * np.exp(-(log_mass - np.log10(m_c))**2 / (2 * (sigma**2)))) 
        dn_dlogm += ((log_mass  > 0) * a * b * mass**(-x))
        
        if log_mode:
            # Number per logarithmic mass range, i.e., dN/dlog(M)
            return dn_dlogm
        else:
            # Number per linear mass range, i.e., dN/dM
            return dn_dlogm / (mass * np.log(10))