*)
Keep[x_]:= PutAppend[Definition[x],"more.m"];
SetAttributes[Keep,HoldFirst];
Fix[x_] := Edit[Definition[x]]
SetAttributes[Fix,HoldFirst];
(Unprotect[In,Out]; Format[In]=SphereIn;
Format[Out]=SphereOut; Protect[In,Out];)
$SpherePrecision:= 16;
Ns[x_]:= N[x,$SpherePrecision];
zeta := 1/(2*ArcTan[2^(1/2)/5])//Ns;
pt := (-Pi/3 + 4*ArcTan[2^(1/2)/5])//Ns;
doct := ((Pi - 4*ArcTan[2^(1/2)/5])/(2*2^(1/2)))//Ns;
chi = -(x1*x4^2) + x1*x4*x5 + x2*x4*x5 - x2*x5^2 + x1*x4*x6 + x3*x4*x6 +
x2*x5*x6 + x3*x5*x6 - 2*x4*x5*x6 - x3*x6^2
eta[x_,y_,z_]:= Module[{s= (x+y+z)/2}, x y z/4/Sqrt[s (s-x)(s-y)(s-z)]]//Ns;
volR[a_, b_, c_] := Sqrt[a^2*(b^2 - a^2)*(c^2 - b^2)]/6//Ns;
solR[x_, y_, z_] := 2*ArcTan[Sqrt[((z - y)*(y - x))/((z + y)*(x + y))]]//Ns;
dihR[x_, y_, z_] := ArcTan[Sqrt[(z^2 - y^2)/(y^2 - x^2)]]//Ns;
vorR[a_, b_, c_] := 4*(-(doct*volR[a, b, c]) + solR[a, b, c]/3)//Ns;
denR[a_, b_, c_] := solR[a, b, c]/(3*volR[a, b, c])
tauR[x__] := (-vorR[x] + solR[x]*zeta*pt)//Ns;
ttauR[x__] := (vorR[x] - 3.2*solR[x]*zeta*pt)//Ns;
TRUNC:=251/200//Ns;
Dihedral2[y1_,y2_,y3_,y4_,y5_,y6_] := (* Dihedral along edge 2 *)
Dihedral[y2,y1,y3,y5,y4,y6];
Dihedral3[y1_,y2_,y3_,y4_,y5_,y6_] := (* Dihedral along edge 3 *)
Dihedral[y3,y1,y2,y6,y4,y5];
Chi[x__] := Module[{x1, x2, x3, x4, x5, x6},
{x1, x2, x3, x4, x5, x6} = {x}^2;
-(x1*x4^2) + x1*x4*x5 + x2*x4*x5 - x2*x5^2 + x1*x4*x6 + x3*x4*x6 +
x2*x5*x6 + x3*x5*x6 - 2*x4*x5*x6 - x3*x6^2];
tauAnalytic[x__] := Solid[x]*zeta*pt - VorAnalytic[x]
ttauAnalytic[x__]:= VorAnalytic[x] - 3.2 Solid[x] zeta pt;
tau[x__] := Solid[x]*zeta pt - Gamma[x];
ttau[x__] := Gamma[x] - 3.2*Solid[x]*zeta pt;
VorAnalytic[x__]:= Module[{y1,y2,y3,y4,y5,y6,x1,x2,x3,x4,x5,x6,
del,u126,u135,u234,vol},
{y1,y2,y3,y4,y5,y6} = {x};
{x1,x2,x3,x4,x5,x6} = {x}^2;
del = Sqrt[Delta[x]];
u126 = U[x1,x2,x6];
u135 = U[x1,x3,x5];
u234 = U[x2,x3,x4];
vol=1/(48 del) (
(x1 (x2+x6-x1)+ x2 (x1+x6-x2)) Chi[y4,y5,y3,y1,y2,y6]/u126 +
(x2 (x3+x4-x2) + x3 (-x3+x4+x2)) Chi[y6,y5,y1,y3,y2,y4]/u234 +
(x1 (-x1+x3+x5) + x3 (x1-x3+x5)) Chi[y4,y6,y2,y1,y3,y5]/u135 );
4(-doct vol + Solid[x]/3)
];
Delta[y1_,y2_,y3_,y4_,y5_,y6_]:= Module[{x1,x2,x3,x4,x5,x6},
x1= y1*y1; x2=y2*y2; x3=y3*y3;
x4= y4*y4; x5=y5*y5; x6=y6*y6;
(x1*x4*(-x1+x2+x3-x4+x5+x6)+
x2*x5*(x1-x2+x3+x4-x5+x6)
+x3*x6*(x1+x2-x3+x4+x5-x6)
-x2*x3*x4-x1*x3*x5-x1*x2*x6-x4*x5*x6)];
U[x1_,x2_,x6_]:= -x1*x1-x2*x2-x6*x6+2*x1*x6+2*x1*x2+2*x2*x6;
Rho[y1_,y2_,y3_,y4_,y5_,y6_]:= Module[{x1,x2,x3,x4,x5,x6},
x1= y1*y1; x2=y2*y2; x3=y3*y3;
x4= y4*y4; x5=y5*y5; x6=y6*y6;
(-x1*x1*x4*x4-x2*x2*x5*x5-x3*x3*x6*x6+
2*x1*x2*x4*x5+2*x1*x3*x4*x6+2*x2*x3*x5*x6)];
Rad[y1_,y2_,y3_,y4_,y5_,y6_]:=
Sqrt[Rho[y1,y2,y3,y4,y5,y6]/Delta[y1,y2,y3,y4,y5,y6]]/2;
aSolid[y1_,y2_,y3_,y4_,y5_,y6_]:=
y1*y2*y3 + y1*(y2*y2+y3*y3-y4*y4)/2 +
y2*(y1*y1+y3*y3-y5*y5)/2 +
y3*(y1*y1+y2*y2-y6*y6)/2;
Unprotect[Gamma];
Print[" -- Gamma redefined as compression --"];
Attributes[Gamma]= {};
Gamma::usage="Gamma[x] is the compression of the simplex with edge lengths x";
Clear[Gamma];
Gamma[y1_,y2_,y3_,y4_,y5_,y6_]:= Module[{u,t,a1,a2,a3,a4},
u=Delta[y1,y2,y3,y4,y5,y6];
If[u<=0,Return[0]];
t = Sqrt[u]/2;
a1 = aSolid[y1,y2,y3,y4,y5,y6];
a2 = aSolid[y1,y5,y6,y4,y2,y3];
a3 = aSolid[y2,y4,y6,y5,y1,y3];
a4 = aSolid[y4,y5,y3,y1,y2,y6];
(-doct*t/6 + (2/3)*(ArcTan[t/a1]+ArcTan[t/a2]+ArcTan[t/a3]+ArcTan[t/a4]))
];
Protect[Gamma];
Solid[y1_,y2_,y3_,y4_,y5_,y6_]:= Module[{u,t,a1},
u = Delta[y1,y2,y3,y4,y5,y6];
If [N[u]<=0, Return[0]];
t = Sqrt[u]/2;
a1 = aSolid[y1,y2,y3,y4,y5,y6];
2*ArcTan[t/a1]
]//Ns;
Dihedral[y1_,y2_,y3_,y4_,y5_,y6_]:= Module[{ux,x1,x2,x3,x4,x5,x6,t},
x1= y1*y1; x2=y2*y2; x3=y3*y3;
x4= y4*y4; x5=y5*y5; x6=y6*y6;
ux = U[x1,x2,x6]*U[x1,x3,x5];
t=(-2*x1*x4+x1*(x2+x3+x5+x6)+x2*x5+x3*x6-x1*x1 -x2*x3-x5*x6)/Sqrt[ux];
ArcCos[t]
];
Density[x_]:= doct/(1-3 x/(16 Pi));
Norm[x_]:= x.x;
Norm[x_,y_]:= Norm[x-y];
Distance[x_,y_]:= Sqrt[(x-y).(x-y)];
FarFrom[x_,pair_]:= If[Norm[x,pair[[1]]]N[b],Return[0]];
If[N[b]>N[c],Return[0]];
-(a^2 + a*c - 2*c^2)*(c - a)*ArcTan[u]/6 + a*(b^2 - a^2)*u/6 -
(2/3)*c^3*ArcTan[(u*(b - a))/(b + c)]]//Ns;
Qn[y1_,y2_,z_]:= -4 doct (Quoin[y1/2,eta[y1,y2,z],TRUNC] +
Quoin[y2/2,eta[y1,y2,z],TRUNC])//Ns;
(* VorVc is the truncation vor(S,1.255) of the Voronoi cell at t0=1.255 *)
VorVc[y1_,y2_,y3_,y4_,y5_,y6_]:= Module[{h1,h2,h3,a1,a2,a3},
h1 = y1/2; h2 = y2/2; h3 = y3/2;
a1 = Dihedral[y1,y2,y3,y4,y5,y6];
a2 = Dihedral[y2,y3,y1,y5,y6,y4];
a3 = Dihedral[y3,y1,y2,y6,y4,y5];
vorstar[h1,h2,h3,a1,a2,a3]+
Qn[y1,y2,y6]+Qn[y2,y3,y4]+Qn[y1,y3,y5]
];
phi[h_,t_]:= 2 (2 - doct h t (h+t))/3//Ns;
phi0 := phi[TRUNC,TRUNC]//Ns;
(* two function from formulation *)
crown[h_]:= Module[{e}, e = eta0[h];
2Pi (1-h/e) (phi[h,e] - phi0)
];
eta0[h_]:= eta[ 2h,2 TRUNC,2];
anc[y1_,y2_,y6_]:= Module[{h1,h2,b,c}, (* Equation F.4.5 *)
h1 = y1/2;
h2 = y2/2;
b = eta[y1,y2,y6];
c = eta0[h1];
If[N[b]>N[c],Return[0]];
-dihR[h1,b,c] crown[h1]/(2Pi) -solR[h1,b,c] phi0 + vorR[h1,b,c] -
dihR[h2,b,c](1-h2/TRUNC)(phi[h2,TRUNC]-phi0)-
solR[h2,b,c] phi0 + vorR[h2,b,c]
];
(* from SP IV, section 2.5. *)
kappa[y1_,y2_,y3_,y4_,y5_,y6_]:= crown[y1/2] Dihedral[y1,y2,y3,y4,y5,y6]/(2Pi)+
anc[y1,y2,y6]+anc[y1,y3,y5];
vorstar[h1_,h2_,h3_,a1_,a2_,a3_]:= Module[{Mx},
Mx[h_]:= If[N[h]a[[j]],{j,1,Length[a]}];
{f0,Table[D[fs,x[i]]/.sub,{i,1,Length[var]}],
Table[D[fs,x[i],x[j]]/.sub,{i,1,Length[var]},{j,1,Length[var]}]}
];
Print[" -- Sphere Packings initialized -- "];
(*
This page is available for historical purposes only. It is a copy
from www.math.lsa.umich.edu/~hales/countdown. It has not been
maintained since 1998.