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

void tri_solid_angle(){
  gROOT->SetStyle("Plain");
  gStyle->SetTextFont(132);

  Double_t cwidth  = 400;
  Double_t cheight = 400;
  TCanvas* c1 = new TCanvas("c1", "c1", cwidth, cheight);
  c1->SetWindowSize(cwidth + (cwidth - c1->GetWw()), cheight + (cheight - c1->GetWh()));

  c1->SetTheta(22.);
  c1->SetPhi(130.);
  TView *view = TView::CreateView(1);
  Double_t view_size = 1.5;
  view->SetRange(-view_size*1.0,-view_size*1.0,-view_size*1.0,
		  view_size*1.0, view_size*1.0, view_size*1.0);

  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 = 60;
  Double_t theta01 = pi/8.0;
  Double_t theta02 = -pi/5.0;
  Double_t phi0    = pi/4.;
  Double_t psi01   = asin(sqrt(1-pow(cos(theta01)*cos(phi0),2)));
  Double_t psi02   = -asin(sqrt(1-pow(cos(theta02)*cos(phi0),2)));
  
  Double_t r0  =  2.3;
  Double_t x01 =  r0*sin(theta01);
  Double_t x02 =  r0*sin(theta02);
  Double_t y0  =  r0*sin(phi0);

  Double_t gtheta = -c1->GetTheta()*pi/180.;
  Double_t gphi   = -c1->GetPhi()*pi/180.;  
  for (Int_t i=0;i<=n;i++) {
    Double_t r = r0;
    Double_t theta = i*2.*pi/n;
    Double_t x = r*cos(theta);
    Double_t z = r*sin(theta);
    guide[1]->SetNextPoint(x*cos(gphi)+z*sin(gtheta)*sin(gphi),
			   x*sin(gphi)-z*sin(gtheta)*cos(gphi),
			   z*cos(gtheta));
  }

  /*
    [ cos(t)  sin(t) 0] [1   0      0     ] [x]
    [-sin(t)  cos(t) 0] [0   cos(p) sin(p)] [y]
    [0        0      1] [0  -sin(p) cos(p)] [0]

    [ cos(t) sin(t) 0] [x]
    [-sin(t) cos(t) 0] [cos(p)*y]
    [ 0      0      1] [-sin(p)*y]

    [ cos(t)*x +sin(t)*cos(p)*y]
    [-sin(t)*x +cos(t)*cos(p)*y]
    [ sin(p)*y]
   */

  for (Int_t i=0;i<=n;i++) {
    Double_t r = 0.5;
    Double_t theta = -i*(theta01-theta02)/n+pi-theta02;
    guide[2]->SetNextPoint(r*cos(theta),r*sin(theta),0);
  }

  for (Int_t i=0;i<=n;i++) {
    Double_t r1 = 0.5;
    Double_t r2 = 0.60;
    Double_t psi1 = i*psi01/n+pi;
    Double_t psi2 = i*psi02/n+pi;
    Double_t phi1 = -asin(sin(phi0)/
			  sqrt(1-pow(cos(phi0)*cos(theta01),2)));
    Double_t phi2 = asin(sin(phi0)/
			  sqrt(1-pow(cos(phi0)*cos(theta02),2)));
    Double_t x1 = r1*cos(psi1);
    Double_t y1 = r1*sin(psi1);
    Double_t x2 = r2*cos(psi2);
    Double_t y2 = r2*sin(psi2);
    guide[3]->SetNextPoint( x1*cos(theta01)+y1*sin(theta01)*cos(phi1),
			    -x1*sin(theta01)+y1*cos(theta01)*cos(phi1),
			     y1*sin(phi1));
    guide[4]->SetNextPoint( x2*cos(theta02)+y2*sin(theta02)*cos(phi2),
			    -x2*sin(theta02)+y2*cos(theta02)*cos(phi2),
			     y2*sin(phi2));
  }

  for (Int_t i=0;i<=n;i++) {
    Double_t r = r0;
    Double_t theta = -i*(theta01-theta02)/n+pi-theta02;
    guide[7]->SetNextPoint(r*cos(theta),r*sin(theta),0);
  }

  for (Int_t i=0;i<=n;i++) {
    Double_t r = r0;
    Double_t psi1 = i*psi01/n+pi;
    Double_t psi2 = i*psi02/n+pi;
    Double_t phi1 = -asin(sin(phi0)/
			  sqrt(1-pow(cos(phi0)*cos(theta01),2)));
    Double_t phi2 = asin(sin(phi0)/
			  sqrt(1-pow(cos(phi0)*cos(theta02),2)));
    Double_t x1 = r*cos(psi1);
    Double_t y1 = r*sin(psi1);
    Double_t x2 = r*cos(psi2);
    Double_t y2 = r*sin(psi2);
    guide[8]->SetNextPoint( x1*cos(theta01)+y1*sin(theta01)*cos(phi1),
			    -x1*sin(theta01)+y1*cos(theta01)*cos(phi1),
			     y1*sin(phi1));
    guide[9]->SetNextPoint( x2*cos(theta02)+y2*sin(theta02)*cos(phi2),
			    -x2*sin(theta02)+y2*cos(theta02)*cos(phi2),
			     y2*sin(phi2));
  }

  guide[10]->SetNextPoint(-1.5,1.5*tan(theta01),0);
  guide[10]->SetNextPoint(-1.5,0,1.5*tan(phi0));
  guide[10]->SetNextPoint(-1.5,1.5*tan(theta02),0);
  guide[10]->SetNextPoint(-1.5,1.5*tan(theta01),0);

  guide[7]->SetLineWidth(5);
  guide[8]->SetLineWidth(5);
  guide[9]->SetLineWidth(5);
  guide[10]->SetLineWidth(5);

  for (Int_t i=0;i<=n;i++) {
    Double_t r = r0;
    Double_t theta = i*2.*pi/n;
    guide[11]->SetNextPoint(r*cos(theta),r*sin(theta),0);

    Double_t phi1   = -asin(sin(phi0)/
			    sqrt(1-pow(cos(phi0)*cos(theta01),2)));
    Double_t phi2   = asin(sin(phi0)/
			    sqrt(1-pow(cos(phi0)*cos(theta02),2)));
    Double_t x = r*cos(theta);
    Double_t y = r*sin(theta);
    guide[13]->SetNextPoint( x*cos(theta01)+y*sin(theta01)*cos(phi1),
			    -x*sin(theta01)+y*cos(theta01)*cos(phi1),
			     y*sin(phi1));
    guide[14]->SetNextPoint( x*cos(theta02)+y*sin(theta02)*cos(phi2),
			    -x*sin(theta02)+y*cos(theta02)*cos(phi2),
			     y*sin(phi2));
  }
  guide[11]->SetLineStyle(2);
  guide[12]->SetLineStyle(2);
  guide[13]->SetLineStyle(2);
  guide[14]->SetLineStyle(2);

  guide[15]->SetNextPoint(-r0*cos(theta01),x01,0);
  guide[15]->SetNextPoint(0,0,0);
  guide[15]->SetNextPoint(-r0*cos(theta02),x02,0);
  guide[16]->SetNextPoint(-r0*cos(phi0),0,y0);
  guide[16]->SetNextPoint(0,0,0);


//  Double_t virlen = 0.2;
//  guide[13]->SetNextPoint(-r0,virlen,0);
//  guide[13]->SetNextPoint(-r0,virlen,virlen);
//  guide[13]->SetNextPoint(-r0,0,virlen);
//  guide[14]->SetNextPoint(r0,virlen,0);
//  guide[14]->SetNextPoint(r0,virlen,virlen);
//  guide[14]->SetNextPoint(r0,0,virlen);
//
//  guide[15]->SetNextPoint(-1.4,virlen,0);
//  guide[15]->SetNextPoint(-1.4,virlen,virlen);
//  guide[15]->SetNextPoint(-1.4,0,virlen);
//  guide[16]->SetNextPoint(-1.4,virlen,0);
//  guide[16]->SetNextPoint(-1.4+virlen,virlen,0);
//  guide[16]->SetNextPoint(-1.4+virlen,0,0);
//  guide[17]->SetNextPoint(-1.4,0,virlen);
//  guide[17]->SetNextPoint(-1.4+virlen,0,virlen);
//  guide[17]->SetNextPoint(-1.4+virlen,0,0);


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

  Double_t xx[3],xxt[3];
  TLatex Tl;
  Tl.SetTextSize(0.06);
  xx[0]=0.2; xx[1]=0.0; xx[2]=-0.4; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"O");
  xx[0]=-0.8; xx[1]=0.2; xx[2]=0.00; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"#it{a}");
  xx[0]=-0.63; xx[1]=0.2; xx[2]=0.25; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"#it{b}");
  xx[0]=-0.42; xx[1]=-0.3; xx[2]=0.22; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"#it{c}");

  xx[0]=-2.0; xx[1]=0.2; xx[2]=1.85; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"#it{A}");
  xx[0]=-2.05; xx[1]=-1.3; xx[2]=0.0; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"#it{B}");
  xx[0]=-2.0; xx[1]=0.8; xx[2]=-0.4; view->WCtoNDC(xx,xxt);
  Tl.DrawLatex(xxt[0],xxt[1],"#it{C}");
}