docs: def pdf(self, mag_1, mag_2, mag_err_1, mag_err_2, 
              distance_modulus=None, delta_mag=0.03, steps=10000):
          """
          Compute isochrone probability for each catalog object.
   
          ADW: This is a memory intensive function, so try as much as
          possible to keep array types at `float32` or smaller (maybe
          using add.at would be good?)
          ADW: Still a little speed to be gained here (broadcasting)
          ADW: Units? [mag^-2] [per sr?]
  
          Parameters:
          -----------
          mag_1 : magnitude of stars (pdf sample points) in first band
          mag_2 : magnitude of stars (pdf sample points) in second band
          mag_err_1 : magnitude error of stars (pdf sample points) in first band
          mag_err_2 : magnitude error of stars (pdf sample points) in second band
          distance_modulus : distance modulus of isochrone
          delta_mag : magnitude binning for evaluating the pdf
          steps : number of isochrone sample points
  
          Returns:
          --------
          u_color : probability that the star belongs to the isochrone [mag^-2]
          """
          nsigma = 5.0
          #pad = 1. # mag
  
          if distance_modulus is None: 
              distance_modulus = self.distance_modulus
  
          # ADW: HACK TO ADD SYSTEMATIC UNCERTAINTY (0.010 mag)
          mag_err_1 = np.sqrt(mag_err_1**2 + 0.01**2)
          mag_err_2 = np.sqrt(mag_err_2**2 + 0.01**2)
   
          # Binned pdf of the isochrone
          histo_pdf,bins_mag_1,bins_mag_2 = self.histogram2d(distance_modulus,delta_mag,steps)
           
          # Keep only isochrone bins that are within the magnitude
          # space of the sample
          mag_1_mesh, mag_2_mesh = np.meshgrid(bins_mag_2[1:], bins_mag_1[1:])
           
          # pdf contribution only calculated out to nsigma,
          # so padding shouldn't be necessary.
          mag_1_max = np.max(mag_1+nsigma*mag_err_1)# +pad 
          mag_1_min = np.min(mag_1-nsigma*mag_err_1)# -pad 
          mag_2_max = np.max(mag_2+nsigma*mag_err_2)# +pad 
          mag_2_min = np.min(mag_2-nsigma*mag_err_2)# -pad 
           
          in_mag_space = ((mag_1_mesh>=mag_1_min)&(mag_1_mesh<=mag_1_max))
          in_mag_space*= ((mag_2_mesh>=mag_2_min)&(mag_2_mesh<=mag_2_max))
          histo_pdf *= in_mag_space
   
          idx_mag_1, idx_mag_2 = np.nonzero(histo_pdf)
          isochrone_pdf = histo_pdf[idx_mag_1, idx_mag_2]
   
          n_catalog = len(mag_1)
          n_isochrone_bins = len(idx_mag_1)
  
          mag_1 = mag_1.reshape([n_catalog, 1])
          mag_err_1 = mag_err_1.reshape([n_catalog, 1])
          mag_2 = mag_2.reshape([n_catalog, 1])
          mag_err_2 = mag_err_2.reshape([n_catalog, 1])
  
          # Calculate (normalized) distance between each catalog object
          # and isochrone bin. Assume normally distributed photometric
          # uncertainties so that the normalized distance is:
          #   norm_dist = (mag_1 - bins_mag_1)/mag_err_1
  
          # ADW: Creating the dist arrays is memory intensive.
          # Can we cut it down (maybe with add.at)?
          dist_mag_1_hi = (mag_1-bins_mag_1[idx_mag_1])/mag_err_1
          dist_mag_1_lo = (mag_1-bins_mag_1[idx_mag_1+1])/mag_err_1
  
          dist_mag_2_hi = (mag_2-bins_mag_2[idx_mag_2])/mag_err_2
          dist_mag_2_lo = (mag_2-bins_mag_2[idx_mag_2+1])/mag_err_2
           
          # Only calculate the PDF using bins that are < nsigma from the
          # data point (i.e., where it is ~nonzero).
          idx_nonzero_0,idx_nonzero_1 = np.nonzero((dist_mag_1_hi > -nsigma) \
                                                  *(dist_mag_1_lo < nsigma)\
                                                  *(dist_mag_2_hi > -nsigma)\
                                                  *(dist_mag_2_lo < nsigma))
  
          # Now calculate the pdf as the delta of the normalized cdf
          # (more accurate than the point evaluation of the pdf)
          pdf_mag_1 = np.zeros([n_catalog, n_isochrone_bins],dtype=np.float32)
          pdf_mag_1[idx_nonzero_0,idx_nonzero_1] = norm_cdf(dist_mag_1_hi[idx_nonzero_0,idx_nonzero_1]) \
              - norm_cdf(dist_mag_1_lo[idx_nonzero_0,idx_nonzero_1])
  
          pdf_mag_2 = np.zeros([n_catalog, n_isochrone_bins],dtype=np.float32)
          pdf_mag_2[idx_nonzero_0,idx_nonzero_1] = norm_cdf(dist_mag_2_hi[idx_nonzero_0,idx_nonzero_1]) \
              - norm_cdf(dist_mag_2_lo[idx_nonzero_0,idx_nonzero_1])
  
          # Signal "color probability" (as opposed to "spatial
          # probability", but more accurately "isochrone probability")
          # is the product of PDFs for each object-bin pair summed over
          # isochrone bins 
  
          #ADW: Here is where add.at would be good...
          u_color = np.sum(pdf_mag_1 * pdf_mag_2 * isochrone_pdf, axis=1)
   
          # Remove the bin size to convert the pdf to units of mag^-2
          u_color /= delta_mag**2
  
          return u_color.astype(np.float32)