// FlatMap.cc
// Contact person: Freija Descamps <fbdescamps@lbl.gov>
// See FlatMap.hh for more details
//———————————————————————//

// Root stuff
#include <TVector3.h>
#include <TVector2.h>
#include <TGraph.h>
#include <TGraph2D.h>
#include <cmath>
#include <algorithm>
#include <RAT/FlatMap.hh>

namespace RAT {
  // From http://www.rwgrayprojects.com/rbfnotes/polyhed/PolyhedraData/Icosahedron/Icosahedron.pdf
  const double gm = ( 1.0 + sqrt( 5.0 ) ) / 2.0; // "Golden mean"
  const TVector3 V2 = TVector3( gm * gm, 0.0, gm * gm * gm ).Unit();
  const TVector3 V6 = TVector3( -gm * gm, 0.0, gm * gm * gm ).Unit(); // V6
  const TVector3 V12 = TVector3( 0.0, gm * gm * gm, gm * gm ).Unit(); // V12
  const TVector3 V17 = TVector3( 0.0, -gm * gm * gm, gm * gm ).Unit(); // V17
  const TVector3 V27 = TVector3( gm * gm * gm, gm * gm, 0.0 ).Unit(); // V27
  const TVector3 V31 = TVector3( -gm * gm * gm, gm * gm, 0.0 ).Unit(); // V31
  const TVector3 V33 = TVector3( -gm * gm * gm, -gm * gm, 0.0 ).Unit(); // V33
  const TVector3 V37 = TVector3( gm * gm * gm, -gm * gm, 0.0 ).Unit(); // V37
  const TVector3 V46 = TVector3( 0.0, gm * gm * gm, -gm * gm ).Unit(); // V46
  const TVector3 V51 = TVector3( 0.0, -gm * gm * gm, -gm * gm ).Unit(); // V51
  const TVector3 V54 = TVector3( gm * gm, 0.0, -gm * gm * gm ).Unit(); // V54
  const TVector3 V58 = TVector3( -gm * gm, 0.0, -gm * gm * gm ).Unit(); // V58

  TVector2 TransformCoord(
                          const TVector3& Vec1,
                          const TVector3& Vec2,
                          const TVector3& Vec3,
                          const TVector2& A1,
                          const TVector2& A2,
                          const TVector2& A3,
                          const TVector3& P )
  {
    /// Transform the co-ord P relative to V1,2,3 to the
    /// returned value relative to A1,2,3
    TVector3 xV = Vec2 - Vec1;
    TVector3 yV = ( ( Vec3 - Vec1 ) + ( Vec3 - Vec2 ) ) * 0.5;
    TVector3 zV = xV.Cross( yV ).Unit();
    double planeD = Vec1.Dot( zV );
    double t = planeD / P.Dot( zV );
    TVector3 localP = t*P - Vec1;
    TVector2 xA = A2 - A1;
    TVector2 yA = ( ( A3 - A1 ) +( A3 - A2 ) )  * 0.5;
    double convUnits = xA.Mod() / xV.Mag();
    TVector2 result;
    result = localP.Dot( xV.Unit() ) * xA.Unit() * convUnits;
    result += localP.Dot( yV.Unit() ) * yA.Unit() * convUnits + A1;
    return result;
  }

  void SphereToIcosahedron(
                           TVector3& pointOnSphere, // The spherical Position
                           TVector2& resultPosition, // The returned position
                           const double totalLength, // Total size of the 2d grid (default to 1)
                           const double rotation) // Rotation: 0 for "Phil" style, 2.12 for "SNO" style
  {
    pointOnSphere.RotateX(rotation);
    TGraph* gA = new TGraph();
    TGraph2D* gV = new TGraph2D();
    int gAC = 0, gVC = 0;
    gV->SetPoint( gVC, V2.X(), V2.Y(), V2.Z() ); gVC++;
    gV->SetPoint( gVC, V6.X(), V6.Y(), V6.Z() ); gVC++;
    gV->SetPoint( gVC, V12.X(), V12.Y(), V12.Z() ); gVC++;
    gV->SetPoint( gVC, V17.X(), V17.Y(), V17.Z() ); gVC++;
    gV->SetPoint( gVC, V27.X(), V27.Y(), V27.Z() ); gVC++;
    gV->SetPoint( gVC, V31.X(), V31.Y(), V31.Z() ); gVC++;
    gV->SetPoint( gVC, V33.X(), V33.Y(), V33.Z() ); gVC++;
    gV->SetPoint( gVC, V37.X(), V37.Y(), V37.Z() ); gVC++;
    gV->SetPoint( gVC, V46.X(), V46.Y(), V46.Z() ); gVC++;
    gV->SetPoint( gVC, V51.X(), V51.Y(), V51.Z() ); gVC++;
    gV->SetPoint( gVC, V54.X(), V54.Y(), V54.Z() ); gVC++;
    gV->SetPoint( gVC, V58.X(), V58.Y(), V58.Z() ); gVC++;
    // Faces {{ 2, 6, 17}, { 2, 12, 6}, { 2, 17, 37}, { 2, 37, 27}, { 2, 27, 12}, {37, 54, 27},
    // {27, 54, 46}, {27, 46, 12}, {12, 46, 31}, {12, 31, 6}, { 6, 31, 33}, { 6, 33, 17},
    // {17, 33, 51}, {17, 51, 37}, {37, 51, 54}, {58, 54, 51}, {58, 46, 54}, {58, 31, 46},
    // {58, 33, 31}, {58, 51, 33}}
    std::vector<TVector3> IcosahedronCentres;
    IcosahedronCentres.push_back( ( V2 + V6 + V17 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V2 + V12 + V6 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V2 + V17 + V37 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V2 + V37 + V27 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V2 + V27 + V12 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V37 + V54 + V27 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V27 + V54 + V46 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V27 + V46 + V12 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V12 + V46 + V31 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V12 + V31 + V6 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V6 + V31 + V33 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V6 + V33 + V17 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V17 + V33 + V51 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V17 + V51 + V37 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V37 + V51 + V54 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V58 + V54 + V51 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V58 + V46 + V54 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V58 + V31 + V46 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V58 + V33 + V31 ) * ( 1.0 / 3.0 ) );
    IcosahedronCentres.push_back( ( V58 + V51 + V33 ) * ( 1.0 / 3.0 ) );
    std::vector<double> distFromCentre;
    unsigned int uLoop;
    //printf("uLoop is %i\n",uLoop);
    for( uLoop = 0; uLoop < IcosahedronCentres.size(); uLoop++ )
      distFromCentre.push_back( ( IcosahedronCentres[uLoop] - pointOnSphere ).Mag() );
    const int face = min_element( distFromCentre.begin(), distFromCentre.end() ) - distFromCentre.begin() + 1;
    const double a = totalLength / 5.5;
    const double b = a * sqrt( 3.0 ) / 2.0;
    const TVector2 A12a = TVector2( a / 2.0, 0.0 );
    gA->SetPoint( gAC, A12a.X(), A12a.Y() ); gAC++;
    const TVector2 A12b = TVector2( 3.0 * a / 2.0, 0.0 );
    gA->SetPoint( gAC, A12b.X(), A12b.Y() ); gAC++;
    const TVector2 A12c = TVector2( 5.0 * a / 2.0, 0.0 );
    gA->SetPoint( gAC, A12c.X(), A12c.Y() ); gAC++;
    const TVector2 A12d = TVector2( 7.0 *a / 2.0, 0.0 );
    gA->SetPoint( gAC, A12d.X(), A12d.Y() ); gAC++;
    const TVector2 A12e = TVector2( 9.0 * a / 2.0, 0.0 );
    gA->SetPoint( gAC, A12e.X(), A12e.Y() ); gAC++;
    const TVector2 A2a = TVector2( 0.0, b );
    gA->SetPoint( gAC, A2a.X(), A2a.Y() ); gAC++;
    const TVector2 A2b = TVector2( 5.0 * a, b );
    gA->SetPoint( gAC, A2b.X(), A2b.Y() ); gAC++;
    const TVector2 A17a = TVector2( a / 2.0 , 2.0 * b );
    gA->SetPoint( gAC, A17a.X(), A17a.Y() ); gAC++;
    const TVector2 A17b = TVector2( 11.0 * a / 2.0 , 2.0 * b );
    gA->SetPoint( gAC, A17b.X(), A17b.Y() ); gAC++;
    const TVector2 A51a = TVector2( a, 3.0 * b );
    gA->SetPoint( gAC, A51a.X(), A51a.Y() ); gAC++;
    const TVector2 A51b = TVector2( 2.0 * a, 3.0 * b );
    gA->SetPoint( gAC, A51b.X(), A51b.Y() ); gAC++;
    const TVector2 A51c = TVector2( 3.0 * a, 3.0 * b );
    gA->SetPoint( gAC, A51c.X(), A51c.Y() ); gAC++;
    const TVector2 A51d = TVector2( 4.0 * a, 3.0 * b );
    gA->SetPoint( gAC, A51d.X(), A51d.Y() ); gAC++;
    const TVector2 A51e = TVector2( 5.0 * a, 3.0 * b );
    gA->SetPoint( gAC, A51e.X(), A51e.Y() ); gAC++;
    const TVector2 A27 = TVector2( 4.0 * a, b );
    gA->SetPoint( gAC, A27.X(), A27.Y() ); gAC++;
    const TVector2 A46 = TVector2( 3.0 * a, b );
    gA->SetPoint( gAC, A46.X(), A46.Y() ); gAC++;
    const TVector2 A31 = TVector2( 2.0 * a, b );
    gA->SetPoint( gAC, A31.X(), A31.Y() ); gAC++;
    const TVector2 A6 = TVector2( a, b );
    gA->SetPoint( gAC, A6.X(), A6.Y() ); gAC++;
    const TVector2 A37 = TVector2( 9.0 * a / 2.0 , 2.0 * b );
    gA->SetPoint( gAC, A37.X(), A37.Y() ); gAC++;
    const TVector2 A33 = TVector2( 3.0 * a / 2.0 , 2.0 * b );
    gA->SetPoint( gAC, A33.X(), A33.Y() ); gAC++;
    const TVector2 A58 = TVector2( 5.0 * a / 2.0 , 2.0 * b );
    gA->SetPoint( gAC, A58.X(), A58.Y() ); gAC++;
    const TVector2 A54 = TVector2( 7.0 * a / 2.0 , 2.0 * b );
    gA->SetPoint( gAC, A54.X(), A54.Y() ); gAC++;

    TVector3 xV, yV;
    TVector2 xA, yA;
    switch(face)
      {
      case 1://{ 2, 6, 17}
        resultPosition = TransformCoord( V2, V6, V17, A2a, A6, A17a, pointOnSphere );
        break;
      case 2://{ 2, 12, 6}
        resultPosition = TransformCoord( V2, V12, V6, A2a, A12a, A6, pointOnSphere );
        break;
      case 3://{ 2, 17, 37}
        resultPosition = TransformCoord( V2, V17, V37, A2b, A17b, A37, pointOnSphere );
        break;
      case 4://{ 2, 37, 27}
        resultPosition = TransformCoord( V2, V37, V27, A2b, A37, A27, pointOnSphere );
        break;
      case 5://{ 2, 27, 12}
        resultPosition = TransformCoord( V2, V27, V12, A2b, A27, A12e, pointOnSphere );
        break;
      case 6://{37, 54, 27}
        resultPosition = TransformCoord( V37, V54, V27, A37, A54, A27, pointOnSphere );
        break;
      case 7://{27, 54, 46}
        resultPosition = TransformCoord( V27, V54, V46, A27, A54, A46, pointOnSphere );
        break;
      case 8://{27, 46, 12}
        resultPosition = TransformCoord( V27, V46, V12, A27, A46, A12d, pointOnSphere );
        break;
      case 9://{12, 46, 31}
        resultPosition = TransformCoord( V12, V46, V31, A12c, A46, A31, pointOnSphere );
        break;
      case 10://{12, 31, 6}
        resultPosition = TransformCoord( V12, V31, V6, A12b, A31, A6, pointOnSphere );
        break;
      case 11://{ 6, 31, 33}
        resultPosition = TransformCoord( V6, V31, V33, A6, A31, A33, pointOnSphere );
        break;
      case 12://{ 6, 33, 17}
        resultPosition = TransformCoord( V6, V33, V17, A6, A33, A17a, pointOnSphere );
        break;
      case 13://{17, 33, 51}
        resultPosition = TransformCoord( V17, V33, V51, A17a, A33, A51a, pointOnSphere );
        break;
      case 14://{17, 51, 37}
        resultPosition = TransformCoord( V17, V51, V37, A17b, A51e, A37, pointOnSphere );
        break;
      case 15://{37, 51, 54}
        resultPosition = TransformCoord( V37, V51, V54, A37, A51d, A54, pointOnSphere );
        break;
      case 16://{58, 54, 51}
        resultPosition = TransformCoord( V58, V54, V51, A58, A54, A51c, pointOnSphere );
        break;
      case 17://{58, 46, 54}
        resultPosition = TransformCoord( V58, V46, V54, A58, A46, A54, pointOnSphere );
        break;
      case 18://{58, 31, 46}
        resultPosition = TransformCoord( V58, V31, V46, A58, A31, A46, pointOnSphere );
        break;
      case 19://{58, 33, 31}
        resultPosition = TransformCoord( V58, V33, V31, A58, A33, A31, pointOnSphere );
        break;
      case 20://{58, 51, 33}
        resultPosition = TransformCoord( V58, V51, V33, A58, A51b, A33, pointOnSphere );
        break;
      }
    return;
  }
} // namespace RAT