#include <G4Polyhedra.hh>
#include <G4RotationMatrix.hh>
#include <G4DisplacedSolid.hh>
#include <G4Sphere.hh>
#include <G4IntersectionSolid.hh>

#include <RAT/SNOMANBellyPlate.hh>
#include <RAT/Log.hh>
using namespace RAT;
using namespace RAT::geo;

#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Units/PhysicalConstants.h>
#include <string>
using namespace CLHEP;
using namespace std;

G4VSolid*
SNOMANBellyPlate::Construct( const std::string& name,
                       DBLinkPtr table ) const
{
  const double octahedronBase = table->GetD( "octahedron_base" ) * mm;
  const double octahedronHalfHeight = table->GetD( "octahedron_half_height" ) * mm;
  const double thickness = table->GetD( "thickness" ) * mm;
  const double rSphere = table->GetD( "r_sphere" ) * mm;
  const double arc = table->GetD( "arc" ) * deg;
  const double rOctahedron = table->GetD( "r_octahedron" ) * mm;

  return ConstructBellyPlate( name, octahedronBase, octahedronHalfHeight, thickness, rSphere, arc, rOctahedron );
}

G4VSolid*
SNOMANBellyPlate::ConstructBellyPlate( const std::string& name,
                                 const double octahedronBase,
                                 const double octahedronHalfHeight,
                                 const double thickness,
                                 const double rSphere,
                                 const double arc,
                                 const double rOctahedron ) const
{
  // Logic:
  // The belly plate is a spherical tile with two beveled extrusions on the inner and outer sides.
  // These extrusions are truncated pyramids, but note the truncation is also concentric with the curvature of the tile.
  // To build this S. Brice used a octahedron truncated by a spherical tile, this approach is followed here.
  // The Belly plate is built at 0, 0, 0 using the rSphere and rOctahedron to position it.

  // First build the octahedron
  double octahedronRadii[] = { 0.0, octahedronBase / cos( pi / 4.0 ), 0.0 }; // Polyhedra radii is to corner not side (definition is side)
  double octahedronZ[] = { octahedronHalfHeight, 0.0, -octahedronHalfHeight };
  G4Polyhedra* octahedronUnRotated = new G4Polyhedra( name + "_octahedron_unrotated_solid", 0.0, twopi, 4, 3, octahedronRadii, octahedronZ );
  // Rotate the octahedron such that point to point is along the x axis
  G4RotationMatrix* rotation = new G4RotationMatrix();
  rotation->rotateY( pi / 2.0 );
  G4DisplacedSolid* octahedronRotated = new G4DisplacedSolid( name + "_octahedron_rotated_solid", octahedronUnRotated, rotation, G4ThreeVector( 0.0, 0.0, 0.0 ) );
  G4RotationMatrix* rotation2 = new G4RotationMatrix();
  rotation2->rotateX( pi / 4.0 );
  // Finish rotation
  G4DisplacedSolid* octahedron = new G4DisplacedSolid( name + "_octahedron_solid", octahedronRotated, rotation2, G4ThreeVector( 0.0, 0.0, 0.0 ) );
  // Now build the spherical tile
  G4Sphere* sphereTile = new G4Sphere( name + "_spherical_tile_undisplaced_solid",
                                       rSphere - thickness / 2.0,
                                       rSphere + thickness / 2.0,
                                       -arc,
                                       2.0 * arc,
                                       pi / 2.0 - arc,
                                       2.0 * arc );
  // Now intersect with the octahedron
  G4VSolid* intersectedSolid = new G4IntersectionSolid( name + "_intersection_solid", sphereTile, octahedron, NULL, G4ThreeVector( rOctahedron, 0.0, 0.0 ) );
  // Finally move the net solid back to the origin
  return new G4DisplacedSolid( name + "_belly_plate", intersectedSolid, NULL, G4ThreeVector( -rSphere, 0.0, 0.0 ) );
}