#include "TROOT.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TView.h"
#include "TPolyLine3D.h"
#include "TLatex.h"

void DrawFrame3D(){
  Int_t k1[3][3] = {
    {0,1,2},
    {2,0,1},
    {1,2,0}};
  Int_t k2[3][3] = {
    {1,0,0},
    {0,1,0},
    {0,0,1}};
  TPolyLine3D *l1=new TPolyLine3D();
  for (Int_t j=0;j<3;j++){
    Double_t thetamin = 0;
    Double_t thetamax = acos(-1.)/2.;
    Int_t    n = 30;
    Double_t delta = (thetamax-thetamin)/n;
    for (Int_t i=0;i<=n;i++) {
      Double_t theta = i * delta;
      Double_t x[3];
      x[0] = cos(theta);
      x[1] = sin(theta);
      x[2] = 0.;
      l1->SetNextPoint(x[k1[j][0]],x[k1[j][1]],x[k1[j][2]]);
    }
    l1->Draw();
  }
  Double_t xyzsize = 1.2;
  TPolyLine3D *l2[3];
  for (Int_t j=0;j<3;j++){
    l2[j] = new TPolyLine3D();
    l2[j]->SetNextPoint(0,0,0);
    l2[j]->SetNextPoint(k2[j][0]*xyzsize,k2[j][1]*xyzsize,k2[j][2]*xyzsize);
    l2[j]->Draw();
  }
  Double_t arrowsize = 0.1;
  TPolyLine3D *arrow[3];
  for (Int_t j=0;j<3;j++){
    arrow[j] = new TPolyLine3D();
    arrow[j]->SetNextPoint(xyzsize*(k2[j][0]*(1+arrowsize)),
			   xyzsize*(k2[j][1]*(1+arrowsize)),
			   xyzsize*(k2[j][2]*(1+arrowsize)));
    arrow[j]->SetNextPoint(xyzsize*(k2[j][0]+0.3*k2[j][1]*arrowsize),
			   xyzsize*(k2[j][1]+0.3*k2[j][2]*arrowsize),
			   xyzsize*(k2[j][2]+0.3*k2[j][0]*arrowsize));
    arrow[j]->SetNextPoint(xyzsize*(k2[j][0]-0.3*k2[j][1]*arrowsize),
			   xyzsize*(k2[j][1]-0.3*k2[j][2]*arrowsize),
			   xyzsize*(k2[j][2]-0.3*k2[j][0]*arrowsize));
    arrow[j]->SetNextPoint(xyzsize*(k2[j][0]*(1+arrowsize)),
			   xyzsize*(k2[j][1]*(1+arrowsize)),
			   xyzsize*(k2[j][2]*(1+arrowsize)));
    arrow[j]->Draw();
  }
}

void rec_solid_angle_def(){
  gROOT->SetStyle("Plain");
  gStyle->SetTextFont(132);
  TCanvas *c1 = new TCanvas("c1","c1",400,400);
  c1->SetTheta(20.);
  c1->SetPhi(-121.);
  TView *view = TView::CreateView(1);
  Double_t view_size = 0.6;
  view->SetRange(-view_size*0.5,-view_size*0.3,-view_size*0.0,
		  view_size*1.5, view_size*1.7, view_size*2.0);
  DrawFrame3D();

  const Int_t nPolyLine3D = 20;
  TPolyLine3D *guide[nPolyLine3D];
  for (Int_t i = 0; i<nPolyLine3D; i++){
    guide[i]=new TPolyLine3D();
  }

  Double_t pi = acos(-1.);
  Int_t    n = 30;
  Double_t r,theta,phi;
  Double_t x[3];
  Double_t theta0 = pi/6.;
  Double_t phi0   = pi/6.;

  r=1.;
  for (Int_t i=0;i<=n;i++) {
    theta = i*(pi/2)/n;
    x[0] = r*cos(theta);
    x[1] = r*sin(theta)*cos(theta0);
    x[2] = r*sin(theta)*sin(theta0);
    guide[0]->SetNextPoint(x[2],x[0],x[1]);
    x[1] = r*sin(theta)*sin(phi0);
    x[2] = r*sin(theta)*cos(phi0);
    guide[1]->SetNextPoint(x[0],x[1],x[2]);
  }
  guide[0]->SetNextPoint(0,0,0);
  guide[1]->SetNextPoint(0,0,0);

  r = 0.15;
  for (Int_t i=0;i<=n;i++) {
    theta = i*theta0/n;
    phi   = i*phi0/n;
    x[0] = r*sin(theta);
    x[1] = 0.;
    x[2] = r*cos(theta);
    guide[2]->SetNextPoint(x[0],x[1],x[2]);
    x[0] = 0.;
    x[1] = r*sin(phi);
    x[2] = r*cos(phi);
    guide[3]->SetNextPoint(x[0],x[1],x[2]);
  }

  x[0]=tan(theta0)/sqrt(1+pow(tan(theta0),2)+pow(tan(phi0),2));
  x[1]=tan(phi0)/sqrt(1+pow(tan(theta0),2)+pow(tan(phi0),2));
  x[2]=1/sqrt(1+pow(tan(theta0),2)+pow(tan(phi0),2));
  guide[4]->SetNextPoint(x[0],x[1],0);
  guide[4]->SetNextPoint(0,x[1],0);
  guide[4]->SetNextPoint(0,x[1],x[2]);
  guide[4]->SetNextPoint(0,0,x[2]);
  guide[4]->SetNextPoint(x[0],0,x[2]);
  guide[4]->SetNextPoint(x[0],0,0);
  guide[4]->SetNextPoint(x[0],x[1],0);
  guide[4]->SetNextPoint(x[0],x[1],x[2]);
  guide[4]->SetLineStyle(2);
  guide[5]->SetNextPoint(0,x[1],x[2]);
  guide[5]->SetNextPoint(x[0],x[1],x[2]);
  guide[5]->SetNextPoint(x[0],0,x[2]);
  guide[5]->SetLineStyle(2);

  for (Int_t i = 0; i<nPolyLine3D; i++){
    guide[i]->Draw();
  }

  Double_t xx[3],xxt[3];
  TLatex Tl;
  xx[0]=0.00; xx[1]=0.03; xx[2]=0.18; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"#it{#varphi}");
  xx[0]=0.13; xx[1]=0.00; xx[2]=0.20; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"#it{#theta}");
  xx[0]=0.00; xx[1]=-0.02; xx[2]=-0.10; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"O");
  xx[0]=x[0]-0.03; xx[1]=x[1]-0.05; xx[2]=x[2]+0.01; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"P");
  xx[0]=0; xx[1]=0.01; xx[2]=x[2]; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"B");
  xx[0]=0; xx[1]=0.01; xx[2]=1.01; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"B'");
  xx[0]=x[0]+0.13; xx[1]=0; xx[2]=x[2]-0.03; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"A");
  xx[0]=sin(theta)+0.12; xx[1]=0; xx[2]=cos(theta)+0.05; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"A'");
  xx[0]=0; xx[1]=x[1]+0.01; xx[2]=x[2]-0.05; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"C");
  xx[0]=0; xx[1]=sin(phi); xx[2]=cos(phi); view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"C'");
  xx[0]=1.3; xx[1]=0; xx[2]=0.1; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"#it{x}");
  xx[0]=0; xx[1]=1.25; xx[2]=0.1; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"#it{y}");
  xx[0]=0; xx[1]=-0.10; xx[2]=1.22; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"#it{z}");
}