def METHOD_NAME(self, variables):
    compiled1 = sv_compile(self.formula1)
    compiled2 = sv_compile(self.formula2)
    compiled3 = sv_compile(self.formula3)
    if self.output_mode == 'XYZ':
        def out_coordinates(x, y, z):
            return x, y, z
    elif self.output_mode == 'CYL':
        def out_coordinates(rho, phi, z):
            return from_cylindrical_np(rho, phi, z, mode='radians')
    else: # SPH
        def out_coordinates(rho, phi, theta):
            return from_spherical_np(rho, phi, theta, mode='radians')
    def function(t):
        variables.update(dict(t=t))
        v1 = safe_eval_compiled(compiled1, variables, allowed_names = safe_names_np)
        v2 = safe_eval_compiled(compiled2, variables, allowed_names = safe_names_np)
        v3 = safe_eval_compiled(compiled3, variables, allowed_names = safe_names_np)
        if not isinstance(v1, np.ndarray):
            v1 = np.full_like(t, v1)
        if not isinstance(v2, np.ndarray):
            v2 = np.full_like(t, v2)
        if not isinstance(v3, np.ndarray):
            v3 = np.full_like(t, v3)
        r = np.array(out_coordinates(v1, v2, v3)).T
        return r
    return function