From eb317b83c42bf9053f6b7b7d0dc32f8b6b8bebce Mon Sep 17 00:00:00 2001 From: Yi-Ting Tu Date: Sun, 22 Dec 2024 14:26:11 -0500 Subject: [PATCH] Use analytical solution for parabolic mirror --- .../js/sceneObjs/mirror/ParabolicMirror.js | 273 +++++++++++++----- .../mirror/ParabolicMirror/colinear.csv | 28 ++ .../mirror/ParabolicMirror/colinear.json | 57 ++++ .../scenes/mirror/ParabolicMirror/default.csv | 36 +-- .../mirror/ParabolicMirror/parallel.csv | 121 ++++++++ .../mirror/ParabolicMirror/parallel.json | 55 ++++ 6 files changed, 485 insertions(+), 85 deletions(-) create mode 100644 test/scenes/mirror/ParabolicMirror/colinear.csv create mode 100644 test/scenes/mirror/ParabolicMirror/colinear.json create mode 100644 test/scenes/mirror/ParabolicMirror/parallel.csv create mode 100644 test/scenes/mirror/ParabolicMirror/parallel.json diff --git a/src/simulator/js/sceneObjs/mirror/ParabolicMirror.js b/src/simulator/js/sceneObjs/mirror/ParabolicMirror.js index cead6659..2ee04dfe 100644 --- a/src/simulator/js/sceneObjs/mirror/ParabolicMirror.js +++ b/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; \ No newline at end of file diff --git a/test/scenes/mirror/ParabolicMirror/colinear.csv b/test/scenes/mirror/ParabolicMirror/colinear.csv new file mode 100644 index 00000000..f4f675ae --- /dev/null +++ b/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 \ No newline at end of file diff --git a/test/scenes/mirror/ParabolicMirror/colinear.json b/test/scenes/mirror/ParabolicMirror/colinear.json new file mode 100644 index 00000000..d864e60f --- /dev/null +++ b/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 +} \ No newline at end of file diff --git a/test/scenes/mirror/ParabolicMirror/default.csv b/test/scenes/mirror/ParabolicMirror/default.csv index 629d6c90..28eb0462 100644 --- a/test/scenes/mirror/ParabolicMirror/default.csv +++ b/test/scenes/mirror/ParabolicMirror/default.csv @@ -1,8 +1,8 @@ Position,Irradiance 0,-0.1451091334635568 -5,-0.181386416829446 -10,-0.19347884461807574 -15,-0.27812583913848404 +5,-0.16929398904081627 +10,-0.21766370019533526 +15,-0.26603341134985425 20,-0.35068040587026256 25,-0.4595122559679304 30,-0.6167138172201172 @@ -11,18 +11,18 @@ Position,Irradiance 45,-1.2817973455947538 50,-1.221335206651605 55,-1.076226073188048 -60,-1.0520412176107885 -65,-0.8101926618381933 -70,-0.822285089626823 -75,-0.9673942230903801 -80,-1.172965495497086 -85,-2.3701158465714203 -90,-2.2250067131078657 -95,-2.176637001953347 -100,-1.6445701792536458 -105,-1.2092427788629752 -110,-0.8827472285699718 -115,-0.7981002340495635 -120,-0.8464699452040826 -125,-0.2539409835612245 -130,-0.012092427788629734 \ No newline at end of file +60,-1.0278563620335288 +65,-0.8343775174154529 +70,-0.8827472285699718 +75,-1.0036715064562693 +80,-1.1850579232857157 +85,-2.3338385632055316 +90,-2.3096537076282724 +95,-2.152452146376088 +100,-1.6566626070422754 +105,-1.172965495497086 +110,-0.9069320841472314 +115,-0.7618229506836742 +120,-0.9190245119358611 +125,-0.229756127983965 +130,-0.048369711154518935 \ No newline at end of file diff --git a/test/scenes/mirror/ParabolicMirror/parallel.csv b/test/scenes/mirror/ParabolicMirror/parallel.csv new file mode 100644 index 00000000..0c484fb4 --- /dev/null +++ b/test/scenes/mirror/ParabolicMirror/parallel.csv @@ -0,0 +1,121 @@ +Position,Irradiance +0,0 +1,0 +2,0 +3,0 +4,0 +5,0 +6,0 +7,0 +8,0 +9,0 +10,0 +11,0 +12,0 +13,0 +14,0 +15,0 +16,0 +17,0 +18,0 +19,0 +20,0 +21,0 +22,0 +23,0 +24,0 +25,0 +26,0 +27,0 +28,0 +29,0 +30,0 +31,0 +32,0 +33,0 +34,0 +35,0 +36,0 +37,0 +38,0 +39,0 +40,0 +41,0 +42,-0.945954298989415 +43,-1.642973256139509 +44,-1.6927603245073728 +45,-1.642973256139509 +46,-1.6927603245073728 +47,-1.6927603245073728 +48,-1.6927603245073728 +49,-1.6927603245073728 +50,-1.6927603245073728 +51,-1.7425473928752366 +52,-1.6927603245073728 +53,-1.7425473928752366 +54,-1.7425473928752366 +55,-1.7425473928752366 +56,-1.7425473928752366 +57,-1.7425473928752366 +58,-1.7425473928752366 +59,-1.7425473928752366 +60,-1.7425473928752366 +61,-1.7425473928752366 +62,-1.7425473928752366 +63,-1.7425473928752366 +64,-1.6927603245073728 +65,-1.7425473928752366 +66,-1.7425473928752366 +67,-1.6927603245073728 +68,-1.7425473928752366 +69,-1.6927603245073728 +70,-1.6927603245073728 +71,-1.7425473928752366 +72,-1.642973256139509 +73,-1.6927603245073728 +74,-1.6927603245073728 +75,-1.642973256139509 +76,-1.642973256139509 +77,-0.945954298989415 +78,0 +79,0 +80,0 +81,0 +82,0 +83,0 +84,0 +85,0 +86,0 +87,0 +88,0 +89,0 +90,0 +91,0 +92,0 +93,0 +94,0 +95,0 +96,0 +97,0 +98,0 +99,0 +100,0 +101,0 +102,0 +103,0 +104,0 +105,0 +106,0 +107,0 +108,0 +109,0 +110,0 +111,0 +112,0 +113,0 +114,0 +115,0 +116,0 +117,0 +118,0 +119,0 \ No newline at end of file diff --git a/test/scenes/mirror/ParabolicMirror/parallel.json b/test/scenes/mirror/ParabolicMirror/parallel.json new file mode 100644 index 00000000..f8478e36 --- /dev/null +++ b/test/scenes/mirror/ParabolicMirror/parallel.json @@ -0,0 +1,55 @@ +{ + "version": 5, + "name": "parallel", + "objs": [ + { + "type": "ParabolicMirror", + "p1": { + "x": 880, + "y": 360 + }, + "p2": { + "x": 880, + "y": 240 + }, + "p3": { + "x": 884.4444444444443, + "y": 300 + } + }, + { + "type": "Beam", + "p1": { + "x": 820, + "y": 240 + }, + "p2": { + "x": 820, + "y": 360 + } + }, + { + "type": "Detector", + "p1": { + "x": 740, + "y": 240 + }, + "p2": { + "x": 740, + "y": 360 + }, + "irradMap": true + } + ], + "width": 1500, + "height": 800, + "rayModeDensity": 20.085536923187668, + "showGrid": true, + "snapToGrid": true, + "lengthScale": 2, + "origin": { + "x": -945.9999999999998, + "y": -285.00000000000006 + }, + "scale": 2.25 +} \ No newline at end of file