Intersection of a line with a specific surface of revolution As part of an art project, I'm trying to write a program to raytrace eggs, as defined by the model described in this paper.
As part of writing this program, I am trying to calculate the intersection between a line in 3d space, and the surface obtained by revolving the equation from that paper, which is.
$
r = T \times (1 + z)^\frac{1}{1 + \lambda} \times (1-z)^\frac{\lambda}{1+\lambda}
$
Where $r$ is the radius of the egg at height $z$ and $T$ and $\lambda$ are parameters of the model.
Using this, and the following equation for a line in 3d space, passing through a point $(x_0, y_0, z_0)$ in the direction $(x_d, y_d, z_d)$:
$
\frac{x - x_0}{x_d} = \frac{y - y_0}{y_d} = \frac{z - z_0}{z_d} 
$
and the following equation for $r$ in terms of $x$ and $y$
$
r^2 = x^2 + y^2
$
I've been able to substitute out the $x$ and $y$ co-ordinates in terms of $z$, to get the following equality that I think should hold at the intersection.
$
\sqrt{(y_0 + y_d  \frac{z - z_0}{z_d})^2 + (x_0 + x_d * \frac{z - z_0}{z_d})^2} = T \times (1 + z)^\frac{1}{1 + \lambda} \times (1-z)^\frac{\lambda}{1+\lambda})
$
This is where I'm stuck.
If I can solve this, I can find the z coordinates of any intersections, and from that, the full coordinates, and this should be sufficient to render the egg.
 A: The function
$$
  r = T\cdot (1+z)^{1/(1+\lambda)} (1-z)^{\lambda/(1+\lambda)}
$$
seems to be more complicated than necessary.  Another form of two-parameter egg equation is an algebraic curve of total degree three, $x^2 + a_0y^2 + a_1xy^2 = 1$, where
$$
  a_0 = \frac{1 + x_{\textrm{max}}^2}{y_{\textrm{max}}^2}
$$
and
$$
  a_1 = \frac{-2x_{\textrm{max}}}{y_{\textrm{max}}^2}.
$$
The $x_{\textrm{max}}$ parameter is the $x$ value where the egg has its maximum and minimum $y$ value and the $y_{\textrm{max}}$ parameter is this maximum $y$ value.  This equation results from a simple transformation of an ellipse equation (see "From the oval to the egg shape").

Here is an example of this algebraic curve for $x_{\textrm{max}}=-0.4$ and $y_{\textrm{max}} = 0.7$:



The equation of the surface of revolution when this algebraic curve is rotated about the $x$ axis is
$$
  y^2 + z^2 = \frac{1-x^2}{a_0 + a_1x}.
$$
Suppose you define your line equation as $(x,y,z) = \vec{P} + t\vec{v}$, where $\vec{P}$ is a point on the line and $\vec{v}$ is the direction of the line.  Then the intersection of the line and the surface of revolution is a cubic equation in the variable $t$ (which can be solved fairly easily).  This is definitely easier to solve than the equation for the intersection of a line with your original egg equation.  You would need to use some nonlinear equation solver in that case.



Here is a comparison of the algebraic egg curve and the fractional exponent egg curve for $x_{\textrm{max}}=-0.2$ and $y_{\textrm{max}}=0.7$ (a $\lambda$ value of 1.5 gives an $x_{\textrm{max}}$ of -0.2):




The fractional exponent curve above is parameterized by the polar angle $\theta$ to get a better spacing of points near $x=-1$ and $x=1$ so that the shape of the curve is illustrated better at the ends.  You can see that the right end of the fractional exponent egg is too sharp (the curvature is actually infinite) and the left end is too flat.  The paper from "The Auk" didn't consider the goodness of the fit at the ends of the generatrix curve.  The paper only computed a few interior diameters of the generatrix and compared these numbers with width measurements of real eggs.  Here is an image of a real chicken egg to show the shape at the ends:


A: Given $f(x)=r$ as a curve for $-1\le x\le 1$ the revolution surface around the $x$ axis is obtained as $$S\to F(x,y,z) = f(x)-\sqrt{y^2+z^2}=0$$. Defining now a line as
$$
L\to p = p_0 +\mu \vec v,\ \ \ p = (x,y,z),\ \ \ p_0 =(x_0,y_0,z_0),\ \ \ \vec v = (v_x,v_y,v_z),\ \ \ \mu\in\mathbb{R}
$$
the intersection $S\cap L$ is obtained by solving for $\mu$
$$
G(\mu) = F(x_0+\mu v_x,y_0+\mu v_y, z_0+\mu v_z) = 0
$$
Here $G(\mu)$ gives us the clue to find the $\mu^*$ satisfying $G(\mu^*)=0$ when it exists.
Calculating $\mu_m=\arg\max_{\mu}G(\mu)$ if $G(\mu_m) < 0$ then $S,L$ doesn't intersect. When $G(\mu_m) \gt 0$ we have two solutions and we choose the lower valued as the solution. The determination of $\mu_m,\mu^*$ can be done using an iterative procedure like Newton's by solving $G'(\mu)=0,G(\mu)=0$.
Follows a MATHEMATICA script with the needed iterative procedure.
Here is $f(x)$
T = 0.6;
lambda = 0.7;
Plot[T (1 + x)^(1/(1 + lambda)) (1 - x)^(lambda/(1 + lambda)), {x, -1, 1}, PlotRange -> {{-1, 1}, {-1, 1}}, AspectRatio -> 1, PlotStyle -> {Thick, Blue}]


follows the determination of $G(\mu)$ and $\mu_m$
G = T (1 + x)^(1/(1 + lambda)) (1 - x)^(lambda/(1 + lambda)) - Sqrt[y^2 + z^2];
p = {x, y, z};
p0 = {1.2, 1, 1};
v = {-1, -1, -2};
L = p0 + mu v;
Gmu = G /. Thread[p -> L];
solmax = Chop[FindMaximum[Gmu, mu]]

and finally the intersection point determination when $G(\mu_m) \gt 0$.
dGmu = D[Gmu, mu];
deltamu = Gmu/dGmu;
mu0 = 0.5 mu /. solmax[[2]];

For[i = 1, i <= 10, i++,
   deltamu0 = deltamu /. {mu -> mu0};
   mu1 = mu0 - deltamu0;
   If[Abs[deltamu0] < 10^-6, Print[mu1]; Break[]];
   mu0 = mu1;
]

pint = L /. {mu -> mu1};
grL = ParametricPlot3D[L, {mu, -2, 2}];
grpti = Graphics3D[{Red, Sphere[pint, 0.02]}];
gr0 = ContourPlot3D[T (1 + x)^(1/(1 + lambda)) (1 - x)^(lambda/(1 + lambda)) == Sqrt[y^2 + z^2], {x, -1, 1}, {y, -1, 1}, {z, -1, 1}, ContourStyle -> {Yellow, Opacity[0.6]}, Mesh -> None, BoundaryStyle -> None]
Show[gr0, grL, grpti]


Note that for the determination of $\mu_m$ we can use as well an iterative procedure.
