Committing work in process of alternative ellipse algorithm so hadbn …

…can review and help debug. Approach shows promise however results aren't as expected (see Turfjs#2739 (comment)) for details.

packages/turf-ellipse/bench.ts

@@ -1,5 +1,5 @@
 import Benchmark from "benchmark";
-import { ellipse } from "./index.js";
+import { ellipse, ellipseOld, ellipseRiemann } from "./index.js";
 
 /**
  * Benchmark Results
@@ -15,14 +15,25 @@ const xSemiAxis = 50;
 const ySemiAxis = 10;
 
 suite
-  .add("turf-ellipse - 8 steps", () =>
-    ellipse(center, xSemiAxis, ySemiAxis, { steps: 8 })
+  // .add("turf-ellipse - 8 steps", () =>
+  //   ellipse(center, xSemiAxis, ySemiAxis, { steps: 8 })
+  // )
+  // .add("turf-ellipse - 64 steps", () =>
+  //   ellipse(center, xSemiAxis, ySemiAxis, { steps: 64 })
+  // )
+  .add("turf-ellipse - old", () => ellipseOld(center, xSemiAxis, ySemiAxis, {}))
+  .add("turf-ellipse - new q1", () =>
+    ellipse(center, xSemiAxis, ySemiAxis, { accuracy: 1 })
   )
-  .add("turf-ellipse - 64 steps", () =>
-    ellipse(center, xSemiAxis, ySemiAxis, { steps: 64 })
+  .add("turf-ellipse - new q2", () =>
+    ellipse(center, xSemiAxis, ySemiAxis, { accuracy: 2 })
   )
-  .add("turf-ellipse - 256 steps", () =>
-    ellipse(center, xSemiAxis, ySemiAxis, { steps: 256 })
+  .add("turf-ellipse - new q3", () =>
+    ellipse(center, xSemiAxis, ySemiAxis, { accuracy: 3 })
   )
+  .add("turf-ellipse - riemann", () =>
+    ellipseRiemann(center, xSemiAxis, ySemiAxis, {})
+  )
+
   .on("cycle", (e) => console.log(String(e.target)))
   .run();

packages/turf-ellipse/index.ts

@@ -6,12 +6,15 @@ import {
   Coord,
   Units,
   point,
+  lineString,
 } from "@turf/helpers";
+import { rhumbBearing } from "@turf/rhumb-bearing";
+import { rhumbDestination } from "@turf/rhumb-destination";
 import { destination } from "@turf/destination";
 import { distance } from "@turf/distance";
 import { transformRotate } from "@turf/transform-rotate";
 import { getCoord } from "@turf/invariant";
-import { GeoJsonProperties, Feature, Polygon } from "geojson";
+import { GeoJsonProperties, Feature, Polygon, Point, Position } from "geojson";
 
 /**
  * Takes a {@link Point} and calculates the ellipse polygon given two semi-axes expressed in variable units and steps for precision.
@@ -51,7 +54,7 @@ function ellipse(
 ): Feature<Polygon> {
   // Optional params
   options = options || {};
-  const steps = options.steps || 64;
+  let steps = options.steps || 64;
   const units = options.units || "kilometers";
   const angle = options.angle || 0;
   const pivot = options.pivot || center;
@@ -154,5 +157,251 @@ function ellipse(
   return polygon([coordinates], properties);
 }
 
-export { ellipse };
+function ellipseRiemann(
+  center: Coord,
+  xSemiAxis: number,
+  ySemiAxis: number,
+  options: {
+    steps?: number;
+    units?: Units;
+    angle?: number;
+    pivot?: Coord;
+    properties?: GeoJsonProperties;
+  }
+): Feature<Polygon> {
+  // Optional params
+  options = options || {};
+  let steps = options.steps || 64;
+  const units = options.units || "kilometers";
+  const angle = options.angle || 0;
+  const pivot = options.pivot || center;
+  const properties = options.properties || {};
+  // validation
+  if (!center) throw new Error("center is required");
+  if (!xSemiAxis) throw new Error("xSemiAxis is required");
+  if (!ySemiAxis) throw new Error("ySemiAxis is required");
+  if (!isObject(options)) throw new Error("options must be an object");
+  if (!isNumber(steps)) throw new Error("steps must be a number");
+  if (!isNumber(angle)) throw new Error("angle must be a number");
+
+  const centerCoords = getCoord(
+    transformRotate(point(getCoord(center)), angle, { pivot })
+  );
+
+  const coordinates: number[][] = [];
+
+  let r = Math.sqrt(
+    (Math.pow(xSemiAxis, 2) * Math.pow(ySemiAxis, 2)) /
+      (Math.pow(xSemiAxis * Math.cos(degreesToRadians(-angle)), 2) +
+        Math.pow(ySemiAxis * Math.sin(degreesToRadians(-angle)), 2))
+  );
+  let currentCoords = getCoord(
+    destination(centerCoords, r, 0, { units: units })
+  );
+  coordinates.push(currentCoords);
+
+  // Divide steps by 4 for one quadrant
+  steps = steps / 4;
+  steps = 10;
+
+  let quadrant: Position[] = [];
+
+  const a = xSemiAxis;
+  const b = ySemiAxis;
+
+  // Gradient x intersect
+  const c = b;
+
+  // Gradient of line
+  const m = (a - b) / (Math.PI / 2);
+
+  // Area under line
+  const A = ((a + b) * Math.PI) / 4;
+
+  // Weighting function
+  const v = 0.5;
+
+  const k = steps;
+
+  let w = 0;
+  let x = 0;
+
+  for (let i = 0; i <= steps; i++) {
+    x += w;
+
+    if (m === 0) {
+      // It's a circle, so use simplified c*w - A/k == 0
+      w = A / k / c;
+    } else {
+      // Otherwise, use full (v*m)*w^2 + (m*x+c)*w - A/k == 0
+      // Solve as quadratic ax^2 + bx + c = 0
+      w =
+        (-(m * x + c) +
+          Math.sqrt(Math.pow(m * x + c, 2) - 4 * (v * m) * -(A / k))) /
+        (2 * (v * m));
+    }
+
+    // console.log("Aw", v * m * Math.pow(w, 2) + (m * x + c) * w);
+
+    // Scale x (between 0 and PI/2) up to the length of the xSemiAxis
+    const scaleX = (x / (Math.PI / 2)) * a;
+
+    // x^2 / a^2 + y^2 / b^2 = 1, solve for y
+    const y = Math.sqrt(
+      (1 - Math.pow(scaleX, 2) / Math.pow(a, 2)) * Math.pow(b, 2)
+    );
+
+    // console.log("scaleX", scaleX, y);
+    const dist = Math.sqrt(Math.pow(scaleX, 2) + Math.pow(y, 2));
+    const baseAngle =
+      (Math.atan(Math.abs(y) / Math.abs(scaleX)) * 180) / Math.PI;
+    quadrant.push([scaleX, y, baseAngle, dist]);
+
+    console.log(scaleX);
+  }
+
+  // Add final point which is just [a, 0]
+  // const dist = Math.sqrt(Math.pow(a, 2) + Math.pow(0, 2));
+  // const baseAngle = (Math.atan(Math.abs(0) / Math.abs(a)) * 180) / Math.PI;
+  // quadrant.push([a, 0, baseAngle, dist]);
+
+  // console.log(JSON.stringify(quadrant));
+  const geometricPoints: Position[] = [];
+  const geodesicPoints: Position[] = [];
+
+  // Push the quadrant forwards, with x values negated (NW quadrant)
+  for (const pt of quadrant) {
+    const [baseAngle, dist] = [pt[2], pt[3]];
+    const bearing = 0 - (90 - baseAngle);
+    geodesicPoints.push(
+      destination(center, dist, bearing).geometry.coordinates
+    );
+    geometricPoints.push([-pt[0], pt[1]]);
+  }
+
+  // Push the quadrant backwards, with x and y values negated (SW quadrant)
+  quadrant.reverse();
+  for (const pt of quadrant) {
+    const [baseAngle, dist] = [pt[2], pt[3]];
+    const bearing = -90 - baseAngle;
+    geodesicPoints.push(
+      destination(center, dist, bearing).geometry.coordinates
+    );
+    geometricPoints.push([-pt[0], -pt[1]]);
+  }
+
+  // Push the quadrant forwards, with y values negated (SE quadrant)
+  quadrant.reverse();
+  for (const pt of quadrant) {
+    const [baseAngle, dist] = [pt[2], pt[3]];
+    const bearing = 180 - (90 - baseAngle);
+    geodesicPoints.push(
+      destination(center, dist, bearing).geometry.coordinates
+    );
+    geometricPoints.push([pt[0], -pt[1]]);
+  }
+
+  // Push the quadrant backwards, verbatim (NE quadrant)
+  quadrant.reverse();
+  for (const pt of quadrant) {
+    const [baseAngle, dist] = [pt[2], pt[3]];
+    const bearing = 90 - baseAngle;
+    geodesicPoints.push(
+      destination(center, dist, bearing).geometry.coordinates
+    );
+    geometricPoints.push([pt[0], pt[1]]);
+  }
+
+  console.log(JSON.stringify(lineString(geometricPoints)));
+
+  return polygon([geodesicPoints], properties);
+}
+
+function ellipseOld(
+  center: Coord,
+  xSemiAxis: number,
+  ySemiAxis: number,
+  options: {
+    steps?: number;
+    units?: Units;
+    angle?: number;
+    pivot?: Coord;
+    properties?: GeoJsonProperties;
+  }
+): Feature<Polygon> {
+  // Optional params
+  options = options || {};
+  const steps = options.steps || 64;
+  const units = options.units || "kilometers";
+  const angle = options.angle || 0;
+  const pivot = options.pivot || center;
+  const properties = options.properties || {};
+
+  // validation
+  if (!center) throw new Error("center is required");
+  if (!xSemiAxis) throw new Error("xSemiAxis is required");
+  if (!ySemiAxis) throw new Error("ySemiAxis is required");
+  if (!isObject(options)) throw new Error("options must be an object");
+  if (!isNumber(steps)) throw new Error("steps must be a number");
+  if (!isNumber(angle)) throw new Error("angle must be a number");
+
+  const centerCoords = getCoord(center);
+  if (units !== "degrees") {
+    const xDest = rhumbDestination(center, xSemiAxis, 90, { units });
+    const yDest = rhumbDestination(center, ySemiAxis, 0, { units });
+    xSemiAxis = getCoord(xDest)[0] - centerCoords[0];
+    ySemiAxis = getCoord(yDest)[1] - centerCoords[1];
+  }
+
+  const coordinates: number[][] = [];
+  for (let i = 0; i < steps; i += 1) {
+    const stepAngle = (i * -360) / steps;
+    let x =
+      (xSemiAxis * ySemiAxis) /
+      Math.sqrt(
+        Math.pow(ySemiAxis, 2) +
+          Math.pow(xSemiAxis, 2) * Math.pow(getTanDeg(stepAngle), 2)
+      );
+    let y =
+      (xSemiAxis * ySemiAxis) /
+      Math.sqrt(
+        Math.pow(xSemiAxis, 2) +
+          Math.pow(ySemiAxis, 2) / Math.pow(getTanDeg(stepAngle), 2)
+      );
+
+    if (stepAngle < -90 && stepAngle >= -270) x = -x;
+    if (stepAngle < -180 && stepAngle >= -360) y = -y;
+    if (units === "degrees") {
+      const angleRad = degreesToRadians(angle);
+      const newx = x * Math.cos(angleRad) + y * Math.sin(angleRad);
+      const newy = y * Math.cos(angleRad) - x * Math.sin(angleRad);
+      x = newx;
+      y = newy;
+    }
+
+    coordinates.push([x + centerCoords[0], y + centerCoords[1]]);
+  }
+  coordinates.push(coordinates[0]);
+  if (units === "degrees") {
+    return polygon([coordinates], properties);
+  } else {
+    return transformRotate(polygon([coordinates], properties), angle, {
+      pivot,
+    });
+  }
+}
+
+/**
+ * Get Tan Degrees
+ *
+ * @private
+ * @param {number} deg Degrees
+ * @returns {number} Tan Degrees
+ */
+function getTanDeg(deg: number) {
+  const rad = (deg * Math.PI) / 180;
+  return Math.tan(rad);
+}
+
+export { ellipse, ellipseOld, ellipseRiemann };
 export default ellipse;

packages/turf-ellipse/package.json

@@ -57,6 +57,7 @@
     "@turf/bbox-polygon": "workspace:^",
     "@turf/circle": "workspace:^",
     "@turf/intersect": "workspace:^",
+    "@turf/meta": "workspace:^",
     "@turf/truncate": "workspace:^",
     "@types/benchmark": "^2.1.5",
     "@types/tape": "^4.2.32",
@@ -75,6 +76,8 @@
     "@turf/distance": "workspace:^",
     "@turf/helpers": "workspace:^",
     "@turf/invariant": "workspace:^",
+    "@turf/rhumb-bearing": "workspace:^",
+    "@turf/rhumb-destination": "workspace:^",
     "@turf/transform-rotate": "workspace:^",
     "@types/geojson": "^7946.0.10",
     "tslib": "^2.6.2"