OK, I've found what's wrong - the theta = atan2(y, z) parametization is wrong. I plotted the tangent vectors and they weren't tangent to the ellipse, so I saw that the geometric angle the hit point makes with the z axis isn't the same as the theta I want.
It turns out. y / z = b sin(theta) / c cos(theta) = (b/c)tan(theta)
so theta = atan2((c / b) y, z)