diff --git a/CMakeLists.txt b/CMakeLists.txt
index 14f8ffc49..e6e63ec9b 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -123,7 +123,7 @@ pybind11_add_module(${PROJECT_NAME}
     src/simsoptpp/biot_savart_py.cpp
     src/simsoptpp/biot_savart_vjp_py.cpp
     src/simsoptpp/regular_grid_interpolant_3d_py.cpp
-    src/simsoptpp/curve.cpp src/simsoptpp/curverzfourier.cpp src/simsoptpp/curvexyzfourier.cpp
+    src/simsoptpp/curve.cpp src/simsoptpp/curverzfourier.cpp src/simsoptpp/curvexyzfourier.cpp src/simsoptpp/curveplanarfourier.cpp
     src/simsoptpp/surface.cpp src/simsoptpp/surfacerzfourier.cpp src/simsoptpp/surfacexyzfourier.cpp
     src/simsoptpp/integral_BdotN.cpp
     src/simsoptpp/dipole_field.cpp src/simsoptpp/permanent_magnet_optimization.cpp
diff --git a/docs/source/simsopt.geo.rst b/docs/source/simsopt.geo.rst
index 9bb006c3c..0a1bf5c32 100644
--- a/docs/source/simsopt.geo.rst
+++ b/docs/source/simsopt.geo.rst
@@ -60,6 +60,14 @@ simsopt.geo.curverzfourier module
    :undoc-members:
    :show-inheritance:
 
+simsopt.geo.curveplanarfourier module
+-------------------------------------
+
+.. automodule:: simsopt.geo.curveplanarfourier
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 simsopt.geo.curvexyzfourier module
 ----------------------------------
 
diff --git a/examples/2_Intermediate/stage_two_optimization_planar_coils.py b/examples/2_Intermediate/stage_two_optimization_planar_coils.py
new file mode 100755
index 000000000..8c9d62b94
--- /dev/null
+++ b/examples/2_Intermediate/stage_two_optimization_planar_coils.py
@@ -0,0 +1,199 @@
+#!/usr/bin/env python
+r"""
+This coil optimization script is similar to stage_two_optimization.py. However
+in this version, the coils are constrained to be planar, by using the curve type
+CurvePlanarFourier. Also the LinkingNumber objective is used to prevent coils
+from becoming topologically linked with each other.
+
+In this example we solve a FOCUS like Stage II coil optimisation problem: the
+goal is to find coils that generate a specific target normal field on a given
+surface.  In this particular case we consider a vacuum field, so the target is
+just zero.
+
+The objective is given by
+
+    J = (1/2) \int |B dot n|^2 ds
+        + LENGTH_WEIGHT * (sum CurveLength)
+        + DISTANCE_WEIGHT * MininumDistancePenalty(DISTANCE_THRESHOLD)
+        + CURVATURE_WEIGHT * CurvaturePenalty(CURVATURE_THRESHOLD)
+        + MSC_WEIGHT * MeanSquaredCurvaturePenalty(MSC_THRESHOLD)
+        + LinkingNumber
+
+if any of the weights are increased, or the thresholds are tightened, the coils
+are more regular and better separated, but the target normal field may not be
+achieved as well. This example demonstrates the adjustment of weights and
+penalties via the use of the `Weight` class.
+
+The target equilibrium is the QA configuration of arXiv:2108.03711.
+"""
+
+import os
+from pathlib import Path
+import numpy as np
+from scipy.optimize import minimize
+from simsopt.field import BiotSavart, Current, coils_via_symmetries
+from simsopt.geo import (
+    CurveLength, CurveCurveDistance,
+    MeanSquaredCurvature, LpCurveCurvature, CurveSurfaceDistance, LinkingNumber,
+    SurfaceRZFourier, curves_to_vtk, create_equally_spaced_planar_curves,
+)
+from simsopt.objectives import Weight, SquaredFlux, QuadraticPenalty
+from simsopt.util import in_github_actions
+
+# Number of unique coil shapes, i.e. the number of coils per half field period:
+# (Since the configuration has nfp = 2, multiply by 4 to get the total number of coils.)
+ncoils = 4
+
+# Major radius for the initial circular coils:
+R0 = 1.0
+
+# Minor radius for the initial circular coils:
+R1 = 0.5
+
+# Number of Fourier modes describing each Cartesian component of each coil:
+order = 5
+
+# Weight on the curve lengths in the objective function. We use the `Weight`
+# class here to later easily adjust the scalar value and rerun the optimization
+# without having to rebuild the objective.
+LENGTH_WEIGHT = Weight(10)
+
+# Threshold and weight for the coil-to-coil distance penalty in the objective function:
+CC_THRESHOLD = 0.08
+CC_WEIGHT = 1000
+
+# Threshold and weight for the coil-to-surface distance penalty in the objective function:
+CS_THRESHOLD = 0.12
+CS_WEIGHT = 10
+
+# Threshold and weight for the curvature penalty in the objective function:
+CURVATURE_THRESHOLD = 10.
+CURVATURE_WEIGHT = 1e-6
+
+# Threshold and weight for the mean squared curvature penalty in the objective function:
+MSC_THRESHOLD = 10
+MSC_WEIGHT = 1e-6
+
+# Number of iterations to perform:
+MAXITER = 50 if in_github_actions else 400
+
+# File for the desired boundary magnetic surface:
+TEST_DIR = (Path(__file__).parent / ".." / ".." / "tests" / "test_files").resolve()
+filename = TEST_DIR / 'input.LandremanPaul2021_QA'
+
+# Directory for output
+OUT_DIR = "./output/"
+os.makedirs(OUT_DIR, exist_ok=True)
+
+#######################################################
+# End of input parameters.
+#######################################################
+
+# Initialize the boundary magnetic surface:
+nphi = 32
+ntheta = 32
+s = SurfaceRZFourier.from_vmec_input(filename, range="half period", nphi=nphi, ntheta=ntheta)
+
+# Create the initial coils:
+base_curves = create_equally_spaced_planar_curves(ncoils, s.nfp, stellsym=True, R0=R0, R1=R1, order=order)
+base_currents = [Current(1e5) for i in range(ncoils)]
+# Since the target field is zero, one possible solution is just to set all
+# currents to 0. To avoid the minimizer finding that solution, we fix one
+# of the currents:
+base_currents[0].fix_all()
+
+coils = coils_via_symmetries(base_curves, base_currents, s.nfp, True)
+bs = BiotSavart(coils)
+bs.set_points(s.gamma().reshape((-1, 3)))
+
+curves = [c.curve for c in coils]
+curves_to_vtk(curves, OUT_DIR + "curves_init")
+pointData = {"B_N": np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2)[:, :, None]}
+s.to_vtk(OUT_DIR + "surf_init", extra_data=pointData)
+
+# Define the individual terms objective function:
+Jf = SquaredFlux(s, bs)
+Jls = [CurveLength(c) for c in base_curves]
+Jccdist = CurveCurveDistance(curves, CC_THRESHOLD, num_basecurves=ncoils)
+Jcsdist = CurveSurfaceDistance(curves, s, CS_THRESHOLD)
+Jcs = [LpCurveCurvature(c, 2, CURVATURE_THRESHOLD) for c in base_curves]
+Jmscs = [MeanSquaredCurvature(c) for c in base_curves]
+linkNum = LinkingNumber(curves)
+
+# Form the total objective function. To do this, we can exploit the
+# fact that Optimizable objects with J() and dJ() functions can be
+# multiplied by scalars and added:
+JF = Jf \
+    + LENGTH_WEIGHT * QuadraticPenalty(sum(Jls), 2.6*ncoils) \
+    + CC_WEIGHT * Jccdist \
+    + CS_WEIGHT * Jcsdist \
+    + CURVATURE_WEIGHT * sum(Jcs) \
+    + MSC_WEIGHT * sum(QuadraticPenalty(J, MSC_THRESHOLD) for J in Jmscs) \
+    + linkNum
+
+# We don't have a general interface in SIMSOPT for optimisation problems that
+# are not in least-squares form, so we write a little wrapper function that we
+# pass directly to scipy.optimize.minimize
+
+
+def fun(dofs):
+    JF.x = dofs
+    J = JF.J()
+    grad = JF.dJ()
+    jf = Jf.J()
+    BdotN = np.mean(np.abs(np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2)))
+    MaxBdotN = np.max(np.abs(np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2)))
+    mean_AbsB = np.mean(bs.AbsB())
+    outstr = f"J={J:.1e}, Jf={jf:.1e}, ⟨B·n⟩={BdotN:.1e}"
+    cl_string = ", ".join([f"{J.J():.1f}" for J in Jls])
+    kap_string = ", ".join(f"{np.max(c.kappa()):.1f}" for c in base_curves)
+    msc_string = ", ".join(f"{J.J():.1f}" for J in Jmscs)
+    outstr += f", Len=sum([{cl_string}])={sum(J.J() for J in Jls):.1f}, ϰ=[{kap_string}], ∫ϰ²/L=[{msc_string}]"
+    outstr += f", C-C-Sep={Jccdist.shortest_distance():.2f}, C-S-Sep={Jcsdist.shortest_distance():.2f}"
+    outstr += f", ║∇J║={np.linalg.norm(grad):.1e}"
+    outstr += f", ⟨B·n⟩/|B|={BdotN/mean_AbsB:.1e}"
+    outstr += f", (Max B·n)/|B|={MaxBdotN/mean_AbsB:.1e}"
+    outstr += f", Link Number = {linkNum.J()}"
+    print(outstr)
+    return J, grad
+
+
+print("""
+################################################################################
+### Perform a Taylor test ######################################################
+################################################################################
+""")
+f = fun
+dofs = JF.x
+np.random.seed(1)
+h = np.random.uniform(size=dofs.shape)
+J0, dJ0 = f(dofs)
+dJh = sum(dJ0 * h)
+for eps in [1e-3, 1e-4, 1e-5, 1e-6, 1e-7]:
+    J1, _ = f(dofs + eps*h)
+    J2, _ = f(dofs - eps*h)
+    print("err", (J1-J2)/(2*eps) - dJh)
+
+print("""
+################################################################################
+### Run the optimisation #######################################################
+################################################################################
+""")
+res = minimize(fun, dofs, jac=True, method='L-BFGS-B', options={'maxiter': MAXITER, 'maxcor': 300}, tol=1e-15)
+curves_to_vtk(curves, OUT_DIR + "curves_opt_short")
+pointData = {"B_N": np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2)[:, :, None]}
+s.to_vtk(OUT_DIR + "surf_opt_short", extra_data=pointData)
+
+
+# We now use the result from the optimization as the initial guess for a
+# subsequent optimization with reduced penalty for the coil length. This will
+# result in slightly longer coils but smaller `B·n` on the surface.
+dofs = res.x
+LENGTH_WEIGHT *= 0.1
+res = minimize(fun, dofs, jac=True, method='L-BFGS-B', options={'maxiter': MAXITER, 'maxcor': 300}, tol=1e-15)
+curves_to_vtk(curves, OUT_DIR + "curves_opt_long")
+pointData = {"B_N": np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2)[:, :, None]}
+s.to_vtk(OUT_DIR + "surf_opt_long", extra_data=pointData)
+
+# Save the optimized coil shapes and currents so they can be loaded into other scripts for analysis:
+bs.save(OUT_DIR + "biot_savart_opt.json")
diff --git a/examples/run_serial_examples b/examples/run_serial_examples
index 163ff6ad4..b87f64702 100755
--- a/examples/run_serial_examples
+++ b/examples/run_serial_examples
@@ -16,6 +16,7 @@ set -ex
 ./2_Intermediate/QSC.py
 ./2_Intermediate/boozer.py
 ./2_Intermediate/stage_two_optimization.py
+./2_Intermediate/stage_two_optimization_planar_coils.py
 ./2_Intermediate/stage_two_optimization_stochastic.py
 ./2_Intermediate/stage_two_optimization_finite_beta.py
 ./2_Intermediate/strain_optimization.py
diff --git a/src/simsopt/geo/__init__.py b/src/simsopt/geo/__init__.py
index bc6d0efe2..beb5e57b6 100644
--- a/src/simsopt/geo/__init__.py
+++ b/src/simsopt/geo/__init__.py
@@ -8,6 +8,7 @@
 from .curvexyzfourier import *
 from .curveperturbed import *
 from .curveobjectives import *
+from .curveplanarfourier import *
 from .framedcurve import *
 from .finitebuild import *
 from .plotting import *
@@ -28,6 +29,7 @@
 __all__ = (curve.__all__ + curvehelical.__all__ +
            curverzfourier.__all__ + curvexyzfourier.__all__ +
            curveperturbed.__all__ + curveobjectives.__all__ +
+           curveplanarfourier.__all__ +
            finitebuild.__all__ + plotting.__all__ +
            boozersurface.__all__ + qfmsurface.__all__ +
            surface.__all__ +
diff --git a/src/simsopt/geo/curve.py b/src/simsopt/geo/curve.py
index a19fc92a2..bda3a5e61 100644
--- a/src/simsopt/geo/curve.py
+++ b/src/simsopt/geo/curve.py
@@ -11,7 +11,7 @@
 from .jit import jit
 from .plotting import fix_matplotlib_3d
 
-__all__ = ['Curve', 'RotatedCurve', 'curves_to_vtk', 'create_equally_spaced_curves']
+__all__ = ['Curve', 'RotatedCurve', 'curves_to_vtk', 'create_equally_spaced_curves', 'create_equally_spaced_planar_curves']
 
 
 @jit
@@ -878,3 +878,46 @@ def create_equally_spaced_curves(ncurves, nfp, stellsym, R0=1.0, R1=0.5, order=6
         curve.x = curve.x  # need to do this to transfer data to C++
         curves.append(curve)
     return curves
+
+
+def create_equally_spaced_planar_curves(ncurves, nfp, stellsym, R0=1.0, R1=0.5, order=6, numquadpoints=None):
+    """
+    Create ``ncurves`` curves of type
+    :obj:`~simsopt.geo.curveplanarfourier.CurvePlanarFourier` of order
+    ``order`` that will result in circular equally spaced coils (major
+    radius ``R0`` and minor radius ``R1``) after applying
+    :obj:`~simsopt.field.coil.coils_via_symmetries`.
+    """
+
+    if numquadpoints is None:
+        numquadpoints = 15 * order
+    curves = []
+    from simsopt.geo.curveplanarfourier import CurvePlanarFourier
+    for k in range(ncurves):
+        angle = (k+0.5)*(2*np.pi) / ((1+int(stellsym))*nfp*ncurves)
+        curve = CurvePlanarFourier(numquadpoints, order, nfp, stellsym)
+
+        rcCoeffs = np.zeros(order+1)
+        rcCoeffs[0] = R1
+        rsCoeffs = np.zeros(order)
+        center = [R0 * cos(angle), R0 * sin(angle), 0]
+        rotation = [1, -cos(angle), -sin(angle), 0]
+        dofs = np.zeros(len(curve.get_dofs()))
+
+        j = 0
+        for i in rcCoeffs:
+            dofs[j] = i
+            j += 1
+        for i in rsCoeffs:
+            dofs[j] = i
+            j += 1
+        for i in rotation:
+            dofs[j] = i
+            j += 1
+        for i in center:
+            dofs[j] = i
+            j += 1
+
+        curve.set_dofs(dofs)
+        curves.append(curve)
+    return curves
\ No newline at end of file
diff --git a/src/simsopt/geo/curveobjectives.py b/src/simsopt/geo/curveobjectives.py
index baf41d598..bb81ec318 100644
--- a/src/simsopt/geo/curveobjectives.py
+++ b/src/simsopt/geo/curveobjectives.py
@@ -520,7 +520,7 @@ def J(self):
                     integrals = sopp.linkNumber(R1, R2, dR1, dR2) * dS * dT
                     linkNum[i-1][j-1] = 1/(4*np.pi) * (integrals)
         linkNumSum = sum(sum(abs(linkNum)))
-        return linkNumSum
+        return round(linkNumSum)
 
     @derivative_dec
     def dJ(self):
diff --git a/src/simsopt/geo/curveplanarfourier.py b/src/simsopt/geo/curveplanarfourier.py
new file mode 100644
index 000000000..4880d05d6
--- /dev/null
+++ b/src/simsopt/geo/curveplanarfourier.py
@@ -0,0 +1,75 @@
+import numpy as np
+
+import simsoptpp as sopp
+from .curve import Curve
+
+__all__ = ['CurvePlanarFourier']
+
+
+class CurvePlanarFourier(sopp.CurvePlanarFourier, Curve):
+    r"""
+    ``CurvePlanarFourier`` is a curve that is restricted to lie in a plane. The
+    shape of the curve within the plane is represented by a Fourier series in
+    polar coordinates. The resulting planar curve is then rotated in three
+    dimensions using a quaternion, and finally a translation is applied. The
+    Fourier series in polar coordinates is
+
+    .. math::
+
+       r(\phi) = \sum_{m=0}^{\text{order}} r_{c,m}\cos(m \phi) + \sum_{m=1}^{\text{order}} r_{s,m}\sin(m \phi).
+
+    The rotation quaternion is
+
+    .. math::
+
+       \bf{q} &= [q_0,q_i,q_j,q_k]
+
+       &= [\cos(\theta / 2), \hat{x}\sin(\theta / 2), \hat{y}\sin(\theta / 2), \hat{z}\sin(\theta / 2)]
+
+    where :math:`\theta` is the counterclockwise rotation angle about a unit axis
+    :math:`(\hat{x},\hat{y},\hat{z})`. Details of the quaternion rotation can be
+    found for example in pages 575-576 of
+    https://www.cis.upenn.edu/~cis5150/ws-book-Ib.pdf.
+
+
+    A quaternion is used for rotation rather than other methods for rotation to
+    prevent gimbal locking during optimization. The quaternion is normalized
+    before being applied to prevent scaling of the curve. The dofs themselves are not normalized. This
+    results in a redundancy in the optimization, where several different sets of 
+    dofs may correspond to the same normalized quaternion. Normalizing the dofs 
+    directly would create a dependence between the quaternion dofs, which may cause 
+    issues during optimization.
+
+    The dofs are stored in the order
+
+    .. math::
+       [r_{c,0}, \cdots, r_{c,\text{order}}, r_{s,1}, \cdots, r_{s,\text{order}}, q_0, q_i, q_j, q_k, x_{\text{center}}, y_{\text{center}}, z_{\text{center}}]
+
+
+    """
+
+    def __init__(self, quadpoints, order, nfp, stellsym, dofs=None):
+        if isinstance(quadpoints, int):
+            quadpoints = list(np.linspace(0, 1., quadpoints, endpoint=False))
+        elif isinstance(quadpoints, np.ndarray):
+            quadpoints = list(quadpoints)
+        sopp.CurvePlanarFourier.__init__(self, quadpoints, order, nfp, stellsym)
+        if dofs is None:
+            Curve.__init__(self, external_dof_setter=CurvePlanarFourier.set_dofs_impl,
+                           x0=self.get_dofs())
+        else:
+            Curve.__init__(self, external_dof_setter=CurvePlanarFourier.set_dofs_impl,
+                           dofs=dofs)
+
+    def get_dofs(self):
+        """
+        This function returns the dofs associated to this object.
+        """
+        return np.asarray(sopp.CurvePlanarFourier.get_dofs(self))
+
+    def set_dofs(self, dofs):
+        """
+        This function sets the dofs associated to this object.
+        """
+        self.local_x = dofs
+        sopp.CurvePlanarFourier.set_dofs(self, dofs)
diff --git a/src/simsoptpp/curveplanarfourier.cpp b/src/simsoptpp/curveplanarfourier.cpp
new file mode 100644
index 000000000..9e11bc3e8
--- /dev/null
+++ b/src/simsoptpp/curveplanarfourier.cpp
@@ -0,0 +1,1598 @@
+#include "curveplanarfourier.h"
+
+template<class Array>
+double CurvePlanarFourier<Array>::inv_magnitude() {
+    double s = q[0] * q[0] + q[1] * q[1] + q[2] * q[2] + q[3] * q[3];
+    if(s != 0) {
+        return 1 / std::sqrt(s);
+    } else {
+        return 1;
+    }
+
+}
+
+template<class Array>
+void CurvePlanarFourier<Array>::gamma_impl(Array& data, Array& quadpoints) {
+    double sinphi, cosphi, siniphi, cosiphi;
+    int numquadpoints = quadpoints.size();
+    data *= 0;
+
+    Array q_norm = q * inv_magnitude();
+
+
+    for (int k = 0; k < numquadpoints; ++k) {
+        double phi = 2 * M_PI * quadpoints[k];
+        double cosphi = cos(phi);
+        double sinphi = sin(phi);
+        data(k, 0) = rc[0] * cosphi;
+        data(k, 1) = rc[0] * sinphi;
+        for (int i = 1; i < order+1; ++i) {
+            cosiphi = cos(i*phi);
+            siniphi = sin(i*phi);
+            data(k, 0) += (rc[i] * cosiphi + rs[i-1] * siniphi) * cosphi;
+            data(k, 1) += (rc[i] * cosiphi + rs[i-1] * siniphi) * sinphi;
+        }
+    }
+    for (int m = 0; m < numquadpoints; ++m) {
+        double i = data(m, 0);
+        double j = data(m, 1);
+        double k = data(m, 2);
+
+        /* Performs quaternion based rotation, see https://www.cis.upenn.edu/~cis5150/ws-book-Ib.pdf page 575, 576 for details regarding this rotation*/
+        data(m, 0) = (i - 2 * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * i + 2 * (q_norm[1] * q_norm[2] - q_norm[3] * q_norm[0]) * j + 2 * (q_norm[1] * q_norm[3] + q_norm[2] * q_norm[0]) * k) + center[0];
+        data(m, 1) = (2 * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * i + j - 2 * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * j + 2 * (q_norm[2] * q_norm[3] - q_norm[1] * q_norm[0]) * k) + center[1];
+        data(m, 2) = (2 * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * i + 2 * (q_norm[2] * q_norm[3] + q_norm[1] * q_norm[0]) * j + k - 2 * (q_norm[1] * q_norm[1] + q_norm[2] * q_norm[2]) * k) + center[2];
+    }
+}
+
+template<class Array>
+void CurvePlanarFourier<Array>::gammadash_impl(Array& data) {
+    data *= 0;
+
+    double inv_sqrt_s = inv_magnitude();
+    Array q_norm = q * inv_sqrt_s;
+
+    double cosiphi, siniphi;
+    for (int k = 0; k < numquadpoints; ++k) {
+        double phi = 2 * M_PI * quadpoints[k];
+        double cosphi = cos(phi);
+        double sinphi = sin(phi);
+        data(k, 0) = rc[0] * (-sinphi);
+        data(k, 1) = rc[0] * (cosphi);
+        for (int i = 1; i < order+1; ++i) {
+            cosiphi = cos(i*phi);
+            siniphi = sin(i*phi);
+            data(k, 0) += rc[i] * ( -(i) * siniphi * cosphi - cosiphi * sinphi) 
+                + rs[i-1] * ( (i) * cosiphi * cosphi - siniphi * sinphi);
+            data(k, 1) += rc[i] * ( -(i) * siniphi * sinphi + cosiphi * cosphi)
+                + rs[i-1] * ( (i) * cosiphi * sinphi + siniphi * cosphi);
+        }
+    }
+        
+    data *= (2*M_PI);
+    for (int m = 0; m < numquadpoints; ++m) {
+        double i = data(m, 0);
+        double j = data(m, 1);
+        double k = data(m, 2);
+
+        /* Performs quaternion based rotation, see https://www.cis.upenn.edu/~cis5150/ws-book-Ib.pdf page 575, 576 for details regarding this rotation*/
+        data(m, 0) = (i - 2 * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * i + 2 * (q_norm[1] * q_norm[2] - q_norm[3] * q_norm[0]) * j + 2 * (q_norm[1] * q_norm[3] + q_norm[2] * q_norm[0]) * k);
+        data(m, 1) = (2 * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * i + j - 2 * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * j + 2 * (q_norm[2] * q_norm[3] - q_norm[1] * q_norm[0]) * k);
+        data(m, 2) = (2 * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * i + 2 * (q_norm[2] * q_norm[3] + q_norm[1] * q_norm[0]) * j + k - 2 * (q_norm[1] * q_norm[1] + q_norm[2] * q_norm[2]) * k);
+    }
+    
+}
+
+template<class Array>
+void CurvePlanarFourier<Array>::gammadashdash_impl(Array& data) {
+    data *= 0;
+
+    Array q_norm = q * inv_magnitude();
+
+    double cosiphi, siniphi;
+    for (int k = 0; k < numquadpoints; ++k) {
+        double phi = 2 * M_PI * quadpoints[k];
+        double cosphi = cos(phi);
+        double sinphi = sin(phi);
+        data(k, 0) = rc[0] * (-cosphi);
+        data(k, 1) = rc[0] * (-sinphi);
+        for (int i = 1; i < order+1; ++i) {
+            cosiphi = cos(i*phi);
+            siniphi = sin(i*phi);
+            data(k, 0) += rc[i] * (+2*(i)*siniphi*sinphi-(i*i+1)*cosiphi*cosphi)
+                + rs[i-1] * (-(i*i+1)*siniphi*cosphi - 2*(i)*cosiphi*sinphi);
+            data(k, 1) += rc[i] * (-2*(i)*siniphi*cosphi-(i*i+1)*cosiphi*sinphi)
+                + rs[i-1] * (-(i*i+1)*siniphi*sinphi + 2*(i)*cosiphi*cosphi);
+        }
+    }
+    data *= 2*M_PI*2*M_PI;
+    for (int m = 0; m < numquadpoints; ++m) {
+        double i = data(m, 0);
+        double j = data(m, 1);
+        double k = data(m, 2);
+
+        /* Performs quaternion based rotation*/
+        data(m, 0) = (i - 2 * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * i + 2 * (q_norm[1] * q_norm[2] - q_norm[3] * q_norm[0]) * j + 2 * (q_norm[1] * q_norm[3] + q_norm[2] * q_norm[0]) * k);
+        data(m, 1) = (2 * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * i + j - 2 * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * j + 2 * (q_norm[2] * q_norm[3] - q_norm[1] * q_norm[0]) * k);
+        data(m, 2) = (2 * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * i + 2 * (q_norm[2] * q_norm[3] + q_norm[1] * q_norm[0]) * j + k - 2 * (q_norm[1] * q_norm[1] + q_norm[2] * q_norm[2]) * k);
+
+    }
+}
+
+template<class Array>
+void CurvePlanarFourier<Array>::gammadashdashdash_impl(Array& data) {
+    data *= 0;
+
+    Array q_norm = q * inv_magnitude();
+
+    double cosiphi, siniphi;
+    for (int k = 0; k < numquadpoints; ++k) {
+        double phi = 2 * M_PI * quadpoints[k];
+        double cosphi = cos(phi);
+        double sinphi = sin(phi);
+
+        data(k, 0) = rc[0]*(+sinphi);
+        data(k, 1) = rc[0]*(-cosphi);
+        for (int i = 1; i < order+1; ++i) {
+            cosiphi = cos(i*phi);
+            siniphi = sin(i*phi);
+            data(k, 0) += rc[i]*(
+                    +(3*i*i + 1)*cosiphi*sinphi
+                    +(i*i + 3)*(i)*siniphi*cosphi
+                    ) + rs[i-1]*(
+                    -(i*i+3) * (i) * cosiphi*cosphi
+                    +(3*i*i+1) * siniphi*sinphi
+                    );
+            data(k, 1) += rc[i]*(
+                    +(i*i + 3)*(i)*siniphi*sinphi
+                    -(3*i*i + 1)*cosiphi*cosphi
+                    ) + rs[i-1]*(
+                    -(i*i+3)*(i)*cosiphi*sinphi 
+                    -(3*i*i+1)*siniphi*cosphi 
+                    );
+        }
+    }
+    data *= 2*M_PI*2*M_PI*2*M_PI;
+    for (int m = 0; m < numquadpoints; ++m) {
+        double i = data(m, 0);
+        double j = data(m, 1);
+        double k = data(m, 2);
+
+        /* Performs quaternion based rotation*/
+        data(m, 0) = (i - 2 * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * i + 2 * (q_norm[1] * q_norm[2] - q_norm[3] * q_norm[0]) * j + 2 * (q_norm[1] * q_norm[3] + q_norm[2] * q_norm[0]) * k);
+        data(m, 1) = (2 * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * i + j - 2 * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * j + 2 * (q_norm[2] * q_norm[3] - q_norm[1] * q_norm[0]) * k);
+        data(m, 2) = (2 * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * i + 2 * (q_norm[2] * q_norm[3] + q_norm[1] * q_norm[0]) * j + k - 2 * (q_norm[1] * q_norm[1] + q_norm[2] * q_norm[2]) * k);
+
+    }
+}
+
+template<class Array>
+void CurvePlanarFourier<Array>::dgamma_by_dcoeff_impl(Array& data) {
+    data *= 0;
+    
+    Array q_norm = q * inv_magnitude();
+
+    double cosnphi, sinnphi;
+    for (int m = 0; m < numquadpoints; ++m) {
+        double phi = 2 * M_PI * quadpoints[m];
+        int counter = 0;
+        double i;
+        double j;
+        double k;
+
+        double cosphi = cos(phi);
+        double sinphi = sin(phi);
+        double cosnphi = 0;
+        double sinnphi = 0;
+
+        for (int n = 0; n < order+1; ++n) {
+            cosnphi = cos(n * phi);
+            i = cosnphi * cosphi;
+            j = cosnphi * sinphi;
+            k = 0;
+
+            data(m, 0, counter) = (i - 2 * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * i + 2 * (q_norm[1] * q_norm[2] - q_norm[3] * q_norm[0]) * j + 2 * (q_norm[1] * q_norm[3] + q_norm[2] * q_norm[0]) * k);
+            data(m, 1, counter) = (2 * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * i + j - 2 * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * j + 2 * (q_norm[2] * q_norm[3] - q_norm[1] * q_norm[0]) * k);
+            data(m, 2, counter) = (2 * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * i + 2 * (q_norm[2] * q_norm[3] + q_norm[1] * q_norm[0]) * j + k - 2 * (q_norm[1] * q_norm[1] + q_norm[2] * q_norm[2]) * k);
+
+            counter++;
+        }
+        
+        for (int n = 1; n < order+1; ++n) {
+            sinnphi = sin(n * phi);
+            i = sinnphi * cosphi;
+            j = sinnphi * sinphi;
+            k = 0;
+
+            data(m, 0, counter) = (i - 2 * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * i + 2 * (q_norm[1] * q_norm[2] - q_norm[3] * q_norm[0]) * j + 2 * (q_norm[1] * q_norm[3] + q_norm[2] * q_norm[0]) * k);
+            data(m, 1, counter) = (2 * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * i + j - 2 * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * j + 2 * (q_norm[2] * q_norm[3] - q_norm[1] * q_norm[0]) * k);
+            data(m, 2, counter) = (2 * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * i + 2 * (q_norm[2] * q_norm[3] + q_norm[1] * q_norm[0]) * j + k - 2 * (q_norm[1] * q_norm[1] + q_norm[2] * q_norm[2]) * k);
+
+            counter++;
+        }
+
+        i = rc[0] * cosphi;
+        j = rc[0] * sinphi;
+        k = 0;
+
+        for (int n = 1; n < order+1; ++n) {
+            cosnphi = cos(n * phi);
+            sinnphi = sin(n * phi);
+            i += (rc[n] * cosnphi + rs[n-1] * sinnphi) * cosphi;
+            j += (rc[n] * cosnphi + rs[n-1] * sinnphi) * sinphi;
+        }
+        double inv_sqrt_s = inv_magnitude();
+        
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (inv_sqrt_s - q[0] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (- q[1] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (- q[2] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0]) 
+                            * (- q[3] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (inv_sqrt_s - q[0] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (- q[1] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (- q[2] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (- q[3] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (inv_sqrt_s - q[0] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (- q[1] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (- q[2] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (- q[3] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0] 
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (- q[0] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (inv_sqrt_s - q[1] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (- q[2] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0])
+                            * (- q[3] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (- q[0] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (inv_sqrt_s - q[1] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (- q[2] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (- q[3] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (- q[0] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (inv_sqrt_s - q[1] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (- q[2] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (- q[3] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0] 
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (- q[0] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            + 
+                            (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (- q[1] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (inv_sqrt_s - q[2] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0]) 
+                            * (- q[3] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (- q[0] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (- q[1] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (inv_sqrt_s - q[2] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (- q[3] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (- q[0] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (- q[1] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (inv_sqrt_s - q[2] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (- q[3] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0] 
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (- q[0] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            + (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (- q[1] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (- q[2] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +                            
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0]) 
+                            * (inv_sqrt_s - q[3] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (- q[0] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (- q[1] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (- q[2] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (inv_sqrt_s - q[3] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (- q[0] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (- q[1] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (- q[2] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (inv_sqrt_s - q[3] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+
+        for (int i = 0; i < 3; ++i) {
+            data(m, 0, counter) = 0;
+            data(m, 1, counter) = 0;
+            data(m, 2, counter) = 0;
+            data(m, i, counter) = 1;
+            counter++;
+        }
+    }
+}
+
+template<class Array>
+void CurvePlanarFourier<Array>::dgammadash_by_dcoeff_impl(Array& data) {
+    data *= 0;
+
+    Array q_norm = q * inv_magnitude();
+
+    double cosnphi, sinnphi;
+    for (int m = 0; m < numquadpoints; ++m) {
+        double phi = 2 * M_PI * quadpoints[m];
+        double cosphi = cos(phi);
+        double sinphi = sin(phi);
+        int counter = 0;
+        double i;
+        double j;
+        double k;
+        for (int n = 0; n < order+1; ++n) {
+            cosnphi = cos(n*phi);
+            sinnphi = sin(n*phi);
+            i = ( -(n) * sinnphi * cosphi - cosnphi * sinphi);
+            j = ( -(n) * sinnphi * sinphi + cosnphi * cosphi);
+            k = 0;
+
+            i *= (2*M_PI);
+            j *= (2*M_PI);
+
+            data(m, 0, counter) = (i - 2 * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * i + 2 * (q_norm[1] * q_norm[2] - q_norm[3] * q_norm[0]) * j + 2 * (q_norm[1] * q_norm[3] + q_norm[2] * q_norm[0]) * k);
+            data(m, 1, counter) = (2 * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * i + j - 2 * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * j + 2 * (q_norm[2] * q_norm[3] - q_norm[1] * q_norm[0]) * k);
+            data(m, 2, counter) = (2 * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * i + 2 * (q_norm[2] * q_norm[3] + q_norm[1] * q_norm[0]) * j + k - 2 * (q_norm[1] * q_norm[1] + q_norm[2] * q_norm[2]) * k);
+            
+
+            counter++;
+        }
+        
+        for (int n = 1; n < order+1; ++n) {
+            cosnphi = cos(n*phi);
+            sinnphi = sin(n*phi);
+            i = ( (n) * cosnphi * cosphi - sinnphi * sinphi);
+            j = ( (n) * cosnphi * sinphi + sinnphi * cosphi);
+            k = 0;
+
+            i *= (2*M_PI);
+            j *= (2*M_PI);
+
+            data(m, 0, counter) = (i - 2 * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * i + 2 * (q_norm[1] * q_norm[2] - q_norm[3] * q_norm[0]) * j + 2 * (q_norm[1] * q_norm[3] + q_norm[2] * q_norm[0]) * k);
+            data(m, 1, counter) = (2 * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * i + j - 2 * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * j + 2 * (q_norm[2] * q_norm[3] - q_norm[1] * q_norm[0]) * k);
+            data(m, 2, counter) = (2 * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * i + 2 * (q_norm[2] * q_norm[3] + q_norm[1] * q_norm[0]) * j + k - 2 * (q_norm[1] * q_norm[1] + q_norm[2] * q_norm[2]) * k);
+
+            counter++;
+            
+        }
+        
+        i = rc[0] * (-sinphi);
+        j = rc[0] * (cosphi);
+        k = 0;
+        for (int n = 1; n < order+1; ++n) {
+            cosnphi = cos(n*phi);
+            sinnphi = sin(n*phi);
+            i += rc[n] * ( -(n) * sinnphi * cosphi - cosnphi * sinphi)
+                + rs[n-1] * ( (n) * cosnphi * cosphi - sinnphi * sinphi);
+            j += rc[n] * ( -(n) * sinnphi * sinphi + cosnphi * cosphi)
+                + rs[n-1] * ( (n) * cosnphi * sinphi + sinnphi * cosphi);
+        }
+        i *= (2*M_PI);
+        j *= (2*M_PI);
+        double inv_sqrt_s = inv_magnitude();
+
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (inv_sqrt_s - q[0] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (- q[1] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (- q[2] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0]) 
+                            * (- q[3] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (inv_sqrt_s - q[0] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (- q[1] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (- q[2] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (- q[3] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (inv_sqrt_s - q[0] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (- q[1] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (- q[2] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (- q[3] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0] 
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (- q[0] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (inv_sqrt_s - q[1] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (- q[2] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0])
+                            * (- q[3] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (- q[0] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (inv_sqrt_s - q[1] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (- q[2] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (- q[3] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (- q[0] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (inv_sqrt_s - q[1] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (- q[2] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (- q[3] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0] 
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (- q[0] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            + 
+                            (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (- q[1] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (inv_sqrt_s - q[2] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0]) 
+                            * (- q[3] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (- q[0] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (- q[1] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (inv_sqrt_s - q[2] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (- q[3] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (- q[0] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (- q[1] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (inv_sqrt_s - q[2] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (- q[3] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0] 
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (- q[0] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            + (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (- q[1] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (- q[2] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +                            
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0]) 
+                            * (inv_sqrt_s - q[3] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (- q[0] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (- q[1] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (- q[2] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (inv_sqrt_s - q[3] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (- q[0] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (- q[1] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (- q[2] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (inv_sqrt_s - q[3] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+        for (int i = 0; i < 2; ++i) {
+            data(m, 0, counter) = 0;
+            data(m, 1, counter) = 0;
+            data(m, 2, counter) = 0;
+            
+            counter++;
+        }
+        
+    }
+}
+
+template<class Array>
+void CurvePlanarFourier<Array>::dgammadashdash_by_dcoeff_impl(Array& data) {
+    data *= 0;
+
+    Array q_norm = q * inv_magnitude();
+
+    double cosnphi, sinnphi;
+    for (int m = 0; m < numquadpoints; ++m) {
+        double phi = 2 * M_PI * quadpoints[m];
+        double cosphi = cos(phi);
+        double sinphi = sin(phi);
+        int counter = 0;
+        double i;
+        double j;
+        double k;
+        for (int n = 0; n < order+1; ++n) {
+            cosnphi = cos(n*phi);
+            sinnphi = sin(n*phi);
+            i = (+2*(n)*sinnphi*sinphi-(n*n+1)*cosnphi*cosphi);
+            j = (-2*(n)*sinnphi*cosphi-(n*n+1)*cosnphi*sinphi);
+            k = 0;
+
+            i *= 2*M_PI*2*M_PI;
+            j *= 2*M_PI*2*M_PI;
+            k *= 2*M_PI*2*M_PI;
+
+            data(m, 0, counter) = (i - 2 * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * i + 2 * (q_norm[1] * q_norm[2] - q_norm[3] * q_norm[0]) * j + 2 * (q_norm[1] * q_norm[3] + q_norm[2] * q_norm[0]) * k);
+            data(m, 1, counter) = (2 * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * i + j - 2 * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * j + 2 * (q_norm[2] * q_norm[3] - q_norm[1] * q_norm[0]) * k);
+            data(m, 2, counter) = (2 * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * i + 2 * (q_norm[2] * q_norm[3] + q_norm[1] * q_norm[0]) * j + k - 2 * (q_norm[1] * q_norm[1] + q_norm[2] * q_norm[2]) * k);
+            
+            counter++;
+        }
+        
+        for (int n = 1; n < order+1; ++n) {
+            cosnphi = cos(n*phi);
+            sinnphi = sin(n*phi);
+            i = (-(n*n+1)*sinnphi*cosphi - 2*(n)*cosnphi*sinphi);
+            j = (-(n*n+1)*sinnphi*sinphi + 2*(n)*cosnphi*cosphi);
+            k = 0;
+
+            i *= 2*M_PI*2*M_PI;
+            j *= 2*M_PI*2*M_PI;
+            k *= 2*M_PI*2*M_PI;
+
+            data(m, 0, counter) = (i - 2 * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * i + 2 * (q_norm[1] * q_norm[2] - q_norm[3] * q_norm[0]) * j + 2 * (q_norm[1] * q_norm[3] + q_norm[2] * q_norm[0]) * k);
+            data(m, 1, counter) = (2 * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * i + j - 2 * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * j + 2 * (q_norm[2] * q_norm[3] - q_norm[1] * q_norm[0]) * k);
+            data(m, 2, counter) = (2 * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * i + 2 * (q_norm[2] * q_norm[3] + q_norm[1] * q_norm[0]) * j + k - 2 * (q_norm[1] * q_norm[1] + q_norm[2] * q_norm[2]) * k);
+            
+            counter++;
+        }
+        
+        i = rc[0] * (-cosphi);
+        j = rc[0] * (-sinphi);
+        k = 0;
+        for (int n = 1; n < order+1; ++n) {
+            cosnphi = cos(n*phi);
+            sinnphi = sin(n*phi);
+            i += rc[n] * (+2*(n)*sinnphi*sinphi-(n*n+1)*cosnphi*cosphi)
+                + rs[n-1] * (-(n*n+1)*sinnphi*cosphi - 2*(n)*cosnphi*sinphi);
+            j += rc[n] * (-2*(n)*sinnphi*cosphi-(n*n+1)*cosnphi*sinphi)
+                + rs[n-1] * (-(n*n+1)*sinnphi*sinphi + 2*(n)*cosnphi*cosphi);
+        }
+        i *= 2*M_PI*2*M_PI;
+        j *= 2*M_PI*2*M_PI;
+        k *= 2*M_PI*2*M_PI;
+        double inv_sqrt_s = inv_magnitude();
+        
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (inv_sqrt_s - q[0] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (- q[1] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (- q[2] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0]) 
+                            * (- q[3] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (inv_sqrt_s - q[0] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (- q[1] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (- q[2] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (- q[3] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (inv_sqrt_s - q[0] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (- q[1] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (- q[2] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (- q[3] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0] 
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (- q[0] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (inv_sqrt_s - q[1] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (- q[2] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0])
+                            * (- q[3] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (- q[0] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (inv_sqrt_s - q[1] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (- q[2] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (- q[3] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (- q[0] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (inv_sqrt_s - q[1] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (- q[2] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (- q[3] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0] 
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (- q[0] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            + 
+                            (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (- q[1] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (inv_sqrt_s - q[2] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0]) 
+                            * (- q[3] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (- q[0] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (- q[1] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (inv_sqrt_s - q[2] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (- q[3] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (- q[0] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (- q[1] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (inv_sqrt_s - q[2] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (- q[3] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0] 
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (- q[0] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            + (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (- q[1] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (- q[2] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +                            
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0]) 
+                            * (inv_sqrt_s - q[3] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (- q[0] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (- q[1] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (- q[2] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (inv_sqrt_s - q[3] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (- q[0] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (- q[1] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (- q[2] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (inv_sqrt_s - q[3] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+        for (int i = 0; i < 3; ++i) {
+            data(m, 0, counter) = 0;
+            data(m, 1, counter) = 0;
+            data(m, 2, counter) = 0;
+            
+            counter++;
+        }
+    }
+}
+
+template<class Array>
+void CurvePlanarFourier<Array>::dgammadashdashdash_by_dcoeff_impl(Array& data) {
+    data *= 0;
+
+    Array q_norm = q * inv_magnitude();
+
+    double cosnphi, sinnphi;
+    for (int m = 0; m < numquadpoints; ++m) {
+        double phi = 2 * M_PI * quadpoints[m];
+        double cosphi = cos(phi);
+        double sinphi = sin(phi);
+        int counter = 0;
+        double i;
+        double j;
+        double k;
+        for (int n = 0; n < order+1; ++n) {
+            cosnphi = cos(n*phi);
+            sinnphi = sin(n*phi);
+            i = (
+                    +(3*n*n + 1)*cosnphi*sinphi
+                    +(n*n + 3)*(n)*sinnphi*cosphi
+                    );
+            j = (
+                    +(n*n + 3)*(n)*sinnphi*sinphi
+                    -(3*n*n + 1)*cosnphi*cosphi
+                    );
+            k = 0;
+
+            i *= 2*M_PI*2*M_PI*2*M_PI;
+            j *= 2*M_PI*2*M_PI*2*M_PI;
+            k *= 2*M_PI*2*M_PI*2*M_PI;
+
+            data(m, 0, counter) = (i - 2 * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * i + 2 * (q_norm[1] * q_norm[2] - q_norm[3] * q_norm[0]) * j + 2 * (q_norm[1] * q_norm[3] + q_norm[2] * q_norm[0]) * k);
+            data(m, 1, counter) = (2 * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * i + j - 2 * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * j + 2 * (q_norm[2] * q_norm[3] - q_norm[1] * q_norm[0]) * k);
+            data(m, 2, counter) = (2 * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * i + 2 * (q_norm[2] * q_norm[3] + q_norm[1] * q_norm[0]) * j + k - 2 * (q_norm[1] * q_norm[1] + q_norm[2] * q_norm[2]) * k);
+
+            counter++;
+        }
+        
+        for (int n = 1; n < order+1; ++n) {
+            cosnphi = cos(n*phi);
+            sinnphi = sin(n*phi);
+            i = (
+                    -(n*n+3) * (n) * cosnphi*cosphi
+                    +(3*n*n+1) * sinnphi*sinphi
+                    );
+            j = (
+                    -(n*n+3)*(n)*cosnphi*sinphi 
+                    -(3*n*n+1)*sinnphi*cosphi 
+                    );
+            k = 0;
+
+            i *= 2*M_PI*2*M_PI*2*M_PI;
+            j *= 2*M_PI*2*M_PI*2*M_PI;
+            k *= 2*M_PI*2*M_PI*2*M_PI;
+
+            data(m, 0, counter) = (i - 2 * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * i + 2 * (q_norm[1] * q_norm[2] - q_norm[3] * q_norm[0]) * j + 2 * (q_norm[1] * q_norm[3] + q_norm[2] * q_norm[0]) * k);
+            data(m, 1, counter) = (2 * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * i + j - 2 * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * j + 2 * (q_norm[2] * q_norm[3] - q_norm[1] * q_norm[0]) * k);
+            data(m, 2, counter) = (2 * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * i + 2 * (q_norm[2] * q_norm[3] + q_norm[1] * q_norm[0]) * j + k - 2 * (q_norm[1] * q_norm[1] + q_norm[2] * q_norm[2]) * k);
+
+            counter++;
+        }
+        
+        i = rc[0]*(sinphi);
+        j = rc[0]*(-cosphi);
+        k = 0;
+        for (int n = 1; n < order+1; ++n) {
+            cosnphi = cos(n*phi);
+            sinnphi = sin(n*phi);
+            i += rc[n]*(
+                    +(3*n*n + 1)*cosnphi*sinphi
+                    +(n*n + 3)*(n)*sinnphi*cosphi
+                    ) + rs[n-1]*(
+                    -(n*n+3) * (n) * cosnphi*cosphi
+                    +(3*n*n+1) * sinnphi*sinphi
+                    );
+            j += rc[n]*(
+                    +(n*n + 3)*(n)*sinnphi*sinphi
+                    -(3*n*n + 1)*cosnphi*cosphi
+                    ) + rs[n-1]*(
+                    -(n*n+3)*(n)*cosnphi*sinphi 
+                    -(3*n*n+1)*sinnphi*cosphi 
+                    );
+        }
+        i *= 2*M_PI*2*M_PI*2*M_PI;
+        j *= 2*M_PI*2*M_PI*2*M_PI;
+        k *= 2*M_PI*2*M_PI*2*M_PI;
+        double inv_sqrt_s = inv_magnitude();
+        
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (inv_sqrt_s - q[0] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (- q[1] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (- q[2] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0]) 
+                            * (- q[3] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (inv_sqrt_s - q[0] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (- q[1] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (- q[2] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (- q[3] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (inv_sqrt_s - q[0] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (- q[1] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (- q[2] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (- q[3] * q[0] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0] 
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (- q[0] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (inv_sqrt_s - q[1] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (- q[2] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0])
+                            * (- q[3] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (- q[0] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (inv_sqrt_s - q[1] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (- q[2] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (- q[3] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (- q[0] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (inv_sqrt_s - q[1] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (- q[2] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (- q[3] * q[1] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0] 
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (- q[0] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            + 
+                            (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (- q[1] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (inv_sqrt_s - q[2] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0]) 
+                            * (- q[3] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (- q[0] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (- q[1] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (inv_sqrt_s - q[2] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (- q[3] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (- q[0] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (- q[1] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (inv_sqrt_s - q[2] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (- q[3] * q[2] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+        data(m, 0, counter) = (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[0] 
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[0]
+                            - 2 * j * q_norm[3]) 
+                            * (- q[0] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            + (4 * i * (q_norm[2] * q_norm[2] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[2]) 
+                            * (- q[1] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[0] * q_norm[0] + q_norm[1] * q_norm[1]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[1]) 
+                            * (- q[2] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +                            
+                            (- 4 * i * (q_norm[1] * q_norm[1] + q_norm[0] * q_norm[0]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[2] - q_norm[0] * q_norm[3]) * q_norm[3]
+                            - 2 * j * q_norm[0]) 
+                            * (inv_sqrt_s - q[3] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 1, counter) = (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[0]
+                            + 2 * i * q_norm[3]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[0])
+                            * (- q[0] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[2]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[1]
+                            - 4 * j * q_norm[1])
+                            * (- q[1] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[2]
+                            + 2 * i * q_norm[1]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[2])
+                            * (- q[2] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[2] + q_norm[3] * q_norm[0]) * q_norm[3]
+                            + 2 * i * q_norm[0]
+                            + 4 * j * (q_norm[1] * q_norm[1] + q_norm[3] * q_norm[3]) * q_norm[3]
+                            - 4 * j * q_norm[3])
+                            * (inv_sqrt_s - q[3] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+
+
+        data(m, 2, counter) = (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[0]
+                            - 2 * i * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[0]
+                            + 2 * j * q_norm[1])
+                            * (- q[0] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 4 * i * (q_norm[1] * q_norm[3] - q_norm[2] * q_norm[0]) * q_norm[1]
+                            + 2 * i * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[1]
+                            + 2 * j * q_norm[0])
+                            * (- q[1] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (- 2 * i * q_norm[0]
+                            - 4 * i * (q_norm[1] * q_norm[3] - q_norm[0] * q_norm[2]) * q_norm[2]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[2]
+                            + 2 * j * q_norm[3])
+                            * (- q[2] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s)
+                            +
+                            (2 * i * q_norm[1]
+                            + 4 * i * (q_norm[2] * q_norm[0] - q_norm[1] * q_norm[3]) * q_norm[3]
+                            - 4 * j * (q_norm[1] * q_norm[0] + q_norm[2] * q_norm[3]) * q_norm[3]
+                            + 2 * j * q_norm[2])
+                            * (inv_sqrt_s - q[3] * q[3] * inv_sqrt_s * inv_sqrt_s * inv_sqrt_s);
+        counter++;
+
+        for (int i = 0; i < 3; ++i) {
+            data(m, 0, counter) = 0;
+            data(m, 1, counter) = 0;
+            data(m, 2, counter) = 0;
+
+            counter++;
+        }
+    }
+}
+
+
+
+#include "xtensor-python/pyarray.hpp"     // Numpy bindings
+typedef xt::pyarray<double> Array;
+template class CurvePlanarFourier<Array>;
diff --git a/src/simsoptpp/curveplanarfourier.h b/src/simsoptpp/curveplanarfourier.h
new file mode 100644
index 000000000..901bde313
--- /dev/null
+++ b/src/simsoptpp/curveplanarfourier.h
@@ -0,0 +1,121 @@
+#pragma once
+
+#include "curve.h"
+
+template<class Array>
+class CurvePlanarFourier : public Curve<Array> {
+    /*
+        CurvePlanarFourier is a curve that is restricted to lie in a plane. In
+        the plane, the curve is represented using a Fourier series in plane polar coordinates:
+
+           r(phi) = \sum_{n=0}^{order} x_{c,n}cos(n*nfp*phi) + \sum_{n=1}^order x_{s,n}sin(n*nfp*phi)
+        
+        The plane is rotated using a quarternion
+
+            q = [q_0, q_1, q_2, q_3] = [cos(\theta/2), x * sin(\theta/2), y * sin(\theta/2), z * sin(\theta/2)]
+
+        Details of the quaternion rotation can be found for example in pages
+        575-576 of https://www.cis.upenn.edu/~cis5150/ws-book-Ib.pdf.
+
+        The simsopt dofs for the quaternion need not generally have unit norm. The
+        quaternion is normalized before being applied to the curve to prevent scaling the curve.
+
+        A translation vector is used to specify the location of the center of the curve:
+
+            c = [c_x, c_y, c_z]
+
+        The dofs are stored in the order
+
+           [r_{c,0},...,r_{c,order},r_{s,1},...,r_{s,order}, q_0, q_1, q_2, q_3, c_x, c_y, c_z]
+
+       */
+    public:
+        const int order;
+        const int nfp;
+        const bool stellsym;
+        using Curve<Array>::quadpoints;
+        using Curve<Array>::numquadpoints;
+        using Curve<Array>::check_the_persistent_cache;
+
+        Array rc;
+        Array rs;
+        Array q;
+        Array center;
+
+        CurvePlanarFourier(int _numquadpoints, int _order, int _nfp, bool _stellsym) : Curve<Array>(_numquadpoints), order(_order), nfp(_nfp), stellsym(_stellsym) {
+            rc = xt::zeros<double>({order + 1});
+            rs = xt::zeros<double>({order});
+            q = xt::zeros<double>({4});
+            center = xt::zeros<double>({3});
+        }
+
+        CurvePlanarFourier(vector<double> _quadpoints, int _order, int _nfp, bool _stellsym) : Curve<Array>(_quadpoints), order(_order), nfp(_nfp), stellsym(_stellsym) {
+            rc = xt::zeros<double>({order + 1});
+            rs = xt::zeros<double>({order});
+            q = xt::zeros<double>({4});
+            center = xt::zeros<double>({3});
+        }
+
+        CurvePlanarFourier(Array _quadpoints, int _order, int _nfp, bool _stellsym) : Curve<Array>(_quadpoints), order(_order), nfp(_nfp), stellsym(_stellsym) {
+            rc = xt::zeros<double>({order + 1});
+            rs = xt::zeros<double>({order});
+            q = xt::zeros<double>({4});
+            center = xt::zeros<double>({3});
+        }
+
+        inline int num_dofs() override {
+            return (2*order+1)+7;
+        }
+
+        void set_dofs_impl(const vector<double>& dofs) override {
+            int counter = 0;
+            for (int i = 0; i < order + 1; ++i)
+                rc.data()[i] = dofs[counter++];
+            for (int i = 0; i < order; ++i)
+                rs.data()[i] = dofs[counter++];
+            for (int i = 0; i < 4; ++i)
+                q.data()[i] = dofs[counter++];
+            for (int i = 0; i < 3; ++i)
+                center.data()[i] = dofs[counter++];
+        }
+
+        vector<double> get_dofs() override {
+            auto res = vector<double>(num_dofs(), 0.);
+            int counter = 0;
+            for (int i = 0; i < order + 1; ++i)
+                res[counter++] = rc[i];
+            for (int i = 0; i < order; ++i)
+                res[counter++] = rs[i];
+            for (int i = 0; i < 4; ++i)
+                res[counter++] = q[i];
+            for (int i = 0; i < 3; ++i)
+                res[counter++] = center[i];
+            return res;
+        }
+
+        Array& dgamma_by_dcoeff() override {
+            return check_the_persistent_cache("dgamma_by_dcoeff", {numquadpoints, 3, num_dofs()}, [this](Array& A) { return dgamma_by_dcoeff_impl(A);});
+        }
+        Array& dgammadash_by_dcoeff() override {
+            return check_the_persistent_cache("dgammadash_by_dcoeff", {numquadpoints, 3, num_dofs()}, [this](Array& A) { return dgammadash_by_dcoeff_impl(A);});
+        }
+        Array& dgammadashdash_by_dcoeff() override {
+            return check_the_persistent_cache("dgammadashdash_by_dcoeff", {numquadpoints, 3, num_dofs()}, [this](Array& A) { return dgammadashdash_by_dcoeff_impl(A);});
+        }
+        Array& dgammadashdashdash_by_dcoeff() override {
+            return check_the_persistent_cache("dgammadashdashdash_by_dcoeff", {numquadpoints, 3, num_dofs()}, [this](Array& A) { return dgammadashdashdash_by_dcoeff_impl(A);});
+        }
+
+        void gamma_impl(Array& data, Array& quadpoints) override;
+        void gammadash_impl(Array& data) override;
+        void gammadashdash_impl(Array& data) override;
+        void gammadashdashdash_impl(Array& data) override;
+        void dgamma_by_dcoeff_impl(Array& data) override;
+        void dgammadash_by_dcoeff_impl(Array& data) override;
+        void dgammadashdash_by_dcoeff_impl(Array& data) override;
+        void dgammadashdashdash_by_dcoeff_impl(Array& data) override;
+
+private:
+        double inv_magnitude();
+
+};
diff --git a/src/simsoptpp/python_curves.cpp b/src/simsoptpp/python_curves.cpp
index 58adbecd3..90760ee4b 100644
--- a/src/simsoptpp/python_curves.cpp
+++ b/src/simsoptpp/python_curves.cpp
@@ -13,6 +13,8 @@ namespace py = pybind11;
 typedef CurveXYZFourier<PyArray> PyCurveXYZFourier;
 #include "curverzfourier.h"
 typedef CurveRZFourier<PyArray> PyCurveRZFourier; 
+#include "curveplanarfourier.h"
+typedef CurvePlanarFourier<PyArray> PyCurvePlanarFourier;
 
 template <class PyCurveXYZFourierBase = PyCurveXYZFourier> class PyCurveXYZFourierTrampoline : public PyCurveTrampoline<PyCurveXYZFourierBase> {
     public:
@@ -55,6 +57,27 @@ template <class PyCurveRZFourierBase = PyCurveRZFourier> class PyCurveRZFourierT
             PyCurveRZFourierBase::gamma_impl(data, quadpoints);
         }
 };
+
+template <class PyCurvePlanarFourierBase = PyCurvePlanarFourier> class PyCurvePlanarFourierTrampoline : public PyCurveTrampoline<PyCurvePlanarFourierBase> {
+    public:
+        using PyCurveTrampoline<PyCurvePlanarFourierBase>::PyCurveTrampoline; // Inherit constructors
+
+        int num_dofs() override {
+            return PyCurvePlanarFourierBase::num_dofs();
+        }
+
+        void set_dofs_impl(const vector<double>& _dofs) override {
+            PyCurvePlanarFourierBase::set_dofs_impl(_dofs);
+        }
+
+        vector<double> get_dofs() override {
+            return PyCurvePlanarFourierBase::get_dofs();
+        }
+
+        void gamma_impl(PyArray& data, PyArray& quadpoints) override {
+            PyCurvePlanarFourierBase::gamma_impl(data, quadpoints);
+        }
+};
 template <typename T, typename S> void register_common_curve_methods(S &c) {
     c.def("gamma", &T::gamma)
      .def("gamma_impl", &T::gamma_impl)
@@ -110,4 +133,15 @@ void init_curves(py::module_ &m) {
         .def_readonly("stellsym", &PyCurveRZFourier::stellsym)
         .def_readonly("nfp", &PyCurveRZFourier::nfp);
     register_common_curve_methods<PyCurveRZFourier>(pycurverzfourier);
+
+    auto pycurveplanarfourier = py::class_<PyCurvePlanarFourier, shared_ptr<PyCurvePlanarFourier>, PyCurvePlanarFourierTrampoline<PyCurvePlanarFourier>, PyCurve>(m, "CurvePlanarFourier")
+        .def(py::init<vector<double>, int, int, bool>())
+        .def_readwrite("rc", &PyCurvePlanarFourier::rc)
+        .def_readwrite("rs", &PyCurvePlanarFourier::rs)
+        .def_readwrite("q", &PyCurvePlanarFourier::q)
+        .def_readwrite("center", &PyCurvePlanarFourier::center)
+        .def_readonly("order", &PyCurvePlanarFourier::order)
+        .def_readonly("stellsym", &PyCurvePlanarFourier::stellsym)
+        .def_readonly("nfp", &PyCurvePlanarFourier::nfp);
+    register_common_curve_methods<PyCurvePlanarFourier>(pycurveplanarfourier);
 }
diff --git a/tests/field/test_coil.py b/tests/field/test_coil.py
index 0b1125b34..47dee47fa 100644
--- a/tests/field/test_coil.py
+++ b/tests/field/test_coil.py
@@ -7,8 +7,8 @@
 from simsopt.geo.curvexyzfourier import CurveXYZFourier, JaxCurveXYZFourier
 from simsopt.geo.curverzfourier import CurveRZFourier
 from simsopt.geo.curvehelical import CurveHelical
-from simsopt.geo.curve import RotatedCurve
-from simsopt.field.coil import Coil, Current, ScaledCurrent, CurrentSum
+from simsopt.geo.curve import RotatedCurve, create_equally_spaced_curves, create_equally_spaced_planar_curves
+from simsopt.field.coil import Coil, Current, ScaledCurrent, CurrentSum, coils_via_symmetries
 from simsopt.field.coil import coils_to_makegrid, coils_to_focus, load_coils_from_makegrid_file
 from simsopt.field.biotsavart import BiotSavart
 from simsopt._core.json import GSONEncoder, GSONDecoder, SIMSON
@@ -187,6 +187,36 @@ def test_load_coils_from_makegrid_file(self):
         np.testing.assert_allclose(B, loaded_B)
         np.testing.assert_allclose(gamma, loaded_gamma)
 
+    def test_equally_spaced_planar_curves(self):
+        ncoils = 4
+        nfp = 4
+        stellsym = False
+        R0 = 2.3
+        R1 = 0.9
+
+        curves = create_equally_spaced_curves(ncoils, nfp, stellsym, R0=R0, R1=R1)
+        currents = [Current(1e5) for i in range(ncoils)]
+
+        curves_planar = create_equally_spaced_planar_curves(ncoils, nfp, stellsym, R0=R0, R1=R1)
+        currents_planar = [Current(1e5) for i in range(ncoils)]
+
+        coils = coils_via_symmetries(curves, currents, nfp, stellsym)
+        coils_planar = coils_via_symmetries(curves_planar, currents_planar, nfp, stellsym)
+        bs = BiotSavart(coils)
+        bs_planar = BiotSavart(coils_planar)
+
+        x1d = np.linspace(R0, R0 + 0.3, 4)
+        y1d = np.linspace(0, 0.2, 3)
+        z1d = np.linspace(-0.2, 0.4, 5)
+        x, y, z = np.meshgrid(x1d, y1d, z1d)
+        points = np.ascontiguousarray(np.array([x.ravel(), y.ravel(), z.ravel()]).T)
+
+        bs.set_points(points)
+        bs_planar.set_points(points)
+
+        np.testing.assert_allclose(bs.B(), bs_planar.B(), atol=1e-16)
+
+
 
 if __name__ == "__main__":
     unittest.main()
diff --git a/tests/geo/test_curve.py b/tests/geo/test_curve.py
index ab92b5175..acf3314bf 100644
--- a/tests/geo/test_curve.py
+++ b/tests/geo/test_curve.py
@@ -9,6 +9,7 @@
 from simsopt._core.json import GSONEncoder, GSONDecoder, SIMSON
 from simsopt.geo.curvexyzfourier import CurveXYZFourier, JaxCurveXYZFourier
 from simsopt.geo.curverzfourier import CurveRZFourier
+from simsopt.geo.curveplanarfourier import CurvePlanarFourier
 from simsopt.geo.curvehelical import CurveHelical
 from simsopt.geo.curve import RotatedCurve, curves_to_vtk
 from simsopt.geo import parameters
@@ -45,7 +46,6 @@ def taylor_test(f, df, x, epsilons=None, direction=None):
         fminuseps = f(x - eps * direction)
         dfest = (fpluseps-fminuseps)/(2*eps)
         err = np.linalg.norm(dfest - dfx)
-        # print(err)
         assert err < 1e-9 or err < 0.3 * err_old
         if err < 1e-9:
             break
@@ -71,6 +71,8 @@ def get_curve(curvetype, rotated, x=np.asarray([0.5])):
         curve = CurveHelical(x, order, 5, 2, 1.0, 0.3)
     elif curvetype == "CurveHelicalInitx0":
         curve = CurveHelical(x, order, 5, 2, 1.0, 0.3, x0=np.ones((2*order,)))
+    elif curvetype == "CurvePlanarFourier":
+        curve = CurvePlanarFourier(x, order, 2, True)
     else:
         assert False
 
@@ -79,7 +81,7 @@ def get_curve(curvetype, rotated, x=np.asarray([0.5])):
         dofs[1] = 1.
         dofs[2*order + 3] = 1.
         dofs[4*order + 3] = 1.
-    elif curvetype in ["CurveRZFourier"]:
+    elif curvetype in ["CurveRZFourier", "CurvePlanarFourier"]:
         dofs[0] = 1.
         dofs[1] = 0.1
         dofs[order+1] = 0.1
@@ -96,8 +98,7 @@ def get_curve(curvetype, rotated, x=np.asarray([0.5])):
 
 class Testing(unittest.TestCase):
 
-    curvetypes = ["CurveXYZFourier", "JaxCurveXYZFourier", "CurveRZFourier", "CurveHelical",
-                  "CurveHelicalInitx0"]
+    curvetypes = ["CurveXYZFourier", "JaxCurveXYZFourier", "CurveRZFourier", "CurvePlanarFourier", "CurveHelical", "CurveHelicalInitx0"]
 
     def test_curve_helical_xyzfourier(self):
         x = np.asarray([0.6])