Poncelet's Porism

Author

Eric W. Weisstein
December 21, 2004

This notebook downloaded from http://mathworld.wolfram.com/notebooks/PlaneGeometry/PonceletsPorism.nb.

For more information, see Eric's MathWorld entry http://mathworld.wolfram.com/PonceletsPorism.html.

©2005 Wolfram Research, Inc. except for portions noted otherwise

Code

(*Poncelet*)(*PonceletClosed[3,d_,r1_:1]:=(r1^2-d^2)/(2r1) PonceletClosed[4,d_,r1_:1]:=(r1^2-d^2)/Sqrt[2(d^2+r1^2)]*)PonceletClosed[n_Integer/;n>2,d_,R_:1.]:=Module[{r,a=1/(R+d),b=1/(R-d),c,lambda,omega,m},c=1/r;
lambda=1+2c^2(a^2-b^2)/a^2/(b^2-c^2);
omega=ArcCosh[lambda];
m=1-Exp[-2omega];
r/.FindRoot[JacobiSC[EllipticK[m]/n,m]==(c Sqrt[b^2-a^2]+b Sqrt[c^2-a^2])/a/(c+b),{r,{.1,.2}}]]

PonceletClosedTransversePlot[n_,x_,t_:0]:=PonceletTransversePlot[n,x,PonceletClosed[n,x],t]

PonceletDiagonalsPlot[n_?EvenQ,d_,t_:0]:=Module[{r=PonceletClosed[n,d],v,lim},lim=LimitingPoints[d,r,1][[-1]];
v=Transpose[Partition[PonceletTransverse[n,d,r,t],n/2]];
{{RGBColor[0,0,1],Line/@v,PointSize[.05],Point[{lim,0}]},PonceletTransversePlot[n,d,r,t]}]

PonceletDiagonalsPlot[n_?OddQ,d_,t_:0]:=Module[{r=PonceletClosed[n,d],pts,v,lim},pts=PonceletTangents[n,d,r,t];
lim=LimitingPoints[d,r,1][[-1]];
v=Transpose[{RotateLeft[#[[1]],(n+1)/2],#[[2]]}&[pts]];
{{RGBColor[0,0,1],Line/@v,PointSize[.05],Point[{lim,0}]},Circle[{0,0},1],Circle[{d,0},r],{RGBColor[1,0,0],ClosedLine[pts[[2]]],PointSize[.05],Point[{Cos[#],Sin[#]}&[t]]}}]

(*rewrite*)
PonceletTangentAngles[n_,x_,r_,t_]:=Module[{p,pl={},pp,tl={t},i,j,tsoln},Off[Solve::ifun];
Do[p=(pp/.Solve[(r+x Cos[pp]-Cos[pp-tl[[-1]]])==0,pp]);
tsoln=Apply[ArcTan,Flatten[Table[Intersections[{Line[{Cos[tl[[-1]]],Sin[tl[[-1]]]}],Circle[x+r Cos[p[[i]]],r Sin[p[[i]]]]},1],{i,2}],1],1];
(*make sure all angles are positive so that ordering will work*)p=Join[p,p+2Pi,p+4Pi];
AppendTo[pl,If[j==1,Min[Select[p,#-t>.01&]],(*Else*)(*Print["p:",j,p];*)pl[[-1]]+Min[Select[p-pl[[-1]],#>.01&]]]];
(*Print["p:",j,pl];*)tsoln=Join[tsoln,tsoln+2Pi,tsoln+4Pi];
AppendTo[tl,If[j==1,Min[Select[tsoln,#-t>.01&]],(*Else*)(*Print["t:",j,tsoln];*)tl[[-1]]+Min[Select[tsoln-tl[[-1]],#>.01&]]]];,{j,n}];
{pl,Rest[tl]}]

PonceletTangents[n_,x_,r_,t_]:=Module[{a=PonceletTangentAngles[n,x,r,t]},{({x,0}+r{Cos[#],Sin[#]})&/@a[[1]],{Cos[#],Sin[#]}&/@a[[2]]}]

PonceletTransverseAngles[n_,x_,r_,t_]:=Module[{p,pp,tl={t},i,j,tsoln},Off[Solve::ifun];
Do[p=(pp/.Solve[(r+x Cos[pp]-Cos[pp-tl[[-1]]])==0,pp]);
tsoln=Apply[ArcTan,Flatten[Table[Intersections[Line[{Cos[tl[[-1]]],Sin[tl[[-1]]]}],Circle[{x+r Cos[p[[i]]],r Sin[p[[i]]]},1]],{i,2}],1],1];
Print[tsoln];
(*make sure all angles are positive so that ordering will work*)tsoln=If[#≥0,#,2Pi+#]&/@tsoln;
AppendTo[tl,If[j==1,Select[tsoln,Abs[#-tl[[-1]]]>.01&][[2]],(*Else*)Select[tsoln,Abs[#-tl[[-1]]]>.01&&Abs[#-tl[[-2]]]>.01&][[1]]]],{j,n}];
tl]

PonceletTransverse[n_,x_,d_,t_]:={Cos[#],Sin[#]}&/@PonceletTransverseAngles[n,x,d,t]

PonceletClosed::usage="PonceletClosed[n,d,R] gives the radius of the smaller of two circles offset by a distance d from the larger (with default radius R=1) so that the bicentric n-gon will close."

PonceletClosedTransversePlot::usage="PonceletClosedTransversePlot[n,d,theta] gives a plot of Poncelet's porism for a closed n-gon with offset d and starting angle theta."

PonceletDiagonalsPlot::usage="PonceletDiagonalsPlot[n,d,<ang>] plots Poncelet's closed transversal for an offset d and offset angle ang, together with diagonals."

PonceletTangent::usage="PonceletTangent[n,d,R,t] returns a list of tangent points of the first n Poncelet transverses of two circles whose centers are offset by d and the larger of which has radius R."

PonceletTangentAngles::usage="PonceletTangentAngles[n,d,R,t] ???."

PonceletTangents::usage="PonceletTangents[n,d,R,t] ???."

PonceletTransverse::usage="PonceletTransverse[n,d,R,t] returns a list of the first n Poncelet transverses of two circles, one a unit circle centered at (0,0), and the other a circle of radius R centered at (d,0).  The angle t gives the starting point around the unit circle."

PonceletTransverseLines::usage="PonceletTransverseLines [n,d,R,t] gives the lines of the Poncelet transverse."

Example

Show[Graphics[PonceletClosedTransversePlot[4,.2,.5],AspectRatio->Automatic]]

Solution

Solve[r1^2-2r1 r2==d^2,r2]

PowerExpand[FullSimplify[Solve[1/(r1-d)^2+1/(r1+d)^2==1/r2^2,r2]]]

?PonceletClosedTransversePlot

3 snapshots for 3- and 4-gons

Show[GraphicsArray[{{
Show[PonceletClosedTransversePlot[3,.3,0.],AspectRatio->Automatic,DisplayFunction->Identity] ,
Show[PonceletClosedTransversePlot[3,.3,1.1],AspectRatio->Automatic,DisplayFunction->Identity] ,
Show[PonceletClosedTransversePlot[3,.3,2.9],AspectRatio->Automatic,DisplayFunction->Identity] 
},
{
Show[PonceletClosedTransversePlot[4,.3,0.],AspectRatio->Automatic,DisplayFunction->Identity] ,
Show[PonceletClosedTransversePlot[4,.3,1.5],AspectRatio->Automatic,DisplayFunction->Identity] ,
Show[PonceletClosedTransversePlot[4,.3,1.8],
AspectRatio->Automatic,DisplayFunction->Identity] 
}}]]

Inverting about limit points

Show[GraphicsArray[{
Graphics[{
PonceletDiagonalsPlot[3,.3,0],
PointSize[.05],
RGBColor[0,1,0],Point[{LimitingPoints[.3,PonceletClosed[3,.3],1][[1]],0}]
},AspectRatio->Automatic],
Graphics[Invert[{{LimitingPoints[.3,PonceletClosed[3,.3],1][[1]],0},1},{PonceletDiagonalsPlot[3,.3,0],
PointSize[.05],
RGBColor[0,1,0],Point[{LimitingPoints[.3,PonceletClosed[3,.3],1][[1]],0}]}
],AspectRatio->Automatic],
Graphics[Invert[{{LimitingPoints[.3,PonceletClosed[3,.3],1][[2]],0},1},{PonceletDiagonalsPlot[3,.3,0],
PointSize[.05],
RGBColor[0,1,0],Point[{LimitingPoints[.3,PonceletClosed[3,.3],1][[1]],0}]}
],AspectRatio->Automatic]
}]]

Animation

Export[gifs<>"poncelet.gif",
Table[
Show[GraphicsArray[{
Show[Graphics[PonceletClosedTransversePlot[3,.2,t]],AspectRatio->Automatic,DisplayFunction->Identity] ,
Show[Graphics[PonceletClosedTransversePlot[4,.2,t]],AspectRatio->Automatic,DisplayFunction->Identity ],
Show[Graphics[PonceletClosedTransversePlot[5,.2,t]],AspectRatio->Automatic,DisplayFunction->Identity ]
}]],
{t,0.,2Pi-2Pi/30,2Pi/30}]
]

Diagonals animation

Export[gifs<>"ponclim.gif",
Table[Show[GraphicsArray[{
Graphics[{PonceletDiagonalsPlot[3,.2,t]},AspectRatio->Automatic],
Graphics[{PonceletDiagonalsPlot[4,.2,t]},AspectRatio->Automatic],
Graphics[{PonceletDiagonalsPlot[5,.2,t]},AspectRatio->Automatic],
Graphics[{PonceletDiagonalsPlot[6,.2,t]},AspectRatio->Automatic]
}]],{t,0.,2Pi-2Pi/20,2Pi/20}],
ImageSize->300,
ConversionOptions->{"AnimationDisplayTime"->0.1}
]

Special Cases

From Kerawala

a=1/(R+d);
b=1/(R-d);
c=1/r;
e0=a+b+c;
e1=-a+b+c;
e2=a-b+c;
e3=a+b-c;
E1=-a^2+b^2+c^2;
E2=a^2-b^2+c^2;
E3=a^2+b^2-c^2;
F1=-E2 E3+E3 E1+E1 E2;
F2= E2 E3-E3 E1+E1 E2;
F3= E2 E3+E3 E1-E1 E2;
F0= E2 E3+E3 E1+E1 E2 (* e0e1e2e3*);
g0= E1 E2 E3+2a b E1 E2+2b c E2 E3+2c a E3 E1;
g1= E1 E2 E3-2a b E1 E2+2b c E2 E3-2c a E3 E1;
g2= E1 E2 E3-2a b E1 E2-2b c E2 E3+2c a E3 E1;
g3= E1 E2 E3+2a b E1 E2-2b c E2 E3-2c a E3 E1;

n=3

Together[a+b-c]

Expand[(Sqrt[R-d-r]+Sqrt[R+d-r])^2-Sqrt[2R]^2]

Expand[(-d-r+R)(d-r+R)-r^2]

n=4

FullSimplify[Expand[(a+b)(b+c)(c+a)-a^3-b^3-c^3]]

FullSimplify[Expand[(a+b+c)^3-4(a^3+b^3+c^3)]]

FullSimplify[Expand[(-a+b+c)(a-b+c)(a+b-c)+4a b c]]

FullSimplify[Expand[Det[{{e0,e3,e2},{e3,e0,e1},{e2,e1,e0}}]]]

det[l,{x1,x2,x3}]:=Module[{
			a=Coefficient[l,x1^2],
			b=Coefficient[l,x2^2],
			c=Coefficient[l,x3^2],
			f=Coefficient[l,x2 x3]/2,
			g=Coefficient[l,x1 x3]/2,
			h=Coefficient[l,x1 x2]/2
			},
		Det[{{a,h,g},{h,b,f},{g,f,c}}]
]

dd=x1^2+x2^2-x3^2 R^2;
cc=x1^2+x2^2-2x1 x3 d+x3^2(d^2-r^2);
t dd+cc

Det[{{a,h,g},{h,b,f},{g,f,c}}]

Simplify[Det[{{a ,0 , g},{ 0,a ,f},{ g, f, c}}]]

x=Together[Expand[det[t  dd+cc,{x1,x2,x3}]/det[cc,{x1,x2,x3}]]]

Together[CoefficientList[Series[Sqrt[x],{t,0,4}],t]]

A[n_]:=CoefficientList[%454,t][[n+1]]

deta[n_Integer?OddQ]:=Module[{m=(n-1)/2,i},
	Det[Table[A[i+j],{i,m},{j,m}]]
]

deta[n_Integer?EvenQ]:=Module[{m=n/2,i},
		Det[Table[A[i+j+1],{i,m-1},{j,m-1}]]
]

deta[3]

R^2-x^2-2Rr

Simplify[deta[3]]

Expand[2r^2(R^2+d^2)-(R^2-d^2)^2]

n=6

eqn=Numerator[Factor[Simplify[1/E1+1/E2-1/E3]]]

eqn=-FullSimplify[-3 (d^8)+4 (d^6) (r^2)+12 (d^6) (R^2)-4 (d^4) (r^2) (R^2)+16 (d^2) (r^4) (R^2)-18 (d^4) (R^4)-4 (d^2) (r^2) (R^4)+12 (d^2) (R^6)+4 (r^2) (R^6)-3 (R^8)]

eqn2=(3(R^2-d^2)^4-4r^2(R^2+d^2)(R^2-d^2)^2-16r^4 d^2 R^2);

Expand[eqn-eqn2]

Weill Approximation

Solve[{d==2r/3,R==4r/3},{d,r}]

eqn2 /. {d->R/2,r->3/4R}

n=8

Kerawla

FullSimplify[Numerator[Factor[Simplify[1/E1^2+1/E2^2-1/E3^2]]]]

eqn1=Expand[((d^16)-8 (d^14) ((r^2)+(R^2))+(R^12) (8 (r^4)-8 (r^2) (R^2)+(R^4))+4 (d^12) (2 (r^4)+10 (r^2) (R^2)+7 (R^4))+2 (d^8) (R^2) (2 (r^2)+(R^2)) (-32 (r^4)-50 (r^2) (R^2)+35 (R^4))+8 (d^10) (6 (r^4) (R^2)-9 (r^2) (R^4)-7 (R^6))+4 (d^4) (R^6) (32 (r^6)-66 (r^4) (R^2)-18 (r^2) (R^4)+7 (R^6))+8 (d^6) (R^2) (16 (r^8)+16 (r^6) (R^2)+52 (r^4) (R^4)+5 (r^2) (R^6)-7 (R^8))-8 (d^2) (R^6) (-16 (r^8)+16 (r^6) (R^2)-6 (r^4) (R^4)-5 (r^2) (R^6)+(R^8)))]

Jacobi

p:=(R+d)/r
q:=(R-d)/r

eqn2=Factor[Numerator[Together[Expand[16p^4q^4(p^2-1)(q^2-1)-(p^2+q^2-p^2q^2)^4]]]]

eqn1+eqn2

16p^4q^4(p^2-1)(q^2-1)-(p^2+q^2-p^2q^2)^4

eqn3=Factor[Numerator[Together[16p^4q^4(p^2-1)(q^2-1)-(p^2+q^2-p^2q^2)^4]]]

Steiner (Wrong)

Expand[8 r^2((R^2-d^2)^2-r^2(R^2+d^2))((R^2+d^2)((R^2-d^2)^4+4r^4d^2R^2)-8r^2d^2R^2(R^2-d^2)^2)-((R^2-d^2)^4-4r^4d^2R^2)^2]

Expand[8 r^2((R^2-d^2)^2-r^2(R^2+d^2))((R^2+d^2)((R^2-d^2)^4+nb r^4d^2R^2)-cc r^2d^2R^2(R^2-d^2)^2)-((R^2-d^2)^4+dd r^4d^2R^2)^2]

Weill Approximation

Solve[{d==r/4,R==15r/4},{d,r}]

eqn3 /. {d->R/15,r->4R/15}

N[%]

n=10

Weill Approximation

Solve[{r==9d/4,R==d Sqrt[10]},{r,d}]

N[%]

n=14

Kerawala

g1==0

FullSimplify[Numerator[Together[Simplify[g1]]]]

Weill Approximation

α=((d+r)^2-R^2)/r^2;
β=((d-r)^2-R^2)/r^2;

either

α==(α-((α-1)/(β-1))^2)^2/(1-β((α-1)/(β-1))^2)^2

or

α==-1/beta

Exact Form

Poncelet[n_,{a_,b_,c_}]:=Module[{k2=1-Exp[-2ArcCosh[1+((2(c^2)((a^2)-(b^2)))/((a^2)((b^2)-(c^2))))]]},JacobiSC[((EllipticK[k2])/n),k2]-((c(Sqrt[(b^2)-(a^2)])+b(Sqrt[(c^2)-(a^2)]))/(a(b+c)))
]

Condition on a and b

0<a<b

Table[{a,b,Poncelet[3,{a,b,a+b}]}/.{a->Random[Real,{0,10},30],b->Random[Real,{0,10},30]}//Chop,{10}]

Table[Poncelet[3,{a,b,a+b}]/.{a->Random[Real,{0,5},50],b->Random[Real,{5,10},50]},{200}]//Chop//Union

n=3

Poncelet[3,{a,b,c}]

Poncelet[3,{a,b,a+b}]

Refine[FullSimplify[FunctionExpand[1-(E^(-2 ArcCosh[1+((2 (((a+b))^2) ((a^2)-(b^2)))/((a^2) ((b^2)-(((a+b))^2))))]))],0<a<b],0<a<b]

Poncelet[3,{a,b,c}]/.{a->1.2,b->2.2,c->3.4}//Chop

n=4

Poncelet[4,{a,b,c}]

Poncelet[4,{a,b,c}]/.{a->1.2,b->2.2,c->2.505992817228334`}//Chop