---
title: Plotting the brightest 200 stars
format:
    html:
        copy-code: true
        embed-resources: true
fig-responsive: true
fig-width: 8
---

Grabbing our libraries
```{python}
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
plt.style.use('seaborn-v0_8-darkgrid')
```

## Solar System
Generating our data, wherein all the density is in the initial bin (the Sun)
```{python}
aus = np.arange(0.1, 5, 0.1)
m_sun = 2E30
d_init = m_sun / (4/3 * np.pi * 0.1**3)
dens = np.zeros(shape=aus.shape)
dens[0] = d_init

ss_data = pd.DataFrame({'R_au': aus, 'density': dens})
ss_data.head()

ss_data.to_csv('L19_solar.csv', index=False)
```

Plotting it up:
```{python}
plt.bar(ss_data.R_au, ss_data.density, width=-0.1, align='edge')
plt.show()
```

Computing the mass and then cumulative mass:
```{python}
ss_data['mass'] = 4/3 * np.pi * ss_data.R_au**3 * ss_data.density
ss_data['cmass'] = ss_data.mass.cumsum()
ss_data.head()
```

Computing the velocity:
```{python}
G = 6.602E-11
au = 1.496E11
ss_data['vel'] = np.sqrt(G * ss_data['cmass'] / (ss_data.R_au * au))
ss_data.head()
```

Graphing it up:
```{python}
plt.plot(ss_data.R_au, ss_data.vel, 'o')
plt.xlabel('Distance away (AU)')
plt.ylabel('Velocity (m/s)')
plt.show()
```

Sanity check: how long for Earth to go around?
```{python}
dist = 2*np.pi * 1 * au
time = dist / ss_data[ss_data.R_au == 1].vel
time / 86400
```


## Galaxy
Generating our data, where the bulge will have a constant density that will then drop off
```{python}
rs = np.arange(10,20000,10) # lyrs
dens_init = 4E30 / (4/3 * np.pi * 10**3)
dens = dens_init/(1+np.exp((rs - 3000)/500))

gal_data = pd.DataFrame({'R_ly': rs, 'density': dens})
gal_data.head()
gal_data.to_csv('L19_galaxy.csv', index=False)
```

Plotting it up?
```{python}
plt.bar(gal_data.R_ly, gal_data.density, width=-10, align='edge')
plt.show()
```

Computing the mass and then cumulative mass:
```{python}
gal_data['mass'] = 4/3 * np.pi * gal_data.R_ly**3 * gal_data.density
gal_data['cmass'] = gal_data.mass.cumsum()
gal_data.head()
```

Computing the velocity:
```{python}
G = 6.602E-11
ly = 9.461E15
gal_data['vel'] = np.sqrt(G * gal_data['cmass'] / (gal_data.R_ly * ly))
gal_data.head()
```

Graphing it up:
```{python}
plt.plot(gal_data.R_ly, gal_data.vel, 'o')
plt.xlabel('Distance away (lyrs)')
plt.ylabel('Velocity (m/s)')
plt.show()
```