An easier way to find the properties of an ellipse from its equation?
a = 0.00228797 ; b = 0.00138781/2; c = 0.00209387; d = -0.281261; e = -0.832702;
f = -1.43328;
ellipse1 = a x^2 + 2 b x y + c y^2 + d x + e y + f;
ContourPlot[ellipse1 == 0, {x, -500, 500}, {y, -500, 500}, ImageSize -> 200]
uv =First@Solve[ForAll[{x, y},
a (x - u)^2 + 2 b (x - u) (y - v) + c (y - v)^2 + f2 == ellipse1], {u, v, f2}, Reals]
ellipse2 = ellipse1 /. {x -> x + u, y -> y + v} /. uv // Simplify // Chop
ContourPlot[{ellipse1 == 0, ellipse2 == 0}, {x, -500, 500}, {y, -500, 500},
ImageSize -> 200, ContourStyle -> {Red, Blue}, Axes -> True, AxesOrigin -> {0, 0}]
solphi = First@Solve[Tan[2 phi] == 2 b/(a - c), phi] // Quiet
rotate = RotationMatrix[phi /. solphi].{x, y};
ellipse3 = ellipse2 /. {x -> First@rotate, y -> Last@rotate} // Simplify // Chop
ContourPlot[{ellipse1 == 0, ellipse2 == 0, ellipse3 == 0}, {x, -500, 500}, {y, -500,
500}, ImageSize -> 200, ContourStyle -> {Red, Blue, Green}, Axes -> True,
AxesOrigin -> {0, 0}]
Edit 1 Way 2
Clear["Global`*"]
ellipse1 = -1.43328 - 0.281261 x + 0.00228797 x^2 - 0.832702 y +
0.00138781 x y + 0.00209387 y^2;
ellipse2 = Collect[#, {x, y}] &@(c (x^2/a^2 + y^2/b^2 - 1) /. {x ->
x Cos[\[Alpha]] - y Sin[\[Alpha]], y -> x Sin[\[Alpha]] + y Cos[\[Alpha]]} /.
{x -> x - h, y -> y - k});
It must has a “c” in this function,because if f1==0 and f2==0 then f1 == c f2 not f1==f2.
func = Flatten@MapThread[Equal, CoefficientList[#, {x, y}] & /@
{ellipse1, ellipse2}, 2]// DeleteCases[#, True] & // Simplify
There are six functions for six variables(a,b,c,h,k,[Alpha]).It is too difficult for NSolve to solve this six functions,so use FindRoot.
FindRoot[func, {a, 40}, {b, 50}, {c, 30}, {h, 2}, {k, 300}, {\[Alpha],Pi/4}]
I post this for completeness (accepting that it was somewhat opaque).
In the series of questions it was known that the quadric form was an ellipse (this could be established using determinants as a previous poster has used).
Accepting an ellipse quadratic form,e.g. the test one provided
ellipse = -1.43328 - 0.281261 x + 0.00228797 x^2 - 0.832702 y +
0.00138781 x y + 0.00209387 y^2
The following:
cf[pol_] := Module[{c, h, k, cons},
c = CoefficientRules[pol];
{h, k} =
0.5 Inverse[{{-({2, 0} /. c), -0.5 ({1, 1} /.
c)}, {-0.5 ({1, 1} /. c), -({0, 2} /. c)}}].{{1, 0} /.
cr, {0, 1} /. c};
cons = {({0, 0} /.
c), ({2, 0} /. c) h^2, ({0, 2} /. c) k^2, ({1, 1} /.
c) h k}.{1, -1, -1, -1};
{h, k, cons};
Expand[({2, 0} /. c) (x)^2 + ({0, 2} /. c) (y)^2 + ({1, 1} /.
c) x y + cons]]
prax[pol_] := Module[{c, mat},
c = CoefficientRules[pol];
mat = {{{2, 0} /. c,
0.5 ({1, 1} /. c)}, {0.5 ({1, 1} /. c), {0, 2} /. c}};
Eigenvalues[mat].{x^2, y^2} + ({0, 0} /. c)
]
cf
centers the ellipse. The ugly code is just completing the square and subtracting the constant to correct.
prax
: just uses eigenvalues of matrix of quadric form to eliminate the cross term (x y).
Applying to test case:
e2 = cf[ellipse];
e3 = prax[e2];
Column[TraditionalForm /@ {ellipse, e2, e3}]
This is consistent (reassuringly) with Chenminqis answer.
Visualizing:
ContourPlot[{ellipse == 0, e2 == 0, e3 == 0}, {x, -500,
500}, {y, -500, 500}, ContourStyle -> {Red, Green, Blue}]
Something little bit different.
a = 0.00228797; b =
0.00138781/
2; c = 0.00209387; d = -0.281261; e = -0.832702; f = -1.43328;
ellipse1 = a x^2 + 2 b x y + c y^2 + d x + e y + f;
point1 = {x, y} /. Maximize[{y, ellipse1 == 0}, {x, y}][[2]];
point2 = {x, y} /. Maximize[{-y, ellipse1 == 0}, {x, y}][[2]];
center = (point1 + point2)/2;
angle = phi /. Solve[Tan[2 phi] == 2 b/(a - c), phi][[1]] // Quiet;
ellipse2 = ellipse1 /. Thread[{x, y} -> {x, y} + center];
r = RotationTransform[angle];
ellipse3 = ellipse2 /. Thread[{x, y} -> r[{x, y}]];
ContourPlot[{ellipse1 == 0, ellipse2 == 0, ellipse3 == 0}, {x, -300,
300}, {y, -300, 500}, ImageSize -> 300, Axes -> True,
AspectRatio -> Automatic]