Repository: sympy/sympy

Issue:
Constraint solver not used for auxiliary equations in KanesMethod
As can be seen in [KanesMethod.kanes_equations](https://github.com/sympy/sympy/blob/cba1e7c6379eac88a469c98b8d86e324c06d1831/sympy/physics/mechanics/kane.py#L725) the constraint solver argument passed to `KanesMethod.__init__` is not used when forming the auxiliary equations.

To solve this issue I would propose to add a `self._constraint_solver` attribute, such that we can retrieve the passed constraint solver in `KanesMethod.kanes_equations`.

Discussion:

Patch:
diff --git a/.mailmap b/.mailmap
index 444c559f4a..9769de1266 100644
--- a/.mailmap
+++ b/.mailmap
@@ -1174,7 +1174,6 @@ Remco de Boer <29308176+redeboer@users.noreply.github.com> Remco de Boer <redebo
 Renato Coutinho <renato.coutinho@gmail.com>
 Renato Orsino <renato.orsino@gmail.com>
 Rhea Parekh <rheaparekh12@gmail.com>
-Riccardo Di Girolamo <riccardodigirolamo01@gmail.com>
 Riccardo Gori <goriccardo@gmail.com>
 Riccardo Magliocchetti <riccardo.magliocchetti@gmail.com>
 Rich LaSota <rjlasota@gmail.com>
diff --git a/doc/src/modules/physics/mechanics/api/deprecated_classes.rst b/doc/src/modules/physics/mechanics/api/deprecated_classes.rst
deleted file mode 100644
index 346385bda8..0000000000
--- a/doc/src/modules/physics/mechanics/api/deprecated_classes.rst
+++ /dev/null
@@ -1,17 +0,0 @@
-===============================
-Deprecated Classes (Docstrings)
-===============================
-
-.. deprecated:: 1.13
-    :class:`~.Body` and :class:`~.JointsMethod` have been deprecated. The
-    functionality of :class:`~.Body` is fully captured by :class:`~.RigidBody`
-    and :class:`~.Particle` and the functionality of :class:`~.JointsMethod` is
-    fully captured by :class:`~.System`.
-
-.. autoclass:: sympy.physics.mechanics.body.Body
-    :members:
-    :inherited-members:
-
-.. autoclass:: sympy.physics.mechanics.jointsmethod.JointsMethod
-    :members:
-    :inherited-members:
diff --git a/doc/src/modules/physics/mechanics/api/index.rst b/doc/src/modules/physics/mechanics/api/index.rst
index 63608ffb26..08583c2c06 100644
--- a/doc/src/modules/physics/mechanics/api/index.rst
+++ b/doc/src/modules/physics/mechanics/api/index.rst
@@ -14,4 +14,3 @@ Mechanics API Reference
    pathway.rst
    actuator.rst
    wrapping_geometry.rst
-   deprecated_classes.rst
diff --git a/sympy/geometry/ellipse.py b/sympy/geometry/ellipse.py
index bd4902e5f5..f096556f4b 100644
--- a/sympy/geometry/ellipse.py
+++ b/sympy/geometry/ellipse.py
@@ -748,14 +748,27 @@ def is_tangent(self, o):
                 return True
             # might return None if it can't decide
             return hit[0].equals(hit[1])
-        elif isinstance(o, (Segment2D, Ray2D)):
+        elif isinstance(o, Ray2D):
             intersect = self.intersection(o)
             if len(intersect) == 1:
-                return o in self.tangent_lines(intersect[0])[0]
+                return intersect[0] != o.source and not self.encloses_point(o.source)
             else:
                 return False
-        elif isinstance(o, Polygon):
-            return all(self.is_tangent(s) for s in o.sides)
+        elif isinstance(o, (Segment2D, Polygon)):
+            all_tangents = False
+            segments = o.sides if isinstance(o, Polygon) else [o]
+            for segment in segments:
+                intersect = self.intersection(segment)
+                if len(intersect) == 1:
+                    if not any(intersect[0] in i for i in segment.points) \
+                        and not any(self.encloses_point(i) for i in segment.points):
+                        all_tangents = True
+                        continue
+                    else:
+                        return False
+                else:
+                    return False
+            return all_tangents
         elif isinstance(o, (LinearEntity3D, Point3D)):
             raise TypeError('Entity must be two dimensional, not three dimensional')
         else:
diff --git a/sympy/integrals/heurisch.py b/sympy/integrals/heurisch.py
index 2a1b61c27d..344edf250a 100644
--- a/sympy/integrals/heurisch.py
+++ b/sympy/integrals/heurisch.py
@@ -1,8 +1,7 @@
 from __future__ import annotations
 
-from collections import defaultdict
-from functools import reduce
 from itertools import permutations
+from functools import reduce
 
 from sympy.core.add import Add
 from sympy.core.basic import Basic
@@ -504,16 +503,7 @@ def heurisch(f, x, rewrite=False, hints=None, mappings=None, retries=3,
         # optimizing the number of permutations of mapping              #
         assert mapping[-1][0] == x # if not, find it and correct this comment
         unnecessary_permutations = [mapping.pop(-1)]
-        # only permute types of objects and let the ordering
-        # of types take care of the order of replacement
-        types = defaultdict(list)
-        for i in mapping:
-            types[type(i)].append(i)
-        mapping = [types[i] for i in types]
-        def _iter_mappings():
-            for i in permutations(mapping):
-                yield [j for i in i for j in i]
-        mappings = _iter_mappings()
+        mappings = permutations(mapping)
     else:
         unnecessary_permutations = unnecessary_permutations or []
 
diff --git a/sympy/parsing/sympy_parser.py b/sympy/parsing/sympy_parser.py
index f32dcc4f2c..c9fba7634b 100644
--- a/sympy/parsing/sympy_parser.py
+++ b/sympy/parsing/sympy_parser.py
@@ -1138,7 +1138,9 @@ def visit_Compare(self, node):
             new_node = ast.Call(
                 func=ast.Name(id=sympy_class, ctx=ast.Load()),
                 args=[left, right],
-                keywords=[ast.keyword(arg='evaluate', value=ast.Constant(value=False))]
+                keywords=[ast.keyword(arg='evaluate', value=ast.Constant(value=False, ctx=ast.Load()))],
+                starargs=None,
+                kwargs=None
             )
             return new_node
         return node
@@ -1169,7 +1171,9 @@ def visit_BinOp(self, node):
                 right = ast.Call(
                     func=ast.Name(id='Mul', ctx=ast.Load()),
                     args=[ast.UnaryOp(op=ast.USub(), operand=ast.Constant(1)), right],
-                    keywords=[ast.keyword(arg='evaluate', value=ast.Constant(value=False))]
+                    keywords=[ast.keyword(arg='evaluate', value=ast.Constant(value=False, ctx=ast.Load()))],
+                    starargs=None,
+                    kwargs=None
                 )
             elif isinstance(node.op, ast.Div):
                 if isinstance(node.left, ast.UnaryOp):
@@ -1178,13 +1182,17 @@ def visit_BinOp(self, node):
                     left = ast.Call(
                     func=ast.Name(id='Pow', ctx=ast.Load()),
                     args=[left, ast.UnaryOp(op=ast.USub(), operand=ast.Constant(1))],
-                    keywords=[ast.keyword(arg='evaluate', value=ast.Constant(value=False))]
+                    keywords=[ast.keyword(arg='evaluate', value=ast.Constant(value=False, ctx=ast.Load()))],
+                    starargs=None,
+                    kwargs=None
                 )
                 else:
                     right = ast.Call(
                     func=ast.Name(id='Pow', ctx=ast.Load()),
                     args=[right, ast.UnaryOp(op=ast.USub(), operand=ast.Constant(1))],
-                    keywords=[ast.keyword(arg='evaluate', value=ast.Constant(value=False))]
+                    keywords=[ast.keyword(arg='evaluate', value=ast.Constant(value=False, ctx=ast.Load()))],
+                    starargs=None,
+                    kwargs=None
                 )
 
             if rev:  # undo reversal
@@ -1192,7 +1200,9 @@ def visit_BinOp(self, node):
             new_node = ast.Call(
                 func=ast.Name(id=sympy_class, ctx=ast.Load()),
                 args=[left, right],
-                keywords=[ast.keyword(arg='evaluate', value=ast.Constant(value=False))]
+                keywords=[ast.keyword(arg='evaluate', value=ast.Constant(value=False, ctx=ast.Load()))],
+                starargs=None,
+                kwargs=None
             )
 
             if sympy_class in ('Add', 'Mul'):
@@ -1205,7 +1215,7 @@ def visit_BinOp(self, node):
     def visit_Call(self, node):
         new_node = self.generic_visit(node)
         if isinstance(node.func, ast.Name) and node.func.id in self.functions:
-            new_node.keywords.append(ast.keyword(arg='evaluate', value=ast.Constant(value=False)))
+            new_node.keywords.append(ast.keyword(arg='evaluate', value=ast.Constant(value=False, ctx=ast.Load())))
         return new_node
 
 
diff --git a/sympy/physics/mechanics/body.py b/sympy/physics/mechanics/body.py
index efc367158b..310610d059 100644
--- a/sympy/physics/mechanics/body.py
+++ b/sympy/physics/mechanics/body.py
@@ -16,10 +16,6 @@ class Body(RigidBody, Particle):  # type: ignore
     passed in and central_inertia is left as None, the Particle object is
     created. Otherwise a RigidBody object will be created.
 
-    .. deprecated:: 1.13
-        The Body class is deprecated. Its functionality is captured by
-        :class:`~.RigidBody` and :class:`~.Particle`.
-
     Explanation
     ===========
 
@@ -69,20 +65,12 @@ class Body(RigidBody, Particle):  # type: ignore
     Examples
     ========
 
-    As Body has been deprecated, the following examples are for illustrative
-    purposes only. The functionality of Body is fully captured by
-    :class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
-    warning we can use the ignore_warnings context manager.
-
-        >>> from sympy.utilities.exceptions import ignore_warnings
-
     Default behaviour. This results in the creation of a RigidBody object for
     which the mass, mass center, frame and inertia attributes are given default
     values. ::
 
         >>> from sympy.physics.mechanics import Body
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     body = Body('name_of_body')
+        >>> body = Body('name_of_body')
 
     This next example demonstrates the code required to specify all of the
     values of the Body object. Note this will also create a RigidBody version of
@@ -96,8 +84,7 @@ class Body(RigidBody, Particle):  # type: ignore
         >>> frame = ReferenceFrame('frame')
         >>> ixx = Symbol('ixx')
         >>> body_inertia = inertia(frame, ixx, 0, 0)
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     body = Body('name_of_body', masscenter, mass, frame, body_inertia)
+        >>> body = Body('name_of_body', masscenter, mass, frame, body_inertia)
 
     The minimal code required to create a Particle version of the Body object
     involves simply passing in a name and a mass. ::
@@ -105,8 +92,7 @@ class Body(RigidBody, Particle):  # type: ignore
         >>> from sympy import Symbol
         >>> from sympy.physics.mechanics import Body
         >>> mass = Symbol('mass')
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     body = Body('name_of_body', mass=mass)
+        >>> body = Body('name_of_body', mass=mass)
 
     The Particle version of the Body object can also receive a masscenter point
     and a reference frame, just not an inertia.
@@ -218,19 +204,12 @@ def kinetic_energy(self, frame):
         Examples
         ========
 
-        As Body has been deprecated, the following examples are for illustrative
-        purposes only. The functionality of Body is fully captured by
-        :class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
-        warning we can use the ignore_warnings context manager.
-
-        >>> from sympy.utilities.exceptions import ignore_warnings
         >>> from sympy.physics.mechanics import Body, ReferenceFrame, Point
         >>> from sympy import symbols
         >>> m, v, r, omega = symbols('m v r omega')
         >>> N = ReferenceFrame('N')
         >>> O = Point('O')
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     P = Body('P', masscenter=O, mass=m)
+        >>> P = Body('P', masscenter=O, mass=m)
         >>> P.masscenter.set_vel(N, v * N.y)
         >>> P.kinetic_energy(N)
         m*v**2/2
@@ -240,8 +219,7 @@ def kinetic_energy(self, frame):
         >>> b.set_ang_vel(N, omega * b.x)
         >>> P = Point('P')
         >>> P.set_vel(N, v * N.x)
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     B = Body('B', masscenter=P, frame=b)
+        >>> B = Body('B', masscenter=P, frame=b)
         >>> B.kinetic_energy(N)
         B_ixx*omega**2/2 + B_mass*v**2/2
 
@@ -286,17 +264,10 @@ def apply_force(self, force, point=None, reaction_body=None, reaction_point=None
         Example
         =======
 
-        As Body has been deprecated, the following examples are for illustrative
-        purposes only. The functionality of Body is fully captured by
-        :class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
-        warning we can use the ignore_warnings context manager.
-
-        >>> from sympy.utilities.exceptions import ignore_warnings
         >>> from sympy import symbols
         >>> from sympy.physics.mechanics import Body, Point, dynamicsymbols
         >>> m, g = symbols('m g')
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     B = Body('B')
+        >>> B = Body('B')
         >>> force1 = m*g*B.z
         >>> B.apply_force(force1) #Applying force on B's masscenter
         >>> B.loads
@@ -321,12 +292,10 @@ def apply_force(self, force, point=None, reaction_body=None, reaction_point=None
         consider two bodies connected through a spring.
 
         >>> from sympy.physics.mechanics import Body, dynamicsymbols
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     N = Body('N') #Newtonion Frame
+        >>> N = Body('N') #Newtonion Frame
         >>> x = dynamicsymbols('x')
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     B1 = Body('B1')
-        ...     B2 = Body('B2')
+        >>> B1 = Body('B1')
+        >>> B2 = Body('B2')
         >>> spring_force = x*N.x
 
         Now let's apply equal and opposite spring force to the bodies.
@@ -393,17 +362,10 @@ def apply_torque(self, torque, reaction_body=None):
         Example
         =======
 
-        As Body has been deprecated, the following examples are for illustrative
-        purposes only. The functionality of Body is fully captured by
-        :class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
-        warning we can use the ignore_warnings context manager.
-
-        >>> from sympy.utilities.exceptions import ignore_warnings
         >>> from sympy import symbols
         >>> from sympy.physics.mechanics import Body, dynamicsymbols
         >>> t = symbols('t')
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     B = Body('B')
+        >>> B = Body('B')
         >>> torque1 = t*B.z
         >>> B.apply_torque(torque1)
         >>> B.loads
@@ -427,10 +389,9 @@ def apply_torque(self, torque, reaction_body=None):
         a torque `T` is acting on one body, and `-T` on the other.
 
         >>> from sympy.physics.mechanics import Body, dynamicsymbols
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     N = Body('N') #Newtonion frame
-        ...     B1 = Body('B1')
-        ...     B2 = Body('B2')
+        >>> N = Body('N') #Newtonion frame
+        >>> B1 = Body('B1')
+        >>> B2 = Body('B2')
         >>> v = dynamicsymbols('v')
         >>> T = v*N.y #Torque
 
@@ -473,15 +434,8 @@ def clear_loads(self):
         Example
         =======
 
-        As Body has been deprecated, the following examples are for illustrative
-        purposes only. The functionality of Body is fully captured by
-        :class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
-        warning we can use the ignore_warnings context manager.
-
-        >>> from sympy.utilities.exceptions import ignore_warnings
         >>> from sympy.physics.mechanics import Body
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     B = Body('B')
+        >>> B = Body('B')
         >>> force = B.x + B.y
         >>> B.apply_force(force)
         >>> B.loads
@@ -510,15 +464,8 @@ def remove_load(self, about=None):
         Example
         =======
 
-        As Body has been deprecated, the following examples are for illustrative
-        purposes only. The functionality of Body is fully captured by
-        :class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
-        warning we can use the ignore_warnings context manager.
-
-        >>> from sympy.utilities.exceptions import ignore_warnings
         >>> from sympy.physics.mechanics import Body, Point
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     B = Body('B')
+        >>> B = Body('B')
         >>> P = Point('P')
         >>> f1 = B.x
         >>> f2 = B.y
@@ -558,16 +505,9 @@ def masscenter_vel(self, body):
         Example
         =======
 
-        As Body has been deprecated, the following examples are for illustrative
-        purposes only. The functionality of Body is fully captured by
-        :class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
-        warning we can use the ignore_warnings context manager.
-
-        >>> from sympy.utilities.exceptions import ignore_warnings
         >>> from sympy.physics.mechanics import Body
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     A = Body('A')
-        ...     B = Body('B')
+        >>> A = Body('A')
+        >>> B = Body('B')
         >>> A.masscenter.set_vel(B.frame, 5*B.frame.x)
         >>> A.masscenter_vel(B)
         5*B_frame.x
@@ -596,18 +536,10 @@ def ang_vel_in(self, body):
         Example
         =======
 
-        As Body has been deprecated, the following examples are for illustrative
-        purposes only. The functionality of Body is fully captured by
-        :class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
-        warning we can use the ignore_warnings context manager.
-
-        >>> from sympy.utilities.exceptions import ignore_warnings
         >>> from sympy.physics.mechanics import Body, ReferenceFrame
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     A = Body('A')
+        >>> A = Body('A')
         >>> N = ReferenceFrame('N')
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     B = Body('B', frame=N)
+        >>> B = Body('B', frame=N)
         >>> A.frame.set_ang_vel(N, 5*N.x)
         >>> A.ang_vel_in(B)
         5*N.x
@@ -636,16 +568,9 @@ def dcm(self, body):
         Example
         =======
 
-        As Body has been deprecated, the following examples are for illustrative
-        purposes only. The functionality of Body is fully captured by
-        :class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
-        warning we can use the ignore_warnings context manager.
-
-        >>> from sympy.utilities.exceptions import ignore_warnings
         >>> from sympy.physics.mechanics import Body
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     A = Body('A')
-        ...     B = Body('B')
+        >>> A = Body('A')
+        >>> B = Body('B')
         >>> A.frame.orient_axis(B.frame, B.frame.x, 5)
         >>> A.dcm(B)
         Matrix([
@@ -688,15 +613,8 @@ def parallel_axis(self, point, frame=None):
         Example
         =======
 
-        As Body has been deprecated, the following examples are for illustrative
-        purposes only. The functionality of Body is fully captured by
-        :class:`~.RigidBody` and :class:`~.Particle`. To ignore the deprecation
-        warning we can use the ignore_warnings context manager.
-
-        >>> from sympy.utilities.exceptions import ignore_warnings
         >>> from sympy.physics.mechanics import Body
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     A = Body('A')
+        >>> A = Body('A')
         >>> P = A.masscenter.locatenew('point', 3 * A.x + 5 * A.y)
         >>> A.parallel_axis(P).to_matrix(A.frame)
         Matrix([
diff --git a/sympy/physics/mechanics/jointsmethod.py b/sympy/physics/mechanics/jointsmethod.py
index df7bd56360..7cc312366f 100644
--- a/sympy/physics/mechanics/jointsmethod.py
+++ b/sympy/physics/mechanics/jointsmethod.py
@@ -11,10 +11,6 @@
 class JointsMethod(_Methods):
     """Method for formulating the equations of motion using a set of interconnected bodies with joints.
 
-    .. deprecated:: 1.13
-        The JointsMethod class is deprecated. Its functionality has been
-        replaced by the new :class:`~.System` class.
-
     Parameters
     ==========
 
@@ -49,14 +45,6 @@ class JointsMethod(_Methods):
     Examples
     ========
 
-    As Body and JointsMethod have been deprecated, the following examples are
-    for illustrative purposes only. The functionality of Body is fully captured
-    by :class:`~.RigidBody` and :class:`~.Particle` and the functionality of
-    JointsMethod is fully captured by :class:`~.System`. To ignore the
-    deprecation warning we can use the ignore_warnings context manager.
-
-    >>> from sympy.utilities.exceptions import ignore_warnings
-
     This is a simple example for a one degree of freedom translational
     spring-mass-damper.
 
@@ -65,14 +53,12 @@ class JointsMethod(_Methods):
     >>> from sympy.physics.vector import dynamicsymbols
     >>> c, k = symbols('c k')
     >>> x, v = dynamicsymbols('x v')
-    >>> with ignore_warnings(DeprecationWarning):
-    ...     wall = Body('W')
-    ...     body = Body('B')
+    >>> wall = Body('W')
+    >>> body = Body('B')
     >>> J = PrismaticJoint('J', wall, body, coordinates=x, speeds=v)
     >>> wall.apply_force(c*v*wall.x, reaction_body=body)
     >>> wall.apply_force(k*x*wall.x, reaction_body=body)
-    >>> with ignore_warnings(DeprecationWarning):
-    ...     method = JointsMethod(wall, J)
+    >>> method = JointsMethod(wall, J)
     >>> method.form_eoms()
     Matrix([[-B_mass*Derivative(v(t), t) - c*v(t) - k*x(t)]])
     >>> M = method.mass_matrix_full
@@ -240,15 +226,6 @@ def form_eoms(self, method=KanesMethod):
         Examples
         ========
 
-        As Body and JointsMethod have been deprecated, the following examples
-        are for illustrative purposes only. The functionality of Body is fully
-        captured by :class:`~.RigidBody` and :class:`~.Particle` and the
-        functionality of JointsMethod is fully captured by :class:`~.System`. To
-        ignore the deprecation warning we can use the ignore_warnings context
-        manager.
-
-        >>> from sympy.utilities.exceptions import ignore_warnings
-
         This is a simple example for a one degree of freedom translational
         spring-mass-damper.
 
@@ -258,14 +235,12 @@ def form_eoms(self, method=KanesMethod):
         >>> q = dynamicsymbols('q')
         >>> qd = dynamicsymbols('q', 1)
         >>> m, k, b = symbols('m k b')
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     wall = Body('W')
-        ...     part = Body('P', mass=m)
+        >>> wall = Body('W')
+        >>> part = Body('P', mass=m)
         >>> part.potential_energy = k * q**2 / S(2)
         >>> J = PrismaticJoint('J', wall, part, coordinates=q, speeds=qd)
         >>> wall.apply_force(b * qd * wall.x, reaction_body=part)
-        >>> with ignore_warnings(DeprecationWarning):
-        ...     method = JointsMethod(wall, J)
+        >>> method = JointsMethod(wall, J)
         >>> method.form_eoms(LagrangesMethod)
         Matrix([[b*Derivative(q(t), t) + k*q(t) + m*Derivative(q(t), (t, 2))]])
 
diff --git a/sympy/physics/mechanics/kane.py b/sympy/physics/mechanics/kane.py
index d8266d7300..3ee73d3260 100644
--- a/sympy/physics/mechanics/kane.py
+++ b/sympy/physics/mechanics/kane.py
@@ -227,7 +227,7 @@ def __init__(self, frame, q_ind, u_ind, kd_eqs=None, q_dependent=None,
         self._bodylist = bodies
 
         self.explicit_kinematics = explicit_kinematics
-
+        self._constraint_solver = constraint_solver
         self._initialize_vectors(q_ind, q_dependent, u_ind, u_dependent,
                 u_auxiliary)
         _validate_coordinates(self.q, self.u)
@@ -720,14 +720,15 @@ def kanes_equations(self, bodies=None, loads=None):
         if self._uaux:
             if not self._udep:
                 km = KanesMethod(self._inertial, self.q, self._uaux,
-                             u_auxiliary=self._uaux)
+                             u_auxiliary=self._uaux, constraint_solver=self._constraint_solver)
             else:
                 km = KanesMethod(self._inertial, self.q, self._uaux,
                         u_auxiliary=self._uaux, u_dependent=self._udep,
                         velocity_constraints=(self._k_nh * self.u +
                         self._f_nh),
                         acceleration_constraints=(self._k_dnh * self._udot +
-                        self._f_dnh)
+                        self._f_dnh),
+                        constraint_solver=self._constraint_solver
                         )
             km._qdot_u_map = self._qdot_u_map
             self._km = km
diff --git a/sympy/physics/vector/functions.py b/sympy/physics/vector/functions.py
index 95d9efb72d..b0291cfde0 100644
--- a/sympy/physics/vector/functions.py
+++ b/sympy/physics/vector/functions.py
@@ -1,7 +1,7 @@
 from functools import reduce
 
 from sympy import (sympify, diff, sin, cos, Matrix, symbols,
-                                Function, S, Symbol, linear_eq_to_matrix)
+                                Function, S, Symbol)
 from sympy.integrals.integrals import integrate
 from sympy.simplify.trigsimp import trigsimp
 from .vector import Vector, _check_vector
@@ -573,16 +573,10 @@ def partial_velocity(vel_vecs, gen_speeds, frame):
         raise TypeError('Generalized speeds must be contained in an iterable')
 
     vec_partials = []
-    gen_speeds = list(gen_speeds)
-    for vel in vel_vecs:
-        partials = [Vector(0) for _ in gen_speeds]
-        for components, ref in vel.args:
-            mat, _ = linear_eq_to_matrix(components, gen_speeds)
-            for i in range(len(gen_speeds)):
-                for dim, direction in enumerate(ref):
-                    if mat[dim, i] != 0:
-                        partials[i] += direction * mat[dim, i]
-
+    for vec in vel_vecs:
+        partials = []
+        for speed in gen_speeds:
+            partials.append(vec.diff(speed, frame, var_in_dcm=False))
         vec_partials.append(partials)
 
     return vec_partials
diff --git a/sympy/printing/c.py b/sympy/printing/c.py
index 21d015bde1..5a1e677496 100644
--- a/sympy/printing/c.py
+++ b/sympy/printing/c.py
@@ -13,7 +13,6 @@
 
 from __future__ import annotations
 from typing import Any
-import re
 
 from functools import wraps
 from itertools import chain
@@ -89,7 +88,6 @@
 
 reserved_words_c99 = ['inline', 'restrict']
 
-
 def get_math_macros():
     """ Returns a dictionary with math-related macros from math.h/cmath
 
@@ -269,13 +267,6 @@ def _traverse_matrix_indices(self, mat):
         rows, cols = mat.shape
         return ((i, j) for i in range(rows) for j in range(cols))
 
-    @staticmethod
-    def _check_valid_c_variable_name(var_name):
-        pattern = re.compile(r"^[a-zA-Z_][a-zA-Z0-9_]{0,30}$")
-        msg = 'SymPy symbol {} not printable in C.'
-        if not pattern.search(var_name):
-            raise ValueError(msg.format(var_name))
-
     @_as_macro_if_defined
     def _print_Mul(self, expr, **kwargs):
         return super()._print_Mul(expr, **kwargs)
@@ -396,8 +387,6 @@ def _print_MatrixElement(self, expr):
 
     def _print_Symbol(self, expr):
         name = super()._print_Symbol(expr)
-        # if name is an indexed, e.g. a[0], only check the variable name
-        self._check_valid_c_variable_name(name.split('[')[0])
         if expr in self._settings['dereference']:
             return '(*{})'.format(name)
         else:


Test patch:
diff --git a/sympy/geometry/tests/test_ellipse.py b/sympy/geometry/tests/test_ellipse.py
index a79eba8c35..e0e5b61180 100644
--- a/sympy/geometry/tests/test_ellipse.py
+++ b/sympy/geometry/tests/test_ellipse.py
@@ -452,8 +452,6 @@ def test_is_tangent():
     assert c1.is_tangent(Ray((-3, -2), (-15, -20))) is False
     assert c1.is_tangent(Ray((-3, -22), (15, 20))) is False
     assert c1.is_tangent(Ray((9, 20), (9, -20))) is True
-    assert c1.is_tangent(Ray((2, 5), (9, 5))) is True
-    assert c1.is_tangent(Segment((2, 5), (9, 5))) is True
     assert e1.is_tangent(Segment((2, 2), (-7, 7))) is False
     assert e1.is_tangent(Segment((0, 0), (1, 2))) is False
     assert c1.is_tangent(Segment((0, 0), (-5, -2))) is False
diff --git a/sympy/integrals/tests/test_heurisch.py b/sympy/integrals/tests/test_heurisch.py
index e220132f58..2b4ffa0684 100644
--- a/sympy/integrals/tests/test_heurisch.py
+++ b/sympy/integrals/tests/test_heurisch.py
@@ -2,7 +2,7 @@
 from sympy.core.add import Add
 from sympy.core.function import (Derivative, Function, diff)
 from sympy.core.numbers import (I, Rational, pi)
-from sympy.core.relational import Eq, Ne
+from sympy.core.relational import Ne
 from sympy.core.symbol import (Symbol, symbols)
 from sympy.functions.elementary.exponential import (LambertW, exp, log)
 from sympy.functions.elementary.hyperbolic import (asinh, cosh, sinh, tanh)
@@ -12,8 +12,6 @@
 from sympy.functions.special.bessel import (besselj, besselk, bessely, jn)
 from sympy.functions.special.error_functions import erf
 from sympy.integrals.integrals import Integral
-from sympy.logic.boolalg import And
-from sympy.matrices import Matrix
 from sympy.simplify.ratsimp import ratsimp
 from sympy.simplify.simplify import simplify
 from sympy.integrals.heurisch import components, heurisch, heurisch_wrapper
@@ -367,33 +365,3 @@ def f(x):
     Uz = integrate(f(z), z)
     Ut = integrate(f(t), t)
     assert Ut == Uz.subs(z, t)
-
-
-def test_heurisch_complex_erf_issue_26338():
-    r = symbols('r', real=True)
-    a = exp(-r**2/(2*(2 - I)**2))
-    assert heurisch(a, r, hints=[]) is None  # None, not a wrong soln
-    a = sqrt(pi)*erf((1 + I)/2)/2
-    assert integrate(exp(-I*r**2/2), (r, 0, 1)) == a - I*a
-
-    a = exp(-x**2/(2*(2 - I)**2))
-    assert heurisch(a, x, hints=[]) is None  # None, not a wrong soln
-    a = sqrt(pi)*erf((1 + I)/2)/2
-    assert integrate(exp(-I*x**2/2), (x, 0, 1)) == a - I*a
-
-
-def test_issue_15498():
-    Z0 = Function('Z0')
-    k01, k10, t, s= symbols('k01 k10 t s', real=True, positive=True)
-    m = Matrix([[exp(-k10*t)]])
-    _83 = Rational(83, 100)  # 0.83 works, too
-    [a, b, c, d, e, f, g] = [100, 0.5, _83, 50, 0.6, 2, 120]
-    AIF_btf = a*(d*e*(1 - exp(-(t - b)/e)) + f*g*(1 - exp(-(t - b)/g)))
-    AIF_atf = a*(d*e*exp(-(t - b)/e)*(exp((c - b)/e) - 1
-        ) + f*g*exp(-(t - b)/g)*(exp((c - b)/g) - 1))
-    AIF_sym = Piecewise((0, t < b), (AIF_btf, And(b <= t, t < c)), (AIF_atf, c <= t))
-    aif_eq = Eq(Z0(t), AIF_sym)
-    f_vec = Matrix([[k01*Z0(t)]])
-    integrand = m*m.subs(t, s)**-1*f_vec.subs(aif_eq.lhs, aif_eq.rhs).subs(t, s)
-    solution = integrate(integrand[0], (s, 0, t))
-    assert solution is not None  # does not hang and takes less than 10 s
diff --git a/sympy/integrals/tests/test_integrals.py b/sympy/integrals/tests/test_integrals.py
index 8605808968..6fa978f3bf 100644
--- a/sympy/integrals/tests/test_integrals.py
+++ b/sympy/integrals/tests/test_integrals.py
@@ -1148,7 +1148,7 @@ def test_issue_3940():
     assert integrate(exp(-x**2 + I*c*x), x) == \
         -sqrt(pi)*exp(-c**2/4)*erf(I*c/2 - x)/2
     assert integrate(exp(a*x**2 + b*x + c), x) == \
-        sqrt(pi)*exp(c - b**2/(4*a))*erfi((2*a*x + b)/(2*sqrt(a)))/(2*sqrt(a))
+        sqrt(pi)*exp(c)*exp(-b**2/(4*a))*erfi(sqrt(a)*x + b/(2*sqrt(a)))/(2*sqrt(a))
 
     from sympy.core.function import expand_mul
     from sympy.abc import k
diff --git a/sympy/printing/tests/test_c.py b/sympy/printing/tests/test_c.py
index 7ba03cb7a2..11836539f0 100644
--- a/sympy/printing/tests/test_c.py
+++ b/sympy/printing/tests/test_c.py
@@ -30,19 +30,6 @@
 x, y, z = symbols('x,y,z')
 
 
-def test_invalid_variable_names():
-    long_var = 'w'*32
-    syms = symbols(r'a_{\delta}[0], I_{x}, a*, 9t, bus#, %s' % long_var)
-    for s in syms:
-        with raises(ValueError):
-            ccode(s)
-    # these should not raise
-    short_var = 'w'*31
-    syms = symbols('funny_var[5], b[1][2], a[idx], out[i*N+j], %s' % short_var)
-    for s in syms:
-        ccode(s)
-
-
 def test_printmethod():
     class fabs(Abs):
         def _ccode(self, printer):