diff --git a/python/wrap/core.cpp b/python/wrap/core.cpp
index f473799..e61c995 100644
--- a/python/wrap/core.cpp
+++ b/python/wrap/core.cpp
@@ -1,117 +1,120 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "logger.hh"
 #include "statistics.hh"
 #include "wrap.hh"
 #include <pybind11/stl.h>
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 /* -------------------------------------------------------------------------- */
 namespace wrap {
 
 using namespace py::literals;
 
 /* -------------------------------------------------------------------------- */
 template <UInt dim>
 void wrapStatistics(py::module& mod) {
   auto name = makeDimensionName("Statistics", dim);
   py::class_<Statistics<dim>>(mod, name.c_str())
       .def_static("computePowerSpectrum",
                   &Statistics<dim>::computePowerSpectrum,
                   py::return_value_policy::move, "Compute PSD of surface")
       .def_static(
           "computeAutocorrelation", &Statistics<dim>::computeAutocorrelation,
           py::return_value_policy::move, "Compute autocorrelation of surface")
       .def_static("computeMoments", &Statistics<dim>::computeMoments,
                   "Compute spectral moments")
       .def_static("computeSpectralRMSSlope",
                   &Statistics<dim>::computeSpectralRMSSlope,
                   "Compute hrms' in Fourier space")
       .def_static("computeRMSHeights", &Statistics<dim>::computeRMSHeights,
                   "Compute hrms")
       .def_static("contact", &Statistics<dim>::contact, "tractions"_a,
                   "perimeter"_a = 0,
                   "Compute the (corrected) contact area. Permieter is the "
-                  "total contact perimeter in number of segments.");
+                  "total contact perimeter in number of segments.")
+      .def_static("graphArea", &Statistics<dim>::graphArea, "zdisplacement"_a,
+                  "Compute area defined by function graph.");
 }
 
 /* -------------------------------------------------------------------------- */
 void wrapCore(py::module& mod) {
   // Exposing logging facility
   py::enum_<LogLevel>(mod, "LogLevel")
       .value("debug", LogLevel::debug)
       .value("info", LogLevel::info)
       .value("warning", LogLevel::warning)
       .value("error", LogLevel::error);
   mod.def("set_log_level", [](LogLevel level) { Logger::setLevel(level); });
   mod.def("get_log_level", Logger::getCurrentLevel);
 
   py::class_<Logger>(mod, "Logger")
       .def(py::init<>())
-      .def("get",
-           [](Logger& logger, LogLevel level) -> Logger& {
-             logger.get(level);
-             return logger;
-           },
-           "Get a logger object for a log level")
-      .def("__lshift__",
-           [](Logger& logger, std::string msg) -> Logger& {
-             auto level = logger.getWishLevel();
-             if (level >= Logger::getCurrentLevel() and
-                 not(mpi::rank() != 0 and level == LogLevel::info)) {
-               py::print(msg,
-                         "file"_a = py::module::import("sys").attr("stderr"));
-             }
-             return logger;
-           },
-           "Log a message");
+      .def(
+          "get",
+          [](Logger& logger, LogLevel level) -> Logger& {
+            logger.get(level);
+            return logger;
+          },
+          "Get a logger object for a log level")
+      .def(
+          "__lshift__",
+          [](Logger& logger, std::string msg) -> Logger& {
+            auto level = logger.getWishLevel();
+            if (level >= Logger::getCurrentLevel() and
+                not(mpi::rank() != 0 and level == LogLevel::info)) {
+              py::print(msg,
+                        "file"_a = py::module::import("sys").attr("stderr"));
+            }
+            return logger;
+          },
+          "Log a message");
 
   wrapStatistics<1>(mod);
   wrapStatistics<2>(mod);
 
-  mod.def("to_voigt",
-          [](const Grid<Real, 3>& field) {
-            TAMAAS_DEPRECATE("tamaas.to_voigt()", "tamaas.compute.to_voigt()");
+  mod.def(
+      "to_voigt",
+      [](const Grid<Real, 3>& field) {
+        TAMAAS_DEPRECATE("tamaas.to_voigt()", "tamaas.compute.to_voigt()");
 
-            if (field.getNbComponents() == 9) {
-              Grid<Real, 3> voigt(field.sizes(), 6);
-              Loop::loop(
-                  [] CUDA_LAMBDA(MatrixProxy<const Real, 3, 3> in,
-                                 SymMatrixProxy<Real, 3> out) {
-                    out.symmetrize(in);
-                  },
-                  range<MatrixProxy<const Real, 3, 3>>(field),
-                  range<SymMatrixProxy<Real, 3>>(voigt));
-              return voigt;
-            }
+        if (field.getNbComponents() == 9) {
+          Grid<Real, 3> voigt(field.sizes(), 6);
+          Loop::loop([] CUDA_LAMBDA(
+                         MatrixProxy<const Real, 3, 3> in,
+                         SymMatrixProxy<Real, 3> out) { out.symmetrize(in); },
+                     range<MatrixProxy<const Real, 3, 3>>(field),
+                     range<SymMatrixProxy<Real, 3>>(voigt));
+          return voigt;
+        }
 
-            TAMAAS_EXCEPTION("Wrong number of components to symmetrize");
-          },
-          "Convert a 3D tensor field to voigt notation",
-          py::return_value_policy::copy);
+        TAMAAS_EXCEPTION("Wrong number of components to symmetrize");
+      },
+      "Convert a 3D tensor field to voigt notation",
+      py::return_value_policy::copy);
 }
 
 }  // namespace wrap
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
diff --git a/src/core/statistics.cpp b/src/core/statistics.cpp
index 5cc8bfc..e466dc5 100644
--- a/src/core/statistics.cpp
+++ b/src/core/statistics.cpp
@@ -1,211 +1,249 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "statistics.hh"
 #include "fft_engine.hh"
 #include "loop.hh"
 #include "static_types.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 template <UInt dim>
 Real Statistics<dim>::computeRMSHeights(Grid<Real, dim>& surface) {
   return std::sqrt(surface.var());
 }
 
 template <UInt dim>
 Real Statistics<dim>::computeSpectralRMSSlope(Grid<Real, dim>& surface) {
   const auto h_size =
       GridHermitian<Real, dim>::hermitianDimensions(surface.sizes());
   auto wavevectors =
       FFTEngine::template computeFrequencies<Real, dim, true>(h_size);
   wavevectors *= 2 * M_PI;  // need q for slopes
   const auto psd = computePowerSpectrum(surface);
 
   const Real rms_slope_mean = Loop::reduce<operation::plus>(
       [] CUDA_LAMBDA(VectorProxy<const Real, dim> q, const Complex& psd_val) {
         // Checking if we're in the zone that does not have hermitian symmetry
         if (std::abs(q.back()) < 1e-15)
           return q.l2squared() * psd_val.real();
         return 2 * q.l2squared() * psd_val.real();
       },
       range<VectorProxy<const Real, dim>>(wavevectors), psd);
 
   return std::sqrt(rms_slope_mean);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 GridHermitian<Real, dim>
 Statistics<dim>::computePowerSpectrum(Grid<Real, dim>& surface) {
   const auto h_size =
       GridHermitian<Real, dim>::hermitianDimensions(surface.sizes());
   GridHermitian<Real, dim> psd(h_size, surface.getNbComponents());
 
   FFTEngine::makeEngine()->forward(surface, psd);
 
   Real factor = 1. / surface.getGlobalNbPoints();
 
   // Squaring the fourier transform of surface and normalizing
   Loop::loop(
       [factor] CUDA_LAMBDA(Complex & c) {
         c *= factor;
         c *= conj(c);
       },
       psd);
 
   return psd;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt dim>
 Grid<Real, dim>
 Statistics<dim>::computeAutocorrelation(Grid<Real, dim>& surface) {
   Grid<Real, dim> acf(surface.sizes(), surface.getNbComponents());
   auto psd = computePowerSpectrum(surface);
 
   FFTEngine::makeEngine()->backward(acf, psd);
   acf *= acf.getGlobalNbPoints();
   return acf;
 }
 
 /* -------------------------------------------------------------------------- */
 template <UInt dim>
 Real Statistics<dim>::contact(const Grid<Real, dim>& tractions,
                               UInt perimeter) {
   Real points = 0;
   UInt nc = tractions.getNbComponents();
 
   switch (nc) {
   case 1:
     points = Loop::reduce(
         [] CUDA_LAMBDA(const Real& t) -> Real { return t > 0; }, tractions);
     break;
   case 2:
     points = Loop::reduce(
         [] CUDA_LAMBDA(VectorProxy<const Real, 2> t) -> Real {
           return t.back() > 0;
         },
         range<VectorProxy<const Real, 2>>(tractions));
     break;
   case 3:
     points = Loop::reduce(
         [] CUDA_LAMBDA(VectorProxy<const Real, 3> t) -> Real {
           return t.back() > 0;
         },
         range<VectorProxy<const Real, 3>>(tractions));
     break;
   default:
     TAMAAS_EXCEPTION("Invalid number of components in traction");
   }
 
   auto area = points / tractions.getGlobalNbPoints();
 
   if (dim == 1)
     perimeter = 0;
 
   // Correction from Yastrebov et al. (Trib. Intl., 2017)
   // 10.1016/j.triboint.2017.04.023
   return area - (M_PI - 1 + std::log(2)) /
                     (24. * tractions.getGlobalNbPoints()) * perimeter;
 }
 
 /* -------------------------------------------------------------------------- */
 namespace {
 template <UInt dim>
 class moment_helper {
 public:
   moment_helper(const std::array<UInt, dim>& exp) : exponent(exp) {}
 
   CUDA_LAMBDA Complex operator()(VectorProxy<Real, dim> q,
                                  const Complex& phi) const {
     Real mul = 1;
     for (UInt i = 0; i < dim; ++i)
       mul *= std::pow(q(i), exponent[i]);
 
     // Do not duplicate everything from hermitian symmetry
     if (std::abs(q.back()) < 1e-15)
       return mul * phi;
     return 2 * mul * phi;
   }
 
 private:
   std::array<UInt, dim> exponent;
 };
 }  // namespace
 
+/* -------------------------------------------------------------------------- */
 template <>
 std::vector<Real> Statistics<1>::computeMoments(Grid<Real, 1>& surface) {
   constexpr UInt dim = 1;
   std::vector<Real> moments(3);
   const auto psd = computePowerSpectrum(surface);
   auto wavevectors =
       FFTEngine::template computeFrequencies<Real, dim, true>(psd.sizes());
   // we don't multiply by 2 pi because moments are computed with k
   moments[0] = Loop::reduce<operation::plus>(moment_helper<dim>{{{0}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   moments[1] = Loop::reduce<operation::plus>(moment_helper<dim>{{{2}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   moments[2] = Loop::reduce<operation::plus>(moment_helper<dim>{{{4}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   return moments;
 }
 
 template <>
 std::vector<Real> Statistics<2>::computeMoments(Grid<Real, 2>& surface) {
   constexpr UInt dim = 2;
   std::vector<Real> moments(3);
   const auto psd = computePowerSpectrum(surface);
   auto wavevectors =
       FFTEngine::template computeFrequencies<Real, dim, true>(psd.sizes());
   // we don't multiply by 2 pi because moments are computed with k
   moments[0] = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 0}}},
                                              range<PVector>(wavevectors), psd)
                    .real();
   auto m02 = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 2}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   auto m20 = Loop::reduce<operation::plus>(moment_helper<dim>{{{2, 0}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   moments[1] = (m02 + m20) / 2;
   auto m22 = Loop::reduce<operation::plus>(moment_helper<dim>{{{2, 2}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   auto m40 = Loop::reduce<operation::plus>(moment_helper<dim>{{{4, 0}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   auto m04 = Loop::reduce<operation::plus>(moment_helper<dim>{{{0, 4}}},
                                            range<PVector>(wavevectors), psd)
                  .real();
   moments[2] = (3 * m22 + m40 + m04) / 3;
   return moments;
 }
 
+/* -------------------------------------------------------------------------- */
+template <UInt dim>
+Real Statistics<dim>::graphArea(const Grid<Real, dim>& zdisplacement) {
+  const auto h_size =
+      GridHermitian<Real, dim>::hermitianDimensions(zdisplacement.sizes());
+  GridHermitian<Real, dim> fourier_disp(h_size, 1),
+      fourier_gradient(h_size, dim);
+  Grid<Real, dim> gradient(zdisplacement.sizes(), dim);
+
+  // Compute gradient in Fourier domain
+  FFTEngine::makeEngine()->forward(zdisplacement, fourier_disp);
+
+  auto wavevectors =
+      FFTEngine::template computeFrequencies<Real, dim, true>(h_size);
+  wavevectors *= 2 * M_PI;
+
+  Loop::loop(
+      [] CUDA_LAMBDA(VectorProxy<const Real, dim> q,
+                     VectorProxy<Complex, dim> grad, const Complex& f) {
+        grad = q * f;
+        grad *= Complex{0, 1};
+      },
+      range<VectorProxy<const Real, dim>>(wavevectors),
+      range<VectorProxy<Complex, dim>>(fourier_gradient), fourier_disp);
+
+  FFTEngine::makeEngine()->backward(gradient, fourier_gradient);
+
+  // Integrate area of graph formula
+  Real factor = Loop::reduce<operation::plus>(
+      [] CUDA_LAMBDA(VectorProxy<const Real, dim> grad) {
+        return std::sqrt(1 + grad.l2squared());
+      },
+      range<VectorProxy<const Real, dim>>(gradient));
+  return factor / gradient.getGlobalNbPoints();
+  // return 2 * M_PI * factor / gradient.getGlobalNbPoints();
+}
+
 template struct Statistics<1>;
 template struct Statistics<2>;
 
 }  // namespace tamaas
diff --git a/src/core/statistics.hh b/src/core/statistics.hh
index 3cf584f..93afcba 100644
--- a/src/core/statistics.hh
+++ b/src/core/statistics.hh
@@ -1,55 +1,57 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program is distributed in the hope that it will be useful,
  *  but WITHOUT ANY WARRANTY; without even the implied warranty of
  *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  *  GNU Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef STATISTICS_HH
 #define STATISTICS_HH
 /* -------------------------------------------------------------------------- */
 #include "fft_engine.hh"
 #include "grid.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 /**
  * @brief Suitcase class for all statistics related functions
  */
 template <UInt dim>
 struct Statistics {
   /// Compute hrms
   static Real computeRMSHeights(Grid<Real, dim>& surface);
   /// Compute hrms' in fourier space
   static Real computeSpectralRMSSlope(Grid<Real, dim>& surface);
   /// Compute PSD of surface
   static GridHermitian<Real, dim>
   computePowerSpectrum(Grid<Real, dim>& surface);
   /// Compute autocorrelation
   static Grid<Real, dim> computeAutocorrelation(Grid<Real, dim>& surface);
   /// Compute spectral moments
   static std::vector<Real> computeMoments(Grid<Real, dim>& surface);
   /// Compute (corrected) contact area fraction
   static Real contact(const Grid<Real, dim>& tractions, UInt perimeter = 0);
+  /// Compute the area scaling factor of a periodic graph
+  static Real graphArea(const Grid<Real, dim>& displacement);
 
   using PVector = VectorProxy<Real, dim>;
 };
 
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
 #endif  // STATISTICS_HH
diff --git a/tests/test_surface.py b/tests/test_surface.py
index 6aec50a..368ad91 100644
--- a/tests/test_surface.py
+++ b/tests/test_surface.py
@@ -1,133 +1,146 @@
 # -*- coding: utf-8 -*-
 #
 # Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
 # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Affero General Public License as published
 # by the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # This program is distributed in the hope that it will be useful,
 # but WITHOUT ANY WARRANTY; without even the implied warranty of
 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 # GNU Affero General Public License for more details.
 #
 # You should have received a copy of the GNU Affero General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 from __future__ import print_function, division
 
 import numpy as np
 import tamaas as tm
 import pytest
 import scipy.integrate as integrate
-from scipy.special import jv
+from scipy.special import jv, ellipe as E
 from numpy import pi
 from scipy.interpolate import interp2d
 
 
 def spectrum(r, isopowerlaw):
     kl, kr, ks = isopowerlaw.q0, isopowerlaw.q1, isopowerlaw.q2
     H = isopowerlaw.hurst
     C = 1.
 
     if r > ks:
         return 0
     elif r > kr:
-        return C * (r/float(kr))**(-2.*(H+1))
+        return C * (r / float(kr))**(-2. * (H + 1))
     elif r > kl:
         return C
     else:
         return 0
 
 
 def integrate_bessel(x, spectrum, n):
-    return 2 * pi * integrate.quad(lambda y: y*jv(0, y*x) * spectrum(n*y),
-                                   0, np.inf, limit=200)[0]
+    return 2 * pi * integrate.quad(
+        lambda y: y * jv(0, y * x) * spectrum(n * y), 0, np.inf, limit=200)[0]
 
 
 def radial_average(field):
     x = np.linspace(-1, 1, field.shape[0])
     y = np.linspace(-1, 1, field.shape[1])
     interpolation = interp2d(x, y, field)
 
     n = field.shape[0] // 2
     radii = np.linspace(0, 1, n)
     thetas = np.linspace(0, 2 * pi, endpoint=False)
 
     sum = np.zeros_like(radii)
 
     for Θ in thetas:
         for i, r in enumerate(radii):
             sum[i] += interpolation(r * np.cos(Θ), r * np.sin(Θ))
     return sum / thetas.size
 
 
 integrate_bessel = np.vectorize(integrate_bessel, otypes=[tm.dtype])
 
 
-@pytest.fixture(scope="module",
-                params=[tm.SurfaceGeneratorFilter2D,
-                        tm.SurfaceGeneratorRandomPhase2D])
+@pytest.fixture(
+    scope="module",
+    params=[tm.SurfaceGeneratorFilter2D, tm.SurfaceGeneratorRandomPhase2D])
 def stats(request):
     iso = tm.Isopowerlaw2D()
 
     iso.q0 = 16
     iso.q1 = 32
     iso.q2 = 64
     iso.hurst = 0.8
 
     N = 243
 
     sg = request.param()
     sg.spectrum = iso
     sg.shape = [N, N]
 
     analytical_stats = {
         "rms_heights": iso.rmsHeights(),
         "rms_slopes_spectral": iso.rmsSlopes(),
         "moments": iso.moments(),
         "alpha": iso.alpha()
     }
 
     stats = {k: [] for k in analytical_stats}
 
     acf = np.zeros(sg.shape, dtype=tm.dtype)
 
     realizations = 1000
 
     for i in range(realizations):
         sg.random_seed = i
         s = sg.buildSurface()
-        stats['rms_heights'].append(
-            tm.Statistics2D.computeRMSHeights(s))
+        stats['rms_heights'].append(tm.Statistics2D.computeRMSHeights(s))
         stats['rms_slopes_spectral'].append(
             tm.Statistics2D.computeSpectralRMSSlope(s))
         stats['moments'].append(tm.Statistics2D.computeMoments(s))
         acf += tm.Statistics2D.computeAutocorrelation(s) / realizations
 
     L = 1 * sg.spectrum.q0
     n = sg.shape[0] // 2
-    x = np.linspace(0, L/2, n)
+    x = np.linspace(0, L / 2, n)
 
     true_acf = (L / (2 * pi))**2 \
         * integrate_bessel(x, lambda x: spectrum(x, sg.spectrum), L / (2 * pi))
 
     return analytical_stats, stats, acf, true_acf
 
 
 @pytest.mark.parametrize('key', ['rms_heights', 'rms_slopes_spectral'])
 def test_statistics(stats, key):
     analytical_stats, stats, _, _ = stats
     mu = np.mean(stats[key])
     error = np.abs(mu - analytical_stats[key]) / analytical_stats[key]
     assert error < 3e-3
 
 
 def test_autocorrelation(stats):
     _, _, acf, true_acf = stats
     acf = radial_average(np.fft.fftshift(acf))
-    acf_error = np.sum((acf-true_acf)**2) / np.sum(true_acf**2)
+    acf_error = np.sum((acf - true_acf)**2) / np.sum(true_acf**2)
 
     assert acf_error < 2e-3
+
+
+def test_graph_area():
+    N = 100
+    x = np.linspace(0, 1, N, endpoint=False)
+    f = np.sin(2 * np.pi * x)
+    df = 2 * np.pi * np.cos(2 * np.pi * x)
+
+    num_area = np.sqrt(1 + df**2).mean()
+    area = 2 * np.sqrt(1 + 4 * np.pi**2) / np.pi * E(4 * np.pi**2 /
+                                                     (1 + 4 * np.pi**2))
+    fourier_area = tm.Statistics1D.graphArea(f)
+    assert np.abs(num_area - fourier_area) / num_area < 1e-10
+    assert np.abs(area - fourier_area) / area < 2e-10