CoCalc Shared FilesInellipse of convex quadrilateral.sagews.html
Author: Dominique Laurain
Views : 10

Author Dominique Laurain
Date 2016-12-03T07:56:51

Inellipse

The word “inellipse” is used in place of “ellipse inscribed”.

In the worsheet we simply compute the drawing of a quadrilateral inellipse.

Some of the key points in Horvitz’s paper :

• using Marden’s theorem

• center of quadrilateral inellipse is lying in the segment linking midpoints of diagonals

References

#%hideimport syssys.path.append('./modules')import DrawingConstants as dc%load ./modules/RationalTrigonometry.sage%load ./modules/EuclideanGeometry.sage%load ./modules/InversiveGeometry.sage%load ./modules/PinwheelGeometry.sage%load ./modules/TrilinearCoordinates.sage%load ./modules/MatrixRepresentationConics.sage# Alias from dc module[nd] = [dc.nd][c0,c1,c2,c3,c4,c5,c6,c7,c8] = [dc.c0,dc.c1,dc.c2,dc.c3,dc.c4,dc.c5,dc.c6,dc.c7,dc.c8][ls] = [dc.ls][th1,th2,th3,th4] = [dc.th1,dc.th2,dc.th3,dc.th4][fs1,fs2,fs3,fs4,fs5] = [dc.fs1,dc.fs2,dc.fs3,dc.fs4,dc.fs5]
# Drawing using Horwitz notation# Initialize graphicsg = Graphics()# Parameterss = 4t = 2# Set convex quadrilateral ABDCA = [0,0] ; B = [1,0] ; C = [0,1] ; D = [s,t]pA = vector(A) ; tA = text("$A$",pA,fontsize=fs2,color=c1); g += tApB = vector(B) ; tB = text("$B$",pB,fontsize=fs2,color=c1); g += tBpC = vector(C) ; tC = text("$C$",pC,fontsize=fs2,color=c1); g += tCpD = vector(D) ; tD = text("$D$",pD,fontsize=fs2,color=c1); g += tDg += line([pA,pB,pD,pC,pA],color=c1,thickness=th2)# Intersecting L2 and L3 : x = 0 and y = -t/(s-1)E = [0,-t/(s-1)]pE = vector(E) ; tE = text("$E$",pE,fontsize=fs2,color=c1); g += tEg += line([pA,pE,pB],color=c1,thickness=th1,linestyle=ls)# Set midpoints of diagonalsM1 = [1/2,1/2] ; M2 = [1/2*s,1/2*t]pM1 = vector(M1) ; tM1 = text("$M_1$",pM1,fontsize=fs2,color=c1); g += tM1pM2 = vector(M2) ; tM2 = text("$M_2$",pM2,fontsize=fs2,color=c1); g += tM2g += line([pA,pD],color=c1,thickness=th1,linestyle=ls)g += line([pB,pC],color=c1,thickness=th1,linestyle=ls)g += line([pM1,pM2],color=c1,thickness=th1)  # segment Z, joining diagonals midpoints# equation of line M1M2 : y = L(x) = 1/2*(s - t + 2*x(t - 1))/(s - 1)  with x in [1/2,1/2*s]x = var('x')L(x) =  1/2*(s - t + 2*x*(t - 1))/(s - 1)# Set center of the inellipseh = 4/3 ; k = L(h) #; print "L(h) = k = ",k  # 7/9G = [h,k]pG = vector(G) ; tG = text("$G$",pG,fontsize=fs2,color=c1); g += tG# From h,k deducing other parameters for ellipse and pointst1 = 2*h - 1 - 2*k*(s-1)/t ; t2 = 1 - 2*h ; t3 = 1 - t1 - t2z1 = 0 ; z2 = 1 ; z3 = -t/(s - 1)*Ie1 = (t2*z3 + t3*z2)/(t2 + t3) ; e2 = (t3*z1 + t1*z3)/(t3 + t1) ; e3 = (t1*z2 + t2*z1)/(t1 + t2)E1 = [e1.real(),e1.imag()]; pE1 = vector(E1) ; tE1 = text("$E_1$",pE1,fontsize=fs2,color=c1); g += tE1E2 = [e2.real(),e2.imag()]; pE2 = vector(E2) ; tE2 = text("$E_2$",pE2,fontsize=fs2,color=c1); g += tE2E3 = [e3.real(),e3.imag()]; pE3 = vector(E3) ; tE3 = text("$E_3$",pE3,fontsize=fs2,color=c1); g += tE3# Drawing symetric points if we want five points on the ellipse, without computing fociipE4 = Symetric(pG,pE1); tE4 = text("$E_4$",pE4,fontsize=fs2,color=c1); g += tE4pE5 = Symetric(pG,pE2); tE5 = text("$E_5$",pE5,fontsize=fs2,color=c1); g += tE5# Compute fociiw = var('w')focii(w) = ((s - 1)*w^2 + (I*t*(t1 + t2) + (s - 1)*(t2 - 1))*w - I*t1*t).factor() #; print "focii(w) = ",focii(w)res = solve([focii == 0],w) #; print "res = ",res# Converting to numeric values, roots of equationw1 = res[0].rhs().n() #; print "w1 = ",w1w2 = res[1].rhs().n() #; print "w2 = ",w2W1 = [w1.real(),w1.imag()]; pW1 = vector(W1) #; tW1 = text("$W_1$",pW1,fontsize=fs2,color=c1); g += tW1W2 = [w2.real(),w2.imag()]; pW2 = vector(W2) #; tW2 = text("$W_2$",pW2,fontsize=fs2,color=c1); g += tW2# Draw the ellipsis after major-axis and minor-axis lengths computationf = RT_Distance(pG,pW1).n()p = ((RT_Distance(pW1,pE1) + RT_Distance(pW2,pE1))/2).n()  # choosing any point (here E1) on ellipseq = sqrt(p^2 - f^2).n()if pW1[0] != pW2[0] :    theta = atan((pW1[1]-W2[1])/(pW1[0]-pW2[0]))else:    theta = pi/2g += ellipse(center=pG,r1=p,r2=q,angle=theta,color=c2)# Check formula for squared area of ellipseSquaredArea_Ellipse = (pi/(2*(s-1)))^2*(2*h - 1)*(s + 2*h*(t - 1))*(s - 2*h)print "SquaredArea_Ellipse = ",SquaredArea_Ellipse.n()," vs (pi*p*q)^2 = ",((pi*p*q)^2).n()# Check formula area for in ellipse triangle CDE using areal coordinates of center G of ellipseArea_CDE = TriangleArea(pC,pD,pE)Area_CDG = TriangleArea(pC,pD,pG) ; k = Area_CDG / Area_CDEArea_DEG = TriangleArea(pD,pE,pG) ; l = Area_DEG / Area_CDEArea_ECG = TriangleArea(pE,pC,pG) ; m = Area_ECG / Area_CDESquaredArea = (pi*Area_CDE)^2*(1 - 2*k)*(1 - 2*l)*(1 - 2*m)print "SquaredArea = ",SquaredArea.n()# Computing maximum for squared areahh = var('hh')f0(hh) = (pi/(2*(s-1)))^2*(2*hh - 1)*(s + 2*hh*(t - 1))*(s - 2*hh); print "Solving max of f0(hh) = ",f0print "f0(1/2) = ",f0(1/2).n()," vs f0(s/2) = ",f0(s/2).n()f1(hh) = f0.derivative(hh).factor(); print "f1 = ",f1res = solve([f1],hh)hh0 = res[0].rhs().n(); print "root for extrema hh0 = ",hh0 # hh0 = 4/3hh1 = res[1].rhs().n(); print "root for extrema hh1 = ",hh1 # hh1 = -1print "f0(hh0) = ",f0(hh0).n()# Show graphicsshow(g,axes=False,aspect_ratio=True)#save(g,'QuadrilateralInellipse.png',axes=False,aspect_ratio=True)