Greetings. I'm trying to do an optical simulation and wrote a ray, ellipsoid intersection function for Matlab/Octave. I fire a bundle of rays from 1 focus, and reflect them off the ellipse (see pictures).

In ellipse2.png, I choose the focus to be (0, -sqrt(1.2^2 - 1)) as it should be. I compute |hit - foci0| + |hit - foci1| and they're a constant 2.4 as expected, but the rays don't converge.

In ellipse1.png, I set the focus to (0, -.46). The rays now converge, but |hit - foci0| + |hit - foci1| no longer seem to be constant.

So I'm guessing the problem is with my intersection routine. Can anyone point out what's wrong?

(my code is a 2D version of the 3D case (YZ plane), so theta = 0 means straight up, instead of left)

1. First I find ray, ellipse intersection by expressing the ray as parametric equations (simple quadratic formula)

2. Find theta parameter associated with intersection: parametric form of ellipse is: y = b * sin(theta) z = c * cos(theta)

theta = atan2(hit.y, hit.z)

3. compute tangent vector dp_dtheta = [b * cos(theta), -c * sin(theta)]

4. rotate tangent vector -90 degrees to get surface normal at intersection

5. reflect the incident vector about the normal vector:

generate a rotation matrix with vector perpendicular to both incident & normal vector as axis. apply the rotation matrix to the incident vector to get the outgoing vector.

My main uncertainty is #3. I do this based on my knowledge of the PBRT raytracer. They compute the surface normal as CrossProduct(dp_du, dp_dv), where u & v, the surface parametization is in [0, 1].

update:

I've estimated the surface tangent by evaluating P(theta + 0.001) - P(theta). The estimated tangent vector matches the 1 computed above to 3 figures.

So now it's really a mystery. Is the reflection incorrect or am I using the wrong ellipse focus formula?