0. small programs

Clear["*Global`*"]

(*t=τ(golden ratio). For "show", just add //N, then all numbers turn to N*)
t=(1+Sqrt[5])/2//N;

length[x_,y_]:=FullSimplify[Sqrt[Total[(x-y)^2]]];
Nlength[x_,y_]:=Sqrt[Total[(N[x]-N[y])^2]];
myrule={1.618033988749895`->"τ",1.6180339887498951`->"τ",-1.618033988749895`->"-τ",
-1.6180339887498951`->"-τ",0.6180339887498948`->"(τ^(-1))",-0.6180339887498948`->"-(τ^(-1))",1.0000000000000002`->"1",-1.0000000000000002`->"-1",2/(1+Sqrt[5])->"(τ^(-1))",(1+Sqrt[5])/2->"τ",-2/(1+Sqrt[5])->"-(τ^(-1))",-(1+Sqrt[5])/2->"-τ",(-1+Sqrt[5])/2->"(τ^(-1))",(1-Sqrt[5])/2->"-(τ^(-1))"};


nearbys[x_]:=Select[vecs24,Nlength[#,x]^2-(-1+Sqrt[5])^2<0.000001&]//Quiet;

(*aとbを τ:1に内分する点*)
Clear[a,b];
naibun[a_,b_]=1/(1+t) a+t/(1+t) b;

(*xを頂点とし，「正三角形」yを底辺とする正四面体が存在するかチェックする述語関数*)
NtetraQ[x_,y_List]:=(Nlength[x,y[[1]]]^2+2Sqrt[5]-6)^2+(Nlength[x,y[[2]]]^2+2Sqrt[5]-6)^2+(Nlength[x,y[[3]]]^2+2Sqrt[5]-6)^2<0.000001;

tetraQ[x_,y_List]:=(length[x,y[[1]]]^2+2Sqrt[5]-6)^2+(length[x,y[[2]]]^2+2Sqrt[5]-6)^2+(length[x,y[[3]]]^2+2Sqrt[5]-6)^2==0;

icosahedronQ[x_,y_,z_,w_]:=((length[x,y]^2+2Sqrt[5]-6)^2==0)&&((length[x,z]^2+2Sqrt[5]-6)^2==0)&&((length[x,w]^2+2Sqrt[5]-6)^2==0)&&((length[y,z]^2+2Sqrt[5]-6)^2==0)&&((length[y,w]^2+2Sqrt[5]-6)^2==0)&&((length[z,w]^2+2Sqrt[5]-6)^2==0);

icosahedronQ[x_,{y_,z_,w_}]:=icosahedronQ[x,y,z,w]
icosahedronQ[{x_,y_,z_,w_}]:=icosahedronQ[x,y,z,w]

findfuta[x_]:=Select[futa,length[#,x]^2==(Sqrt[5]-1)^2&]

proj[mode_Integer,ushift_,x_List]:=Switch[mode,
1,Table[Drop[x[[k]],{4}],{k,1,Length[x]}],
2,Table[Drop[x[[k]],{3}],{k,1,Length[x]}],
3,Table[1/(x[[k,4]]+ushift) * Drop[x[[k]],{4}],{k,1,Length[x]}]
];

1. 正24胞体の胞(正8面体)の辺を1:τに内分する点を使い,  頂点の座標を作成するプログラム．(c1111などは各正8面体から作られる12個の頂点. cells24はその総和．また cells24//Union//Flatten = vecs24)

正24胞体の座標を作成してから　600胞体の座標を作成するプログラム.（直接作成するプログラムは 「 点と辺makingOld.nb」参照)

(*vertices of Octas(24 cells) tetrasC の作成に使う.*)
(*VofOctas=Module[{base,p,q},
base={{p,q,0,0},{p,0,q,0},{p,0,0,q},{0,p,q,0},{0,p,0,q},{0,0,p,q}};
t*Flatten[{base/.{p->1,q->1},base/.{p->1,q->-1},base/.{p->-1,q->1},base/.{p->-1,q->-1}},1]];*)

cellOctaA[x_List]:=(
mixed:=Partition[Riffle[x,{0,0,0,0}],2];
combi={{1,1,2,2},{2,2,1,1},{1,2,1,2},{2,1,2,1},{1,2,2,1},{2,1,1,2}};
Table[Table[mixed[[i,combi[[k,i]]]],{i,1,4}],{k,1,6}])

cellOctaB[x_List]:=(pos:=Position[x,0];
posc:=If[MemberQ[x,1],Position[x,1],Position[x,-1]];cellOctaB0:=Table[{ReplacePart[x,k->1],ReplacePart[x,k->-1]},{k,pos}];
If[EvenQ[posc[[1,1]]],Flatten[cellOctaB0[[{1,3,2}]],1],Flatten[cellOctaB0,1]]);
(*600胞体の座標をうまく作るために 偶奇で場合わけ*)

o1111=cellOctaA[{1,1,1,1}];o1112=cellOctaA[{1,1,1,-1}];o1121=cellOctaA[{1,1,-1,1}];
o1122=cellOctaA[{1,1,-1,-1}];o1211=cellOctaA[{1,-1,1,1}];o1212=cellOctaA[{1,-1,1,-1}];
o1221=cellOctaA[{1,-1,-1,1}];o1222=cellOctaA[{1,-1,-1,-1}];o2111=cellOctaA[{-1,1,1,1}];
o2112=cellOctaA[{-1,1,1,-1}];o2121=cellOctaA[{-1,1,-1,1}];o2122=cellOctaA[{-1,1,-1,-1}];
o2211=cellOctaA[{-1,-1,1,1}];o2212=cellOctaA[{-1,-1,1,-1}];o2221=cellOctaA[{-1,-1,-1,1}];
o2222=cellOctaA[{-1,-1,-1,-1}];
o1000=cellOctaB[{1,0,0,0}];o0100=cellOctaB[{0,1,0,0}];o0010=cellOctaB[{0,0,1,0}];o0001=cellOctaB[{0,0,0,1}];o2000=cellOctaB[{-1,0,0,0}];o0200=cellOctaB[{0,-1,0,0}];o0020=cellOctaB[{0,0,-1,0}];o0002=cellOctaB[{0,0,0,-1}];

octas={o1111,o1112,o1121,o1122,o1211,o1212,o1221,o1222,
o2111,o2112,o2121,o2122,o2211,o2212,o2221,o2222,o1000,o0100,o0010,o0001,o2000,o0200,o0020,o0002};

Clear[a,b];
naibun[a_,b_]=1/(1+t) a+t/(1+t) b;

creatCell[x_List]:=(
parted:=t*Partition[x,2];
FullSimplify[{
naibun[parted[[3,2]],parted[[2,2]]],
naibun[parted[[2,2]],parted[[1,1]]],
naibun[parted[[3,1]],parted[[2,2]]],
naibun[parted[[2,2]],parted[[1,2]]],
naibun[parted[[1,2]],parted[[3,2]]],
naibun[parted[[1,1]],parted[[3,2]]],
naibun[parted[[2,1]],parted[[1,1]]],
naibun[parted[[1,1]],parted[[3,1]]],
naibun[parted[[1,2]],parted[[3,1]]],
naibun[parted[[2,1]],parted[[1,2]]],
naibun[parted[[3,2]],parted[[2,1]]],
naibun[parted[[3,1]],parted[[2,1]]]}])

(*cells24 は 正20面体上の頂点の集合．(ダブルカウントあり). futa はそれ以外の頂点の集合.*)
cells24=Map[creatCell,octas];

{c1111,c1112,c1121,c1122,c1211,c1212,c1221,c1222,
c2111,c2112,c2121,c2122,c2211,c2212,c2221,c2222,c1000,c0100,c0010,c0001,c2000,c0200,c0020,c0002}=cells24;


futa1={{1,1,1,1},{1,1,1,-1},{1,1,-1,1},{1,1,-1,-1},
{1,-1,1,1},{1,-1,1,-1},{1,-1,-1,1},{1,-1,-1,-1},
{-1,1,1,1},{-1,1,1,-1},{-1,1,-1,1},{-1,1,-1,-1},
{-1,-1,1,1},{-1,-1,1,-1},{-1,-1,-1,1},{-1,-1,-1,-1}
};

futa2={{2,0,0,0},{0,2,0,0},{0,0,2,0},{0,0,0,2},{-2,0,0,0},{0,-2,0,0},{0,0,-2,0},{0,0,0,-2}};
futa=Flatten[{futa1,futa2},1];

vecs24=Union[Flatten[cells24,1]];

2. [TetrasAの作成] 正24胞体の正8面体胞 を使い,  正20面体の表面を作成．（c1111上の三角形が t1111．その総和がtriangles24. )  さらに，蓋をして 正四面体の胞: tetrasAを 24*20=480室作る.

triangleQ[x_List]:=(Abs[Nlength[x[[1]],x[[2]]]^2+2Sqrt[5]-6]<0.0001&&Abs[Nlength[x[[2]],x[[3]]]^2+2Sqrt[5]-6]<0.0001&&Abs[Nlength[x[[3]],x[[1]]]^2+2Sqrt[5]-6]<0.0001);

icosahedron[x_List]:=(
north=x[[1]];
grp1=x[[{2,3,4,5,6,2}]];
grp2=x[[{7,8,9,10,11,7}]];
south=x[[12]];
paired1=Partition[grp1,2,1];
paired2=Partition[grp2,2,1];
triangles1= Map[Prepend[#,north]&,paired1];
triangles4= Map[Prepend[#,south]&,paired2];
triangles2=Table[Append[paired1[[i]],grp2[[i+1]]],{i,1,5}];
triangles3=Table[Append[paired2[[i]],grp1[[i]]],{i,1,5}];
Return[Flatten[{triangles1,triangles2,triangles3,triangles4},1]])


triangles24=Map[icosahedron,cells24];
{t1111,t1112,t1121,t1122,t1211,t1212,t1221,t1222,
t2111,t2112,t2121,t2122,t2211,t2212,t2221,t2222,t1000,t0100,t0010,t0001,t2000,t0200,t0020,t0002}=triangles24;(*definition of t1111...t0004*)

tetrasA0:=Table[Map[Append[#,futa[[k]]]&,triangles24[[k]]],{k,1,24}];
tetrasA:=Flatten[tetrasA0,1];

3. [TetrasBの作成]  正20面体の表面の三角形と，他の正20面体の頂点から作られる正四面体胞tetrasB (96個)を作成．

[方法1] 正攻法.線分PQを20面体Vの△ABCに一番近い正8面体上線分とする.　PQを共有するV以外の正8面体V1,V2は，PQから最も遠い点Rを共有する．最後にPRをτ:1に内分する点が△ABCに対応する正四面体の頂点.

Clear[x,y,z];
findOcta[x_,y_,z_]:=Intersection[Select[octas,MemberQ[#,x]& ],Select[octas,MemberQ[#,y]&],Select[octas,MemberQ[#,z]&]]
findAGpOutofseg[{{index11_,index12_}},{{index21_,index22_}}]:=Complement[{1,2,3},{index11},{index21}][[1]];
change1And2[{{index1_,index2_}}]={{index1,3-index2}};


findTetra[closest_,othersideofseg_,currentOcta_]:=(
OctaParted:=Partition[currentOcta,2];
cgp:=findAGpOutofseg[Position[OctaParted,closest],Position[OctaParted,othersideofseg]];
outofseg1=Extract[OctaParted,{cgp,1}];
outofseg2=Extract[OctaParted,{cgp,2}];
nearbyOcta:=Complement[findOcta[closest,othersideofseg,outofseg1],{currentOcta}][[1]];
nearbyPaired:=Partition[nearbyOcta,2];
oppositePtPos:=Position[nearbyPaired,outofseg1]//change1And2;
oppositeOfoutofseg=Extract[nearbyPaired,oppositePtPos][[1]];
V4:=naibun[oppositeOfoutofseg,closest];
V1:=naibun[outofseg1,closest];
V2:=naibun[outofseg2,closest];
V3:=naibun[closest,othersideofseg];
t*{V1,V2,V3,V4}//FullSimplify
)


tetrasB0=Table[{
findTetra[octas[[k,1]],octas[[k,5]],octas[[k]]],
findTetra[octas[[k,1]],octas[[k,6]],octas[[k]]],
findTetra[octas[[k,2]],octas[[k,5]],octas[[k]]],
findTetra[octas[[k,2]],octas[[k,6]],octas[[k]]],
findTetra[octas[[k,3]],octas[[k,1]],octas[[k]]],
findTetra[octas[[k,3]],octas[[k,2]],octas[[k]]],
findTetra[octas[[k,4]],octas[[k,1]],octas[[k]]],
findTetra[octas[[k,4]],octas[[k,2]],octas[[k]]],
findTetra[octas[[k,5]],octas[[k,3]],octas[[k]]],
findTetra[octas[[k,5]],octas[[k,4]],octas[[k]]],
findTetra[octas[[k,6]],octas[[k,3]],octas[[k]]],
findTetra[octas[[k,6]],octas[[k,4]],octas[[k]]]},
{k,1,24}];

tetrasB1=Flatten[tetrasB0,1];
tetrasB=Union[Table[Union[tetrasB1[[k]]],{k,1,Length[tetrasB1]}]];

4. [TetrasCの作成] 正20面体A,B,C,D,E,Fからそれぞれ1つの辺(2点)をえらんで作られる正四面体胞tetrasC(24個)を作成.

(*vertices of Octas(24 cells)*)
VofOctas=Module[{base,p,q},
base={{p,q,0,0},{p,0,q,0},{p,0,0,q},{0,p,q,0},{0,p,0,q},{0,0,p,q}};
t*Flatten[{base/.{p->1,q->1},base/.{p->1,q->-1},base/.{p->-1,q->1},base/.{p->-1,q->-1}},1]];
(*length between each vertice of tetrahedron and each vertice of 24 cells is sqrt(3-sqrt(5)).*)
tetrasC=Table[Select[vecs24,Abs[Nlength[#,VofOctas[[k]]]^2-3+Sqrt[5]]<0.001&],{k,1,24}];
%/.myrule;

5. 投影図.

skelton[x_List]:={Line[{x[[2]],x[[3]],x[[4]],x[[2]]}],Line[{x[[1]],x[[2]]}],Line[{x[[1]],x[[3]]}],Line[{x[[1]],x[[4]]}]};
body[x_List]:={Polygon[{x[[2]],x[[3]],x[[4]]}],Polygon[{x[[1]],x[[2]],x[[3]]}],
Polygon[{x[[1]],x[[2]],x[[4]]}],Polygon[{x[[1]],x[[3]],x[[4]]}]};

tetras=Join[tetrasA,tetrasB,tetrasC];

Proj[]の第1引数は mode. mode=1はu=0への直投影, mode=2は xyu空間への直投影, mode=3は中心投影.
即ち 「(τ+ushift)だけ u軸方向に平行移動した600胞体」を u=1の上に Oから中心投影している.
Skeltonは骨組みのみ，Bodyは胞も表示. kは濃い色で塗られる 胞の番号. 
Skelton は opacity ≥ 0.5で. BodyはOpacity =   0.1 ぐらいが面白い．

Manipulate[
mode:=Switch[prjmode,"中心投影",3,"直投影//u",1,"直投影//z",2];
tetraP:=Map[proj[mode,t+ushift,#]&,tetras];
Skelton:=Map[skelton,tetraP];
Body:=Map[body,tetraP];
If[style=="Body",Graphics3D[{Opacity[0.7],Body,Body[[k]]},Boxed->False],
Graphics3D[{Opacity[0.7],Skelton,Body[[k]]},Boxed->False]]
,
(*controllers*)
{{prjmode,"中心投影",Style["mode",Blue,Bold]},{"中心投影","直投影 //u","直投影 //z"}},
{{style,"Body",Style["style",Blue,Bold]},{"Body","Skelton"}},
{{ushift,0.5,Style["ushift",Blue,Bold]},0.1,2},
{{k,1,Style["cell#",Blue,Bold]},1,600}
,ControlPlacement->Left
]