def METHOD_NAME(self):
    # Test points:
    rng = np.random.default_rng(42)
    ntest = 4
    Rs = rng.uniform(1, 15, size=ntest) * u.kpc
    phis = rng.uniform(0, 2 * np.pi, size=ntest) * u.radian
    zs = rng.uniform(1, 15, size=ntest) * u.kpc
    cyl = CylindricalRepresentation(Rs, phis, zs)
    xyz = cyl.to_cartesian().xyz
    self.Rs = Rs.to_value(ro)
    self.phis = phis.to_value(u.rad)
    self.zs = zs.to_value(ro)
    self.Rpz_iter = np.array(list(zip(self.Rs, self.phis, self.zs))).copy()
    self.xyz = xyz.copy()
    Jac = np.zeros((len(cyl), 3, 3))
    Jac[:, 0, 0] = xyz[0] / cyl.rho
    Jac[:, 0, 1] = xyz[1] / cyl.rho
    Jac[:, 1, 0] = (-xyz[1] / cyl.rho**2).to_value(1 / ro)
    Jac[:, 1, 1] = (xyz[0] / cyl.rho**2).to_value(1 / ro)
    Jac[:, 2, 2] = 1.0
    self.Jac = Jac