Editing SSC/Structure/BiPolytropes/AnalyzeStepFunction (section)

==Model Amodel2==

Under the "AModel2FD" tab of an Excel spreadsheet titled, "qAndNuMaxAug21", we have constructed a discrete model of the core of a bipolytrope that has the following properties:

<table border="1" align="center" cellpadding="5">
<tr>
  <td align="center" colspan="2">Model A</td>
</tr>

<tr>
  <td align="center" colspan="1"><math>\mu_e/\mu_c</math></td>
  <td align="right">0.31</td>
</tr>

<tr>
  <td align="center" colspan="1"><math>\nu \equiv \frac{M_\mathrm{core}}{M_\mathrm{tot}}</math></td>
  <td align="right">0.337217</td>
</tr>

<tr>
  <td align="center" colspan="1"><math>q \equiv \frac{r_\mathrm{core}}{R_\mathrm{tot}}</math></td>
  <td align="right">0.0755023</td>
</tr>

<tr>
  <td align="center" colspan="2">Core Properties (at interface)</td>
</tr>

<tr>
  <td align="center" colspan="1"><math>\xi_i</math></td>
  <td align="right">9.0149598</td>
</tr>

<tr>
  <td align="center" colspan="1"><math>\theta_i</math></td>
  <td align="right">0.1886798</td>
</tr>

<tr>
  <td align="center" colspan="1"><math>\frac{d\theta}{d\xi}\biggr|_i</math></td>
  <td align="right">-0.0201845</td>
</tr>

<tr>
  <td align="center" colspan="1"><math>M_\mathrm{core} = 
\biggl( \frac{2\cdot 3}{\pi } \biggr)^{1/2} \biggl[ \xi^3 \biggl( 1 + \frac{1}{3}\xi^2 \biggr)^{-3/2} \biggr]_i</math></td>
  <td align="right">6.8009303</td>
</tr>

<tr>
  <td align="center" colspan="1"><math>r_\mathrm{core} = \biggl(\frac{3}{2\pi}\biggr)^{1 / 2} \xi_i</math></td>
  <td align="right">6.2292317</td>
</tr>
</table>
We have divided this model into 100 equally spaced radial zones (center and surface included), that is, each with the following radial-shell thickness:  &nbsp; <math>\delta \xi \equiv \xi_i/99 = 0.0910602</math>, and <math>\delta r = r_\mathrm{core}/99 = 0.0629215</math>.  Analytic expressions providing the physical properties of our <math>(n_c, n_e) = (5, 1)</math> bipolytropes at each radial shell have been drawn from [[SSC/Structure/BiPolytropes/Analytic51#BiPolytrope_with_nc_=_5_and_ne_=_1|an accompanying chapter]] &#8212; for example &hellip;
<table border="0" cellpadding="8" align="center">

<tr>
  <td align="right"><math>r^*</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>\biggl(\frac{3}{2\pi}\biggr)^{1 / 2}\xi \, ,</math>
  </td>
</tr>

<tr>
  <td align="right"><math>\rho^*</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>\biggl(1 + \frac{\xi^2}{3}\biggr)^{-5 / 2} \, ,</math>
  </td>
</tr>

<tr>
  <td align="right"><math>P^*</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>\biggl(1 + \frac{\xi^2}{3}\biggr)^{-3} \, ,</math>
  </td>
</tr>

<tr>
  <td align="right"><math>m</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>

We note as well that, in equilibrium,
<table border="0" cellpadding="8" align="center">

<tr>
  <td align="right"><math>g_0 \equiv \frac{m}{(r^*)^2}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
- 4\pi (r^*)^2 \cdot \frac{dP^*}{dm}
=
- \frac{1}{\rho^*} \cdot \frac{dP^*}{dr^*}
\, .
</math>
  </td>
</tr>
</table>
Hence, in equilibrium we have,
<table border="0" cellpadding="8" align="center">

<tr>
  <td align="right"><math>\frac{dP^*}{dm}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\frac{1}{4\pi (r^*)^2\rho^*} \biggl( \frac{\xi}{r^*} \biggr) \frac{dP^*}{d\xi}
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\frac{\xi}{4\pi (r^*)^3\rho^*} \cdot \biggl[-3\biggl(1 + \frac{\xi^2}{3}\biggr)^{-4} \frac{2\xi}{3}\biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
\frac{1}{4\pi} \biggl[ \biggl(\frac{3}{2\pi}\biggr)^{1 / 2}\xi \biggr]^{-3} 
\biggl[ \biggl(1 + \frac{\xi^2}{3}\biggr)^{-5 / 2} \biggr]^{-1} 
\biggl[-2 \xi^2 \biggl(1 + \frac{\xi^2}{3}\biggr)^{-4} \biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
-\frac{1}{2\pi} \biggl[ \biggl(\frac{3}{2\pi}\biggr)^{-3 / 2}\xi^{-3} \biggr] 
\biggl[ \biggl(1 + \frac{\xi^2}{3}\biggr)^{5 / 2} \biggr] 
\biggl[\xi^2 \biggl(1 + \frac{\xi^2}{3}\biggr)^{-4} \biggr]
</math>
  </td>
</tr>

<tr>
  <td align="right">&nbsp;</td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
-\biggl(\frac{2^3 \pi^3}{3^3} \cdot \frac{1}{2^2 \pi^2}\biggr)^{1 / 2} 
\biggl(1 + \frac{\xi^2}{3}\biggr)^{-3 / 2}  
\xi^{-1}
</math>
  </td>
</tr>

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

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

<tr>
  <td align="right"><math>4\pi (r^*)^2 \cdot \frac{dP^*}{dm}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
- 2^2\pi \biggl[ \biggl(\frac{3}{2\pi}\biggr)^{1 / 2}\xi \biggr]^2  
\biggl(\frac{2 \pi}{3^3} \biggr)^{1 / 2} 
\biggl(1 + \frac{\xi^2}{3}\biggr)^{-3 / 2}  
\xi^{-1}
</math>
  </td>
</tr>

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

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

<tr>
  <td align="right"><math>4\pi (r^*)^4 \cdot \frac{dP^*}{dm}</math></td>
  <td align="center"><math>=</math></td>
  <td align="left">
<math>
- \biggl[ \biggl(\frac{3}{2\pi}\biggr)^{1 / 2}\xi \biggr]^2 
2 \biggl(\frac{2 \pi}{3} \biggr)^{1 / 2} 
\biggl(1 + \frac{\xi^2}{3}\biggr)^{-3 / 2}  
\xi
</math>
  </td>
</tr>

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

The following table provides a sample of variable values in the central region (shells 0 - 4) and near the interface (shells 95 - 99) of the equilibrium configuration.

<table border="1" align="center" cellpadding="5">
<tr>
  <td align="center" colspan="8">Analytically Determined Physical Properties at Various Radial Shells</td>
</tr>

<tr>
  <td align="center">Shell</td>
  <td align="center"><math>\xi</math></td>
  <td align="center"><math>m = -4\pi (r^*)^4\frac{dP^*}{dm} </math></td>
  <td align="center"><math>r^*</math></td>
  <td align="center"><math>\rho^*</math></td>
  <td align="center"><math>P^*</math></td>
  <td align="center"><math>-~\frac{dP^*}{dm}</math></td>
  <td align="center"><math>-4\pi (r^*)^2\frac{dP^*}{dm} = g_0</math></td>
</tr>

<tr>
  <td align="center">0</td>
  <td align="right">0.000000</td>
  <td align="right">0.000000</td>
  <td align="right">0.000000</td>
  <td align="right">1.000000</td>
  <td align="right">1.000000</td>
  <td align="center"><math>\infty</math></td>
  <td align="right">0.000000</td>
</tr>

<tr>
  <td align="center">1</td>
  <td align="right">0.091060</td>
  <td align="right">0.001039</td>
  <td align="right">0.062922</td>
  <td align="right">0.993123</td>
  <td align="right">0.991754</td>
  <td align="right">5.275715</td>
  <td align="right">0.262476</td>
</tr>

<tr>
  <td align="center">2</td>
  <td align="right">0.182120</td>
  <td align="right">0.008211</td>
  <td align="right">0.125843</td>
  <td align="right">0.972886</td>
  <td align="right">0.967552</td>
  <td align="right">2.605474</td>
  <td align="right">0.518508</td>
</tr>

<tr>
  <td align="center">3</td>
  <td align="right">0.273181</td>
  <td align="right">0.027155</td>
  <td align="right">0.188765</td>
  <td align="right">0.940420</td>
  <td align="right">0.928937</td>
  <td align="right">1.701968</td>
  <td align="right">0.762083</td>
</tr>

<tr>
  <td align="center">4</td>
  <td align="right">0.364241</td>
  <td align="right">0.062586</td>
  <td align="right">0.251686</td>
  <td align="right">0.897462</td>
  <td align="right">0.878252</td>
  <td align="right">1.241164</td>
  <td align="right">0.988001</td>
</tr>

<tr>
  <td align="center" colspan="8"><font size="+1"><b>&vellip;</b></td>
</tr>

<tr>
  <td align="center">95</td>
  <td align="right">8.650719</td>
  <td align="right">6.769823</td>
  <td align="right">5.977546</td>
  <td align="right">0.000292</td>
  <td align="right">0.000057</td>
  <td align="right">0.000422</td>
  <td align="right">0.189466</td>
</tr>

<tr>
  <td align="center">96</td>
  <td align="right">8.741779</td>
  <td align="right">6.777942</td>
  <td align="right">6.040467</td>
  <td align="right">0.000277</td>
  <td align="right">0.000054</td>
  <td align="right">0.000405</td>
  <td align="right">0.185762</td>
</tr>

<tr>
  <td align="center">97</td>
  <td align="right">8.832839</td>
  <td align="right">6.785828</td>
  <td align="right">6.103389</td>
  <td align="right">0.000264</td>
  <td align="right">0.000051</td>
  <td align="right">0.000389</td>
  <td align="right">0.182163</td>
</tr>

<tr>
  <td align="center">98</td>
  <td align="right">8.923900</td>
  <td align="right">6.793488</td>
  <td align="right">6.166310</td>
  <td align="right">0.000251</td>
  <td align="right">0.000048</td>
  <td align="right">0.000374</td>
  <td align="right">0.178666</td>
</tr>

<tr>
  <td align="center">99</td>
  <td align="right">9.014960</td>
  <td align="right">6.800930</td>
  <td align="right">6.229232</td>
  <td align="right">0.000239</td>
  <td align="right">0.000045</td>
  <td align="right">0.000359</td>
  <td align="right">0.175267</td>
</tr>
</table>

This last expression is precisely the same (in magnitude, but opposite in sign) as the expression that we have presented for <math>m</math>.  If we therefore return to <font color="red>STEP 6</font>, we appreciate that the right-hand side of the "eigenvector" expression,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>+\biggl[ \frac{\omega^2}{G\rho_c} \biggr] \biggl(\frac{\delta r^*}{r^*}\biggr) (r^*)</math>
  </td>
  <td align="center">
<math>\approx</math>
  </td>
  <td align="left">
<math>
\biggl[ 4\pi (r^*)^2 \cdot \frac{dP^*}{dm} + \frac{m}{(r^*)^2}\biggr] \, ,
</math>
  </td>
</tr>
</table>
goes to zero at every radial shell if the configuration is in equilibrium.

<table border="1" align="center" cellpadding="10" width="80%"><tr><td align="left">
<font color="red">ASIDE:</font>&nbsp; Next, we will focus on building a discrete, finite-difference representation of each equilibrium model, as well as of this "eigenvector" expression.  In doing so, we will immediately find that the two terms on the right-hand-side ''do not'' exactly sum to zero, even for an equilibrium configuration. We should nevertheless try to construct a finite-difference representation for which the two terms cancel to a relatively high degree of precision.   We have considered rewriting this "eigenvector" expression in the form,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>+\biggl[ \frac{\omega^2}{G\rho_c} \biggr] \biggl(\frac{\delta r^*}{r^*}\biggr) (r^*)^3</math>
  </td>
  <td align="center">
<math>\approx</math>
  </td>
  <td align="left">
<math>
\biggl[ 4\pi (r^*)^4 \cdot \frac{dP^*}{dm} + m \biggr] \, ,
</math>
  </td>
</tr>
</table>

because this form isolates the Lagrangian mass, <math>m</math>, which is time-invariant.  But as the numbers in the third column of the above table illustrate, the terms on the right-hand-side vary by several orders of magnitude as we move from shell 1 to shell 100.  If it is written instead in the form that we have initially suggested, namely,
<table border="0" cellpadding="5" align="center">

<tr>
  <td align="right">
<math>+\biggl[ \frac{\omega^2}{G\rho_c} \biggr] \biggl(\frac{\delta r^*}{r^*}\biggr) (r^*)</math>
  </td>
  <td align="center">
<math>\approx</math>
  </td>
  <td align="left">
<math>
\biggl[ 4\pi (r^*)^2 \cdot \frac{dP^*}{dm} + \frac{m}{(r^*)^2}\biggr] \, ,
</math>
  </td>
</tr>
</table>
then, as is illustrated by the numbers in the last column of the table and by the following figure, the terms on the right-hand-side vary by no more than one order of magnitude as we move from the center to the surface of the configuration.
<div align="center">[[File:G0variationBipolytrope.png|450px|Variation of <math>g_0</math> from center to interface]]</div>  This choice should facilitate cancellation to a higher degree of precision in our finite-difference-based model.
</td></tr></table>