From d9a6572831a6e489d248fae0eea919dee4aba054 Mon Sep 17 00:00:00 2001
From: Florian Wechsung <f.wechsung@gmail.com>
Date: Tue, 19 Apr 2022 11:34:17 -0400
Subject: [PATCH 01/14] copied over multifilament code from fw/multifil

---
 .../stage_two_optimization_finitebuild.py     | 151 +++++++++++
 src/simsopt/field/coil.py                     |  37 ++-
 src/simsopt/geo/curvefilament.py              | 251 ++++++++++++++++++
 src/simsopt/util/zoo.py                       |   1 +
 tests/geo/test_curvefilament.py               |  76 ++++++
 5 files changed, 505 insertions(+), 11 deletions(-)
 create mode 100755 examples/2_Intermediate/stage_two_optimization_finitebuild.py
 create mode 100644 src/simsopt/geo/curvefilament.py
 create mode 100644 tests/geo/test_curvefilament.py

diff --git a/examples/2_Intermediate/stage_two_optimization_finitebuild.py b/examples/2_Intermediate/stage_two_optimization_finitebuild.py
new file mode 100755
index 000000000..689ca4090
--- /dev/null
+++ b/examples/2_Intermediate/stage_two_optimization_finitebuild.py
@@ -0,0 +1,151 @@
+#!/usr/bin/env python
+r"""
+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 + alpha * (sum CurveLength) + beta * MininumDistancePenalty
+
+if alpha or beta are increased, the coils are more regular and better
+separated, but the target normal field may not be achieved as well.
+
+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.geo.surfacerzfourier import SurfaceRZFourier
+from simsopt.objectives.fluxobjective import SquaredFlux
+from simsopt.geo.curve import curves_to_vtk, create_equally_spaced_curves
+from simsopt.field.biotsavart import BiotSavart
+from simsopt.field.coil import Current, coils_via_symmetries
+from simsopt.geo.curveobjectives import CurveLength, MinimumDistance
+from simsopt.geo.curvefilament import create_multifilament_grid
+
+# 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:
+ALPHA = 1e-6
+
+# Threshhold for the coil-to-coil distance penalty in the objective function:
+MIN_DIST = 0.1
+
+# Weight on the coil-to-coil distance penalty term in the objective function:
+BETA = 10
+
+# Number of iterations to perform:
+ci = "CI" in os.environ and os.environ['CI'].lower() in ['1', 'true']
+MAXITER = 50 if ci 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:
+numfilaments_n = 1
+numfilaments_b = 5
+gapsize_n = 0.01
+gapsize_b = 0.01
+base_curves = create_equally_spaced_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()
+
+
+base_curves_finite_build = sum([
+    create_multifilament_grid(c, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b, rotation_order=3) for c in base_curves], [])
+base_currents_finite_build = sum([[c]*numfilaments_n*numfilaments_b for c in base_currents], [])
+
+
+coils = coils_via_symmetries(base_curves_finite_build, base_currents_finite_build, 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 objective function:
+Jf = SquaredFlux(s, bs)
+Jls = [CurveLength(c) for c in base_curves]
+curves_for_dist = [c.curve for c in coils[::(numfilaments_n*numfilaments_b)]]
+Jdist = MinimumDistance(curves_for_dist, MIN_DIST)
+
+# 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 + ALPHA * sum(Jls) + BETA * Jdist
+
+# 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()
+    cl_string = ", ".join([f"{J.J():.3f}" for J in Jls])
+    mean_AbsB = np.mean(bs.AbsB())
+    jf = Jf.J()
+    print(f"J={J:.3e}, Jflux={jf:.3e}, sqrt(Jflux)/Mean(|B|)={np.sqrt(jf)/mean_AbsB:.3e}, CoilLengths=[{cl_string}], ||∇J||={np.linalg.norm(grad):.3e}")
+    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': 400}, tol=1e-15)
+curves_to_vtk(curves, OUT_DIR + "curves_opt")
+pointData = {"B_N": np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2)[:, :, None]}
+s.to_vtk(OUT_DIR + "surf_opt", extra_data=pointData)
diff --git a/src/simsopt/field/coil.py b/src/simsopt/field/coil.py
index 1684de283..264bfec9c 100644
--- a/src/simsopt/field/coil.py
+++ b/src/simsopt/field/coil.py
@@ -72,6 +72,29 @@ def __neg__(self):
         return ScaledCurrent(self, -1.)
 
 
+def apply_symmetries_to_curves(base_curves, nfp, stellsym):
+    flip_list = [False, True] if stellsym else [False]
+    curves = []
+    for k in range(0, nfp):
+        for flip in flip_list:
+            for i in range(len(base_curves)):
+                if k == 0 and not flip:
+                    curves.append(base_curves[i])
+                else:
+                    rotcurve = RotatedCurve(base_curves[i], 2*pi*k/nfp, flip)
+                    curves.append(rotcurve)
+    return curves
+
+def apply_symmetries_to_currents(base_currents, nfp, stellsym):
+    flip_list = [False, True] if stellsym else [False]
+    currents = []
+    for k in range(0, nfp):
+        for flip in flip_list:
+            for i in range(len(base_currents)):
+                current = ScaledCurrent(base_currents[i], -1.) if flip else base_currents[i]
+                currents.append(current)
+    return currents
+
 def coils_via_symmetries(curves, currents, nfp, stellsym):
     """
     Take a list of ``n`` curves and return ``n * nfp * (1+int(stellsym))``
@@ -80,15 +103,7 @@ def coils_via_symmetries(curves, currents, nfp, stellsym):
     """
 
     assert len(curves) == len(currents)
-    flip_list = [False, True] if stellsym else [False]
-    coils = []
-    for k in range(0, nfp):
-        for flip in flip_list:
-            for i in range(len(curves)):
-                if k == 0 and not flip:
-                    coils.append(Coil(curves[i], currents[i]))
-                else:
-                    rotcurve = RotatedCurve(curves[i], 2*pi*k/nfp, flip)
-                    current = ScaledCurrent(currents[i], -1.) if flip else currents[i]
-                    coils.append(Coil(rotcurve, current))
+    curves = apply_symmetries_to_curves(curves, nfp, stellsym)
+    currents = apply_symmetries_to_currents(currents, nfp, stellsym)
+    coils = [Coil(curv, curr) for (curv, curr) in zip(curves, currents)]
     return coils
diff --git a/src/simsopt/geo/curvefilament.py b/src/simsopt/geo/curvefilament.py
new file mode 100644
index 000000000..0ed9d2797
--- /dev/null
+++ b/src/simsopt/geo/curvefilament.py
@@ -0,0 +1,251 @@
+import simsoptpp as sopp
+from .._core.graph_optimizable import Optimizable
+from .._core.derivative import Derivative
+from .curve import Curve
+from .jit import jit
+import numpy as np
+import jax.numpy as jnp
+from jax import vjp, jvp
+
+
+
+def create_multifilament_grid(curve, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b, rotation_order=None):
+    if numfilaments_n % 2 == 1:
+        shifts_n = np.arange(numfilaments_n) - numfilaments_n//2
+    else:
+        shifts_n = np.arange(numfilaments_n) - numfilaments_n/2 + 0.5
+    shifts_n = shifts_n * gapsize_n
+    if numfilaments_b % 2 == 1:
+        shifts_b = np.arange(numfilaments_b) - numfilaments_b//2
+    else:
+        shifts_b = np.arange(numfilaments_b) - numfilaments_b/2 + 0.5
+    shifts_b = shifts_b * gapsize_b
+
+    if rotation_order is None:
+        rotation = ZeroRotation(curve.quadpoints)
+    else:
+        rotation = FilamentRotation(curve.quadpoints, rotation_order)
+    filaments = []
+    for i in range(numfilaments_n):
+        for j in range(numfilaments_b):
+            filaments.append(CurveFilament(curve, shifts_n[i], shifts_b[j], rotation))
+    return filaments
+
+
+
+
+        
+
+class CurveFilament(sopp.Curve, Curve):
+
+    def __init__(self, curve, dn, db, rotation=None):
+        """
+        Implementation of the centroid frame introduced in
+        doi:10.1017/S0022377820000756. Given a curve, one defines a normal and
+        binormal vector and then creates a grid of curves by shifting along the
+        normal and binormal vector. In addition, we specify an angle along the
+        curve that allows us to optimise for the rotation of the winding pack.
+
+        The idea is explained well in Figure 1 in the reference above.
+
+        Args:
+            curve: the underlying curve
+            dn: how far to move in normal direction
+            db: how far to move in binormal direction
+            rotation: angle along the curve to rotate the frame.
+        """
+        self.curve = curve
+        sopp.Curve.__init__(self, curve.quadpoints)
+        deps = [curve]
+        if rotation is not None:
+            deps.append(rotation)
+        Curve.__init__(self, depends_on=deps)
+        self.curve = curve
+        self.dn = dn
+        self.db = db
+        if rotation is None:
+            rotation = ZeroRotation(curve.quadpoints)
+        self.rotation = rotation
+
+    def recompute_bell(self, parent=None):
+        self.invalidate_cache()
+
+    def gamma_impl(self, gamma, quadpoints):
+        assert quadpoints.shape[0] == self.curve.quadpoints.shape[0]
+        assert np.linalg.norm(quadpoints - self.curve.quadpoints) < 1e-15
+        c = self.curve
+        t, n, b = rotated_centroid_frame(c.gamma(), c.gammadash(), self.rotation.alpha(c.quadpoints))
+        gamma[:] = self.curve.gamma() + self.dn * n + self.db * b
+
+    def gammadash_impl(self, gammadash):
+        c = self.curve
+        td, nd, bd = rotated_centroid_frame_dash(
+            c.gamma(), c.gammadash(), c.gammadashdash(),
+            self.rotation.alpha(c.quadpoints), self.rotation.alphadash(c.quadpoints)
+        )
+        gammadash[:] = self.curve.gammadash() + self.dn * nd + self.db * bd
+
+    def dgamma_by_dcoeff_vjp(self, v):
+        g = self.curve.gamma()
+        gd = self.curve.gammadash()
+        a = self.rotation.alpha(self.curve.quadpoints)
+        zero = np.zeros_like(v)
+        vg = rotated_centroid_frame_dcoeff_vjp0(g, gd, a, (zero, self.dn*v, self.db*v))
+        vgd = rotated_centroid_frame_dcoeff_vjp1(g, gd, a, (zero, self.dn*v, self.db*v))
+        va = rotated_centroid_frame_dcoeff_vjp2(g, gd, a, (zero, self.dn*v, self.db*v))
+        return self.curve.dgamma_by_dcoeff_vjp(v + vg) \
+            + self.curve.dgammadash_by_dcoeff_vjp(vgd) \
+            + self.rotation.dalpha_by_dcoeff_vjp(self.curve.quadpoints, va)
+
+    def dgammadash_by_dcoeff_vjp(self, v):
+        g = self.curve.gamma()
+        gd = self.curve.gammadash()
+        gdd = self.curve.gammadashdash()
+        a = self.rotation.alpha(self.curve.quadpoints)
+        ad = self.rotation.alphadash(self.curve.quadpoints)
+        zero = np.zeros_like(v)
+
+        vg = rotated_centroid_frame_dash_dcoeff_vjp0(g, gd, gdd, a, ad, (zero, self.dn*v, self.db*v))
+        vgd = rotated_centroid_frame_dash_dcoeff_vjp1(g, gd, gdd, a, ad, (zero, self.dn*v, self.db*v))
+        vgdd = rotated_centroid_frame_dash_dcoeff_vjp2(g, gd, gdd, a, ad, (zero, self.dn*v, self.db*v))
+        va = rotated_centroid_frame_dash_dcoeff_vjp3(g, gd, gdd, a, ad, (zero, self.dn*v, self.db*v))
+        vad = rotated_centroid_frame_dash_dcoeff_vjp4(g, gd, gdd, a, ad, (zero, self.dn*v, self.db*v))
+        return self.curve.dgamma_by_dcoeff_vjp(vg) \
+            + self.curve.dgammadash_by_dcoeff_vjp(v+vgd) \
+            + self.curve.dgammadashdash_by_dcoeff_vjp(vgdd) \
+            + self.rotation.dalpha_by_dcoeff_vjp(self.curve.quadpoints, va) \
+            + self.rotation.dalphadash_by_dcoeff_vjp(self.curve.quadpoints, vad)
+
+
+class FilamentRotation(Optimizable):
+
+    def __init__(self, quadpoints, order):
+        """
+        The rotation of the multifilament pack; alpha in Figure 1 of
+        doi:10.1017/S0022377820000756
+        """
+        self.order = order
+        Optimizable.__init__(self, x0=np.zeros((2*order+1, )))
+        self.quadpoints = quadpoints
+        self.jac = rotation_dcoeff(quadpoints, order)
+        self.jacdash = rotationdash_dcoeff(quadpoints, order)
+        self.jax_alpha = jit(lambda dofs, points: jaxrotation_pure(dofs, points, self.order))
+        self.jax_alphadash = jit(lambda dofs, points: jaxrotationdash_pure(dofs, points, self.order))
+
+    def alpha(self, quadpoints):
+        return self.jax_alpha(self._dofs.full_x, quadpoints)
+
+    def alphadash(self, quadpoints):
+        return self.jax_alphadash(self._dofs.full_x, quadpoints)
+
+    def dalpha_by_dcoeff_vjp(self, quadpoints, v):
+        return Derivative({self: sopp.vjp(v, self.jac)})
+
+    def dalphadash_by_dcoeff_vjp(self, quadpoints, v):
+        return Derivative({self: sopp.vjp(v, self.jacdash)})
+
+
+class ZeroRotation():
+
+    def __init__(self, quadpoints):
+        self.zero = np.zeros((quadpoints.size, ))
+
+    def alpha(self, quadpoints):
+        return self.zero
+
+    def alphadash(self, quadpoints):
+        return self.zero
+
+    def dalpha_by_dcoeff_vjp(self, quadpoints, v):
+        return Derivative({})
+
+    def dalphadash_by_dcoeff_vjp(self, quadpoints, v):
+        return Derivative({})
+
+
+@jit
+def rotated_centroid_frame(gamma, gammadash, alpha):
+    t = gammadash
+    t *= 1./jnp.linalg.norm(gammadash, axis=1)[:, None]
+    R = jnp.mean(gamma, axis=0)  # centroid
+    delta = gamma - R[None, :]
+    n = delta - jnp.sum(delta * t, axis=1)[:, None] * t
+    n *= 1./jnp.linalg.norm(n, axis=1)[:, None]
+    b = jnp.cross(t, n, axis=1)
+
+    # now rotate the frame by alpha
+    nn = jnp.cos(alpha)[:, None] * n - jnp.sin(alpha)[:, None] * b
+    bb = jnp.sin(alpha)[:, None] * n + jnp.cos(alpha)[:, None] * b
+    return t, nn, bb
+
+
+rotated_centroid_frame_dash = jit(
+    lambda gamma, gammadash, gammadashdash, alpha, alphadash: jvp(rotated_centroid_frame,
+                                                                  (gamma, gammadash, alpha),
+                                                                  (gammadash, gammadashdash, alphadash))[1])
+
+rotated_centroid_frame_dcoeff_vjp0 = jit(
+    lambda gamma, gammadash, alpha, v: vjp(
+        lambda g: rotated_centroid_frame(g, gammadash, alpha), gamma)[1](v)[0])
+
+rotated_centroid_frame_dcoeff_vjp1 = jit(
+    lambda gamma, gammadash, alpha, v: vjp(
+        lambda gd: rotated_centroid_frame(gamma, gd, alpha), gammadash)[1](v)[0])
+
+rotated_centroid_frame_dcoeff_vjp2 = jit(
+    lambda gamma, gammadash, alpha, v: vjp(
+        lambda a: rotated_centroid_frame(gamma, gammadash, a), alpha)[1](v)[0])
+
+rotated_centroid_frame_dash_dcoeff_vjp0 = jit(
+    lambda gamma, gammadash, gammadashdash, alpha, alphadash, v: vjp(
+        lambda g: rotated_centroid_frame_dash(g, gammadash, gammadashdash, alpha, alphadash), gamma)[1](v)[0])
+
+rotated_centroid_frame_dash_dcoeff_vjp1 = jit(
+    lambda gamma, gammadash, gammadashdash, alpha, alphadash, v: vjp(
+        lambda gd: rotated_centroid_frame_dash(gamma, gd, gammadashdash, alpha, alphadash), gammadash)[1](v)[0])
+
+rotated_centroid_frame_dash_dcoeff_vjp2 = jit(
+    lambda gamma, gammadash, gammadashdash, alpha, alphadash, v: vjp(
+        lambda gdd: rotated_centroid_frame_dash(gamma, gammadash, gdd, alpha, alphadash), gammadashdash)[1](v)[0])
+
+rotated_centroid_frame_dash_dcoeff_vjp3 = jit(
+    lambda gamma, gammadash, gammadashdash, alpha, alphadash, v: vjp(
+        lambda a: rotated_centroid_frame_dash(gamma, gammadash, gammadashdash, a, alphadash), alpha)[1](v)[0])
+
+rotated_centroid_frame_dash_dcoeff_vjp4 = jit(
+    lambda gamma, gammadash, gammadashdash, alpha, alphadash, v: vjp(
+        lambda ad: rotated_centroid_frame_dash(gamma, gammadash, gammadashdash, alpha, ad), alphadash)[1](v)[0])
+
+
+def jaxrotation_pure(dofs, points, order):
+    rotation = jnp.zeros((len(points), ))
+    rotation += dofs[0]
+    for j in range(1, order+1):
+        rotation += dofs[2*j-1] * jnp.sin(2*np.pi*j*points)
+        rotation += dofs[2*j] * jnp.cos(2*np.pi*j*points)
+    return rotation
+
+
+def jaxrotationdash_pure(dofs, points, order):
+    rotation = jnp.zeros((len(points), ))
+    for j in range(1, order+1):
+        rotation += dofs[2*j-1] * 2*np.pi*j*jnp.cos(2*np.pi*j*points)
+        rotation -= dofs[2*j] * 2*np.pi*j*jnp.sin(2*np.pi*j*points)
+    return rotation
+
+
+def rotation_dcoeff(points, order):
+    jac = np.zeros((len(points), 2*order+1))
+    jac[:, 0] = 1
+    for j in range(1, order+1):
+        jac[:, 2*j-1] = np.sin(2*np.pi*j*points)
+        jac[:, 2*j+0] = np.cos(2*np.pi*j*points)
+    return jac
+
+
+def rotationdash_dcoeff(points, order):
+    jac = np.zeros((len(points), 2*order+1))
+    for j in range(1, order+1):
+        jac[:, 2*j-1] = +2*np.pi*j*np.cos(2*np.pi*j*points)
+        jac[:, 2*j+0] = -2*np.pi*j*np.sin(2*np.pi*j*points)
+    return jac
diff --git a/src/simsopt/util/zoo.py b/src/simsopt/util/zoo.py
index c570e0c00..493a517d7 100644
--- a/src/simsopt/util/zoo.py
+++ b/src/simsopt/util/zoo.py
@@ -42,4 +42,5 @@ def get_ncsx_data(Nt_coils=25, Nt_ma=10, ppp=10):
     ma = CurveRZFourier(numpoints, Nt_ma, nfp, True)
     ma.rc[:] = cR[0:(Nt_ma+1)]
     ma.zs[:] = sZ[0:Nt_ma]
+    ma.x = ma.get_dofs()
     return (curves, currents, ma)
diff --git a/tests/geo/test_curvefilament.py b/tests/geo/test_curvefilament.py
new file mode 100644
index 000000000..6bdd25d98
--- /dev/null
+++ b/tests/geo/test_curvefilament.py
@@ -0,0 +1,76 @@
+import unittest
+from simsopt.util.zoo import get_ncsx_data
+from simsopt.geo.curvefilament import CurveFilament, FilamentRotation
+import numpy as np
+
+
+class MultifilamentTesting(unittest.TestCase):
+
+    def test_multifilament_gammadash(self):
+        curves, currents, ma = get_ncsx_data(Nt_coils=6, ppp=80)
+        c = curves[0]
+
+        rotation = FilamentRotation(c.quadpoints, 1)
+        rotation.x = np.array([0, 0.1, 0.3])
+        c = CurveFilament(c, 0.01, 0.01, rotation)
+        g = c.gamma()
+        gd = c.gammadash()
+        idx = 16
+
+        dphi = c.quadpoints[1]
+        weights = [1/280, -4/105, 1/5, -4/5, 0, 4/5, -1/5, 4/105, -1/280]
+        est = 0
+        for j in range(-4, 5):
+            est += weights[j+4] * g[idx+j, :]
+        est *= 1./dphi
+        print(est)
+        print(gd[idx])
+        assert np.all(np.abs(est - gd[idx]) < 1e-10)
+
+    def test_multifilament_coefficient_derivative(self):
+
+        curves, currents, ma = get_ncsx_data(Nt_coils=4, ppp=10)
+        c = curves[0]
+
+        rotation = FilamentRotation(c.quadpoints, 1)
+        rotation.x = np.array([0, 0.1, 0.1])
+
+        c = CurveFilament(c, 0.02, 0.02, rotation)
+
+        dofs = c.x
+
+        g = c.gamma()
+        v = np.ones_like(g)
+        np.random.seed(1)
+
+        v = np.random.standard_normal(size=g.shape)
+        h = np.random.standard_normal(size=dofs.shape)
+        df = np.sum(c.dgamma_by_dcoeff_vjp(v)(c)*h)
+        dg = np.sum(c.dgammadash_by_dcoeff_vjp(v)(c)*h)
+
+        errf_old = 1e10
+        errg_old = 1e10
+
+        for i in range(12, 17):
+            eps = 0.5**i
+            c.x = dofs + eps*h
+            f1 = np.sum(c.gamma()*v)
+            c.x = dofs - eps*h
+            f2 = np.sum(c.gamma()*v)
+            errf = (f1-f2)/(2*eps) - df
+            print(errf)
+            assert errf < 0.3 * errf_old
+            errf_old = errf
+
+        print("==============")
+        for i in range(10, 17):
+            eps = 0.5**i
+            c.x = dofs + eps*h
+            g1 = np.sum(c.gammadash()*v)
+            c.x = dofs - eps*h
+            g2 = np.sum(c.gammadash()*v)
+            errg = (g1-g2)/(2*eps) - dg
+            # errg = (g1-g0)/(eps) - dg
+            print(errg)
+            assert errg < 0.3 * errg_old
+            errg_old = errg

From 2b26a5efd9ccec1d04672455a58281b93cd65e50 Mon Sep 17 00:00:00 2001
From: Florian Wechsung <f.wechsung@gmail.com>
Date: Thu, 21 Apr 2022 21:35:58 -0400
Subject: [PATCH 02/14] docs and scaling

---
 src/simsopt/geo/curvefilament.py | 39 ++++++++++++++++++++++----------
 1 file changed, 27 insertions(+), 12 deletions(-)

diff --git a/src/simsopt/geo/curvefilament.py b/src/simsopt/geo/curvefilament.py
index 0ed9d2797..f2509ee4b 100644
--- a/src/simsopt/geo/curvefilament.py
+++ b/src/simsopt/geo/curvefilament.py
@@ -9,7 +9,22 @@
 
 
 
-def create_multifilament_grid(curve, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b, rotation_order=None):
+def create_multifilament_grid(curve, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b, rotation_order=None, rotation_scaling=None):
+    """
+    Create a regular grid of `numfilaments_n * numfilaments_b` many filaments to approximate a finite-build coil.
+    Args:
+        curve: the underlying curve.
+        numfilaments_n: number of filaments in normal direction.
+        numfilaments_b: number of filaments in bi-normal direction.
+        gapsize_n: gap between filaments in normal direction.
+        gapsize_b: gap between filaments in bi-normal direction.
+        rotation_order: Fourier order to use in the expression for the rotation
+                        of the filament pack. `None` means that the rotation is not optimized.
+        rotation_scaling: scaling for the rotation degrees of freedom. good
+                           scaling improves the convergence of first order optimization
+                           algorithms. If None, then the default of `1/max(gapsize_n, gapsize_b)`
+                           is used.
+    """
     if numfilaments_n % 2 == 1:
         shifts_n = np.arange(numfilaments_n) - numfilaments_n//2
     else:
@@ -21,10 +36,12 @@ def create_multifilament_grid(curve, numfilaments_n, numfilaments_b, gapsize_n,
         shifts_b = np.arange(numfilaments_b) - numfilaments_b/2 + 0.5
     shifts_b = shifts_b * gapsize_b
 
+    if rotation_scaling is None:
+        rotation_scaling = 1/max(gapsize_n, gapsize_b)
     if rotation_order is None:
         rotation = ZeroRotation(curve.quadpoints)
     else:
-        rotation = FilamentRotation(curve.quadpoints, rotation_order)
+        rotation = FilamentRotation(curve.quadpoints, rotation_order, scale=rotation_scaling)
     filaments = []
     for i in range(numfilaments_n):
         for j in range(numfilaments_b):
@@ -32,10 +49,6 @@ def create_multifilament_grid(curve, numfilaments_n, numfilaments_b, gapsize_n,
     return filaments
 
 
-
-
-        
-
 class CurveFilament(sopp.Curve, Curve):
 
     def __init__(self, curve, dn, db, rotation=None):
@@ -119,7 +132,7 @@ def dgammadash_by_dcoeff_vjp(self, v):
 
 class FilamentRotation(Optimizable):
 
-    def __init__(self, quadpoints, order):
+    def __init__(self, quadpoints, order, scale=1.):
         """
         The rotation of the multifilament pack; alpha in Figure 1 of
         doi:10.1017/S0022377820000756
@@ -127,27 +140,29 @@ def __init__(self, quadpoints, order):
         self.order = order
         Optimizable.__init__(self, x0=np.zeros((2*order+1, )))
         self.quadpoints = quadpoints
+        self.scale = scale
         self.jac = rotation_dcoeff(quadpoints, order)
         self.jacdash = rotationdash_dcoeff(quadpoints, order)
         self.jax_alpha = jit(lambda dofs, points: jaxrotation_pure(dofs, points, self.order))
         self.jax_alphadash = jit(lambda dofs, points: jaxrotationdash_pure(dofs, points, self.order))
 
     def alpha(self, quadpoints):
-        return self.jax_alpha(self._dofs.full_x, quadpoints)
+        return self.scale * self.jax_alpha(self._dofs.full_x, quadpoints)
 
     def alphadash(self, quadpoints):
-        return self.jax_alphadash(self._dofs.full_x, quadpoints)
+        return self.scale * self.jax_alphadash(self._dofs.full_x, quadpoints)
 
     def dalpha_by_dcoeff_vjp(self, quadpoints, v):
-        return Derivative({self: sopp.vjp(v, self.jac)})
+        return Derivative({self: self.scale * sopp.vjp(v, self.jac)})
 
     def dalphadash_by_dcoeff_vjp(self, quadpoints, v):
-        return Derivative({self: sopp.vjp(v, self.jacdash)})
+        return Derivative({self: self.scale * sopp.vjp(v, self.jacdash)})
 
 
-class ZeroRotation():
+class ZeroRotation(Optimizable):
 
     def __init__(self, quadpoints):
+        Optimizable.__init__(self, x0=[])
         self.zero = np.zeros((quadpoints.size, ))
 
     def alpha(self, quadpoints):

From 151d5a166e6890453b9c7924f4dfc88ec702eab1 Mon Sep 17 00:00:00 2001
From: Florian Wechsung <f.wechsung@gmail.com>
Date: Thu, 21 Apr 2022 21:41:19 -0400
Subject: [PATCH 03/14] linting

---
 src/simsopt/field/coil.py        | 2 ++
 src/simsopt/geo/curvefilament.py | 1 -
 2 files changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/simsopt/field/coil.py b/src/simsopt/field/coil.py
index 264bfec9c..462ab05fd 100644
--- a/src/simsopt/field/coil.py
+++ b/src/simsopt/field/coil.py
@@ -85,6 +85,7 @@ def apply_symmetries_to_curves(base_curves, nfp, stellsym):
                     curves.append(rotcurve)
     return curves
 
+
 def apply_symmetries_to_currents(base_currents, nfp, stellsym):
     flip_list = [False, True] if stellsym else [False]
     currents = []
@@ -95,6 +96,7 @@ def apply_symmetries_to_currents(base_currents, nfp, stellsym):
                 currents.append(current)
     return currents
 
+
 def coils_via_symmetries(curves, currents, nfp, stellsym):
     """
     Take a list of ``n`` curves and return ``n * nfp * (1+int(stellsym))``
diff --git a/src/simsopt/geo/curvefilament.py b/src/simsopt/geo/curvefilament.py
index f2509ee4b..e3899f055 100644
--- a/src/simsopt/geo/curvefilament.py
+++ b/src/simsopt/geo/curvefilament.py
@@ -8,7 +8,6 @@
 from jax import vjp, jvp
 
 
-
 def create_multifilament_grid(curve, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b, rotation_order=None, rotation_scaling=None):
     """
     Create a regular grid of `numfilaments_n * numfilaments_b` many filaments to approximate a finite-build coil.

From acc04724d95a90bc835e3c162677da5619bfa641 Mon Sep 17 00:00:00 2001
From: Florian Wechsung <f.wechsung@gmail.com>
Date: Fri, 22 Apr 2022 17:36:47 -0400
Subject: [PATCH 04/14] add more testing

---
 tests/geo/test_curvefilament.py | 42 ++++++++++++++++++++++++++++++++-
 1 file changed, 41 insertions(+), 1 deletion(-)

diff --git a/tests/geo/test_curvefilament.py b/tests/geo/test_curvefilament.py
index 6bdd25d98..5ddcff04e 100644
--- a/tests/geo/test_curvefilament.py
+++ b/tests/geo/test_curvefilament.py
@@ -1,6 +1,6 @@
 import unittest
 from simsopt.util.zoo import get_ncsx_data
-from simsopt.geo.curvefilament import CurveFilament, FilamentRotation
+from simsopt.geo.curvefilament import CurveFilament, FilamentRotation, create_multifilament_grid
 import numpy as np
 
 
@@ -74,3 +74,43 @@ def test_multifilament_coefficient_derivative(self):
             print(errg)
             assert errg < 0.3 * errg_old
             errg_old = errg
+
+    def test_filamentpack(self):
+        curves, currents, ma = get_ncsx_data(Nt_coils=6, ppp=80)
+        c = curves[0]
+
+        numfilaments_n = 2
+        numfilaments_b = 3
+        gapsize_n = 0.01
+        gapsize_b = 0.02
+        fils = create_multifilament_grid(
+            c, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b,
+            rotation_order=None, rotation_scaling=None)
+
+        def check(fils, c):
+            assert len(fils) == numfilaments_n * numfilaments_b
+            dists = np.linalg.norm(fils[0].gamma()-fils[-1].gamma(), axis=1)
+            # check that filaments are equidistant everywhere
+            assert np.var(dists) < 1e-16
+            # check that first and last filament are on opossing corners of filament pack and have correct distance
+            assert abs(dists[0] - (((numfilaments_n-1)*gapsize_n)**2+((numfilaments_b-1)*gapsize_b)**2)**0.5) < 1e-13
+            # check that the coil pack is centered around the underlying curve
+            assert np.linalg.norm(np.mean([f.gamma() for f in fils], axis=0)-c.gamma()) < 1e-13
+
+        check(fils, c)
+
+        fils = create_multifilament_grid(
+            c, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b,
+            rotation_order=3, rotation_scaling=None)
+
+        xr = fils[0].rotation.x
+        fils[0].rotation.x = xr + 1e-2*np.random.standard_normal(size=xr.shape)
+
+        check(fils, c)
+
+        fils = create_multifilament_grid(
+            c, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b,
+            rotation_order=3, rotation_scaling=None)
+
+        xr = fils[0].rotation.x
+        fils[0].rotation.x = xr + 1e-2*np.random.standard_normal(size=xr.shape)

From 0bc3b71323b892c3a1303c88c929d6038cfa0be4 Mon Sep 17 00:00:00 2001
From: Florian Wechsung <f.wechsung@gmail.com>
Date: Mon, 25 Apr 2022 11:20:17 -0400
Subject: [PATCH 05/14] more tests

---
 tests/geo/test_curvefilament.py | 58 +++++++++++++++++++++------------
 1 file changed, 38 insertions(+), 20 deletions(-)

diff --git a/tests/geo/test_curvefilament.py b/tests/geo/test_curvefilament.py
index 5ddcff04e..3482be29d 100644
--- a/tests/geo/test_curvefilament.py
+++ b/tests/geo/test_curvefilament.py
@@ -1,17 +1,28 @@
 import unittest
 from simsopt.util.zoo import get_ncsx_data
-from simsopt.geo.curvefilament import CurveFilament, FilamentRotation, create_multifilament_grid
+from simsopt.geo.curvefilament import CurveFilament, FilamentRotation, \
+    create_multifilament_grid, ZeroRotation
 import numpy as np
 
 
 class MultifilamentTesting(unittest.TestCase):
 
     def test_multifilament_gammadash(self):
+        for order in [None, 1]:
+            with self.subTest(order=order):
+                self.subtest_multifilament_gammadash(order)
+
+    def subtest_multifilament_gammadash(self, order):
+        assert order in [1, None]
         curves, currents, ma = get_ncsx_data(Nt_coils=6, ppp=80)
         c = curves[0]
 
-        rotation = FilamentRotation(c.quadpoints, 1)
-        rotation.x = np.array([0, 0.1, 0.3])
+        if order == 1:
+            rotation = FilamentRotation(c.quadpoints, order)
+            rotation.x = np.array([0, 0.1, 0.3])
+        else:
+            rotation = ZeroRotation(c.quadpoints)
+
         c = CurveFilament(c, 0.01, 0.01, rotation)
         g = c.gamma()
         gd = c.gammadash()
@@ -28,12 +39,21 @@ def test_multifilament_gammadash(self):
         assert np.all(np.abs(est - gd[idx]) < 1e-10)
 
     def test_multifilament_coefficient_derivative(self):
+        for order in [None, 1]:
+            with self.subTest(order=order):
+                self.subtest_multifilament_coefficient_derivative(order)
+
+    def subtest_multifilament_coefficient_derivative(self, order):
+        assert order in [1, None]
 
         curves, currents, ma = get_ncsx_data(Nt_coils=4, ppp=10)
         c = curves[0]
 
-        rotation = FilamentRotation(c.quadpoints, 1)
-        rotation.x = np.array([0, 0.1, 0.1])
+        if order == 1:
+            rotation = FilamentRotation(c.quadpoints, order)
+            rotation.x = np.array([0, 0.1, 0.3])
+        else:
+            rotation = ZeroRotation(c.quadpoints)
 
         c = CurveFilament(c, 0.02, 0.02, rotation)
 
@@ -79,38 +99,36 @@ def test_filamentpack(self):
         curves, currents, ma = get_ncsx_data(Nt_coils=6, ppp=80)
         c = curves[0]
 
-        numfilaments_n = 2
-        numfilaments_b = 3
         gapsize_n = 0.01
         gapsize_b = 0.02
-        fils = create_multifilament_grid(
-            c, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b,
-            rotation_order=None, rotation_scaling=None)
 
-        def check(fils, c):
+        def check(fils, c, numfilaments_n, numfilaments_b):
             assert len(fils) == numfilaments_n * numfilaments_b
             dists = np.linalg.norm(fils[0].gamma()-fils[-1].gamma(), axis=1)
             # check that filaments are equidistant everywhere
             assert np.var(dists) < 1e-16
             # check that first and last filament are on opossing corners of filament pack and have correct distance
-            assert abs(dists[0] - (((numfilaments_n-1)*gapsize_n)**2+((numfilaments_b-1)*gapsize_b)**2)**0.5) < 1e-13
+            assert abs(dists[0] - (((numfilaments_n-1)*gapsize_n) ** 2+((numfilaments_b-1)*gapsize_b)**2)**0.5) < 1e-13
             # check that the coil pack is centered around the underlying curve
             assert np.linalg.norm(np.mean([f.gamma() for f in fils], axis=0)-c.gamma()) < 1e-13
 
-        check(fils, c)
-
+        numfilaments_n = 2
+        numfilaments_b = 3
         fils = create_multifilament_grid(
             c, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b,
-            rotation_order=3, rotation_scaling=None)
-
-        xr = fils[0].rotation.x
-        fils[0].rotation.x = xr + 1e-2*np.random.standard_normal(size=xr.shape)
+            rotation_order=None, rotation_scaling=None)
+        check(fils, c, numfilaments_n, numfilaments_b)
 
-        check(fils, c)
+        numfilaments_n = 3
+        numfilaments_b = 2
+        fils = create_multifilament_grid(
+            c, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b,
+            rotation_order=None, rotation_scaling=None)
+        check(fils, c, numfilaments_n, numfilaments_b)
 
         fils = create_multifilament_grid(
             c, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b,
             rotation_order=3, rotation_scaling=None)
-
         xr = fils[0].rotation.x
         fils[0].rotation.x = xr + 1e-2*np.random.standard_normal(size=xr.shape)
+        check(fils, c, numfilaments_n, numfilaments_b)

From d8ec332491566ff60c10198c6ba101eeeb033869 Mon Sep 17 00:00:00 2001
From: Florian Wechsung <f.wechsung@gmail.com>
Date: Mon, 25 Apr 2022 14:16:21 -0400
Subject: [PATCH 06/14] clean up example

---
 .../stage_two_optimization_finitebuild.py     | 135 +++++++++++-------
 examples/run_serial_examples                  |   1 +
 2 files changed, 82 insertions(+), 54 deletions(-)
 rename examples/{2_Intermediate => 3_Advanced}/stage_two_optimization_finitebuild.py (50%)

diff --git a/examples/2_Intermediate/stage_two_optimization_finitebuild.py b/examples/3_Advanced/stage_two_optimization_finitebuild.py
similarity index 50%
rename from examples/2_Intermediate/stage_two_optimization_finitebuild.py
rename to examples/3_Advanced/stage_two_optimization_finitebuild.py
index 689ca4090..aad55e2e9 100755
--- a/examples/2_Intermediate/stage_two_optimization_finitebuild.py
+++ b/examples/3_Advanced/stage_two_optimization_finitebuild.py
@@ -1,58 +1,76 @@
 #!/usr/bin/env python
 r"""
-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.
+In this example we solve a FOCUS like Stage II coil optimisation problem for finite build coils.
+We approximate each finite build coil using a multifilament approach. To model
+the multilament pack we follow the approach of
+
+    Optimization of finite-build stellarator coils,
+    Singh, Luquant, et al.  Journal of Plasma Physics 86.4 (2020).
+
+This means, that in addition to the degrees of freedom for the shape of the
+coils, we have additional degrees of freedom for the rotation of the coil pack.
 
 The objective is given by
 
-    J = (1/2) \int |B dot n|^2 ds + alpha * (sum CurveLength) + beta * MininumDistancePenalty
+    J = (1/2) ∫ |B_{BiotSavart}·n - B_{External}·n|^2 ds
+        + LENGTH_PEN * Σ ½(CurveLength - L0)^2
+        + DIST_PEN * PairwiseDistancePenalty
 
-if alpha or beta are increased, the coils are more regular and better
-separated, but the target normal field may not be achieved as well.
+The target equilibrium is the QA configuration of
+
+    Magnetic fields with precise quasisymmetry for plasma confinement,
+    Landreman, M., & Paul, E. (2022), Physical Review Letters, 128(3), 035001.
 
-The target equilibrium is the QA configuration of arXiv:2108.03711.
 """
 
 import os
-from pathlib import Path
 import numpy as np
+from pathlib import Path
 from scipy.optimize import minimize
-from simsopt.geo.surfacerzfourier import SurfaceRZFourier
-from simsopt.objectives.fluxobjective import SquaredFlux
-from simsopt.geo.curve import curves_to_vtk, create_equally_spaced_curves
 from simsopt.field.biotsavart import BiotSavart
-from simsopt.field.coil import Current, coils_via_symmetries
-from simsopt.geo.curveobjectives import CurveLength, MinimumDistance
+from simsopt.field.coil import Current, ScaledCurrent, Coil, apply_symmetries_to_curves, apply_symmetries_to_currents
+from simsopt.geo.curve import curves_to_vtk, create_equally_spaced_curves
 from simsopt.geo.curvefilament import create_multifilament_grid
+from simsopt.geo.curveobjectives import CurveLength, CurveCurveDistance, LpCurveCurvature
+from simsopt.geo.surfacerzfourier import SurfaceRZFourier
+from simsopt.objectives.fluxobjective import SquaredFlux
+from simsopt.objectives.utilities import QuadraticPenalty
 
 # 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
+R0 = 1.00
 
 # Minor radius for the initial circular coils:
-R1 = 0.5
+R1 = 0.70
 
 # Number of Fourier modes describing each Cartesian component of each coil:
 order = 5
 
-# Weight on the curve lengths in the objective function:
-ALPHA = 1e-6
+# Weight on the curve length penalty in the objective function:
+LENGTH_PEN = 1e-2
 
-# Threshhold for the coil-to-coil distance penalty in the objective function:
-MIN_DIST = 0.1
+# Threshhold and weight for the coil-to-coil distance penalty in the objective function:
+DIST_MIN = 0.1
+DIST_PEN = 10
 
-# Weight on the coil-to-coil distance penalty term in the objective function:
-BETA = 10
+# Settings for multifilament approximation
+numfilaments_n = 2  # number of filaments in normal direction
+numfilaments_b = 3  # number of filaments in bi-normal direction
+gapsize_n = 0.02  # gap between filaments in normal direction
+gapsize_b = 0.04  # gap between filaments in bi-normal direction
+rot_order = 1  # order of the Fourier expression for the rotation of the filament pack
 
 # Number of iterations to perform:
 ci = "CI" in os.environ and os.environ['CI'].lower() in ['1', 'true']
 MAXITER = 50 if ci else 400
 
+#######################################################
+# End of input parameters.
+#######################################################
+
 # File for the desired boundary magnetic surface:
 TEST_DIR = (Path(__file__).parent / ".." / ".." / "tests" / "test_files").resolve()
 filename = TEST_DIR / 'input.LandremanPaul2021_QA'
@@ -61,52 +79,58 @@
 OUT_DIR = "./output/"
 os.makedirs(OUT_DIR, exist_ok=True)
 
-#######################################################
-# End of input parameters.
-#######################################################
+config_str = f"rot_order_{rot_order}_nfn_{numfilaments_n}_nfb_{numfilaments_b}"
 
 # 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:
-numfilaments_n = 1
-numfilaments_b = 5
-gapsize_n = 0.01
-gapsize_b = 0.01
+nfil = numfilaments_n * numfilaments_b
 base_curves = create_equally_spaced_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()
-
-
+base_currents = []
+for i in range(ncoils):
+    curr = Current(1.)
+    # 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
+    if i == 0:
+        curr.fix_all()
+    base_currents.append(ScaledCurrent(curr, 1e5/nfil))
+
+# use sum here to concatenate lists
 base_curves_finite_build = sum([
-    create_multifilament_grid(c, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b, rotation_order=3) for c in base_curves], [])
-base_currents_finite_build = sum([[c]*numfilaments_n*numfilaments_b for c in base_currents], [])
-
-
-coils = coils_via_symmetries(base_curves_finite_build, base_currents_finite_build, s.nfp, True)
-bs = BiotSavart(coils)
+    create_multifilament_grid(c, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b, rotation_order=rot_order) for c in base_curves], start=[])
+base_currents_finite_build = sum([[c]*nfil for c in base_currents], start=[])
+
+# apply stellarator and rotation symmetries
+curves_fb = apply_symmetries_to_curves(base_curves_finite_build, s.nfp, True)
+currents_fb = apply_symmetries_to_currents(base_currents_finite_build, s.nfp, True)
+# also apply symmetries to the underlying base curves, as we use those in the
+# curve-curve distance penalty
+curves = apply_symmetries_to_curves(base_curves, s.nfp, True)
+
+coils_fb = [Coil(c, curr) for (c, curr) in zip(curves_fb, currents_fb)]
+bs = BiotSavart(coils_fb)
 bs.set_points(s.gamma().reshape((-1, 3)))
 
-curves = [c.curve for c in coils]
 curves_to_vtk(curves, OUT_DIR + "curves_init")
+curves_to_vtk(curves_fb, OUT_DIR + f"curves_init_fb_{config_str}")
+
 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)
+s.to_vtk(OUT_DIR + f"surf_init_fb_{config_str}", extra_data=pointData)
 
 # Define the objective function:
 Jf = SquaredFlux(s, bs)
 Jls = [CurveLength(c) for c in base_curves]
-curves_for_dist = [c.curve for c in coils[::(numfilaments_n*numfilaments_b)]]
-Jdist = MinimumDistance(curves_for_dist, MIN_DIST)
+Jdist = CurveCurveDistance(curves, DIST_MIN)
 
 # 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 + ALPHA * sum(Jls) + BETA * Jdist
+JF = Jf \
+    + LENGTH_PEN * sum(QuadraticPenalty(Jls[i], Jls[i].J()) for i in range(len(base_curves))) \
+    + DIST_PEN * Jdist
 
 # 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
@@ -120,8 +144,9 @@ def fun(dofs):
     cl_string = ", ".join([f"{J.J():.3f}" for J in Jls])
     mean_AbsB = np.mean(bs.AbsB())
     jf = Jf.J()
-    print(f"J={J:.3e}, Jflux={jf:.3e}, sqrt(Jflux)/Mean(|B|)={np.sqrt(jf)/mean_AbsB:.3e}, CoilLengths=[{cl_string}], ||∇J||={np.linalg.norm(grad):.3e}")
-    return J, grad
+    kap_string = ", ".join(f"{np.max(c.kappa()):.1f}" for c in base_curves)
+    print(f"J={J:.3e}, Jflux={jf:.3e}, sqrt(Jflux)/Mean(|B|)={np.sqrt(jf)/mean_AbsB:.3e}, CoilLengths=[{cl_string}], [{kap_string}], ||∇J||={np.linalg.norm(grad):.3e}")
+    return 1e-4*J, 1e-4*grad
 
 
 print("""
@@ -135,7 +160,7 @@ def fun(dofs):
 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]:
+for eps in [1e-3, 1e-4, 1e-5, 1e-6, 1e-7, 1e-8]:
     J1, _ = f(dofs + eps*h)
     J2, _ = f(dofs - eps*h)
     print("err", (J1-J2)/(2*eps) - dJh)
@@ -145,7 +170,9 @@ def fun(dofs):
 ### Run the optimisation #######################################################
 ################################################################################
 """)
-res = minimize(fun, dofs, jac=True, method='L-BFGS-B', options={'maxiter': MAXITER, 'maxcor': 400}, tol=1e-15)
-curves_to_vtk(curves, OUT_DIR + "curves_opt")
+
+res = minimize(fun, dofs, jac=True, method='L-BFGS-B', options={'maxiter': MAXITER, 'maxcor': 400, 'gtol': 1e-20, 'ftol': 1e-20}, tol=1e-20)
+
+curves_to_vtk(curves_fb, OUT_DIR + f"curves_opt_fb_{config_str}")
 pointData = {"B_N": np.sum(bs.B().reshape((nphi, ntheta, 3)) * s.unitnormal(), axis=2)[:, :, None]}
-s.to_vtk(OUT_DIR + "surf_opt", extra_data=pointData)
+s.to_vtk(OUT_DIR + f"surf_opt_fb_{config_str}", extra_data=pointData)
diff --git a/examples/run_serial_examples b/examples/run_serial_examples
index 9dbb3447d..c80260fa2 100755
--- a/examples/run_serial_examples
+++ b/examples/run_serial_examples
@@ -14,3 +14,4 @@ set -ex
 ./2_Intermediate/boozer.py
 ./2_Intermediate/stage_two_optimization.py
 ./2_Intermediate/stage_two_optimization_stochastic.py
+./3_Intermediate/stage_two_optimization_finitebuild.py

From 20cbd0159b9a009728367f979ab2be0b5e5f2224 Mon Sep 17 00:00:00 2001
From: Florian Wechsung <f.wechsung@gmail.com>
Date: Mon, 25 Apr 2022 14:21:09 -0400
Subject: [PATCH 07/14] docs and renaming

---
 src/simsopt/field/coil.py                            | 11 +++++++++++
 src/simsopt/geo/{curvefilament.py => finitebuild.py} |  4 ++++
 2 files changed, 15 insertions(+)
 rename src/simsopt/geo/{curvefilament.py => finitebuild.py} (98%)

diff --git a/src/simsopt/field/coil.py b/src/simsopt/field/coil.py
index 4bc593248..7afbff0d2 100644
--- a/src/simsopt/field/coil.py
+++ b/src/simsopt/field/coil.py
@@ -73,6 +73,12 @@ def __neg__(self):
 
 
 def apply_symmetries_to_curves(base_curves, nfp, stellsym):
+    """
+    Take a list of ``n`` :mod:`simsopt.geo.curve.Curve`s and return ``n * nfp *
+    (1+int(stellsym))`` :mod:`simsopt.geo.curve.Curve` objects obtained by
+    applying rotations and flipping corresponding to ``nfp`` fold rotational
+    symmetry and optionally stellarator symmetry.
+    """
     flip_list = [False, True] if stellsym else [False]
     curves = []
     for k in range(0, nfp):
@@ -87,6 +93,11 @@ def apply_symmetries_to_curves(base_curves, nfp, stellsym):
 
 
 def apply_symmetries_to_currents(base_currents, nfp, stellsym):
+    """
+    Take a list of ``n`` :mod:`Current`s and return ``n * nfp * (1+int(stellsym))``
+    :mod:`Current` objects obtained by copying (for ``nfp`` rotations) and
+    sign-flipping (optionally for stellarator symmetry).
+    """
     flip_list = [False, True] if stellsym else [False]
     currents = []
     for k in range(0, nfp):
diff --git a/src/simsopt/geo/curvefilament.py b/src/simsopt/geo/finitebuild.py
similarity index 98%
rename from src/simsopt/geo/curvefilament.py
rename to src/simsopt/geo/finitebuild.py
index e3899f055..218700a4b 100644
--- a/src/simsopt/geo/curvefilament.py
+++ b/src/simsopt/geo/finitebuild.py
@@ -6,6 +6,10 @@
 import numpy as np
 import jax.numpy as jnp
 from jax import vjp, jvp
+"""
+The functions and classes in this model are used to deal with multifilament
+approximation of finite build coils.
+"""
 
 
 def create_multifilament_grid(curve, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b, rotation_order=None, rotation_scaling=None):

From 2a0cdb57f034166f61182ecddbf0f40861956592 Mon Sep 17 00:00:00 2001
From: Florian Wechsung <f.wechsung@gmail.com>
Date: Mon, 25 Apr 2022 16:33:03 -0400
Subject: [PATCH 08/14] fix renaming

---
 examples/3_Advanced/stage_two_optimization_finitebuild.py | 2 +-
 tests/geo/{test_curvefilament.py => test_finitebuild.py}  | 2 +-
 2 files changed, 2 insertions(+), 2 deletions(-)
 rename tests/geo/{test_curvefilament.py => test_finitebuild.py} (98%)

diff --git a/examples/3_Advanced/stage_two_optimization_finitebuild.py b/examples/3_Advanced/stage_two_optimization_finitebuild.py
index aad55e2e9..c9097faf2 100755
--- a/examples/3_Advanced/stage_two_optimization_finitebuild.py
+++ b/examples/3_Advanced/stage_two_optimization_finitebuild.py
@@ -30,7 +30,7 @@
 from simsopt.field.biotsavart import BiotSavart
 from simsopt.field.coil import Current, ScaledCurrent, Coil, apply_symmetries_to_curves, apply_symmetries_to_currents
 from simsopt.geo.curve import curves_to_vtk, create_equally_spaced_curves
-from simsopt.geo.curvefilament import create_multifilament_grid
+from simsopt.geo.finitebuild import create_multifilament_grid
 from simsopt.geo.curveobjectives import CurveLength, CurveCurveDistance, LpCurveCurvature
 from simsopt.geo.surfacerzfourier import SurfaceRZFourier
 from simsopt.objectives.fluxobjective import SquaredFlux
diff --git a/tests/geo/test_curvefilament.py b/tests/geo/test_finitebuild.py
similarity index 98%
rename from tests/geo/test_curvefilament.py
rename to tests/geo/test_finitebuild.py
index 3482be29d..99e0a5044 100644
--- a/tests/geo/test_curvefilament.py
+++ b/tests/geo/test_finitebuild.py
@@ -1,6 +1,6 @@
 import unittest
 from simsopt.util.zoo import get_ncsx_data
-from simsopt.geo.curvefilament import CurveFilament, FilamentRotation, \
+from simsopt.geo.finitebuild import CurveFilament, FilamentRotation, \
     create_multifilament_grid, ZeroRotation
 import numpy as np
 

From 651a428172718d0330a6d8593f4cacff3bdc175f Mon Sep 17 00:00:00 2001
From: Florian Wechsung <f.wechsung@gmail.com>
Date: Mon, 25 Apr 2022 19:59:13 -0400
Subject: [PATCH 09/14] fix

---
 examples/run_serial_examples | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/examples/run_serial_examples b/examples/run_serial_examples
index c80260fa2..43f8172d0 100755
--- a/examples/run_serial_examples
+++ b/examples/run_serial_examples
@@ -14,4 +14,4 @@ set -ex
 ./2_Intermediate/boozer.py
 ./2_Intermediate/stage_two_optimization.py
 ./2_Intermediate/stage_two_optimization_stochastic.py
-./3_Intermediate/stage_two_optimization_finitebuild.py
+./3_Advanced/stage_two_optimization_finitebuild.py

From 14825c3bbc8f10891cf2108c8e637eb115a73458 Mon Sep 17 00:00:00 2001
From: Florian Wechsung <f.wechsung@gmail.com>
Date: Mon, 25 Apr 2022 20:56:42 -0400
Subject: [PATCH 10/14] fix

---
 examples/3_Advanced/stage_two_optimization_finitebuild.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/examples/3_Advanced/stage_two_optimization_finitebuild.py b/examples/3_Advanced/stage_two_optimization_finitebuild.py
index c9097faf2..836f6d045 100755
--- a/examples/3_Advanced/stage_two_optimization_finitebuild.py
+++ b/examples/3_Advanced/stage_two_optimization_finitebuild.py
@@ -100,8 +100,8 @@
 
 # use sum here to concatenate lists
 base_curves_finite_build = sum([
-    create_multifilament_grid(c, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b, rotation_order=rot_order) for c in base_curves], start=[])
-base_currents_finite_build = sum([[c]*nfil for c in base_currents], start=[])
+    create_multifilament_grid(c, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b, rotation_order=rot_order) for c in base_curves], [])
+base_currents_finite_build = sum([[c]*nfil for c in base_currents], [])
 
 # apply stellarator and rotation symmetries
 curves_fb = apply_symmetries_to_curves(base_curves_finite_build, s.nfp, True)

From 07a6571574311cedd4607a8224d4a6577c5ed4df Mon Sep 17 00:00:00 2001
From: Florian Wechsung <f.wechsung@gmail.com>
Date: Tue, 26 Apr 2022 13:50:16 -0400
Subject: [PATCH 11/14] more complicated test

---
 tests/geo/test_finitebuild.py | 56 +++++++++++++++++++++++++++++++++++
 1 file changed, 56 insertions(+)

diff --git a/tests/geo/test_finitebuild.py b/tests/geo/test_finitebuild.py
index 99e0a5044..07ba4dc0b 100644
--- a/tests/geo/test_finitebuild.py
+++ b/tests/geo/test_finitebuild.py
@@ -132,3 +132,59 @@ def check(fils, c, numfilaments_n, numfilaments_b):
         xr = fils[0].rotation.x
         fils[0].rotation.x = xr + 1e-2*np.random.standard_normal(size=xr.shape)
         check(fils, c, numfilaments_n, numfilaments_b)
+
+    def test_biotsavart_with_symmetries(self):
+        from .surface_test_helpers import get_surface, get_exact_surface
+        from simsopt.field.biotsavart import BiotSavart
+        from simsopt.field.coil import Coil, apply_symmetries_to_curves, apply_symmetries_to_currents
+        from simsopt.geo.qfmsurface import QfmSurface
+        from simsopt.util.zoo import get_ncsx_data
+        from simsopt.geo.curveobjectives import CurveLength, CurveCurveDistance
+        from simsopt.objectives.fluxobjective import SquaredFlux
+        from simsopt.objectives.utilities import QuadraticPenalty
+
+        np.random.seed(1)
+        base_curves, base_currents, ma = get_ncsx_data(Nt_coils=5)
+        base_curves_finite_build = sum(
+            [create_multifilament_grid(c, 2, 2, 0.01, 0.01, rotation_order=1) for c in base_curves], [])
+        base_currents_finite_build = sum([[c]*4 for c in base_currents], [])
+
+        nfp = 3
+
+        curves = apply_symmetries_to_curves(base_curves, nfp, True)
+        curves_fb = apply_symmetries_to_curves(base_curves_finite_build, nfp, True)
+        currents_fb = apply_symmetries_to_currents(base_currents_finite_build, nfp, True)
+
+        coils_fb = [Coil(c, curr) for (c, curr) in zip(curves_fb, currents_fb)]
+
+        bs = BiotSavart(coils_fb)
+        s = get_surface("SurfaceXYZFourier", True)
+        s.fit_to_curve(ma, 0.1)
+        Jf = SquaredFlux(s, bs)
+        Jls = [CurveLength(c) for c in base_curves]
+        Jdist = CurveCurveDistance(curves, 0.5)
+        LENGTH_PEN = 1e-2
+        DIST_PEN = 1e-2
+        JF = Jf \
+            + LENGTH_PEN * sum(QuadraticPenalty(Jls[i], Jls[i].J()) for i in range(len(base_curves))) \
+            + DIST_PEN * Jdist
+
+        def fun(dofs, grad=True):
+            JF.x = dofs
+            return (JF.J(), JF.dJ()) if grad else JF.J()
+
+        dofs = JF.x
+        dofs += 1e-2 * np.random.standard_normal(size=dofs.shape)
+        np.random.seed(1)
+        h = np.random.uniform(size=dofs.shape)
+        J0, dJ0 = fun(dofs)
+        dJh = sum(dJ0 * h)
+        err = 1e6
+        for i in range(10, 15):
+            eps = 0.5**i
+            J1 = fun(dofs + eps*h, grad=False)
+            J2 = fun(dofs - eps*h, grad=False)
+            err_new = abs((J1-J2)/(2*eps) - dJh)
+            assert err_new < 0.55**2 * err
+            err = err_new
+            print("err", err)

From 149ec5ac7e62f7c7105789ea9fbded0ec2106fe4 Mon Sep 17 00:00:00 2001
From: Florian Wechsung <f.wechsung@gmail.com>
Date: Tue, 26 Apr 2022 13:53:15 -0400
Subject: [PATCH 12/14] doc

---
 tests/geo/test_finitebuild.py | 25 +++++++++++++++----------
 1 file changed, 15 insertions(+), 10 deletions(-)

diff --git a/tests/geo/test_finitebuild.py b/tests/geo/test_finitebuild.py
index 07ba4dc0b..de335fac5 100644
--- a/tests/geo/test_finitebuild.py
+++ b/tests/geo/test_finitebuild.py
@@ -1,7 +1,16 @@
 import unittest
-from simsopt.util.zoo import get_ncsx_data
+from .surface_test_helpers import get_surface, get_exact_surface
+from simsopt.field.biotsavart import BiotSavart
+from simsopt.field.coil import Coil, apply_symmetries_to_curves, apply_symmetries_to_currents
+from simsopt.geo.curveobjectives import CurveLength, CurveCurveDistance
 from simsopt.geo.finitebuild import CurveFilament, FilamentRotation, \
     create_multifilament_grid, ZeroRotation
+from simsopt.geo.qfmsurface import QfmSurface
+from simsopt.objectives.fluxobjective import SquaredFlux
+from simsopt.objectives.utilities import QuadraticPenalty
+from simsopt.util.zoo import get_ncsx_data
+from simsopt.util.zoo import get_ncsx_data
+
 import numpy as np
 
 
@@ -134,15 +143,11 @@ def check(fils, c, numfilaments_n, numfilaments_b):
         check(fils, c, numfilaments_n, numfilaments_b)
 
     def test_biotsavart_with_symmetries(self):
-        from .surface_test_helpers import get_surface, get_exact_surface
-        from simsopt.field.biotsavart import BiotSavart
-        from simsopt.field.coil import Coil, apply_symmetries_to_curves, apply_symmetries_to_currents
-        from simsopt.geo.qfmsurface import QfmSurface
-        from simsopt.util.zoo import get_ncsx_data
-        from simsopt.geo.curveobjectives import CurveLength, CurveCurveDistance
-        from simsopt.objectives.fluxobjective import SquaredFlux
-        from simsopt.objectives.utilities import QuadraticPenalty
-
+        """
+        More involved test that checks whether the multifilament code interacts
+        properly with symmetries, biot savart, and objectives that only depend
+        on the underlying curve (not the finite build filaments)
+        """
         np.random.seed(1)
         base_curves, base_currents, ma = get_ncsx_data(Nt_coils=5)
         base_curves_finite_build = sum(

From 1c285bebf8b667d22f186ffe55c8fd309203f382 Mon Sep 17 00:00:00 2001
From: Florian Wechsung <f.wechsung@gmail.com>
Date: Tue, 26 Apr 2022 13:58:38 -0400
Subject: [PATCH 13/14] doc

---
 src/simsopt/geo/finitebuild.py | 8 ++++++++
 1 file changed, 8 insertions(+)

diff --git a/src/simsopt/geo/finitebuild.py b/src/simsopt/geo/finitebuild.py
index 218700a4b..78d18e34c 100644
--- a/src/simsopt/geo/finitebuild.py
+++ b/src/simsopt/geo/finitebuild.py
@@ -165,6 +165,14 @@ def dalphadash_by_dcoeff_vjp(self, quadpoints, v):
 class ZeroRotation(Optimizable):
 
     def __init__(self, quadpoints):
+        """
+        Dummy class that just returns zero for the rotation angle. Equivalent to using 
+        .. code-block:: python
+
+            rot = FilamentRotation(...)
+            rot.fix_all()
+
+        """
         Optimizable.__init__(self, x0=[])
         self.zero = np.zeros((quadpoints.size, ))
 

From 5d267211f3167bebe1f5182f8f0f273c60ed0eef Mon Sep 17 00:00:00 2001
From: Matt Landreman <mattland@umd.edu>
Date: Wed, 27 Apr 2022 16:29:54 -0400
Subject: [PATCH 14/14] Minor additions to and formatting of comments

---
 docs/source/simsopt.geo.rst                   |  8 ++++++
 .../stage_two_optimization_finitebuild.py     |  7 +++--
 src/simsopt/geo/finitebuild.py                | 28 +++++++++++++++----
 3 files changed, 35 insertions(+), 8 deletions(-)

diff --git a/docs/source/simsopt.geo.rst b/docs/source/simsopt.geo.rst
index b4733a4eb..9417ad506 100644
--- a/docs/source/simsopt.geo.rst
+++ b/docs/source/simsopt.geo.rst
@@ -68,6 +68,14 @@ simsopt.geo.curvexyzfourier module
    :undoc-members:
    :show-inheritance:
 
+simsopt.geo.finitebuild module
+------------------------------
+
+.. automodule:: simsopt.geo.finitebuild
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
 simsopt.geo.jit module
 ----------------------
 
diff --git a/examples/3_Advanced/stage_two_optimization_finitebuild.py b/examples/3_Advanced/stage_two_optimization_finitebuild.py
index 836f6d045..02094b882 100755
--- a/examples/3_Advanced/stage_two_optimization_finitebuild.py
+++ b/examples/3_Advanced/stage_two_optimization_finitebuild.py
@@ -56,12 +56,15 @@
 DIST_MIN = 0.1
 DIST_PEN = 10
 
-# Settings for multifilament approximation
+# Settings for multifilament approximation.  In the following
+# parameters, note that "normal" and "binormal" refer not to the
+# Frenet frame but rather to the "coil centroid frame" defined by
+# Singh et al., before rotation.
 numfilaments_n = 2  # number of filaments in normal direction
 numfilaments_b = 3  # number of filaments in bi-normal direction
 gapsize_n = 0.02  # gap between filaments in normal direction
 gapsize_b = 0.04  # gap between filaments in bi-normal direction
-rot_order = 1  # order of the Fourier expression for the rotation of the filament pack
+rot_order = 1  # order of the Fourier expression for the rotation of the filament pack, i.e. maximum Fourier mode number
 
 # Number of iterations to perform:
 ci = "CI" in os.environ and os.environ['CI'].lower() in ['1', 'true']
diff --git a/src/simsopt/geo/finitebuild.py b/src/simsopt/geo/finitebuild.py
index 78d18e34c..e20e177fe 100644
--- a/src/simsopt/geo/finitebuild.py
+++ b/src/simsopt/geo/finitebuild.py
@@ -6,6 +6,7 @@
 import numpy as np
 import jax.numpy as jnp
 from jax import vjp, jvp
+
 """
 The functions and classes in this model are used to deal with multifilament
 approximation of finite build coils.
@@ -14,18 +15,24 @@
 
 def create_multifilament_grid(curve, numfilaments_n, numfilaments_b, gapsize_n, gapsize_b, rotation_order=None, rotation_scaling=None):
     """
-    Create a regular grid of `numfilaments_n * numfilaments_b` many filaments to approximate a finite-build coil.
+    Create a regular grid of ``numfilaments_n * numfilaments_b`` many
+    filaments to approximate a finite-build coil.
+
+    Note that "normal" and "binormal" in the function arguments here
+    refer not to the Frenet frame but rather to the "coil centroid
+    frame" defined by Singh et al., before rotation.
+
     Args:
-        curve: the underlying curve.
+        curve: The underlying curve.
         numfilaments_n: number of filaments in normal direction.
         numfilaments_b: number of filaments in bi-normal direction.
         gapsize_n: gap between filaments in normal direction.
         gapsize_b: gap between filaments in bi-normal direction.
-        rotation_order: Fourier order to use in the expression for the rotation
-                        of the filament pack. `None` means that the rotation is not optimized.
+        rotation_order: Fourier order (maximum mode number) to use in the expression for the rotation
+                        of the filament pack. ``None`` means that the rotation is not optimized.
         rotation_scaling: scaling for the rotation degrees of freedom. good
                            scaling improves the convergence of first order optimization
-                           algorithms. If None, then the default of `1/max(gapsize_n, gapsize_b)`
+                           algorithms. If ``None``, then the default of ``1 / max(gapsize_n, gapsize_b)``
                            is used.
     """
     if numfilaments_n % 2 == 1:
@@ -57,6 +64,8 @@ class CurveFilament(sopp.Curve, Curve):
     def __init__(self, curve, dn, db, rotation=None):
         """
         Implementation of the centroid frame introduced in
+        Singh et al, "Optimization of finite-build stellarator coils",
+        Journal of Plasma Physics 86 (2020),
         doi:10.1017/S0022377820000756. Given a curve, one defines a normal and
         binormal vector and then creates a grid of curves by shifting along the
         normal and binormal vector. In addition, we specify an angle along the
@@ -64,6 +73,10 @@ def __init__(self, curve, dn, db, rotation=None):
 
         The idea is explained well in Figure 1 in the reference above.
 
+        Note that "normal" and "binormal" in the function arguments here
+        refer not to the Frenet frame but rather to the "coil centroid
+        frame" defined by Singh et al., before rotation.
+
         Args:
             curve: the underlying curve
             dn: how far to move in normal direction
@@ -138,6 +151,8 @@ class FilamentRotation(Optimizable):
     def __init__(self, quadpoints, order, scale=1.):
         """
         The rotation of the multifilament pack; alpha in Figure 1 of
+        Singh et al, "Optimization of finite-build stellarator coils",
+        Journal of Plasma Physics 86 (2020),
         doi:10.1017/S0022377820000756
         """
         self.order = order
@@ -166,7 +181,8 @@ class ZeroRotation(Optimizable):
 
     def __init__(self, quadpoints):
         """
-        Dummy class that just returns zero for the rotation angle. Equivalent to using 
+        Dummy class that just returns zero for the rotation angle. Equivalent to using
+
         .. code-block:: python
 
             rot = FilamentRotation(...)