docs: def mean_size(self, p, q):
          '''        
          >>> ds = 1E-6*np.array([240, 360, 450, 562.5, 703, 878, 1097, 1371, 1713, 2141, 2676, 3345, 4181, 5226, 6532])
          >>> numbers = [65, 119, 232, 410, 629, 849, 990, 981, 825, 579, 297, 111, 21, 1]
          >>> psd = ParticleSizeDistribution(ds=ds, fractions=numbers, order=0)
          >>> psd.mean_size(3, 2)
          0.002269321031745045
          '''
          if p != q:
              # Note: D(p, q) = D(q, p); in ISO and proven experimentally
              numerator = sum(self.di_power(i=i, power=p)*self.number_fractions[i] for i in range(self.N))
              denominator = sum(self.di_power(i=i, power=q)*self.number_fractions[i] for i in range(self.N))
              return (numerator/denominator)**(1.0/(p-q))
          else:
              numerator = sum(log(self.di_power(i=i, power=1))*self.di_power(i=i, power=p)*self.number_fractions[i] for i in range(self.N))
              denominator = sum(self.di_power(i=i, power=q)*self.number_fractions[i] for i in range(self.N))
              return exp(numerator/denominator)