Editing SSC/Structure/BiPolytropes/AnalyzeStepFunction (section)

==Even Simpler Core==
<font color="green"><b>STEP 0:</b></font> &nbsp; We construct a finite-difference representation of the initial (unperturbed) equilibrium configuration by dividing the model into 99 radial zones that are equally spaced in <math>0 \le \xi \le \xi_i</math>.  The initial radial coordinate of each zone and the corresponding initial enclosed mass are given, respectively, by the expressions &hellip; 

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>[r_j]_0</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl( \frac{3}{2\pi}\biggr)^{1 / 2} \xi_j \, ,
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; and, &nbsp; &nbsp; &nbsp; 
  <td align="right">
<math>m_j</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl( \frac{2\cdot 3}{\pi } \biggr)^{1/2} \biggl[ \xi^3 \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2} \biggr] \, .
</math>
  </td>
</tr>
</table>

The mass, <math>m_j</math>, will serve as our Lagrangian coordinate; that is, we will perturb the model by modifying the radial location of each shell while fixing the enclosed mass.

<font color="green"><b>STEP 1:</b></font> &nbsp; Guess the eigenvector, <math>{\delta r}_i</math>, remembering that a reasonably good trial eigenfunction for the core is one that has a "parabolic dependence on the radius," namely,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\frac{x}{\alpha_\mathrm{scale}}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
1 - \frac{\xi^2}{15} \, ,
</math>
  </td>
<td align="center">&nbsp; &nbsp; &nbsp; where, &nbsp; &nbsp; &nbsp; 
  <td align="right">
<math>\xi</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl(\frac{2\pi}{3}\biggr)^{1 / 2} r_0 \, .
</math>
  </td>
</tr>
</table>

This means,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>x(r_0)</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\alpha_\mathrm{scale}\biggl[1 - \biggl(\frac{2\pi}{45}\biggr)r_0^2  \biggr]
\, ,
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~ \frac{dx}{dr_0}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
-\alpha_\mathrm{scale}\biggl(\frac{4\pi}{45}\biggr)r_0 
\, .
</math>
  </td>
</tr>
</table>

Once a "guess" for the fractional displacement vector, <math>(\delta r)_j = [x \cdot r_0]_j\, ,</math> has been specified, we recognize that the perturbed location of each radial shell is given by the expression,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>r_j</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
(r_0)_j [ 1 + x_j] \, .
</math>
  </td>
</tr>
</table>

<font color="green"><b>STEP 2:</b></font> &nbsp; Our finite-difference representation of the mass-density at each radial shell in the equilibrium configuration (subscript "0") is,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\biggl[{\bar\rho}_{j-1/2} \biggr]_0</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl[m_j - m_{j-1}\biggr]_0\biggl[\frac{4\pi}{3} (r_j^3 - r_{j-1}^3)  \biggr]_0^{-1} \, .
</math>
  </td>
</tr>
</table>
After perturbing the radial location of each shell &#8212; that is, after setting <math>r_j = (r_0)_j + [x \cdot r_0]_j</math> &#8212; the resulting finite-difference representation of the perturbed mass density of each shell is given by the expression,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>{\bar\rho}_{j-1/2} </math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl[m_j - m_{j-1}\biggr]_0\biggl[\frac{4\pi}{3} (r_j^3 - r_{j-1}^3)  \biggr]^{-1} \, .
</math>
  </td>
</tr>
</table>
(Note that we retain a subscript "0" on the mass, <math>m_j</math>, because it serves as our Lagrangian identifier for each shell.)

<font color="green"><b>STEP 3:</b></font> &nbsp; The pressure can be determined in the equilibrium configuration (subscript "0") and after the perturbation from knowledge of the density and the chosen adiabatic index, <math>\gamma_c = 6/5</math>, via knowledge of the (fixed) specific entropy, namely,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\biggl[ P_{j-1/2} \biggr]_0</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl[ {\bar\rho}_{j-1/2}\biggr]_0^{\gamma_c} \cdot \exp\biggl[ \frac{\mu (\gamma-1)s}{\mathfrak{R}} \biggr]
\, ,
</math>
  </td>
<td align="center">&nbsp;  &nbsp; &nbsp; and, &nbsp; &nbsp; &nbsp;</td>
  <td align="right">
<math>P_{j-1/2}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
\biggl[ {\bar\rho}_{j-1/2}\biggr]^{\gamma_c} \cdot \exp\biggl[ \frac{\mu (\gamma-1)s}{\mathfrak{R}} \biggr]
\, .
</math>
  </td>
</tr>
</table>
The pressure perturbation can therefore be obtained from the simple difference,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\biggl[ \delta P \biggr]_{j-1/2}</math>
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>
P_{j - 1 / 2} - \biggl[ P_{j - 1 / 2}\biggr]_0 \, .
</math>
  </td>
</tr>
</table>

[[File:DisplacementFunctionAModel2.png|350px|right|Displacement Function]]<font color="green"><b>STEP 4:</b></font> &nbsp; If, in the context of our [[#Step_6|above discussion of the perturbed "Euler + Poisson Equations"]], we set LHS = RHS, we obtain,

<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>\biggl[\frac{\omega_i^2}{G\rho_c}\biggr] x r_0</math>
  </td>
  <td align="center">
<math>\approx</math>
  </td>
  <td align="left">
<math>\biggl\{
4\pi r_0^2\biggl[ \frac{d(P_0 + \delta P)}{dm}\biggr]
+ g_0 \biggr\}
- 4xg_0
</math>
  </td>
</tr>

<tr>
  <td align="right">
&nbsp;
  </td>
  <td align="center">
<math>=</math>
  </td>
  <td align="left">
<math>\biggl\{
4\pi r_0^2\biggl[ \frac{d(\delta P)}{dm}\biggr]\biggr\}
- 4xg_0
</math>
  </td>
</tr>

<tr>
  <td align="right">
<math>\Rightarrow ~~~ x</math>
  </td>
  <td align="center">
<math>\approx</math>
  </td>
  <td align="left">
<math>\biggl\{
4\pi r_0^2\biggl[ \frac{d(\delta P)}{dm}\biggr]\biggr\}\biggl\{ \biggl[\frac{\omega_i^2}{G\rho_c}\biggr] r_0 + 4g_0 \biggr\}^{-1} \, .
</math>
  </td>
</tr>
</table>
So for a specified value of the square of the oscillation frequency, <math>[\omega^2/(G\rho_c)]</math> &#8212; same value for all shells &#8212; the strategy should be to (a) guess <math>x_j</math>; (b) evaluate the right-hand-side of this last expression; (c) if the RHS does not equal the "guessed" eigenvector, <math>x_j</math>, then you need to guess a new eigenvector; (d) repeat!

In the figure shown here on the right, the black curve displays the variation with mass, <math>m</math>, of the (parabolic-shaped) displacement function, <math>x/\alpha_\mathrm{scale}</math>, that served as our initial "guess;" while the red dots show how the right-hand side of this last expression varies with <math>m</math> for the case, <math>\omega_i^2 = 0</math>.  They lie almost exactly on top of one another, as hoped/expected.