# Radio Tutorial 2: Density models and radio burst speeds

### The goal for this tutorial

The goal will be to plot a variety of 1 dimensional electron density models of the solar atmosphere, both theoretical and empirical. The empirical models that we will use (Newkirk and Saito) are often used for estimates of radio burst altitude and speed in the corona. We'll be plotting electrons density as a function of height $n(r)$.

N.B. When I use the term height, altitude or distance I mean heliocnetric distance i.e., 1 R$_{\odot}$ is the solar surface.

### First let's start with the hydrostatic plane parallel model

Start with plotting the following equation from 1-10 solar radii
\begin{equation}
n(r) = n_0\mathrm{exp}\bigg(\frac{-r}{H}\bigg)
\end{equation}

where $n_0$ is the density at the base of the corona and $r$ is the distance into the solar atmosphere. $H$ is known as the scale height. And is given by

\begin{equation}
H = \frac{m_p g_{\odot}}{k_B T}
\end{equation}

where $k_B$ is Boltzmann's constant, $T$ is the temperature, $m_p$ is the proton mass, and $g_{\odot}$ is solar gravity. What is the significance of the scale height?


```python
import numpy as np
import matplotlib.pyplot as plt
```

### First define your constants here


```python
n0 = None       # Electron density at 1 Rsun (cm^-3)
T = None        # MK
kb = None       # Boltzman constant J/K
mp = None       # Proton mass kg
Msun = None     # Solar mass kg
rsun = None     # Solar radius m
G = None        # Universal gravitational constant (m3 kg-1 s-2)
gsun = None     # Solar gravity (m/s/s)
```

### First build a 1D distance grid


```python
# Build a 1D array of distance values in units of Rsun.
# Choose a reasonable resolution e.g., <0.1 Rsun. 
```

### Let's look at a plane-parallel hydrostatic model


```python
# Start with defining the scale height.
# Print its value in Mm. Does this make sense? What is the meaning of this value?
```


```python
# Now use your distance grid and scale height to get the density values, using Equation 1 above. 
```


```python
# And plot the result...
'''
fig, ax = plt.subplots()
ax.plot(radius, npp)
ax.set_xlabel('Heliocentric Distance (R$_{\odot}$)')
ax.set_ylabel('Electron density (cm$^{-3})$')
plt.show()
'''
```




    "\nfig, ax = plt.subplots()\nax.plot(radius, npp)\nax.set_xlabel('Heliocentric Distance (R$_{\\odot}$)')\nax.set_ylabel('Electron density (cm$^{-3})$')\nplt.show()\n"



### Increasing the temperature of the atmosphere

Now let's see what happens when we increase the temperature of the atmosphere. Make a scale height for an atmosphere of 1 MK, 5 MK, and 10 MK. What would you expect to happen? 


```python
# The most convenient way to do this is to build a function that takes your distance grid 
# and temperature as input plots the density model. Fill in the function template below.


fig, ax = plt.subplots()
def plothydrostat(radius, temp, label=None):
    
    # Inside here calulate scale height H
    # Then get the density and plot the result
    
    #ax.plot(radius, n, label=label)
    
    return #ax

'''
plothydrostat(radius, T1, label='1 MK')
plothydrostat(radius, T2, label='5 MK')
ax = plothydrostat(radius, T3, label='10 MK')
ax.set_yscale('log')
ax.set_xlabel('Heliocentric Distance (R$_{\odot}$)')
ax.set_ylabel('Electron density (cm$^{-3})$')
plt.legend()
plt.show()
'''
```

### Now let's plot the spherically symmetric solution


Now we'll do something slightly more realistic. Gravity in the solar atmosphere is of course not a constant, but decreases with increasing altitude. The hydrostatic model that's accounts for this is similar:

\begin{equation}
n = n_0\mathrm{exp}\bigg( \frac{1}{H}\bigg\{\frac{1}{r}-1\bigg\} \bigg)
\end{equation}

where $H$ takes on the slightly different form

\begin{equation}
H = \frac{2 k_B T R_{\odot}^2}{G M_{\odot}m_p} 
\end{equation}


```python
# Define these equations here. Take care with your units, units of Rsun or in meters etc.


# the plot up the result, and also compare to the plane parallel solution you calculated previously.
'''
plt.plot(radius, nss, label='Plane parallel')
plt.plot(radius, npp, label='Spherically symmetric')
plt.yscale('log')
plt.ylim(1e4,1e10)
plt.xlabel('Heliocentric Distance (R$_{\odot}$)')
plt.ylabel('Electron density (cm$^{-3})$')
plt.legend()
plt.show()
'''
```




    "\nplt.plot(radius, nss, label='Plane parallel')\nplt.plot(radius, npp, label='Spherically symmetric')\nplt.yscale('log')\nplt.ylim(1e4,1e10)\nplt.xlabel('Heliocentric Distance (R$_{\\odot}$)')\nplt.ylabel('Electron density (cm$^{-3})$')\nplt.legend()\nplt.show()\n"



### What does this look like in terms of plasma frequency?


```python
# Use the plasma frequency formula from the lectures to convert from density to frequency in MHz
# Then plot as a function of radius.

# What is the radius at which you expect 100 MHz?

'''
plt.plot(radius, fplasma_ss)
plt.xlabel('Heliocentric Distance (R$_{\odot}$)')
plt.ylabel('Plasma Frequency (MHz)')
plt.yscale('log')
plt.ylim(1,300)
plt.show()
'''

```




    "\nplt.plot(radius, fplasma_ss)\nplt.xlabel('Heliocentric Distance (R$_{\\odot}$)')\nplt.ylabel('Plasma Frequency (MHz)')\nplt.yscale('log')\nplt.ylim(1,300)\nplt.show()\n"



### Now let's plot some empirical density models

Radio bursts in solar physics are most often analysed with 1D empircal density models. There are many: Newkirk, Saito, Leblanc, Baumbach Allen, Mann. This list goes on. Each are suitable to different parts of the corona. For example, the Newkirk and Saito can be used for the low corona, in the early stages of eruption typically below a few solar radii (although the Newkirk model is known to be slightly higher than epxected densities). The Leblanc and Mann models is for distances further into the heliosphere.

Here we'll take a look at the Newkirk and Saito models. Firtly the newkirk model

\begin{equation}
n(r) = 4.2\times10^4 \times 10^{4.32/r}
\end{equation}

where $r$ is units of solar radii. Plot this model along side the spherically symmetric model from above and compare them.


```python
# Plot the the Newkirk and spherically symmetric model
  

```

Now plot the  Saito model using the following formula

\begin{equation}
n(r) = c_1r^{-d_1} + c_2r^{-d_2}  
\end{equation}

where $r$ is units of solar radii and $c1=1.38\times10^6$, $c2=1.68\times10^8$, $d1=2.14$, $d2=6.13$.



```python
# Use this formula to plot the Saito, overplotting Newkirk, and the hydrostatic model. 
```

### Now let's get the speed of the shock that caused the type II


```python
# Firstly write a function of frequency to density

def freq2dens(freq):
    ne = None
    return ne # cm^-3
```


```python
import numpy as np

result = np.load('/Users/eoincarley/Desktop/WP4-T1/STELLAR_SSW_tutorials/day2_radio_session/ft_coords.npz', allow_pickle=True)

# Read in the frequency time points from Pearse's tutorial.

#freq = np.array([150e6, 100e6, 40e6, 30e6, 20e6])
#tsec = np.array([0, 120, 200, 360, 520])

# Convert the frequencies to densities.
# density = freq2dens(freq)
list(result.keys())
time = result['arr_0']
```


```python
# Now use one of the density models to map density to radius.
# You can rearrange the Newkirk model to get the density directly. 
# However, with a polynomial like Saito, the most convient would 
# be to interpolate between density and radius
# Hint. Use numpy.interp.

# The results should be an array of heliocentric distances.
```


```python
# Plot the results of radii vs time and fit a line to the data points
# Use scipy.stats.linregress
```

### Now estimate the speed using a Saito model.

Follow the same procedure for the Saito model. Compare the distance values that you get with Saito. Compare the the speed given by the Newkirk model and Saito. Which one would you trust more?


```python

```

0fc6f3d29f9277fa4278919bbaa0ab75cf71499f
ipynb
Jupyter Notebook
day2_radio_session/.ipynb_checkpoints/radio2_questions-checkpoint.ipynb
eoincarley/STELLAR_SSW_tutorials
2940f0452234f6e53feb70dcd99823ae6eaceb40
BSD-2-Clause
day2_radio_session/.ipynb_checkpoints/radio2_questions-checkpoint.ipynb
eoincarley/STELLAR_SSW_tutorials
2940f0452234f6e53feb70dcd99823ae6eaceb40
BSD-2-Clause
day2_radio_session/.ipynb_checkpoints/radio2_questions-checkpoint.ipynb
eoincarley/STELLAR_SSW_tutorials
2940f0452234f6e53feb70dcd99823ae6eaceb40
BSD-2-Clause
Qwen/Qwen-72B
1. YES
2. YES