-
-
Save se5a/1135132f0f2d8d84770ca960405b9e41 to your computer and use it in GitHub Desktop.
//taken from: https://doi.org/10.1093/comjnl/14.1.81 | |
public static (int x, int y)[] EllipseArrayFromPaper(double semiMaj, double semiMin, double tilt, int numPoints) | |
{ | |
(int x, int y)[] points = new (int x, int y)[numPoints]; | |
var n = numPoints; | |
var a = semiMaj; | |
var b = semiMin; | |
var xc = 0; //this is the center position of the ellipse; | |
var yc = 0; | |
var deltaPhi = 2 * Math.PI / (n - 1); | |
var ct = Math.Cos(tilt); | |
var st = Math.Sin(tilt); | |
var cdp = Math.Cos(deltaPhi); | |
var sdp = Math.Sin(deltaPhi); | |
var alpha = cdp + sdp * st * ct * (a / b - b / a); | |
var bravo = - sdp * ( (b * st) * (b * st) + (a * ct) * (a * ct)) / (a * b); | |
var charli = sdp * ( (b * ct) * (b * ct) + (a * st) * (a * st)) / (a * b); | |
var delta = cdp + sdp * st * ct * (b / a - a / b); | |
delta = delta - (charli * bravo) / alpha; | |
charli = charli / alpha; | |
var x = a * ct; | |
var y = a * st; | |
for (int i = 0; i < numPoints; i++) | |
{ | |
var xn = xc + x; | |
var yn = yc + y; | |
points[i] = ((int)xn, (int)yn); | |
x = alpha * x + bravo * y; | |
y = charli * x + delta * y; | |
} | |
return points; | |
} |
//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;
}
There is a bug in line #34. It uses the modified
x
in the calculation fory
. Fix this using tuplesin place of lines 33-34.