How to Plot Ellipse Given a General Equation in R

Obtain vertices of the ellipse on an ellipse covariance plot (created by `car::ellipse`)

For practical purposes, @Tensibai's answer is probably good enough. Just use a large enough value for the segments argument so that the points give a good approximation to the true vertices.

If you want something a bit more rigorous, you can solve for the location along the ellipse that maximises/minimises the distance from the center, parametrised by the angle. This is more complex than just taking angle={0, pi/2, pi, 3pi/2} because of the presence of the shape matrix. But it's not too difficult:

# location along the ellipse
# linear algebra lifted from the code for ellipse()
ellipse.loc <- function(theta, center, shape, radius)
{
    vert <- cbind(cos(theta), sin(theta))
    Q <- chol(shape, pivot=TRUE)
    ord <- order(attr(Q, "pivot"))
    t(center + radius*t(vert %*% Q[, ord]))
}

# distance from this location on the ellipse to the center 
ellipse.rad <- function(theta, center, shape, radius)
{
    loc <- ellipse.loc(theta, center, shape, radius)
    (loc[,1] - center[1])^2 + (loc[,2] - center[2])^2
}

# ellipse parameters
center <- c(-0.05, 0.09)
A <- matrix(c(20.43, -8.59, -8.59, 24.03), nrow=2)
radius <- 1.44

# solve for the maximum distance in one hemisphere (hemi-ellipse?)
t1 <- optimize(ellipse.rad, c(0, pi - 1e-5), center=center, shape=A, radius=radius, maximum=TRUE)$m
l1 <- ellipse.loc(t1, center, A, radius)

# solve for the minimum distance
t2 <- optimize(ellipse.rad, c(0, pi - 1e-5), center=center, shape=A, radius=radius)$m
l2 <- ellipse.loc(t2, center, A, radius)

# other points obtained by symmetry
t3 <- pi + t1
l3 <- ellipse.loc(t3, center, A, radius)

t4 <- pi + t2
l4 <- ellipse.loc(t4, center, A, radius)

# plot everything
MASS::eqscplot(center[1], center[2], xlim=c(-7, 7), ylim=c(-7, 7), xlab="", ylab="")
ellipse(center, A, radius, col="red", lty=2)
points(rbind(l1, l2, l3, l4), cex=2, col="blue", lwd=2)

How to draw rotated ellipse in Maxima CAS?

The formula "y: acos(t)*c + bsin(t)*c" does not seem to be correct. Try changing it to "y: acos(t)*s + bsin(t)*c" to see if it works.

Plotting error ellipses in R

A few month ago I wrote a little function to draw ellipses to answer someone else's question. By simplifying it a little we can achieve something useful for your issue I think.

ellipse <- function(a,b,xc,yc,...){
    # a is the length of the axis parallel to the x-axis
    # b is the length of the axis parallel to the y-axis
    # xc and yc are the coordinates of the center of the ellipse
    # ... are any arguments that can be passed to function lines
    t <- seq(0, 2*pi, by=pi/100)
    xt <- xc + a*cos(t)
    yt <- yc + b*sin(t)
    lines(xt,yt,...)
    }

plot(UPb, PbPb, pch=19,
     xlab = expression({}^238*"U/"*{}^206*"Pb"), ylab = expression({}^207*"Pb/"*{}^206*"Pb"), 
     xlim = c(0,2500),ylim = c(0, 1), las=1, cex.lab = 1.5)

apply(upbiso, 1, 
      function(x)ellipse(a=x[2]*x[1]/100, b=x[5]*x[4]/100, 
                               xc=x[1], yc=x[4], col="red"))

How to plot an ellipse on canvas from 2 points on the ellipse, where slope of major axis (rx), and minor axis (ry) length are unknown

In case anyone is interested, here is my solution to the problem, which isn't really "the" solution. If anyone can solve this I would still be happy to know.

Since I can't solve for both slope of the major axis and length of the minor axis, I just take a guess at slope and then test how close it is, and then refine the result by trying in a smaller and smaller region. Since the final ellipse that gets drawn is actually an estimation constructed from bezier curves, I can get close enough in a reasonable amount of time.

function getEllipseFrom2Points (center, u, v, a) {
    function getSemiMinorAxis([x, y], a, t) {
        // equation for rotated ellipse
        // b = a(ycos(t) - xsin(t)) / sqrt(a^2 - x^2cos^2(t) - 2xysin(t)cos(t) - y^2sin^2(t)) and
        // b = a(xsin(t) - ycos(t)) / sqrt(a^2 - x^2cos^2(t) - 2xysin(t)cos(t) - y^2sin^2(t)) 
        // where a^2 !== (xcos(t) + ysin(t))^2
        // and aycos(t) !== axsin(t)

        if (a ** 2 !== (x * Math.cos(t) + y * Math.sin(t)) ** 2 &&
            a * y * Math.cos(t) !== a * x * Math.sin(t)) {
            var b = [],
                q = (Math.sqrt(a ** 2 - x ** 2 * (Math.cos(t)) ** 2 - 2 * x * y * Math.sin(t) * Math.cos(t) - y ** 2 * (Math.sin(t)) ** 2));
            b[0] = (a * (y * Math.cos(t) - x * Math.sin(t))) / q;
            b[1] = (a * (x * Math.sin(t) - y * Math.cos(t))) / q;
            return b;
        }
    }
    function getAngle_radians(point1, point2){
        return Math.atan2(point2[1] - point1[1], point2[0] - point1[0]);
    }
    function getDistance(point1, point2) {
        return Math.sqrt((point2[0] - point1[0]) ** 2 + (point2[1] - point1[1]) ** 2);
    }
    function rotatePoint(point, center, radians) {
        var x = (point[0] - center[0]) * Math.cos(radians) - (point[1] - center[1]) * Math.sin(radians) + center[0];
        var y = (point[1] - center[1]) * Math.cos(radians) + (point[0] - center[0]) * Math.sin(radians) + center[1];
        return [x, y];
    }
    function measure(ellipseRotation, pointOnEllipse, minorAxisLength) {
        var d = getDistance(point, pointOnEllipse);
        if (d < bestDistanceBetweenPointAndEllipse) {
            bestDistanceBetweenPointAndEllipse = d;
            bestEstimationOfB = minorAxisLength;
            bestEstimationOfR = ellipseRotation;
        }
    }
    function getBestEstimate(min, max) {
        var testIncrement = (max - min) / 10;
        for (let r = min; r < max; r = r + testIncrement) {
            if (radPoint1 < r && radPoint2 < r || radPoint1 > r && radPoint2 > r) {//points both on same side of ellipse
                semiMinorAxis = getSemiMinorAxis(v, a, r);
                if (semiMinorAxis) {
                    for (let t = 0; t < circle; t = t + degree) {
                        ellipsePoint1 = [a * Math.cos(t), semiMinorAxis[0] * Math.sin(t)];
                        ellipsePoint2 = [a * Math.cos(t), semiMinorAxis[1] * Math.sin(t)];
                        point = rotatePoint(u, [0, 0], -r);
                        measure(r, ellipsePoint1, semiMinorAxis[0]);
                        measure(r, ellipsePoint2, semiMinorAxis[1]);
                    }
                }
            }
        }
        count++;
        if (new Date().getTime() - startTime < 200 && count < 10) //refine estimate
            getBestEstimate(bestEstimationOfR - testIncrement, bestEstimationOfR + testIncrement);
    }
    if (center instanceof Array &&
        typeof center[0] === "number" &&
        typeof center[1] === "number" &&
        u instanceof Array &&
        typeof u[0] === "number" &&
        typeof u[1] === "number" &&
        v instanceof Array &&
        typeof v[0] === "number" &&
        typeof v[1] === "number" &&
        typeof a === "number") {

        // translate points
        u = [u[0] - center[0], u[1] - center[1]];
        v = [v[0] - center[0], v[1] - center[1]];

        var bestDistanceBetweenPointAndEllipse = a,
            point,
            semiMinorAxis,
            ellipsePoint1,
            ellipsePoint2,
            bestEstimationOfB,
            bestEstimationOfR,
            radPoint1 = getAngle_radians([0, 0], v),
            radPoint2 = getAngle_radians([0, 0], u),
            circle = 2 * Math.PI,
            degree = circle / 360,
            startTime = new Date().getTime(),
            count = 0;

        getBestEstimate(0, circle);

        var ellipseModel = MakerJs.$(new MakerJs.models.Ellipse(a, bestEstimationOfB))
            .rotate(MakerJs.angle.toDegrees(bestEstimationOfR), [0, 0])
            .move(center)
            .originate([0, 0])
            .$result;

        return ellipseModel;
}

I need an equation for equal movement along an ellipse

If you want the real thing

then the planets closer to the primary focus point (center of mass of stellar system ... very close to star) are moving faster so use Kepler's equation here: C++ implementation of mine. Do not forget to check out all the sub-links in that answer you can find there everything you need.

If you want constant speed instead

Then use parametric ellipse equation

x(a)=x0+rx*cos(a)
y(a)=y0+ry*sin(a)

where a is angle <0,2.0*PI> (x0,y0) is the ellipse center and (rx,ry) are the ellipse semi axises (radii).

if a is incremented with constant speed then the area increase is constant so the a is the mean circular angle not the visual on ellipse !!! For more info look here:

• Issue with ellipse angle calculation and point calculation

[edit1] as MartinR pointed out the speed is not constant

so here is approximation with his formula for speed. Ellipse is axis aligned defined by x0,y0,rx,ry (rx>=ry) the perimeter aproximation l:

h=(rx-ry)/(rx+ry); h*=3.0*h; l=M_PI*(rx+ry)*(1.0+(h/(10.0+sqrt(4.0-h))));

if you want to have n chunks of equal sized steps along the perimeter then

l/=n;

initial computations:

double x0,y0,rx,ry,n,l,h;
x0=Form1->ClientWidth>>1;  // center is centered on form
y0=Form1->ClientHeight>>1;
rx=200; // semiaxises rx>=ry !!!
ry=75;
n=40.0; // number of chunks per ellipse (1/speed)
//l=2.0*M_PI*sqrt(0.5*((rx*rx)+(ry*ry)));  // not accurate enough
h=(rx-ry)/(rx+ry); h*=3.0*h; l=M_PI*(rx+ry)*(1.0+(h/(10.0+sqrt(4.0-h)))); // this is more precise
l/=n; // single step size in units,pixels,or whatever

first the slow bruteforce attack (black):

int i;
double a,da,x,y,xx,yy,ll;
a=0.0;
x=x0+rx*cos(a);
y=y0+ry*sin(a);
for (i=n;i>0;i--)
    {
    xx=x; yy=y;
    for (da=a;;)
        {
        a+=0.001;
        x=x0+rx*cos(a);
        y=y0+ry*sin(a);
        ll=sqrt(((xx-x)*(xx-x))+((yy-y)*(yy-y)));
        if (ll>=l) break;
        } da=a-da;
    scr->MoveTo(5.0+50.0*a,5.0);
    scr->LineTo(5.0+50.0*a,5.0+300.0*da);

    scr->MoveTo(x0,y0);
    scr->LineTo(xx,yy);
    scr->LineTo(x ,y );
    ll=sqrt(((xx-x)*(xx-x))+((yy-y)*(yy-y)));
    scr->TextOutA(0.5*(x+xx)+20.0*cos(a),0.5*(y+yy)+20.0*sin(a),floor(ll));
    }

Now the approximation (Blue):

a=0.0; da=0;
x=x0+rx*cos(a);
y=y0+ry*sin(a);
for (i=n;i>0;i--)
    {
    scr->MoveTo(5.0+50.0*a,5.0+300.0*da);
    xx=rx*sin(a);
    yy=ry*cos(a);
    da=l/sqrt((xx*xx)+(yy*yy)); a+=da;
    scr->LineTo(5.0+50.0*a,5.0+300.0*da);

    xx=x; yy=y;
    x=x0+rx*cos(a);
    y=y0+ry*sin(a);

    scr->MoveTo(x0,y0);
    scr->LineTo(xx,yy);
    scr->LineTo(x ,y );
    ll=sqrt(((xx-x)*(xx-x))+((yy-y)*(yy-y)));
    scr->TextOutA(0.5*(x+xx)+40.0*cos(a),0.5*(y+yy)+40.0*sin(a),floor(ll));
    }

This is clean ellipse step (no debug draws)

a=???; // some initial angle
// point on ellipse 
x=x0+rx*cos(a);
y=y0+ry*sin(a);
// next angle by almost constant speed
xx=rx*sin(a);
yy=ry*cos(a);
da=l/sqrt((xx*xx)+(yy*yy)); a+=da;
// next point on ellipse ...
x=x0+rx*cos(a);
y=y0+ry*sin(a);

Here the output of comparison bruteforce and approximation:

[edit2] little precision boost

 a,da=???; // some initial angle and step (last)
 x=x0+rx*cos(a);
 y=y0+ry*sin(a);
 // next angle by almost constant speed
 xx=rx*sin(a+0.5*da); // use half step angle for aproximation ....
 yy=ry*cos(a+0.5*da);
 da=l/sqrt((xx*xx)+(yy*yy)); a+=da;
 // next point on ellipse ...
 x=x0+rx*cos(a);
 y=y0+ry*sin(a);

the half step angle in approximation lead to much closer result to bruteforce attack

---

Related Topics