def test_energy_vs_gradient_sinewave():
    n = 8
    w = 1 / np.pi
    Es = 3. / 4
    cf = SphereCrackFrontERRPenetrationEnergy(n,
                                              w_radius_integral=lambda a, theta: a ** 2 / 2 * w * 2 * np.pi / n,
                                              w_radius=lambda a, theta: a * w * 2 * np.pi / n,
                                              dw_radius=lambda a, theta: w * 2 * np.pi / n, )
    a = 0.65 * np.ones(n)
    # k = np.fft.fftfreq(len(a), 1 / len(a))
    penetration = JKR.penetration(contact_radius=0.6)
    # contact radius perturbation
    da = np.random.uniform(-1, 1, size=n)
    angle = np.arange(n) * 2 * np.pi / n
    da = 0.1 * np.sin(3 * angle)
    epsilons = np.logspace(-3, -9)
    dUel = np.array([cf.energy(a + eps * da, penetration) - cf.energy(a, penetration) for eps in epsilons])

    grad_da = np.vdot(cf.gradient(a, penetration), da)

    import matplotlib.pyplot as plt
    fig, ax = plt.subplots()
    ax.loglog(epsilons,abs(dUel/(grad_da * epsilons)- 1),".")
    # ax.plot(epsilons,)
    #ax.plot(epsilons, dUel / epsilons, ".")
    #ax.axhline(grad_da)

    plt.show(block=True)

    @staticmethod
    def _n_an_2(contact_radius):
        npx = len(contact_radius)
        nq = np.fft.rfftfreq(npx, 1 / npx)
        fourier_scalar_prod_factors = np.ones(npx // 2 + 1) * 2
        fourier_scalar_prod_factors[0] = 1
        if npx % 2 == 0:
            fourier_scalar_prod_factors[-1] == 1
        a_fourier = np.fft.rfft(contact_radius, norm="forward")
        return np.vdot(a_fourier * nq * fourier_scalar_prod_factors, a_fourier).real