Plotting an Ovalloid in Spherical Coordinates

I have to plot a 3D ovalloid using a specific formula in spherical coordinates. The oval can be drawn in 2D polar form, but I am not sure how to draw it in 3D. The full equation for the oval is . How would I convert this equation into 3D form to produce an ovalloid?

function oval(k,P,b)
% k is some eccentricity where 0<k<1
% P is the direction in which the oval's axis points towards
% Normally, P is represented by a vector, but instead I made it a value of
% pi
% b is a value that corresponds with the size and spacing of the oval
x=0:.01:2*pi;
x_dot_P=cos(x).*cos(P)+sin(x).*sin(P); % Find the dot product of x and P, where x and P are both represented by the vectors <cos(x),sin(x)> and <cos(P),sin(P)>

% Spreading the equation of the oval out into multiple lines to make the
% code easier to write

positive_numerator=(b-k^2.*x_dot_P);
positive_root=(b-k^2.*x_dot_P).^2;
norm_P=sqrt((cos(P))^2+(sin(P)).^2);
negative_root=(1-k^2).*(b^2-k^2.*norm_P^2);
negative_numerator=sqrt(positive_root-negative_root);
denominator=1-k^2;
rho_x=(positive_numerator-negative_numerator)./denominator; % This is the final equation for the oval, or rho(x)
clf
polarplot(x,rho_x) % Plot the oval
rlim([0 1.25*b])
hold on

% Now, we will draw a line from the origin to the oval every pi/12 radians. This is not important for the ovalloid

for theta=(P:pi/12:2*pi+P)
    s=theta;
    incident_x_dot_P=cos(s).*cos(P)+sin(s).*sin(P);
    incident_positive_numerator=(b-k^2.*incident_x_dot_P);
    incident_positive_root=(b-k^2.*incident_x_dot_P).^2;
    incident_negative_root=(1-k^2).*(b^2-k^2.*norm_P^2);
    incident_negative_numerator=sqrt(incident_positive_root-incident_negative_root);
    incident_denominator=1-k^2;
    max_incident_rho_x=(incident_positive_numerator-incident_negative_numerator)./incident_denominator;
    incident_rho_x=(0:.01:max_incident_rho_x);
    theta2=s.*(incident_rho_x./incident_rho_x);
    polarplot(theta2,incident_rho_x)
end