import pyvista as pv
import numpy as np
from matplotlib.pyplot import cm
import matplotlib.pyplot as plt

ax_index = 0
ax_label = {0: "x", 1: "y", 2: "z"}[ax_index]
n = 20
l = []

output = pv.MultiBlock()
mesh = pv.Cone()
plotter = pv.Plotter()

rng = np.linspace(mesh.bounds[ax_index * 2] + 1e-8, mesh.bounds[ax_index * 2 + 1], n)
center = mesh.center
col = cm.rainbow(np.linspace(0, 1, n))

for i in range(n):
    center[ax_index] = rng[i]
    clipped, mesh = mesh.clip(normal=ax_label, origin=center, return_clipped=True)

    output.append(clipped)
    l.append(clipped.integrate_data()["Area"])
    plotter.add_mesh(clipped, color=col[i])
    
plotter.show_axes()
plotter.show()
plt.plot(np.cumsum(l))
plt.xlabel(ax_label)
plt.ylabel("Cumulative area integral")
plt.show()