Cut an ellipse in 2 k pieces of equal length?

Posted 13 days ago
84 Views
|
0 Replies
|
0 Total Likes
|
 I need to cut an ellipse given by the inverse FishIM of its weight matrix, and its centre. Here are two examples of matrices, "Good" and "Bad", with various charcteristics (excentricity, length, etc). Centre = {5.\ 180035459148518677402774168851052046478271419194201229075629.\ 825575763955253, 6.765286741741379750483800007963429000445197980403714798821429.\ 825575763955253}; Good = {{0.5879689788535514100187639228053951611316223513920108533286\ 29.524545768291276, 0.271464841220227416810592214350729574581064323382437782577430.\ 137549484531476}, \ {0.271464841220227416810592214350729574581064323382437782577430.\ 137549484531476, 0.152021089466772629675045991875814730649037698172398387773129.\ 56407044570504}}; Bad = {{0.4219874795337851353296839707215259074013647925209521994819\ 29.635459593034348, 0.357303297463283250778605033518552853551896336891054249589430.\ 1919939283046}, \ {0.357303297463283250778605033518552853551896336891054249589430.\ 1919939283046, 0.349704540740865980617302370810882526021971822740217965492229.\ 577952677191043}}; FishIM = Good; Print["Information matrix : ", MatrixForm[FishIM]]; {a, b} = Sqrt[ 1/Eigenvalues@FishIM]; (* racines de VP de \[CapitalSigma]^-1 *) Print[Style["{a,b}=", "Section"], {a, b}]; {Va, Vb} = Eigenvectors@FishIM; \[Theta] = Min[ArcCos[{-1, 0}.Va], ArcCos[{-1, 0}.(-Va)]]; circ = N[\[Pi] (a + b) Hypergeometric2F1[-1/2, -1/2, 1, h]]; (* Gauss-Kummer formula *) Gr = {}; excentricité := Sqrt[1 - Min[a, b]^2/Max[a, b]^2]; (* Mathworld *) h := ((a - b)/(a + b))^2; y[x_] := b/a Sqrt[a^2 - x^2]; Print["Excentricité de l'ellipse: ", N@excentricité]; Print["Circonférence (Gauss-Kummer): ", circ, " Numerical approximation : ", 2 N[ArcLength[{t, y[t]}, {t, -a, a}], Chif]]; Print["surface Théorique: ", N[\[Pi] a b ]]; Print@Graphics[{LightGreen, Ellipsoid[Centre, {a, b}]}, Frame -> True]; Print["surface 1: ", N@RegionMeasure[Ellipsoid[Centre, {a, b}]]]; Print@Graphics[{LightRed, Ellipsoid[Centre, Inverse@FishIM]}, Frame -> True]; Print["surface 2: ", N@RegionMeasure[Ellipsoid[Centre, Inverse@FishIM]]]; EllipsRef = Show[Graphics[{LightGray, Ellipsoid[Centre, Inverse@FishIM]}], ParametricPlot[ Centre + RotationMatrix[\[Theta]].{t, y[t]}, {t, -a, a}, PlotRange -> All, PlotPoints -> 100 , AspectRatio -> Automatic], ParametricPlot[ Centre + RotationMatrix[\[Theta]].{t, -y[t]}, {t, -a, a}, PlotRange -> All, PlotPoints -> 100 , AspectRatio -> Automatic], Graphics[{PointSize@Large, Point[Centre + RotationMatrix[\[Theta]].{-a, 0}]}], Graphics[{PointSize@Large, Yellow, Point[Centre]}], ImageSize -> Large] I wrote this routine to do the job for a centered ellipse: Chif = 30(*\$MachinePrecision*); TourCir[np_, a_, b_, \[Theta]_, circ_] := Block[{X, pred, sol, x, EchEllP, EchEllM, Angles, k, r, eq, u}, eq[u_] := b/a Sqrt[a^2 - u^2]; X = {}; pred = -(a + 0.001); Do[ sol = NMinimize[{(SetPrecision[ ArcLength[ RotationMatrix[\[Theta]].({t, eq[t]}), {t, -a, u}], Chif] - SetPrecision[k circ/(np), Chif])^2, a > u > Max[-a, pred]}, {u}, WorkingPrecision -> Chif]; Print["Verification lengths: ", np ArcLength[ RotationMatrix[\[Theta]].({t, eq[t]}), {t, -a, u /. Last@sol}]/circ, " =?= k=", k]; Print@Chop@sol; pred = u /. Last@sol; AppendTo[X, Chop@pred] , {k, np/2}]; Print["Dim X=", Dimensions@X]; Print["X=", X]; EchEllP = Map[{#, eq[#]} &, X]; EchEllM = Map[-{#, eq[#]} &, X]; Print["Dim EchEll=", Dimensions@EchEllP]; Print@ListPlot[MapAll[Re, Join[EchEllP, EchEllM]], PlotRange -> All, PlotLabel -> "Solutions"]; MapAll[Re, Join[EchEllP, EchEllM]] ]; Parameters: np is half the number of points, a is a eigenvalue of the ellipse a teta is the angle made by the corresponding eigenvector and the vector {-1,0}. From a and teta, we obtain the initial (and arbitrary) point on the ellipse. You already guessed that everything's fine with the Good ellipse: here are the initial point (black), and the centre: Here is the final result (np=4): There is in the routine a verification instruction: the curvilinear abcissa associted with the k th point shoul be k; Print["Verification lengths: ", np ArcLength[ RotationMatrix[\[Theta]].({t, eq[t]}), {t, -a, u /. Last@sol}]/ circ, " =?= k=", k]; This is true for the Good ellipse and false for the Bad one, because the "cut points" are ill-placed: What is wrong with TourCir? I join the complete notebook. Attachments: Answer