docs: 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))