Hello! I established a toy model for spinning star in a paper weaks ago, at undergraduate level, below are prat of this work.
Supposing that a star of mass $M$ consists of incompressible fluid of density $\rho$, and it rotates uniformly around its geologically fixed polar axis at a constant angular velocity of $\omega$, relative to an inertial reference $R$ at rest. To fix a 3D Cartesian coordinate system $R_C[O;\vec{i},\vec{j},\vec{k}]$ and a spherical
one $R_S[O;\vec{e_r},\vec{e_\theta},\vec{e_\varphi}]$, who share the same origin which coincides the geometrical center of the star, and the unit vector $\vec{k}$ forms a right-handed spiral mathematical
structure with the spin. Hence, in equilibrium, an arbitrary mass element with coordinates $P(r,\theta,\varphi)$ receives two body forces:
\begin{equation}
\vec{F_1}=-F_1{\vec{e_r}}
\end{equation}
\begin{equation}
\vec{F_2}=\omega^2r\cos\beta(\cos\beta\vec{e_r}+\sin\beta\vec{e_\theta})
\end{equation}
where $\vec{F_1}$ is gravitation from the remaining mass of the star, $\vec{F_2}$ being the centrifugal inertial force, $\beta$
being the declination, $\beta=\frac{\pi}{2}-\theta$. The equation of equilibrium for the very mass element is$[1]$ :
\begin{equation}
\rho\vec{F}=\nabla p
\end{equation}
Considering that
\begin{equation}
\vec{F}=\vec{F_1}+\vec{F_2}
\end{equation}
and the gradient operator in spherical coordinates is:
\begin{equation}
\nabla=\frac{\partial}{\partial r}\vec{e_r}+\frac{1}{r}\frac{\partial}{\partial\theta}\vec{e_\theta}+\frac{1}{r\sin\theta}\frac{\partial}{\partial\varphi}\vec{e_\varphi}
\end{equation}
Then Eq3 turns into
\begin{subequations}
\begin{equation}
\rho(-F_1+\omega^2r\cos^2\beta)=\frac{\partial p}{\partial r}
\end{equation}
\begin{equation}
\frac{1}{2}\rho\omega^2r\sin2\beta=\frac{1}{r}\frac{\partial
p}{\partial \theta}=-\frac{1}{r}\frac{\partial p}{\partial \beta}
\end{equation}
\begin{equation}
0=\frac{1}{r\sin\theta}\frac{\partial p}{\partial
\varphi}=\frac{\partial p}{\partial \varphi}
\end{equation}
\end{subequations}
These are the componential equations of equilibrium under impressible fluid scenario. Immediately, let's get down to the equation of the stellar surface. \\
The surface is a generalized equipotential one, and the potential energy $U_m$ per unit mass is the resultant of its gravitational potential energy $U_1$ and the centrifugal inertial potential energy $U_2$$[2]$.
\begin{equation}
U_m=U_1+U_2
\end{equation}
Assuming the stellar mass distribution is spherically symmetric, and the gravitational potential energy is set zero for reference in $r=+\infty$, then \emph{approximately}
\begin{equation}
U_1=-\frac{GM}{r}
\end{equation}
As to $U_2$, an alternative approach to $\vec{F_2}$ is:
\begin{equation}
\vec{F_2}=-\nabla U_2=-\frac{\partial U_2}{\partial r}\vec{e_r}-\frac{1}{r}\frac{\partial U_2}{\partial \theta}
\vec {e_\theta}
\end{equation}
Comparison of the coefficients of each orthogonal components in Eq2
and Eq9 lead to
\begin{subequations}
\begin{equation}
\frac{\partial U_2}{\partial r}=-\omega^2r\cos^2\beta
\end{equation}
\begin{equation}
\frac{\partial U_2}{\partial
\theta}=-\frac{1}{2}\omega^2r^2\sin2\beta
\end{equation}
\text{For $\frac{\partial U_2}{\partial \theta}=-\frac{\partial
U_2}{\partial \beta}$, Eq$10b$ also reads}
\begin{equation}
\frac{\partial U_2}{\partial \beta}=\frac{1}{2}\omega^2r^2\sin2\beta
\end{equation}
\end{subequations}
Hence, the proper differential of $U_2$ is
\begin{equation}
\begin{split}
dU_2&=\frac{\partial U_2}{\partial r}dr+\frac{\partial U_2}{\partial
\beta}d\beta\\
&=-\omega^2r\cos^2\beta dr+\frac{1}{2}\omega^2r^2\sin2\beta d\beta\\
&=-d\left(\frac{1}{2}\omega^2r^2\cos^2\beta \right)
\end{split}
\end{equation}
Set the centrifugal potential energy to be zero for reference at
$r=0$ (Noting that $\beta=\frac{\pi}{2}$ corresponds to Riemannian
singularity, and lacks generality), after integral we have
\begin{equation}
U_2=-\frac{1}{2}\omega^2r^2\cos^2\beta
\end{equation}
\begin{equation}
U_m=U_1+U_2=-\frac{GM}{r}-\frac{1}{2}\omega^2r^2\cos^2\beta
\end{equation}
As to Eq13, considering that the zero point for $U_m$, and $U_1$ and $U_2$ are selected individually, the resultant $U_m$ is divergent
for both $r=\infty$ and $r=0$. However, because Eq13 is employed for investigations simply at $r=\bar{R}$ ($\bar{R}$ refers to the mean
value of the stellar's radii) where $U_m$ is conservative, there's no need to re-asset the zero potential point for $U_m$ here. Since
the surface is equipotential, $U_m = \text{constant}=C$, then
\begin{equation}
-\frac{GM}{r}-\frac{1}{2}\omega^2r^2\cos^2\beta=C
\end{equation}
Denote $R_p$ to be the radius between the center and one pole( $R_p$ is observable), then insert the boundary condition $\{\beta=\frac{\pi}{2},r=R_p\}$ into Eq14, we have
\begin{equation}
C=-\frac{GM}{R_p}
\end{equation}
From Eq14 and Eq15, the equation of the surface reads
\begin{equation}
f(r,\beta)=\frac{1}{2}\omega^2r^2\cos^2\beta+\frac{GM}{r}-\frac{GM}{R_p}=0
\end{equation}
or
\begin{equation}
f(r,\theta)=\frac{1}{2}\omega^2r^2\sin^2\theta+\frac{GM}{r}-\frac{GM}{R_p}=0
\end{equation}
The physical background for the geometrical equations of the surface is the equipotential conditions. According to the symmetric
properties, Eq16 or Eq17 refers to the surface of a rotating ellipsoid geometrically, which is generated by the rotation around
the minor axis of the ellipse whose major axis is the diameter of the stellar's equator and minor axis is the diameter between two
opposite poles. In particular, for the equator, $\theta=\frac{\pi}{2}(\beta=0)$, $r=R_E$, then from Eq17 we have
\begin{equation}
\frac{R_E-R_P}{R_P}=\eta=\frac{\omega^2R_E^3}{2GM}
\end{equation}
where $\eta$ is traditionally defined$[3]$ as the flattening of the a rotating ellipsoid. If $\eta$ is slight, the rotating ellipsoid
is also called analogous sphere. Eq18 provides a direct perception of the effect of the horizontal expansion due to spin. The faster
the stellar's spin and the larger the volume($R_E^3$), the more serious the deformation. As for the Earth we live on,
$R_P=6356.766km$, $R_E=6378.160km$, $\omega=\frac{2\pi}{24\times 3600}rad\cdot s^{-1}\approx 7.27\times10^{-5} rad\cdot s^{-1}$,
$G=6.67\times10^{-11}N\cdot m^2\cdot kg^{-2}$, $M=5.796\times 10^{24}kg$$[3]$, hence
\begin{equation}
\eta_e=\frac{\omega^2R_E^3}{2GM}\approx 0.1773\times10^{-2}
\end{equation}
which is about half of the popular value, thus Eq18 is reasonable to a considerable degree. And the flattening of the Earth is small, so it can be treated as a analogous sphere in certain astrophysical problems for accuracy.
You can contact
tianwj1@gmail.com for a texified pdf file with visual formulae.