Space Weather

The Solar Storm of 1859 was caused by a coronal mass ejection that hit the earth. During the storm, Auroras were seen from the arctic all the way down to the Caribbean. Many Telegraph systems, of the time, failed and even shocked their operators. This solar if it happened now would be the cause of a great deal of damage to our power grid, it would damage alot of our satellites, disrupt avionics in aircraft and spacecraft, and disrupt radio communication.
With the ability to predict space weather we can prep our satellites, our energy infrastructure, our aircraft to prevent damage. Through the use of strategically placed satellites and Solar Wind Models like NASA's WSA-ENLIL we can predict space weather up to around 3 days in advanced. WSA-ENLIL is a physics-based prediction model whose roots stem from Parker Wind.

Solar Wind

Solar wind is a stream of charged particles originating from the corona of the sun. The temperature around the corona of the sun is about 1,300,000 K leading the Maxwell-Boltzmann Velocity Distribution.

Picture1
Figure 1: A Plot of the Maxwell-Boltzmann Velocity Distribution of Ions in the Corona Where the typical velocity of solar wind coming out of the Corona is highlighted with a vertical red line

The corona starts at around $2.1\times 10^6$ meters from the center of the sun until $5.36\times 10^9$ meters. With the corresponding escape velocities of $1.12\times 10^7 m/s$ at the start of the corona and $2.23\times 10^5$ at the end of the corona.

Parker Wind

Parker Wind is the Solar Wind model first proposed by Eugene N. Parker in his 1958 paper, "Dynamics of the Interplanetary Gas and Magnetic Fields". On the 4th page of this paper he showed a very interesting differential equation $$NMu\frac{du}{dr}=-\frac{d}{dr}(2NkT)-GNMM_o\frac{1}{r^2}$$
It's derivation is obtained by using the velocity of the solar wind as having only a radial component $\vec{u}=u\hat{r}$. In the equation the pressure is $p$, the density $\rho$, the number of particles $N$, temperature $T$, Gravitational constant $G$, mass of particle $M$, mass of the sun $M_o$, $r$ is the radius, and $k$ is the boltzmann constant.

Derivation

It can be derived from Euler's Fluid Equation $$\frac{D\vec{u}}{Dt}=-\frac{\nabla p}{\rho}+\vec{g}$$ We use the fact that we are only worrying about the radial component and the equation turns to $$\rho u\frac{du}{dr}=-\frac{dp}{dr}+\rho g$$
The density can be calculated using the equation $\rho=NM$, the pressure using $p=2Nkt$, and finally the gravitational acceleration $g = - \frac{GM_o}{r^2}$. By plugging those, you get the above Parker Differential Equation.

Solution

I will be using MATLAB's ODE45 differential equation solver. ODE45 uses an explicit Runge-Kutta (4,5) formula which uses the solution only at the preceding point, $y(t_{n-1})$, to solve the function $y(t_n)$. MATLAB's ODE45 needs the input differenctial equation in the form, for example, for $y''=\frac{A}{B}t_y$ in the form $y_1'=y_2$ and $y_2'=\frac{A}{B}ty_1$.
We use the conservation of particles $\frac{d}{dt}(N(r)u(r)r^2)=0$ with the differential equation to obtain $$Mu\frac{du}{dr}=-ur^2\frac{d}{dr}\big(\frac{2kT}{ur^2}\big)-\frac{GM_oM}{r^2}$$
We use the Ion Thermal Velocity Squared , $v^2=\frac{2kT}{M}$, to further simplify the equation. $$u\frac{du}{dr}=-ur^2\frac{d}{dr}\big(\frac{v^2}{ut^2}\big)-\frac{GM_o}{r^2}$$
We set the system as isothermal so $T$ is constant making $v^2$ constant with respect to $r$ and expanding the middle term. $$u\frac{du}{dr}=\frac{2v^2}{r}+\frac{v^2}{u}\frac{du}{dr}-\frac{GM_o}{r^2}$$ And finally, we solve for $\frac{du}{dr}$. $$\frac{du}{dr}=\bigg(\frac{2v^2}{r}-\frac{GM_o}{r^2}\bigg)\bigg(u-\frac{v^2}{u}\bigg)^{-1}$$
We are now ready to solve using MATLAB, first set constants

k = 1.38064852*10^(-23);  % Boltzmann Constant (m^2kgs^-2K^-1)
T = 1.3*10^6;             % Typical Coronal Temprature (K)
M_p = 1.6726219*10^(-27); % Mass of a Proton (kg)
U = 2*k*T/M_p;            % Ion Thermal Velocity Squared
G = 6.67408*10^(-11);     % Gravitational Constant (m^2kg^-1s^-2)
M_sun = 1.989*10^30;      % Mass of the sun (kg)
r_corona = 5360500000;    % radius from center of sun to outer edge of corona (m)
AU = 149597870700;        % 1 Astronomical Unit in Meters
v_corona = 400000;        % Typical velocity of solar wind at the corona

Then we declare the differential function as an Anonymous Function.

odefun = @(t,y) (2*U/t-G*M_sun/t^2)/(y-U/y);

Finally, we use ode45 to solve

[r,v] = ode45(odefun, [r_corona AU], v_corona);

After Plotting Distance (m) vs Velocity (m/s),
plot
Source Code