def METHOD_NAME(self, order):
    assert order in [1, None]
    curves, currents, ma = get_ncsx_data(Nt_coils=4, ppp=10)
    c = curves[0]
    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)
    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