I am programming a solution to Lamb's problem for a point source in 3-D that is outlined in Richards (1979), but I am confused by how the polynomial
$$ (A - 2P^2)^4 - 16X^2Y^2P^4 $$
is cubic in $P^2$. $\\$
The dimensionless quantities $A, T, X, \ \text{and} \ Y$ are defined as: $\\$
\begin{aligned} A &= \frac{\alpha^2}{\beta^2} (\alpha=\text{P-wave speed}, \ \beta=\text{S-wave speed})\\ T &= \frac{\alpha t}{r} (t=\text{time}, \ r=(x_1^2 + x_2^2)^{1/2} \ (x_1, x_2) \ \text{being the location of the receiver at the free surface.)}\\ X &= (1 - P^2)^{1/2} \ \text{or} \ -i(P^2 - 1)^{1/2} \\ Y &= (A - P^2)^{1/2} \ \text{or} \ -i(P^2 - A)^{1/2} \\ \end{aligned} $\\$
Apparently, $$ (A - 2P^2)^4 - 16X^2Y^2P^4 = 16(1 - A)(P^2 - R_1)(P^2 - R_2)(P^2 - R_3) $$ where $R_1, R_2, \ \text{and} \ R_3 $ are roots of the Rayleigh cubic in $P^2$. $\\$
Richards notes that an effective approach would be to find the largest root with $R_3^{1/2} = \frac{\alpha^2}{\gamma^2}$ ($\gamma = $ Rayleigh wave speed) and factorise $P^2 - R_3$ from the cubic and solve a quadratic for $R_1$ and $R_2$. $\\$
Ultimately I need to be able to solve for the roots for any given $A$, so my thought was to simply put it in standard polynomial form and use an algorithm such as numpy.roots([c1, c2, c3, c4]), where c1, c2, c3, and c4 are the coefficients of the standard form polynomial. This makes sense, but the polynomial is quartic in $P^2$ when I expand it. $\\$
I completely understand the ways of solving for the roots of a polynomial, but I just can't understand how this is a cubic polynomial as Richards describes. Clearly I am misunderstanding what Richards means by 'cubic in $P^2$'. $\\$
If anyone could kindly explain what this means that would be very helpful.