Use analytical solution for parabolic mirror

src/simulator/js/sceneObjs/mirror/ParabolicMirror.js

@@ -255,85 +255,224 @@ class ParabolicMirror extends BaseFilter {
     }
   }
 
-  checkRayIntersects(ray) {
-    if (!this.p3) { return; }
-    if (!this.tmp_points || !this.checkRayIntersectFilter(ray)) return;
-    var i, j;
-    var pts = this.tmp_points;
-    var dir = geometry.distance(this.p2, ray.p1) > geometry.distance(this.p1, ray.p1);
-    var incidentPoint;
-    for (j = 0; j < pts.length - 1; j++) {
-      i = dir ? j : (pts.length - 2 - j);
-      var rp_temp = geometry.linesIntersection(geometry.line(ray.p1, ray.p2), geometry.line(pts[i], pts[i + 1]));
-      var seg = geometry.line(pts[i], pts[i + 1]);
-      // need Simulator.MIN_RAY_SEGMENT_LENGTH check to handle a ray that reflects off mirror multiple times
-      if (geometry.distance(ray.p1, rp_temp) < Simulator.MIN_RAY_SEGMENT_LENGTH * this.scene.lengthScale)
-        continue;
-      if (geometry.intersectionIsOnSegment(rp_temp, seg) && geometry.intersectionIsOnRay(rp_temp, ray)) {
-        if (!incidentPoint || geometry.distance(ray.p1, rp_temp) < geometry.distance(ray.p1, incidentPoint)) {
-          incidentPoint = rp_temp;
-          this.tmp_i = i;
-        }
-      }
-    }
-    if (incidentPoint) return incidentPoint;
+  /**
+   * Transform a point from global coordinates to local parabola coordinates
+   * where the origin is at p3 (vertex), x-axis along p1-p2,
+   * and y-axis perpendicular to p1-p2
+   * @param {Point} point - Point in global coordinates
+   * @returns {Object} Transformed point {x, y}
+   */
+  transformToLocal(point) {
+    const p12d = geometry.distance(this.p1, this.p2);
+    // unit vector from p1 to p2
+    const dir1 = [(this.p2.x - this.p1.x) / p12d, (this.p2.y - this.p1.y) / p12d];
+    // perpendicular direction
+    const dir2 = [dir1[1], -dir1[0]];
+
+    // Use p3 as origin instead of midpoint
+    const dx = point.x - this.p3.x;
+    const dy = point.y - this.p3.y;
+
+    return {
+      x: dx * dir1[0] + dy * dir1[1],
+      y: dx * dir2[0] + dy * dir2[1]
+    };
   }
 
-  onRayIncident(ray, rayIndex, incidentPoint) {
-    var rx = ray.p1.x - incidentPoint.x;
-    var ry = ray.p1.y - incidentPoint.y;
-    var i = this.tmp_i;
-    var pts = this.tmp_points;
-    var seg = geometry.line(pts[i], pts[i + 1]);
-    var mx = seg.p2.x - seg.p1.x;
-    var my = seg.p2.y - seg.p1.y;
-
+  /**
+   * Transform a point from local parabola coordinates back to global coordinates
+   * @param {Object} point - Point in local coordinates {x, y}
+   * @returns {Point} Point in global coordinates
+   */
+  transformToGlobal(point) {
+    const p12d = geometry.distance(this.p1, this.p2);
+    // unit vector from p1 to p2
+    const dir1 = [(this.p2.x - this.p1.x) / p12d, (this.p2.y - this.p1.y) / p12d];
+    // perpendicular direction
+    const dir2 = [dir1[1], -dir1[0]];
+
+    return geometry.point(
+      this.p3.x + point.x * dir1[0] + point.y * dir2[0],
+      this.p3.y + point.x * dir1[1] + point.y * dir2[1]
+    );
+  }
 
-    ray.p1 = incidentPoint;
-    var frac;
-    if (Math.abs(mx) > Math.abs(my)) {
-      frac = (incidentPoint.x - seg.p1.x) / mx;
-    } else {
-      frac = (incidentPoint.y - seg.p1.y) / my;
-    }
+  /**
+   * Get the parabola coefficient 'a' in y = ax²
+   * @returns {number} Coefficient a
+   */
+  getParabolaCoefficient() {
+    const p12d = geometry.distance(this.p1, this.p2);
+    const p1Local = this.transformToLocal(this.p1);
+    // Calculate a from the fact that p1 lies on the parabola
+    return p1Local.y / (p1Local.x * p1Local.x);
+  }
 
-    if ((i == 0 && frac < 0.5) || (i == pts.length - 2 && frac >= 0.5)) {
-      ray.p2 = geometry.point(incidentPoint.x + rx * (my * my - mx * mx) - 2 * ry * mx * my, incidentPoint.y + ry * (mx * mx - my * my) - 2 * rx * mx * my);
-    } else {
-      // Use a simple trick to smooth out the slopes of outgoing rays so that image detection works.
-      // However, a more proper numerical algorithm from the beginning (especially to handle singularities) is still desired.
+  /**
+   * Check if the parabola is degenerate (p1, p2, p3 are collinear)
+   * @returns {boolean} True if degenerate
+   */
+  isDegenerate() {
+    // Calculate area of triangle formed by p1, p2, p3
+    // If area is close to 0, points are collinear
+    const area = Math.abs(
+      (this.p2.x - this.p1.x) * (this.p3.y - this.p1.y) -
+      (this.p3.x - this.p1.x) * (this.p2.y - this.p1.y)
+    );
+    return area < 1e-10;
+  }
 
-      var outx = incidentPoint.x + rx * (my * my - mx * mx) - 2 * ry * mx * my;
-      var outy = incidentPoint.y + ry * (mx * mx - my * my) - 2 * rx * mx * my;
+  checkRayIntersects(ray) {
+    if (!this.p3 || !this.checkRayIntersectFilter(ray)) return;
+
+    // Handle degenerate case (linear mirror)
+    if (this.isDegenerate()) {
+      const intersection = geometry.linesIntersection(
+        geometry.line(ray.p1, ray.p2),
+        geometry.line(this.p1, this.p2)
+      );
+      if (intersection &&
+          geometry.intersectionIsOnSegment(intersection, geometry.line(this.p1, this.p2)) &&
+          geometry.intersectionIsOnRay(intersection, ray) &&
+          geometry.distance(ray.p1, intersection) >= Simulator.MIN_RAY_SEGMENT_LENGTH * this.scene.lengthScale) {
+        return intersection;
+      }
+      return null;
+    }
 
-      var segA;
-      if (frac < 0.5) {
-        segA = geometry.line(pts[i - 1], pts[i]);
-      } else {
-        segA = geometry.line(pts[i + 1], pts[i + 2]);
+    // Normal parabolic case
+    // Transform ray to local coordinates
+    const p1Local = this.transformToLocal(ray.p1);
+    const p2Local = this.transformToLocal(ray.p2);
+
+    // Get normalized direction vector
+    const dx = p2Local.x - p1Local.x;
+    const dy = p2Local.y - p1Local.y;
+    const d = Math.sqrt(dx * dx + dy * dy);
+    const vx = dx / d;
+    const vy = dy / d;
+
+    // Get parabola coefficient
+    const a = this.getParabolaCoefficient();
+
+    // Handle the case where ray is parallel to axis of symmetry
+    // In this case, vx ≈ 0 and the quadratic equation becomes degenerate
+    if (Math.abs(vx) < 1e-10) {
+      // For vertical rays, x stays constant
+      const x = p1Local.x;
+      // Intersection occurs at y = ax²
+      const y = a * x * x;
+      // Calculate time to reach intersection
+      const t = (y - p1Local.y) / vy;
+      if (t > Simulator.MIN_RAY_SEGMENT_LENGTH * this.scene.lengthScale &&
+          Math.abs(x) <= geometry.distance(this.p1, this.p2) / 2) {
+        return this.transformToGlobal({x, y});
       }
-      var rxA = ray.p1.x - incidentPoint.x;
-      var ryA = ray.p1.y - incidentPoint.y;
-      var mxA = segA.p2.x - segA.p1.x;
-      var myA = segA.p2.y - segA.p1.y;
+      return null;
+    }
 
-      var outxA = incidentPoint.x + rxA * (myA * myA - mxA * mxA) - 2 * ryA * mxA * myA;
-      var outyA = incidentPoint.y + ryA * (mxA * mxA - myA * myA) - 2 * rxA * mxA * myA;
+    // Solve quadratic equation: a·dx²·t² + (2a·x₁·dx - dy)·t + (a·x₁² - y₁) = 0
+    const A = a * vx * vx;
+    const B = 2 * a * p1Local.x * vx - vy;
+    const C = a * p1Local.x * p1Local.x - p1Local.y;
+
+    const discriminant = B * B - 4 * A * C;
+    if (discriminant < 0) return null;
+
+    // Get solutions
+    const t1 = (-B + Math.sqrt(discriminant)) / (2 * A);
+    const t2 = (-B - Math.sqrt(discriminant)) / (2 * A);
+
+    // Get intersection points
+    const intersections = [];
+    if (t1 > Simulator.MIN_RAY_SEGMENT_LENGTH * this.scene.lengthScale) {
+      const x = p1Local.x + t1 * vx;
+      const y = p1Local.y + t1 * vy;
+      // Check if point is within the mirror bounds
+      if (Math.abs(x) <= geometry.distance(this.p1, this.p2) / 2) {
+        intersections.push({t: t1, point: this.transformToGlobal({x, y})});
+      }
+    }
+    if (t2 > Simulator.MIN_RAY_SEGMENT_LENGTH * this.scene.lengthScale) {
+      const x = p1Local.x + t2 * vx;
+      const y = p1Local.y + t2 * vy;
+      // Check if point is within the mirror bounds
+      if (Math.abs(x) <= geometry.distance(this.p1, this.p2) / 2) {
+        intersections.push({t: t2, point: this.transformToGlobal({x, y})});
+      }
+    }
 
-      var outxFinal;
-      var outyFinal;
+    // Return closest intersection
+    if (intersections.length === 0) return null;
+    return intersections.reduce((closest, current) => 
+      current.t < closest.t ? current : closest
+    ).point;
+  }
 
-      if (frac < 0.5) {
-        outxFinal = outx * (0.5 + frac) + outxA * (0.5 - frac);
-        outyFinal = outy * (0.5 + frac) + outyA * (0.5 - frac);
-      } else {
-        outxFinal = outxA * (frac - 0.5) + outx * (1.5 - frac);
-        outyFinal = outyA * (frac - 0.5) + outy * (1.5 - frac);
-      }
-      //console.log(frac);
-      ray.p2 = geometry.point(outxFinal, outyFinal);
+  onRayIncident(ray, rayIndex, incidentPoint) {
+    // Handle degenerate case (linear mirror)
+    if (this.isDegenerate()) {
+      const dir = [(this.p2.x - this.p1.x), (this.p2.y - this.p1.y)];
+      const len = Math.sqrt(dir[0] * dir[0] + dir[1] * dir[1]);
+      const nx = -dir[1] / len;  // Normal vector
+      const ny = dir[0] / len;
+
+      // Calculate incident vector
+      const dx = incidentPoint.x - ray.p1.x;
+      const dy = incidentPoint.y - ray.p1.y;
+      const d = Math.sqrt(dx * dx + dy * dy);
+      const vx = dx / d;
+      const vy = dy / d;
+
+      // Calculate reflection using r = v - 2(v·n)n
+      const dot = vx * nx + vy * ny;
+      const rx = vx - 2 * dot * nx;
+      const ry = vy - 2 * dot * ny;
+
+      ray.p1 = incidentPoint;
+      ray.p2 = geometry.point(
+        incidentPoint.x + rx,
+        incidentPoint.y + ry
+      );
+      return;
     }
+
+    // Normal parabolic case
+    // Transform to local coordinates
+    const incidentLocal = this.transformToLocal(incidentPoint);
+    const rayStartLocal = this.transformToLocal(ray.p1);
+
+    // Get incident direction
+    const dx = incidentLocal.x - rayStartLocal.x;
+    const dy = incidentLocal.y - rayStartLocal.y;
+    const d = Math.sqrt(dx * dx + dy * dy);
+    const vx = dx / d;
+    const vy = dy / d;
+
+    // Get parabola coefficient
+    const a = this.getParabolaCoefficient();
+
+    // Calculate normal vector at intersection point (-2ax, 1)
+    const nx = -2 * a * incidentLocal.x;
+    const ny = 1;
+    const nlen = Math.sqrt(nx * nx + ny * ny);
+    const nnx = nx / nlen;
+    const nny = ny / nlen;
+
+    // Calculate reflection direction: r = v - 2(v·n)n
+    const dot = vx * nnx + vy * nny;
+    const rx = vx - 2 * dot * nnx;
+    const ry = vy - 2 * dot * nny;
+
+    // Set new ray endpoint in global coordinates
+    ray.p1 = incidentPoint;
+    const reflectedPoint = this.transformToGlobal({
+      x: incidentLocal.x + rx,
+      y: incidentLocal.y + ry
+    });
+    ray.p2 = reflectedPoint;
   }
+
 };
 
 export default ParabolicMirror;

test/scenes/mirror/ParabolicMirror/colinear.csv

@@ -0,0 +1,28 @@
+Position,Irradiance
+0,-0.09673942230903787
+5,-0.13301670567492707
+10,-0.15720156125218654
+15,-0.21766370019533526
+20,-0.27812583913848404
+25,-0.3869576892361518
+30,-0.532066822699709
+35,-0.7497305228950445
+40,-0.9673942230903801
+45,-1.0036715064562693
+50,-0.7497305228950445
+55,-0.41114254481341134
+60,-0.5683441060655982
+65,-0.6167138172201172
+70,-0.677175956163266
+75,-0.7739153784723041
+80,-0.9673942230903801
+85,-1.2576124900174943
+90,-1.6445701792536458
+95,-1.6445701792536458
+100,-1.2455200622288645
+105,-0.8343775174154529
+110,-0.5804365338542279
+115,-0.4232349726020411
+120,-0.3144031225043733
+125,-0.26603341134985425
+130,-0.024184855577259468

test/scenes/mirror/ParabolicMirror/colinear.json

@@ -0,0 +1,57 @@
+{
+  "version": 5,
+  "name": "colinear",
+  "objs": [
+    {
+      "type": "PointSource",
+      "x": 670,
+      "y": 339,
+      "brightness": 0.2
+    },
+    {
+      "type": "PointSource",
+      "x": 678.6923076923077,
+      "y": 291.61538461538464,
+      "brightness": 0.2
+    },
+    {
+      "type": "Detector",
+      "p1": {
+        "x": 668.8461538461538,
+        "y": 243.00000000000003
+      },
+      "p2": {
+        "x": 655.9230769230769,
+        "y": 372.84615384615387
+      },
+      "irradMap": true,
+      "binSize": 5
+    },
+    {
+      "type": "ParabolicMirror",
+      "p1": {
+        "x": 660,
+        "y": 300
+      },
+      "p2": {
+        "x": 700,
+        "y": 340
+      },
+      "p3": {
+        "x": 680,
+        "y": 320
+      }
+    }
+  ],
+  "width": 1500,
+  "height": 800,
+  "rayModeDensity": 3.307855188319685,
+  "showGrid": true,
+  "snapToGrid": true,
+  "lengthScale": 2,
+  "origin": {
+    "x": -219.8749999999999,
+    "y": -202.87500000000006
+  },
+  "scale": 1.625
+}