Editing Apps/RiemannEllipsoidsCompressible (section)

==Proposed Solution Strategy==
===Preamble:===
Specify the three "polar" boundary locations, <math>a, b,</math> and <math>c</math>; specify the <i>direction</i> but not the magnitude of the rotating frame's angular velocity, for example, <math>(\vec{\omega}/\omega) = \hat{k}</math>; pin the central density to the value <math>\rho_c = 1</math>.  Define the following dimensionless density, velocity vector, angular velocity, enthalpy, gravitational potential, and position vector:
<div align="center">
<math>
\rho^* \equiv \frac{\rho}{\rho_c} ; ~~~~~{\vec{v}}^* \equiv \frac{\vec{v}}{[a^2G\rho_c]^{1/2}} ; ~~~~~\omega^* \equiv \frac{\omega}{[G\rho_c]^{1/2}} ;
</math>

<math>
H^* \equiv \frac{H}{[a^2G\rho_c]} ; ~~~~~\Phi^* \equiv \frac{\Phi}{[a^2G\rho_c]} ; ~~~~~{\vec{x}}^* \equiv \frac{\vec{x}}{a} .
</math>
</div>
From here, on, spatial operators are assumed to be in terms of the dimensionless coordinates.

===Step #1:===
Guess a 3D density distribution &#8212; such as a uniform-density ellipsoid, or one of the converged models from Ou (2006) &#8212; that conforms to a selected set of <i>positional</i> boundary conditions, that is, where the density goes to zero along the three principal axes at <math>x=a</math>, <math>y = b</math>, and <math>z = c</math>.  Solve the Poisson equation in order to derive values for <math>\Phi</math> everywhere inside and on the surface of the 3D configuration:
<div align="center">
<math>
\nabla^2 \Phi^* = 4\pi \rho^* .
</math>
</div>

===Step #2:===
Use the continuity equation and the curl of the Euler equation to numerically derive the <i>structure</i> but not the overall magnitude of the velocity flow-field throughout the 3D configuration.  Take advantage of the fact that the direction, <math>(\vec{\omega}/\omega)</math>, has been specified; and assume that the 3D density distribution is known.  Here are the relevant equations:
<div align="center">
<math>
\nabla\cdot(\rho^* {\vec{v}}^*) = 0 ;
</math>

<math>
\nabla\times \biggl[({\vec{v}}^*\cdot \nabla){\vec{v}}^* +2 {\vec{\omega}}^* \times {\vec{v}}^* \biggr] = 0 .
</math>
</div>

The first of these is a scalar equation; the second is a vector equation and it will presumably provide two useful scalar equations (perhaps constraining the two components of <math>{\vec{v}}^*</math> that are perpendicular to <math>\hat{k}</math> ?).  Since the left-hand-side of the second equation is obviously nonlinear in the velocity, we may have to linearize this set of equations and look for small "corrections" <math>\delta\vec{v}</math> to an initial "guess" for the velocity field, such as the flow field in Riemann S-type ellipsoids, which is also the flow-field adopted by Ou (2006). 

===Step #3:===
Take the divergence of the Euler equation and use it to solve for <math>H</math> throughout the configuration, given the structure of the flow-field obtained in Step #2.  Boundary conditions at the three "poles" of the configuration may suffice to uniquely determine <math>\omega</math>, the overall normalization factor for the flow-field <math>\vec\zeta</math> &#8212; hopefully this is analogous to solving for the vorticity parameter <math>\lambda</math> in Ou (2006) &#8212; and the Bernoulli constant (or something equivalent).  The relevant "Poisson"-like equation is:

<div align="center">
<math>
\nabla^2 \biggl[H^* + \Phi^* -\frac{1}{2}(\omega^*)^2 \biggl(\frac{R}{a}\biggr)^2  \biggr] = - \nabla\cdot [({\vec{v}}^*\cdot \nabla){\vec{v}}^* + 2 {\vec{\omega}}^*\times {\vec{v}}^* ] .
</math>
</div>