// OCCT 8 development API. Build with the Modeling Data and Foundation Classes modules.
#include <Geom_BezierCurve.hxx>
#include <Geom_OffsetCurve.hxx>
#include <Geom_OffsetSurface.hxx>
#include <Geom_SurfaceOfLinearExtrusion.hxx>
#include <Geom_SurfaceOfRevolution.hxx>
#include <GeomGridEval_SurfaceOfExtrusion.hxx>
#include <NCollection_Array1.hxx>
#include <NCollection_Array2.hxx>
#include <Standard_Failure.hxx>
#include <gp_Ax1.hxx>

#include <cassert>
#include <cmath>
#include <cstddef>
#include <iostream>

int main() {
  try {
    NCollection_Array1<gp_Pnt> aPoles(4);
    aPoles.ChangeAt(0) = gp_Pnt(1, 0, 0);
    aPoles.ChangeAt(1) = gp_Pnt(2, 0, 0.5);
    aPoles.ChangeAt(2) = gp_Pnt(1.5, 0, 1.5);
    aPoles.ChangeAt(3) = gp_Pnt(1, 0, 2);
    occ::handle<Geom_Curve> aBasis = new Geom_BezierCurve(aPoles);
    const gp_Dir aDirection(0, 1, 0);
    const double u = 0.4, v = 0.3, d = 0.15;

    Geom_OffsetCurve anOffsetCurve(aBasis, d, aDirection);
    const auto aCurveD1 = anOffsetCurve.EvalD1(u);
    std::cout << "Offset-curve speed: " << aCurveD1.D1.Magnitude() << '\n';

    const occ::handle<Geom_SurfaceOfLinearExtrusion> anExtrusion =
      new Geom_SurfaceOfLinearExtrusion(aBasis, aDirection);
    const auto [aPoint, aDu, aDv] = anExtrusion->EvalD1(u, v);
    assert(aDv.IsEqual(gp_Vec(aDirection), 1.0e-12, 1.0e-12));

    Geom_SurfaceOfRevolution aRevolution(
      aBasis, gp_Ax1(gp_Pnt(0, 0, 0), gp_Dir(0, 0, 1)));
    const auto aRevD2 = aRevolution.EvalD2(0.7, u); // radians, then basis parameter
    std::cout << "Revolution angular speed: " << aRevD2.D1U.Magnitude() << '\n';

    Geom_OffsetSurface anOffsetSurface(anExtrusion, d);
    const gp_Pnt anOffsetPoint = anOffsetSurface.EvalD0(u, v);
    const gp_Vec aNormal = aDu.Crossed(aDv).Normalized(); // regular sample by construction
    assert(anOffsetPoint.Distance(aPoint.Translated(d * aNormal)) < 1.0e-12);

    NCollection_Array1<double> aUParams(21), aVParams(17);
    for (size_t i = 0; i < aUParams.Size(); ++i) {
      aUParams.ChangeAt(i) = double(i) / double(aUParams.Size() - 1);
    }
    for (size_t j = 0; j < aVParams.Size(); ++j) {
      aVParams.ChangeAt(j) = -0.5 + double(j) / double(aVParams.Size() - 1);
    }
    GeomGridEval_SurfaceOfExtrusion anEvaluator(anExtrusion);
    const auto aGrid = anEvaluator.EvaluateGrid(aUParams, aVParams);
    assert(!aGrid.IsEmpty());
    for (size_t i = 0; i < aUParams.Size(); ++i) {
      for (size_t j = 0; j < aVParams.Size(); ++j) {
        const gp_Pnt aDirect = anExtrusion->EvalD0(aUParams.At(i), aVParams.At(j));
        assert(aGrid.At(i, j).Distance(aDirect) < 1.0e-12);
      }
    }
    std::cout << "Compared " << aUParams.Size() * aVParams.Size() << " grid points\n";
  } catch (const Standard_Failure& anError) {
    std::cerr << anError.what() << '\n';
    return 1;
  }
}