se5a/gist:1135132f0f2d8d84770ca960405b9e41

There is a bug in line #34. It uses the modified x in the calculation for y. Fix this using tuples

(x,y) = (alpha * x + bravo * y, charli * x + delta * y);

in place of lines 33-34.

//version using a matrix math to do the roatation around the focal point.
// note it also outputs a vector2 that uses doubles.
public static DVec.Vector2[] EllipseArrayMTX(double semiMaj, double eccentricity, double tilt, int numPoints)
{
DVec.Vector2[] points = new DVec.Vector2[numPoints];

var n = numPoints;
var a = semiMaj;
var b = semiMaj * Math.Sqrt(1 - eccentricity * eccentricity);
var xc = 0; //this is the center position of the ellipse;
var yc = 0;
var linEcc = Math.Sqrt(a * a - b * b);

var deltaPhi = 2 * Math.PI / (n - 1);
var cdp = Math.Cos(deltaPhi);
var sdp = Math.Sin(deltaPhi);

var alpha = cdp;
var bravo = - sdp * (a * a) / (a * b);
var charli = sdp * ( b * b) / (a * b);
var delta = cdp - (charli * bravo) / alpha;
charli = charli / alpha;

double x = a;
double y = 0;

for (int i = 0; i < numPoints; i++)
{
    var xn = xc + x;
    var yn = yc + y;
    points[i] = new DVec.Vector2(xn, yn);
    x = alpha * x + bravo * y;
    y = charli * x + delta * y;
}

var moveMtx = Matrix.IDTranslate(linEcc, 0);
var rotMtx = Matrix.IDRotate(tilt);
var endMtx = moveMtx * rotMtx;
var pnts = endMtx.Transform(points);

return pnts;

}

//this is the un-optomized version from the paper.
public static DVec.Vector2[] EllipsArrayPaper1(double semiMaj, double eccentricity, double tilt, double start, double end, int numPoints)
{
DVec.Vector2[] points = new DVec.Vector2[numPoints];
double deltaPhi = 2 * Math.PI / (numPoints - 1);
double ct = Math.Cos(tilt);
double st = Math.Sin(tilt);
double cdp = Math.Cos(deltaPhi);
double sdp = Math.Sin(deltaPhi);
double cndp = 1.0;
double sndp = 0;
double a = semiMaj;
double b = semiMaj * Math.Sqrt(1 - eccentricity * eccentricity);
double xc = 0;
double yc = 0;

for (int i = 0; i < numPoints; i++)
{
    double x = a * cndp;
    double y = b * sndp;

    double xi = xc + x * ct - y * st;
    double yi = yc + x * st + y * ct;
    points[i] = new DVec.Vector2(xi, yi);

    double temp = cndp * cdp - sndp * sdp;
    sndp = sndp * cdp + cndp * sdp;
    cndp = temp;
}

return points;

}