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)