a=1;(*楕円長軸*)
b=0.641;(*楕円中軸*)
c=0.431;(*楕円短軸*)

ClearAll[fpl20];
num=20;
pl=Table[{
theta1=RandomReal[{0,2 Pi}];
theta2=RandomReal[{0,2 Pi}];
data=Table[{
theta3=RandomReal[{0,2 Pi}];
p=Cos[theta1] Cos[theta3]-Sin[theta1] Cos[theta2] Sin[theta3];
q=Sin[theta2] Sin[theta3];
r=Sin[theta1] Cos[theta3]+Cos[theta1] Cos[theta2] Sin[theta3];
u=Sin[theta1] Sin[theta2];
v=Cos[theta2];
w=-Cos[theta1] Sin[theta2];
d==-Sin[theta1]Cos[theta2] Cos[theta3]-Cos[theta1] Sin[theta3];
e=Sin[theta2] Cos[theta3];
f=Cos[theta1] Cos[theta2] Cos[theta3]-Sin[theta1] Sin[theta3];
k1=u^2 a^2+v^2 b^2+w^2 c^2;
k2=(-2 u p) a^2+(-2 v q)  b^2+(-2 w r)  c^2;
k3=p^2 a^2+q^2 b^2+r^2 c^2;
kb=k1 k3-k2^2/4;
k0=Sqrt[(k1/kb-k3/kb)^2+(k2/kb)^2];
aa=1/Sqrt[(k1/kb+k3/kb-k0)/2];
cc=1/Sqrt[(k1/kb+k3/kb+k0)/2];
area=Pi aa cc;
l={aa,cc,area}
},{i,num}];
max=Flatten[MaximalBy[Flatten[data,1],Last]];
min=Flatten[MinimalBy[Flatten[data,1],Last]];
shape={max[[1]],max[[2]],min[[2]]};
shp={min[[2]]/max[[1]], max[[2]]/max[[1]]}
},{i,100000}];
fpl20=Flatten[pl,1];
ListPlot[fpl20]
fpl20il=Flatten[Take [fpl20, All, -1],1];
List[Length[Select[fpl20il,#≤ b (1.1)&&#≥b (0.9)&]],Length[Select[fpl20il,#≤ b (1.2)&&#≥b (0.8)&]],Length[Select[fpl20il,#≤ b (1.25)&&#≥b (0.75)&]],Length[Select[fpl20il,#≤ b (1.3)&&#≥b (0.7)&]],Length[Select[fpl20il,#≤ b (1.4)&&#≥b (0.6)&]],Length[Select[fpl20il,#≤ b (1.5)&&#≥b (0.5)&]]]
fpl20sl=Flatten[Take [fpl20, All, 1],1];
List[Length[Select[fpl20sl,#≤ c (1.1)&&#≥c (0.9)&]],Length[Select[fpl20sl,#≤ c (1.2)&&#≥c (0.8)&]],Length[Select[fpl20sl,#≤ c (1.25)&&#≥c (0.75)&]],Length[Select[fpl20sl,#≤ c (1.3)&&#≥c (0.7)&]],Length[Select[fpl20sl,#≤ c (1.4)&&#≥c (0.6)&]],Length[Select[fpl20sl,#≤ c (1.5)&&#≥c (0.5)&]]]

DensityHistogram[fpl20,{0.02},PlotRange->{{0,1},{0,1}},ColorFunction->(ColorData["LightTemperatureMap"][Rescale[#,{0,1/0.17102},{0,6}]]&),ImageSize->700,AspectRatio->1,ChartLegends->Placed[Automatic,Right],FrameStyle->Thickness[0.001],Method->{"DistributionAxes"->"Histogram"}]

ClearAll[fpl12];
num=12;
pl=Table[{
theta1=RandomReal[{0,2 Pi}];
theta2=RandomReal[{0,2 Pi}];
data=Table[{
theta3=RandomReal[{0,2 Pi}];
p=Cos[theta1] Cos[theta3]-Sin[theta1] Cos[theta2] Sin[theta3];
q=Sin[theta2] Sin[theta3];
r=Sin[theta1] Cos[theta3]+Cos[theta1] Cos[theta2] Sin[theta3];
u=Sin[theta1] Sin[theta2];
v=Cos[theta2];
w=-Cos[theta1] Sin[theta2];
d==-Sin[theta1]Cos[theta2] Cos[theta3]-Cos[theta1] Sin[theta3];
e=Sin[theta2] Cos[theta3];
f=Cos[theta1] Cos[theta2] Cos[theta3]-Sin[theta1] Sin[theta3];
k1=u^2 a^2+v^2 b^2+w^2 c^2;
k2=(-2 u p) a^2+(-2 v q)  b^2+(-2 w r)  c^2;
k3=p^2 a^2+q^2 b^2+r^2 c^2;
kb=k1 k3-k2^2/4;
k0=Sqrt[(k1/kb-k3/kb)^2+(k2/kb)^2];
aa=1/Sqrt[(k1/kb+k3/kb-k0)/2];
cc=1/Sqrt[(k1/kb+k3/kb+k0)/2];
area=Pi aa cc;
l={aa,cc,area}
},{i,num}];
max=Flatten[MaximalBy[Flatten[data,1],Last]];
min=Flatten[MinimalBy[Flatten[data,1],Last]];
shape={max[[1]],max[[2]],min[[2]]};
shp={min[[2]]/max[[1]], max[[2]]/max[[1]]}
},{i,100000}];
fpl12=Flatten[pl,1];
ListPlot[fpl12]
fpl12il=Flatten[Take [fpl12, All, -1],1];
List[Length[Select[fpl12il,#≤ b (1.1)&&#≥b (0.9)&]],Length[Select[fpl12il,#≤ b (1.2)&&#≥b (0.8)&]],Length[Select[fpl12il,#≤ b (1.25)&&#≥b (0.75)&]],Length[Select[fpl12il,#≤ b (1.3)&&#≥b (0.7)&]],Length[Select[fpl12il,#≤ b (1.4)&&#≥b (0.6)&]],Length[Select[fpl12il,#≤ b (1.5)&&#≥b (0.5)&]]]
fpl12sl=Flatten[Take [fpl12, All, 1],1];
List[Length[Select[fpl12sl,#≤ c (1.1)&&#≥c (0.9)&]],Length[Select[fpl12sl,#≤ c (1.2)&&#≥c (0.8)&]],Length[Select[fpl12sl,#≤ c (1.25)&&#≥c (0.75)&]],Length[Select[fpl12sl,#≤ c (1.3)&&#≥c (0.7)&]],Length[Select[fpl12sl,#≤ c (1.4)&&#≥c (0.6)&]],Length[Select[fpl12sl,#≤ c (1.5)&&#≥c (0.5)&]]]

DensityHistogram[fpl12,{0.02},PlotRange->{{0,1},{0,1}},ColorFunction->(ColorData["LightTemperatureMap"][Rescale[#,{0,1/0.11028},{0,6}]]&),ImageSize->700,AspectRatio->1,ChartLegends->Placed[Automatic,Right],FrameStyle->Thickness[0.001],Method->{"DistributionAxes"->"Histogram"}]

ClearAll[fpl3];
num=3;
pl=Table[{
theta1=RandomReal[{0,2 Pi}];
theta2=RandomReal[{0,2 Pi}];
data=Table[{
theta3=RandomReal[{0,2 Pi}];
p=Cos[theta1] Cos[theta3]-Sin[theta1] Cos[theta2] Sin[theta3];
q=Sin[theta2] Sin[theta3];
r=Sin[theta1] Cos[theta3]+Cos[theta1] Cos[theta2] Sin[theta3];
u=Sin[theta1] Sin[theta2];
v=Cos[theta2];
w=-Cos[theta1] Sin[theta2];
d==-Sin[theta1]Cos[theta2] Cos[theta3]-Cos[theta1] Sin[theta3];
e=Sin[theta2] Cos[theta3];
f=Cos[theta1] Cos[theta2] Cos[theta3]-Sin[theta1] Sin[theta3];
k1=u^2 a^2+v^2 b^2+w^2 c^2;
k2=(-2 u p) a^2+(-2 v q)  b^2+(-2 w r)  c^2;
k3=p^2 a^2+q^2 b^2+r^2 c^2;
kb=k1 k3-k2^2/4;
k0=Sqrt[(k1/kb-k3/kb)^2+(k2/kb)^2];
aa=1/Sqrt[(k1/kb+k3/kb-k0)/2];
cc=1/Sqrt[(k1/kb+k3/kb+k0)/2];
area=Pi aa cc;
l={aa,cc,area}
},{i,num}];
max=Flatten[MaximalBy[Flatten[data,1],Last]];
min=Flatten[MinimalBy[Flatten[data,1],Last]];
shape={max[[1]],max[[2]],min[[2]]};
shp={min[[2]]/max[[1]], max[[2]]/max[[1]]}
},{i,100000}];
fpl3=Flatten[pl,1];
ListPlot[fpl3]
fpl3il=Flatten[Take [fpl3, All, -1],1];
List[Length[Select[fpl3il,#≤ b (1.1)&&#≥b (0.9)&]],Length[Select[fpl3il,#≤ b (1.2)&&#≥b (0.8)&]],Length[Select[fpl3il,#≤ b (1.25)&&#≥b (0.75)&]],Length[Select[fpl3il,#≤ b (1.3)&&#≥b (0.7)&]],Length[Select[fpl3il,#≤ b (1.4)&&#≥b (0.6)&]],Length[Select[fpl3il,#≤ b (1.5)&&#≥b (0.5)&]],Length[Select[fpl3il,#≤ b+0.1&&#≥b-0.1&]]]
fpl3sl=Flatten[Take [fpl3, All, 1],1];
List[Length[Select[fpl3sl,#≤ c (1.1)&&#≥c (0.9)&]],Length[Select[fpl3sl,#≤ c (1.2)&&#≥c (0.8)&]],Length[Select[fpl3sl,#≤ c (1.25)&&#≥c (0.75)&]],Length[Select[fpl3sl,#≤ c (1.3)&&#≥c (0.7)&]],Length[Select[fpl3sl,#≤ c (1.4)&&#≥c (0.6)&]],Length[Select[fpl3sl,#≤ c (1.5)&&#≥c (0.5)&]],Length[Select[fpl3sl,#≤ c+0.1&&#≥c-0.1&]]]

DensityHistogram[fpl3,{0.02},PlotRange->{{0,1},{0,1}},ColorFunction->(ColorData["LightTemperatureMap"][Rescale[#,{0,1/0.02700},{0,6}]]&),ImageSize->700,AspectRatio->1,ChartLegends->Placed[Automatic,Right],FrameStyle->Thickness[0.001],Method->{"DistributionAxes"->"Histogram"}]

ClearAll[fpl6];
num=6;
pl=Table[{
theta1=RandomReal[{0,2 Pi}];
theta2=RandomReal[{0,2 Pi}];
data=Table[{
theta3=RandomReal[{0,2 Pi}];
p=Cos[theta1] Cos[theta3]-Sin[theta1] Cos[theta2] Sin[theta3];
q=Sin[theta2] Sin[theta3];
r=Sin[theta1] Cos[theta3]+Cos[theta1] Cos[theta2] Sin[theta3];
u=Sin[theta1] Sin[theta2];
v=Cos[theta2];
w=-Cos[theta1] Sin[theta2];
d==-Sin[theta1]Cos[theta2] Cos[theta3]-Cos[theta1] Sin[theta3];
e=Sin[theta2] Cos[theta3];
f=Cos[theta1] Cos[theta2] Cos[theta3]-Sin[theta1] Sin[theta3];
k1=u^2 a^2+v^2 b^2+w^2 c^2;
k2=(-2 u p) a^2+(-2 v q)  b^2+(-2 w r)  c^2;
k3=p^2 a^2+q^2 b^2+r^2 c^2;
kb=k1 k3-k2^2/4;
k0=Sqrt[(k1/kb-k3/kb)^2+(k2/kb)^2];
aa=1/Sqrt[(k1/kb+k3/kb-k0)/2];
cc=1/Sqrt[(k1/kb+k3/kb+k0)/2];
area=Pi aa cc;
l={aa,cc,area}
},{i,num}];
max=Flatten[MaximalBy[Flatten[data,1],Last]];
min=Flatten[MinimalBy[Flatten[data,1],Last]];
shape={max[[1]],max[[2]],min[[2]]};
shp={min[[2]]/max[[1]], max[[2]]/max[[1]]}
},{i,100000}];
fpl6=Flatten[pl,1];
ListPlot[fpl6]
fpl6il=Flatten[Take [fpl6, All, -1],1];
List[Length[Select[fpl6il,#≤ b (1.1)&&#≥b (0.9)&]],Length[Select[fpl6il,#≤ b (1.2)&&#≥b (0.8)&]],Length[Select[fpl6il,#≤ b (1.25)&&#≥b (0.75)&]],Length[Select[fpl6il,#≤ b (1.3)&&#≥b (0.7)&]],Length[Select[fpl6il,#≤ b (1.4)&&#≥b (0.6)&]],Length[Select[fpl6il,#≤ b (1.5)&&#≥b (0.5)&]],Length[Select[fpl6il,#≤ b+0.1&&#≥b-0.1&]]]
fpl6sl=Flatten[Take [fpl6, All, 1],1];
List[Length[Select[fpl6sl,#≤ c (1.1)&&#≥c (0.9)&]],Length[Select[fpl6sl,#≤ c (1.2)&&#≥c (0.8)&]],Length[Select[fpl6sl,#≤ c (1.25)&&#≥c (0.75)&]],Length[Select[fpl6sl,#≤ c (1.3)&&#≥c (0.7)&]],Length[Select[fpl6sl,#≤ c (1.4)&&#≥c (0.6)&]],Length[Select[fpl6sl,#≤ c (1.5)&&#≥c (0.5)&]],Length[Select[fpl6sl,#≤ c+0.1&&#≥c-0.1&]]]

DensityHistogram[fpl6,{0.02},PlotRange->{{0,1},{0,1}},ColorFunction->(ColorData["LightTemperatureMap"][Rescale[#,{0,1/0.07383},{0,6}]]&),ImageSize->700,AspectRatio->1,ChartLegends->Placed[Automatic,Right],FrameStyle->Thickness[0.001],Method->{"DistributionAxes"->"Histogram"}]