Implement Symbolic Recursive Integration for Dommaschk (#1496)

Resolves #1467 
Resolves #1278 

Changes implementation of Dommaschk potentials to use `sympy` and
symbolic integration to recursively compute the potentials with symbolic
algebra. I then use the `sympy.lambdify` function to turn this into a
jax-differentiable expression usable with our `ScalarPotentialField`
class.

Author: dpanici

CHANGELOG.md

@@ -14,6 +14,7 @@ New Features
 - ``desc.objectives.Omnigenity`` is now vectorized and able to optimize multiple surfaces at the same time. Previously it was required to use a different objective for each surface.
 - Adds a new objective ``desc.objectives.MirrorRatio`` for targeting a particular mirror ratio on each flux surface, for either an ``Equilibrium`` or ``OmnigenousField``.
 - Adds the output quantities ``wb`` and ``wp`` to ``VMECIO.save``.
+- Change implementation of Dommaschk potentials to use recursive algorithm and symbolic integration.
 - Changes hessian computation to use chunked ``jacfwd`` and ``jacrev``, allowing ``jac_chunk_size`` to now reduce hessian memory usage as well.
 
 Bug Fixes