# AST 341: Stars and Radiation
## homework #1 solutions

I'm using SymPy -- the symbolic math library for python -- to do the algebra here.  The next 2 cells load SymPy for interactive work and define the variables we are going to use (and treat as math symbols)


```python
from sympy import init_session
init_session()
%matplotlib inline
```

    IPython console for SymPy 1.6.2 (Python 3.8.5-64-bit) (ground types: python)
    
    These commands were executed:
    >>> from __future__ import division
    >>> from sympy import *
    >>> x, y, z, t = symbols('x y z t')
    >>> k, m, n = symbols('k m n', integer=True)
    >>> f, g, h = symbols('f g h', cls=Function)
    >>> init_printing()
    
    Documentation can be found at https://docs.sympy.org/1.6.2/
    



```python
rho = symbols('rho', cls=Function)
qc = symbols('q_c')
pc = symbols('p_c')
G = symbols('G')
Mstar, Rstar = symbols('M_\star R_\star')
r = symbols('r')
xi = symbols('xi')
beta = symbols('beta')
```

# 1.

Consider a mass profile in a star:
$$m(r) = 4\pi R^3 \beta \left [ \frac{1}{9} \left (\frac{r}{R} \right )^3 - \frac{1}{18} \left ( \frac{r}{R} \right)^6 \right ]$$
where $\beta$ is a constant

## a. 

What is the density profile, $\rho(r)$, that gives rise to this mass?

First let's plot the function $m(r)$


```python
m = 4*pi*Rstar**3*beta*(Rational(1, 9) * (r/Rstar)**3 - Rational(1, 18) * (r/Rstar)**6)
m
```


```python
# to plot it, we need to make it dimensionless
z = (m/(4*pi*Rstar**3*beta)).subs(r, xi*Rstar)
plot(z, (xi, 0, 1), xlabel=r"$r/R_\star$", ylabel=r"$m(r)/(4\pi R_\star^3\beta)$")
```

Now we get the density as
\begin{equation}
\rho = \frac{1}{4\pi r^2} dM/dr
\end{equation}


```python
rho = m.diff(r, 1)/(4*pi*r**2)
rho = simplify(rho)
rho
```

## b.

What is $\beta$ in terms of the central density?

The central density is just $\rho(0)$


```python
rhoc = rho.subs(r, 0)
rhoc
```

So $\beta = 3 \rho_c$.  We can rewrite the density in terms of $\rho_c$ now


```python
rhoc = symbols("rho_c")
rho = simplify(rho.subs(beta, 3*rhoc))
rho
```

We can also rewrite the mass, $m(r)$, in terms of $\rho_c$:


```python
m = simplify(m.subs(beta, 3*rhoc))
m
```

Notice that the density vanishes at $r = R_\star$.  Now let's plot this.  


```python
z = simplify((rho/rhoc).subs(r, xi*Rstar))
plot(z, (xi, 0, 1), xlabel=r"$r/R_\star$", ylabel=r"$\rho(r)/\rho_c$")
```

## c.

What is the ratio of the central density to the average density, $\rho_c / \bar{\rho}$?

The average density is just the total mass divided by the volume of the star.  First we get the total mass, $M_\star = m(R_\star)$


```python
Mstar = m.subs(r, Rstar)
Mstar
```

and now we compute
\begin{equation}
\bar{\rho} = \frac{M_\star}{(4/3) \pi R_\star^3}
\end{equation}


```python
rhobar = symbols(r"\bar{\rho}")
rhobar = Mstar / (Rational(4,3)*pi*Rstar**3)
rhobar
```

So $\bar{\rho} = \rho_c/2$ and the ratio of $\rho_c/\bar{\rho}$ is 2.

## d.

What is the central pressure of the star for this mass distribution?

Now we can integrate HSE.  We will do this as

$$p = p_c + \int_0^{R_\star} \frac{-Gm(r)}{r^2} \rho dr$$

Here I've written the integral in terms of $r$.
We'll substitute in our expressions for $m(r)$ and $\rho(r)$ in the integrand.


```python
pc = symbols("p_c")
p = pc + integrate(-G*m/r**2*rho, (r, 0, Rstar))
p
```

This was the integral to the surface, so this result is the surface pressure, but we know that $p(R_\star) = 0$, 
so we can enforce that here to find $p_c$


```python
pc = solve(Eq(p, 0), pc)[0]
pc
```

This is $p_c$ in terms of $\rho_c$, but we can substitute in $\rho_c = 2 \bar{\rho}$ to find the central pressure in terms of $M_\star$ and $R_\star$.
(Note to make the code use $M_\star$ and not its equivalent in terms of $\rho_c$ found above, I am defining a new symbol $M$ here that will represent the total mass of the star.)


```python
M = symbols("M")
pc.subs(rhoc, 2*M/(Rational(4,3)*pi*Rstar**3))
```

From this expression, we see that $f = 63/(80\pi)$.

## e.

What is the total gravitational potential energy, $\Omega$, of the star?

We integrate

$$\Omega = -\int_0^{M_\star} \frac{Gm(r)}{r} dm = -\int_0^{R_\star} \frac{G m(r)}{r} 4\pi r^2 \rho dr$$

To do this, we use our expression for $m(r)$ and $\rho(r)$


```python
Omega = symbols("Omega")
Omega = integrate(-G*m*4*pi*r*rho, (r, 0, Rstar))
```


```python
Omega
```

We can rewrite this in terms of mass, by using $\rho_c = 2 \bar{\rho}$ that we found for this model above


```python
Omega = Omega.subs(rhoc, 2*M/(Rational(4,3)*pi*Rstar**3))
Omega
```

So this shows that $\alpha = 153/220$.


```python

```

f4a0fd5353cf22a7d7db72d822939ab3de1c0986
ipynb
Jupyter Notebook
simple_stellar_model/simple-stellar-model.ipynb
zingale/ast341_examples
0a15b9bf0b268b00021c59504eb7a2006b7e5ada
BSD-3-Clause
2020-09-09T15:48:41.000Z
2021-08-09T16:08:51.000Z
simple_stellar_model/simple-stellar-model.ipynb
zingale/ast341_examples
0a15b9bf0b268b00021c59504eb7a2006b7e5ada
BSD-3-Clause
simple_stellar_model/simple-stellar-model.ipynb
zingale/ast341_examples
0a15b9bf0b268b00021c59504eb7a2006b7e5ada
BSD-3-Clause
Qwen/Qwen-72B
1. YES
2. YES