Plot solutions (integral scale) using disk representation
Yang Oct 15 2015

Deterministic sampling on (S^2) using HEALPix

Nside=8;
Vector={};
For[i=1,i<Nside,i++,
For[j=1,j≤4i,j++,
z=1-((i^2)/(3(Nside^2)));
ϕ=(π/(2i))(j-(1/2));
x=Cos[ϕ](Sqrt[1-(z^2)]);
y=Sin[ϕ](Sqrt[1-(z^2)]);
AppendTo[Vector,(({x,y,z})/(Norm[{x,y,z}]))];
]
]
For[i=Nside,i<2Nside,i++,
For[j=1,j≤4Nside,j++,
z=(4/3)-((2i)/(3Nside));
s=Mod[i-Nside+1,2];
ϕ=(π/(2Nside))(j-(s/2));
x=Cos[ϕ](Sqrt[1-(z^2)]);
y=Sin[ϕ](Sqrt[1-(z^2)]);
AppendTo[Vector,(({x,y,z})/(Norm[{x,y,z}]))];
]
]
Vector=Join[Vector,-Vector];
For[j=1,j≤4Nside,j++,
z=0;
s=Mod[Nside+1,2];
ϕ=(π/(2Nside))(j-(s/2));
x=Cos[ϕ](Sqrt[1-(z^2)]);
y=Sin[ϕ](Sqrt[1-(z^2)]);
AppendTo[Vector,(({x,y,z})/(Norm[{x,y,z}]))];
]

Length[Vector]

Plot dots

VecN=Length[Vector];
ZeroPt={0,0,0};
ZAxis=Arrow[{ZeroPt,{0,0,1}}];
Points=Table[{PointSize[Medium],Point[Vector[[i]]]},{i,VecN}];
Show[Graphics3D[Points],Graphics3D[{Red,ZAxis}]]

Convert points to α/β

alphaBeta=Table[Null,{VecN}];
For[i=1,i≤VecN,i++,
z=Vector[[i]];
β=ArcCos[z[[3]]];
If[z[[1]]==0,
(*true*)
Which[z[[2]]≤0,α=((3π)/2),z[[2]]>0,α=(π/2)],
(*false*)
tanα=((z[[2]])/(z[[1]]));
Which[
z[[1]]<0,α=ArcTan[tanα]+π,
z[[1]]>0&&z[[2]]≥0,α=ArcTan[tanα],
z[[1]]>0&&z[[2]]<0,α=ArcTan[tanα]+2π
];
];
alphaBeta[[i]]=N[{α,β},8];
]

alphaBeta[[1]]

Convert data to rotation matrices

Clear[α,β,γ]
ma[α_]:=(({{Cos[α], -Sin[α], 0}, {Sin[α], Cos[α], 0}, {0, 0, 1}}));
mb[β_]:=(({{Cos[β], 0, Sin[β]}, {0, 1, 0}, {-Sin[β], 0, Cos[β]}}));
mc[γ_]:=(({{Cos[γ], -Sin[γ], 0}, {Sin[γ], Cos[γ], 0}, {0, 0, 1}}));
mabc[α_,β_,γ_]:=ma[α].mb[β].mc[γ];
mabc[α,β,γ]//MatrixForm

angleN=48;
rotMList=Flatten[Table[mabc[alphaBeta[[i,1]],alphaBeta[[i,2]],γ],{i,VecN},{γ,0,2π-((2π)/angleN),((2π)/angleN)}],1];

Build polygons

diskR=0.9;
sphereR=15;
dγ=((2π)/angleN);
polygon0={{0,0,sphereR},{Cos[((dγ)/2)]diskR,-Sin[((dγ)/2)]diskR,sphereR},{Cos[((dγ)/2)]diskR,Sin[((dγ)/2)]diskR,sphereR}};
polyList=Table[polygon0.rotMList[[i]]\[Transpose],{i,Length[rotMList]}];

Show[Table[Graphics3D[{EdgeForm[],Polygon[polyList[[i]]]}],{i,Length[rotMList]}]]

Convert rotation matrices to quaternions

quas=Table[Null,{Length[rotMList]}];
For[i=1,i≤Length[rotMList],i++,
rotM=rotMList[[i]];
Which[
(*case 1*)
1+rotM[[1,1]]+rotM[[2,2]]+rotM[[3,3]]>0,
q1=((Sqrt[1+rotM[[1,1]]+rotM[[2,2]]+rotM[[3,3]]])/2);
q2=((rotM[[3,2]]-rotM[[2,3]])/(4q1));
q3=((rotM[[1,3]]-rotM[[3,1]])/(4q1));
q4=((rotM[[2,1]]-rotM[[1,2]])/(4q1)),
(*case 2*)
1+rotM[[1,1]]-rotM[[2,2]]-rotM[[3,3]]>0,
q2=((Sqrt[1+rotM[[1,1]]-rotM[[2,2]]-rotM[[3,3]]])/2);
q1=((rotM[[3,2]]-rotM[[2,3]])/(4q2));
q3=((rotM[[1,2]]+rotM[[2,1]])/(4q2));
q4=((rotM[[1,3]]+rotM[[3,1]])/(4q2)),
(*case 3*)
1-rotM[[1,1]]+rotM[[2,2]]-rotM[[3,3]]>0,
q3=((Sqrt[1-rotM[[1,1]]+rotM[[2,2]]-rotM[[3,3]]])/2);
q1=((rotM[[1,3]]-rotM[[3,1]])/(4q3));
q2=((rotM[[1,2]]+rotM[[2,1]])/(4q3));
q4=((rotM[[2,3]]+rotM[[3,2]])/(4q3)),
(*case 4*)
1-rotM[[1,1]]-rotM[[2,2]]+rotM[[3,3]]>0,
q4=((Sqrt[1-rotM[[1,1]]-rotM[[2,2]]+rotM[[3,3]]])/2);
q1=((rotM[[2,1]]-rotM[[1,2]])/(4q4));
q2=((rotM[[1,3]]+rotM[[3,1]])/(4q4));
q3=((rotM[[2,3]]+rotM[[3,2]])/(4q4))
]
quas[[i]]={q1,q2,q3,q4};
]

Length[quas]

Define Bingham distribution on (S^3)(sol1)

QuaList=quas;
Ndat=Length[quas];

Σ1=(({{1, 0, 0, 0}, {0, 0.7179424849229702`, 0.5975391483281293`, 0.3589865224229396`}, {0, -0.5116735783191485`, 0.09961914423048843`, 0.8527017172796637`}, {0, 0.4734073365676751`, -0.7948943550037976`, 0.37993093283318946`}})).(({{-0.022729940000000018`, -0.24434603`, 0.91056426`, 0.33264058`}, {-0.97860921`, 0.16874238000000003`, 0.05821440000000003`, -0.10227300000000006`}, {-0.1385124800000001`, -0.12402565999999998`, -0.3693558`, 0.91049815`}, {0.15040051999999998`, 0.94680374`, 0.17623496999999994`, 0.22334336`}}))\[Transpose].(({{1.125, 0, 0, 0}, {0, -1.125, 0, 0}, {0, 0, -0.375, 0}, {0, 0, 0, 0.375}})).(({{-0.022729940000000018`, -0.24434603`, 0.91056426`, 0.33264058`}, {-0.97860921`, 0.16874238000000003`, 0.05821440000000003`, -0.10227300000000006`}, {-0.1385124800000001`, -0.12402565999999998`, -0.3693558`, 0.91049815`}, {0.15040051999999998`, 0.94680374`, 0.17623496999999994`, 0.22334336`}})).(({{1, 0, 0, 0}, {0, 0.7179424849229702`, 0.5975391483281293`, 0.3589865224229396`}, {0, -0.5116735783191485`, 0.09961914423048843`, 0.8527017172796637`}, {0, 0.4734073365676751`, -0.7948943550037976`, 0.37993093283318946`}}))\[Transpose];
Σ=Σ1;
Bham[Σ_,q_]:=Exp[q.Σ.q\[Transpose]];

BhamList=Table[Bham[Σ,{QuaList[[i]]}][[1]][[1]],{i,Ndat}];
BhamMax=Max[BhamList]
BhamTotal=Total[BhamList]
AreaProb=0;
j=0;
For[i=1,i≤Ndat,i++,
If[BhamList[[i]]≥0.831BhamMax,
AreaProb=AreaProb+BhamList[[i]];
j++;
];
]
AreaProb=(AreaProb/BhamTotal)
N[(j/Ndat)]

Prob={"Prob",0.05,0.10,0.20,0.40,0.80};
Percent={"Percent",0,0,0,0.097,0.351};
Ratio={"Ratio",0.831,0.736,0.599,0.426,0.143};
Grid[{Prob,Percent,Ratio},Frame -> All, Background->{{LightGray,None},{None,None}}]

ColorScale[x_]:=Blend[{{0.123,White},{0.295,Blue},{0.359,Green},{0.563,Yellow},{0.715,Orange},{0.818,Red}},x];

BhamMax=17.6;

ArL=1.5sphereR;
DiskList=Table[Graphics3D[{EdgeForm[],ColorScale[(1/BhamMax)Bham[Σ,{QuaList[[i]]}][[1]][[1]]],Polygon[polyList[[i]]]}],{i,Length[rotMList]}];
DiskFig=Show[(*Graphics3D[{Black,Text[Style["X",Bold,40],{1.1ArL,0,0}],Text[Style["Y",Bold,40],{0,1.1ArL,0}],Text[Style["Z",Bold,40],{0,0,1.1ArL}]}],*)Graphics3D[{Black,Arrowheads[0.05],Thickness[0.02],Arrow[Tube[{{0,0,0},{-ArL,0,0}}]],Arrow[Tube[{{0,0,0},{0,ArL,0}}]],Arrow[Tube[{{0,0,0},{0,0,ArL}}]]}],Graphics3D[{RGBColor[0.85,0.85,0.85],Sphere[{0,0,0},.99sphereR]}],DiskList,Lighting->{{"Ambient",White}},Boxed->False,ViewPoint->5*{-3,-1,(*-0.5*)0},SphericalRegion->True]

Export["/Users/yangqi/Documents/projects/multi_model/sol2_res_view1_rescale.jpeg",Image[DiskFig,ImageSize->6000]];

Define Bingham distribution on (S^3)(sol2)

QuaList=quas;
Ndat=Length[quas];

(*solution 2: x1*)
Σ2=(({{1, 0, 0, 0}, {0, 0.7179424849229702`, 0.5975391483281293`, 0.3589865224229396`}, {0, -0.5116735783191485`, 0.09961914423048843`, 0.8527017172796637`}, {0, 0.4734073365676751`, -0.7948943550037976`, 0.37993093283318946`}})).(({{-0.10057290732213947`, 0.11747471240125662`, 0.9879471946690229`, -0.0067173510122877445`}, {0.06610566622064512`, 0.08489027737410254`, 0.0033952178126328758`, 0.9941891943673887`}, {-0.0019352845632774063`, 0.9894071279183784`, -0.11841611543842674`, -0.0839488743369324`}, {-0.9927292611180333`, -0.008177284261643214`, -0.09963150428496932`, 0.06704706925781526`}}))\[Transpose].(({{2.875, 0, 0, 0}, {0, -3.875, 0, 0}, {0, 0, -0.625, 0}, {0, 0, 0, 1.625}})).(({{-0.10057290732213947`, 0.11747471240125662`, 0.9879471946690229`, -0.0067173510122877445`}, {0.06610566622064512`, 0.08489027737410254`, 0.0033952178126328758`, 0.9941891943673887`}, {-0.0019352845632774063`, 0.9894071279183784`, -0.11841611543842674`, -0.0839488743369324`}, {-0.9927292611180333`, -0.008177284261643214`, -0.09963150428496932`, 0.06704706925781526`}})).(({{1, 0, 0, 0}, {0, 0.7179424849229702`, 0.5975391483281293`, 0.3589865224229396`}, {0, -0.5116735783191485`, 0.09961914423048843`, 0.8527017172796637`}, {0, 0.4734073365676751`, -0.7948943550037976`, 0.37993093283318946`}}))\[Transpose];
Σ=Σ2;
Bham[Σ_,q_]:=Exp[q.Σ.q\[Transpose]];

BhamList=Table[Bham[Σ,{QuaList[[i]]}][[1]][[1]],{i,Ndat}];
BhamMax=Max[BhamList]
BhamTotal=Total[BhamList]
AreaProb=0;
j=0;
For[i=1,i≤Ndat,i++,
If[BhamList[[i]]≥0.818BhamMax,
AreaProb=AreaProb+BhamList[[i]];
j++;
];
]
AreaProb=(AreaProb/BhamTotal)
N[(j/Ndat)]

Prob={"Prob",0.05,0.10,0.20,0.40,0.80};
Percent={"Percent",0,0,0,0.097,0.351};
Ratio={"Ratio",0.818,0.715,0.563,0.359,0.123};
Grid[{Prob,Percent,Ratio},Frame -> All, Background->{{LightGray,None},{None,None}}]

ColorScale[x_]:=Blend[{{0.123,White},{0.295,Blue},{0.359,Green},{0.563,Yellow},{0.715,Orange},{0.818,Red}},x];

0.123+0.7315(0.359-0.123)

ArL=1.5sphereR;
DiskList=Table[Graphics3D[{EdgeForm[],ColorScale[(1/BhamMax)Bham[Σ,{QuaList[[i]]}][[1]][[1]]],Polygon[polyList[[i]]]}],{i,Length[rotMList]}];
DiskFig=Show[(*Graphics3D[{Black,(*Text[Style["X",Bold,40],{1.1ArL,0,0}],*)Text[Style["Y",Bold,40],{0,1.1ArL,0}],Text[Style["Z",Bold,40],{0,0,1.1ArL}]}],*)Graphics3D[{Black,Arrowheads[0.05],Thickness[0.02],Arrow[Tube[{{0,0,0},{ArL,0,0}}]],Arrow[Tube[{{0,0,0},{0,ArL,0}}]],Arrow[Tube[{{0,0,0},{0,0,ArL}}]]}],Graphics3D[{RGBColor[0.85,0.85,0.85],Sphere[{0,0,0},.99sphereR]}],DiskList,Lighting->{{"Ambient",White}},Boxed->False,ViewPoint->-5*{-3,-1,(*-0.5*)0},SphericalRegion->True]

Export["C:\\Users\\Eninala\\Desktop\\SpAN\\ZLBT-CdpA\\Distribution Calculation\\Figures for presentation\\DiskFigs\\sol2_view2.jpeg",Image[DiskFig,ImageSize->6000]];

ArL=1.5sphereR;
DiskList=Table[Graphics3D[{EdgeForm[],ColorScale[(1/BhamMax)Bham[Σ,{QuaList[[i]]}][[1]][[1]]],Polygon[polyList[[i]]]}],{i,Length[rotMList]}];
DiskFig=Show[(*Graphics3D[{Black,(*Text[Style["X",Bold,40],{1.1ArL,0,0}],*)Text[Style["Y",Bold,40],{0,1.1ArL,0}],Text[Style["Z",Bold,40],{0,0,1.1ArL}]}],*)Graphics3D[{Black,Arrowheads[0.05],Thickness[0.02],Arrow[Tube[{{0,0,0},{-ArL,0,0}}]],Arrow[Tube[{{0,0,0},{0,ArL,0}}]],Arrow[Tube[{{0,0,0},{0,0,ArL}}]]}],Graphics3D[{RGBColor[0.85,0.85,0.85],Sphere[{0,0,0},.99sphereR]}],DiskList,Lighting->{{"Ambient",White}},Boxed->False,ViewPoint->5*{-2,2,(*0.5*)0},SphericalRegion->True]

Export["C:\\Users\\Eninala\\Desktop\\SpAN\\ZLBT-CdpA\\Distribution Calculation\\Figures for presentation\\DiskFigs\\sol2_view3.jpeg",Image[DiskFig,ImageSize->2000]];

ColorScale[x_]:=Blend[{{0.2,White},{0.4,Blue},{0.6,Green},{0.8,Yellow},{0.9,Orange},{0.95,Red}},x];

Nline=6000;
Bar=Graphics[Table[{ColorScale[i/Nline],Line[{{((10i)/Nline),0},{((10i)/Nline),.4}}]},{i,1,Nline}]];
Labels=Graphics[Table[Text[Style[ToString[10i],Bold,80],{i,-.3}],{i,1,9}]];
bBox=Graphics[{EdgeForm[{Directive[Thickness[0.005]],Black}],RGBColor[1,1,1,0],Rectangle[{0,0},{10,0.4}]}];
PopLegend=Show[Bar,bBox,Labels,PlotLabel->Style["Percentile",Black,Bold,80,FontFamily->"Arial"],ImageSize->2000]

Export["C:\\Users\\Eninala\\Desktop\\SpAN\\ZLBT-CdpA\\Distribution Calculation\\Figures for presentation\\DiskFigs\\legend.jpeg",PopLegend];

Clear[α,β,γ]
ma[α_]:=(({{Cos[α], -Sin[α], 0}, {Sin[α], Cos[α], 0}, {0, 0, 1}}));
mb[β_]:=(({{Cos[β], 0, Sin[β]}, {0, 1, 0}, {-Sin[β], 0, Cos[β]}}));
mc[γ_]:=(({{Cos[γ], -Sin[γ], 0}, {Sin[γ], Cos[γ], 0}, {0, 0, 1}}));
mabc[α_,β_,γ_]:=ma[α].mb[β].mc[γ];
mabc[α,β,γ]//MatrixForm

ColorScale[x_]:=Blend[{{0.2,White},{0.4,Blue},{0.6,Green},{0.8,Yellow},{0.9,Orange},{0.95,Red}},x];
Rmat=mabc[((5.5π)/4),(π/3),0];
diskR=1;
sphereR=15;
dγ=((2π)/angleN);
angleN=1000;
rotMList2=Flatten[Table[mabc[((5.5π)/4),(π/3),γ],{i,1},{γ,π,3π-((2π)/angleN),((2π)/angleN)}],1];
polygon0={{0,0,0.75sphereR},{Cos[((dγ)/2)] sphereR,-Sin[((dγ)/2)] sphereR,0.75sphereR},{Cos[((dγ)/2)] sphereR,Sin[((dγ)/2)] sphereR,0.75sphereR}};
polyList=Table[polygon0.rotMList2[[i]]\[Transpose],{i,Length[rotMList2]}];
ExampleDisk=Table[Graphics3D[{EdgeForm[],ColorScale[3PDF[VonMisesDistribution[0,1],2π (i/angleN)-π]],Polygon[polyList[[i]]]}],{i,Length[rotMList2]}];
Show[ExampleDisk,Lighting->{{"Ambient",White}}]

ArL=1.5sphereR;
CoordinateView=Show[ExampleDisk,(*Graphics3D[{Black,Text[Style["(Subscript[X, 1])",Bold,18],{1.2ArL,0,0}],Text[Style["(Subscript[Y, 1])",Bold,18],{0,1.2ArL,0}],Text[Style["(Subscript[Z, 1])",Bold,18],{0,0,1.2ArL}]}],*)Graphics3D[{Red,Arrow[Tube[{{0,0,0},{-1.1ArL,0,0}},0.3]],Arrow[Tube[{{0,0,0},{0,1.1ArL,0}},0.3]],Arrow[Tube[{{0,0,0},{0,0,1.1ArL}},0.3]]}],
(*Graphics3D[{Black,Text[Style["(Subscript[X, 2])",Bold,18],Rmat.{1.1ArL,0,.5ArL}],Text[Style["(Subscript[Y, 2])",Bold,18],Rmat.{0,1.1ArL,.5ArL}],Text[Style["(Subscript[Z, 2])",Bold,18],Rmat.{0,0,1.1ArL}]}],*)Graphics3D[{Blue,Arrow[Tube[{Rmat.{0,0,.5ArL},Rmat.{ArL,0,.5ArL}},0.3]],Arrow[Tube[{Rmat.{0,0,.5ArL},Rmat.{0,ArL,.5ArL}},0.3]],Arrow[Tube[{{0,0,0},Rmat.{0,0,ArL}},0.3]]}],Lighting->{{"Ambient",White}},Boxed->False,ImageSize->Large,ViewPoint->15*{-3,-1,-0.5}]

Image[CoordinateView,ImageSize->800]

ArL=1.5sphereR;
CoordinateView=Image[Show[ExampleDisk,Graphics3D[{Black,Text[Style["(Subscript[X, 1])",Bold,18],{1.2ArL,0,0}],Text[Style["(Subscript[Y, 1])",Bold,18],{0,1.2ArL,0}],Text[Style["(Subscript[Z, 1])",Bold,18],{0,0,1.2ArL}]}],Graphics3D[{Black,Arrow[Tube[{{0,0,0},{1.1ArL,0,0}}]],Arrow[Tube[{{0,0,0},{0,1.1ArL,0}}]],Arrow[Tube[{{0,0,0},{0,0,1.1ArL}}]]}],
Graphics3D[{Black,Text[Style["(Subscript[X, 2])",Bold,18],Rmat.{1.1ArL,0,.5ArL}],Text[Style["(Subscript[Y, 2])",Bold,18],Rmat.{0,1.1ArL,.5ArL}],Text[Style["(Subscript[Z, 2])",Bold,18],Rmat.{0,0,1.1ArL}]}],Graphics3D[{Black,Arrow[Tube[{Rmat.{0,0,.5ArL},Rmat.{ArL,0,.5ArL}}]],Arrow[Tube[{Rmat.{0,0,.5ArL},Rmat.{0,ArL,.5ArL}}]],Arrow[Tube[{{0,0,0},Rmat.{0,0,ArL}}]]}],Lighting->{{"Ambient",White}},Boxed->False,ImageSize->Medium,ViewPoint->{3,-1.8,1.8}]]