#include <TopoDS_Shape.hxx>
#include <TopoDS.hxx>
#include <Geom_BSplineSurface.hxx>
#include <TopoDS_Face.hxx>
#include <TopoDS_Shell.hxx>
#include <TopoDS_Solid.hxx>
#include <BRepBuilderAPI_MakeFace.hxx>
#include <BRepBuilderAPI_MakeShell.hxx>
#include <BRepBuilderAPI_Sewing.hxx>
#include <BRepBuilderAPI_MakeSolid.hxx>
#include <TopExp_Explorer.hxx>
#include <STEPControl_Writer.hxx>
#include <STEPControl_StepModelType.hxx>
#include <gp_Trsf.hxx>
#include <BRepBuilderAPI_Transform.hxx>
#include <BRepAlgoAPI_Cut.hxx>
#include <BRepAlgoAPI_Fuse.hxx>
#include <GC_MakeTrimmedCylinder.hxx>
#include <gce_MakeCylinder.hxx>
#include <BRepPrimAPI_MakeCylinder.hxx>
#include <BRepPrimAPI_MakeSphere.hxx>

#define NCYLS 4

main()
{
  gp_Pnt p[NCYLS];
  int i;
  for (i = 0; i < NCYLS; i++) {
    double angle = 2.0*M_PI*(i/double(NCYLS));
    p[i].SetX(cos(angle));
    p[i].SetY(sin(angle));
    p[i].SetZ(0.0);
  }
  TopoDS_Shape total;
  STEPControl_Writer writer2;
  for (i = 0; i < NCYLS; i++) {
    std::cerr << "Cylinder " << i << "\n";
    gp_Vec vec(p[i], p[(i+1)%NCYLS]);
    gp_Dir dir(vec);
    gp_Dir up(0.0, 0.0, 1.0);
    gp_Ax1 ax1(p[i], dir);
    gp_Ax2 ax2(p[i], dir, up);
    Standard_Real ht = vec.Magnitude();
    BRepPrimAPI_MakeCylinder primcyl(ax2, 0.1, ht);
    TopoDS_Shell shellcyl = primcyl.Shell();
    BRepBuilderAPI_MakeSolid solidcylbuild(shellcyl);
    TopoDS_Solid solidcyl = TopoDS::Solid(solidcylbuild);
    writer2.Transfer(solidcyl, STEPControl_AsIs);
    if (i == 0) {
      total = solidcyl;
    }
    else {
      BRepAlgoAPI_Fuse fuse(total, solidcyl);
      std::cerr << "Fuse error status " << fuse.ErrorStatus() << "\n";
      total = fuse.Shape();
    }

    std::cerr << "Sphere " << i << "\n";
    BRepPrimAPI_MakeSphere primsph(p[i], 0.1);
    TopoDS_Shell shellsph = primsph.Shell();
    BRepBuilderAPI_MakeSolid solidsphbuild(shellsph);
    TopoDS_Solid solidsph = TopoDS::Solid(solidsphbuild);
    writer2.Transfer(solidsph, STEPControl_AsIs);
    BRepAlgoAPI_Fuse fuse(total, solidsph);
    std::cerr << "Fuse error status " << fuse.ErrorStatus() << "\n";
    total = fuse.Shape();

  }
  STEPControl_Writer writer;
  writer.Transfer(total, STEPControl_AsIs);
  writer.Write("STEP_fuse.stp");
  writer2.Write("STEP_fuse2.stp");
}