CoCalc Shared FilesInellipse of convex quadrilateral.sagews
Authors: Dominique Laurain, Harald Schilly
Views : 7
%md

### 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

- 2003 - Horvitz - The locus of centers of ellipses inscribed in quadrilaterals

- Worlfram - Area inellipse given areal coordinates of center - http://mathworld.wolfram.com/Inellipse.html

- Wolfram - Areal coordinates of point with respect to triangle - http://mathworld.wolfram.com/ArealCoordinates.html



### 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

#%hide
import sys
sys.path.append('./modules')
import DrawingConstants as dc

# 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 graphics
g = Graphics()

# Parameters
s = 4
t = 2

A = [0,0] ; B = [1,0] ; C = [0,1] ; D = [s,t]
pA = vector(A) ; tA = text("$A$",pA,fontsize=fs2,color=c1); g += tA
pB = vector(B) ; tB = text("$B$",pB,fontsize=fs2,color=c1); g += tB
pC = vector(C) ; tC = text("$C$",pC,fontsize=fs2,color=c1); g += tC
pD = vector(D) ; tD = text("$D$",pD,fontsize=fs2,color=c1); g += tD
g += 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 += tE
g += line([pA,pE,pB],color=c1,thickness=th1,linestyle=ls)

# Set midpoints of diagonals
M1 = [1/2,1/2] ; M2 = [1/2*s,1/2*t]
pM1 = vector(M1) ; tM1 = text("$M_1$",pM1,fontsize=fs2,color=c1); g += tM1
pM2 = vector(M2) ; tM2 = text("$M_2$",pM2,fontsize=fs2,color=c1); g += tM2
g += 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 inellipse
h = 4/3 ; k = L(h) #; print "L(h) = k = ",k  # 7/9
G = [h,k]
pG = vector(G) ; tG = text("$G$",pG,fontsize=fs2,color=c1); g += tG

# From h,k deducing other parameters for ellipse and points
t1 = 2*h - 1 - 2*k*(s-1)/t ; t2 = 1 - 2*h ; t3 = 1 - t1 - t2
z1 = 0 ; z2 = 1 ; z3 = -t/(s - 1)*I
e1 = (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 += tE1
E2 = [e2.real(),e2.imag()]; pE2 = vector(E2) ; tE2 = text("$E_2$",pE2,fontsize=fs2,color=c1); g += tE2
E3 = [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 focii
pE4 = Symetric(pG,pE1); tE4 = text("$E_4$",pE4,fontsize=fs2,color=c1); g += tE4
pE5 = Symetric(pG,pE2); tE5 = text("$E_5$",pE5,fontsize=fs2,color=c1); g += tE5

# Compute focii
w = 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 equation
w1 = res[0].rhs().n() #; print "w1 = ",w1
w2 = res[1].rhs().n() #; print "w2 = ",w2
W1 = [w1.real(),w1.imag()]; pW1 = vector(W1) #; tW1 = text("$W_1$",pW1,fontsize=fs2,color=c1); g += tW1
W2 = [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 computation
f = RT_Distance(pG,pW1).n()
p = ((RT_Distance(pW1,pE1) + RT_Distance(pW2,pE1))/2).n()  # choosing any point (here E1) on ellipse
q = sqrt(p^2 - f^2).n()
if pW1[0] != pW2[0] :
theta = atan((pW1[1]-W2[1])/(pW1[0]-pW2[0]))
else:
theta = pi/2
g += ellipse(center=pG,r1=p,r2=q,angle=theta,color=c2)

# Check formula for squared area of ellipse
SquaredArea_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 ellipse
Area_CDE = TriangleArea(pC,pD,pE)
Area_CDG = TriangleArea(pC,pD,pG) ; k = Area_CDG / Area_CDE
Area_DEG = TriangleArea(pD,pE,pG) ; l = Area_DEG / Area_CDE
Area_ECG = TriangleArea(pE,pC,pG) ; m = Area_ECG / Area_CDE
SquaredArea = (pi*Area_CDE)^2*(1 - 2*k)*(1 - 2*l)*(1 - 2*m)
print "SquaredArea = ",SquaredArea.n()

# Computing maximum for squared area
hh = 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) = ",f0
print "f0(1/2) = ",f0(1/2).n()," vs f0(s/2) = ",f0(s/2).n()
f1(hh) = f0.derivative(hh).factor(); print "f1 = ",f1
res = solve([f1],hh)
hh0 = res[0].rhs().n(); print "root for extrema hh0 = ",hh0 # hh0 = 4/3
hh1 = res[1].rhs().n(); print "root for extrema hh1 = ",hh1 # hh1 = -1
print "f0(hh0) = ",f0(hh0).n()

# Show graphics
show(g,axes=False,aspect_ratio=True)