CoCalc -- Inellipse of convex quadrilateral.sagews

Project: Pinwheel Geometry and Number Theory

Path: Inellipse of convex quadrilateral.sagews

Views: ⁹¹

%md
## Inellipse of convex quadrilateral

### Inellipse

  The word "inellipse" is used in place of "ellipse inscribed".

  ![Quadrilateral Inellipse](QuadrilateralInellipse.png "Drawing of a quadrilateral inellipse")

  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 of convex quadrilateral

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

#%hide
import sys
sys.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]

DrawingConstants module loading ...
...DrawingConstants module loaded
RationalTrigonometry (RT) module loading ...
...RationalTrigonometry module loaded
EuclideanGeometry (EG) module loading ...
...EuclideanGeometry module loaded
InversiveGeometry(IG) module loading ...
...EuclideanGeometry module loaded
Pinwheel Geometry(PWG) module loading ...
...Pinwheel Geometry module loaded
TrilinearCoordinates (TC) module loading ...
...Trilinear coordinates module loaded
Matrix representation of Conics (MR) module loading ...
... Matrix representation of Conics module loaded


# Drawing using Horwitz notation


# Initialize graphics
g = Graphics()

# Parameters
s = 4
t = 2

# Set convex quadrilateral ABDC
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)
#save(g,'QuadrilateralInellipse.png',axes=False,aspect_ratio=True)

SquaredArea_Ellipse =  4.06156559715611  vs (pi*p*q)^2 =  4.06156559715612
SquaredArea =  4.06156559715611
Solving max of f0(hh) =  hh |--> -1/9*pi^2*(2*hh - 1)*(hh + 2)*(hh - 2)
f0(1/2) =  0.000000000000000  vs f0(s/2) =  0.000000000000000
f1 =  hh |--> -2/9*pi^2*(3*hh - 4)*(hh + 1)
root for extrema hh0 =  1.33333333333333
root for extrema hh1 =  -1.00000000000000
f0(hh0) =  4.06156559715611