```python
%matplotlib inline 
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets
```


```python
N = 1024
```


```python
T = np.logspace(1, 8, N)
```


```python
T_eV = T / 11605.
log_T_eV = np.log(T_eV)
```


```python
k1 = np.exp(-32.71396786375
          + 13.53655609057*log_T_eV
          - 5.739328757388*log_T_eV**2 
          + 1.563154982022*log_T_eV**3
          - 0.2877056004391*log_T_eV**4
          + 0.03482559773736999*log_T_eV**5
          - 0.00263197617559*log_T_eV**6
          + 0.0001119543953861*log_T_eV**7
          - 2.039149852002e-6*log_T_eV**8)
k2 = 4.881357e-6*T**(-1.5)* (1.+1.14813e2 * T**(-0.407))**(-2.242)
```


```python
def sampling_rate(N):
    plt.loglog(T, k1, label='k1')
    plt.loglog(T[::N], k1[::N], label='k1s')
    plt.loglog(T, k2, label='k2')
    plt.loglog(T[::N], k2[::N], label='k2s')
    plt.ylim(1e-19, 1e-7)
    plt.legend(loc="lower left")
```


```python
sampling_rate(16)
```

$dS/dt$

So for atomic hydrogen, we have destruction from:

* H + de
* H + H
* H + He
 
For atomic hydrogen, we have formation from:

* H+ + de


```python
import sympy
```


```python
HI, HII, HeI, HeII, HeIII, de = sympy.sympify(
"HI, HII, HeI, HeII, HeIII, de")
k1, k2, k3, k4, k5, k6, k7, k8 = sympy.sympify(
"k1, k2, k3, k4, k5, k6, k7, k8")
```


```python
r1 = (HI + de), (HII + de + de), k1
r2 = (HII + de), (HI), k2
r3 = (HeI + de), (HeII + de + de), k3
r4 = (HeII + de), (HeI), k4
r5 = (HeII + de), (HeIII + de + de), k5
r6 = (HeIII + de), (HeII), k6
r7 = (HI + HI), (HII + de + HI), k7
r8 = (HI + HeI), (HII + de + HeI), k8
```


```python
all_reactions = [r1, r2, r3, r4, r5, r6, r7, r8]
all_reactions = [r1, r2]
```


```python
def find_formation(species):
    rxns = []
    for r in all_reactions:
        if species in r[1].atoms():
            rxns.append(r)
    return rxns

def find_destruction(species):
    rxns = []
    for r in all_reactions:
        if species in r[0].atoms():
            rxns.append(r)
    return rxns
```


```python
find_formation(HI)
```




    [(HII + de, HI, k2)]




```python
find_destruction(HI)
```




    [(HI + de, HII + 2*de, k1)]




```python

```


```python
def get_rhs(species):
    dSdt = 0
    for lhs, rhs, coeff in find_formation(species):
        term = coeff
        for atom in list(lhs.atoms()):
            term *= atom
        dSdt += term
    for lhs, rhs, coeff in find_destruction(species):
        term = -coeff
        for atom in list(lhs.atoms()):
            term *= atom
        dSdt += term
    return dSdt
```


```python
get_rhs(HI)
```




    -HI*de*k1 + HII*de*k2




```python
get_rhs(HII)
```




    HI*de*k1 - HII*de*k2




```python
sympy.solve(get_rhs(HI))
```




    [{HI: HII*k2/k1}, {de: 0}]




```python
plt.loglog(T, k1/k2, "-")
plt.ylim(1e-7, 10.0)
```


```python

```

0daab527ab208bcb4e18df815ef8e5563513ad6a
ipynb
Jupyter Notebook
week12/examples_week12_2.ipynb
MatthewTurk/astr496-spr2018
382ed6910a81a22b6f30967a075e15cd5d6f7adf
BSD-3-Clause
2018-01-16T17:27:28.000Z
2018-03-01T18:26:32.000Z
week12/examples_week12_2.ipynb
MatthewTurk/astr496-spr2018
382ed6910a81a22b6f30967a075e15cd5d6f7adf
BSD-3-Clause
2018-02-19T21:32:07.000Z
2018-02-20T21:03:45.000Z
week12/examples_week12_2.ipynb
MatthewTurk/astr496-spr2018
382ed6910a81a22b6f30967a075e15cd5d6f7adf
BSD-3-Clause
2018-02-12T01:18:41.000Z
2018-02-19T19:01:14.000Z
Qwen/Qwen-72B
1. YES
2. YES