1#!/usr/bin/env python3
2"""
3Parker Solar Wind Model
4
5This module implements the Parker solar wind solution, which explains how
6the solar corona expands to create the solar wind that fills the heliosphere.
7
8Key features:
9- Isothermal Parker solution with critical point
10- Transonic flow: subsonic → supersonic transition at r_c
11- Multiple solution branches: breeze, wind, accretion
12- Comparison with observations
13- Extension to non-isothermal case
14
15The critical point occurs where the flow velocity equals the local sound speed.
16
17Author: MHD Course Examples
18License: MIT
19"""
20
21import numpy as np
22import matplotlib.pyplot as plt
23from scipy.integrate import odeint
24from scipy.optimize import fsolve
25
26
27class ParkerSolarWind:
28 """
29 Parker solar wind model solver.
30
31 Attributes:
32 M_sun (float): Solar mass (kg)
33 R_sun (float): Solar radius (m)
34 T_corona (float): Coronal temperature (K)
35 n_base (float): Base density (m^-3)
36 gamma (float): Polytropic index
37 """
38
39 def __init__(self, M_sun=1.989e30, R_sun=6.96e8,
40 T_corona=1.5e6, n_base=1e14, gamma=1.0):
41 """
42 Initialize Parker solar wind model.
43
44 Parameters:
45 M_sun (float): Solar mass (kg)
46 R_sun (float): Solar radius (m)
47 T_corona (float): Coronal temperature (K)
48 n_base (float): Base number density (m^-3)
49 gamma (float): Polytropic index (1.0 = isothermal)
50 """
51 self.M_sun = M_sun
52 self.R_sun = R_sun
53 self.T_corona = T_corona
54 self.n_base = n_base
55 self.gamma = gamma
56
57 # Constants
58 self.G = 6.674e-11 # Gravitational constant
59 self.k_B = 1.38e-23 # Boltzmann constant
60 self.m_p = 1.67e-27 # Proton mass
61
62 # Sound speed (isothermal)
63 self.c_s = np.sqrt(2 * self.k_B * T_corona / self.m_p)
64
65 # Critical radius (sonic point)
66 self.r_c = self.G * M_sun / (2 * self.c_s**2)
67
68 # Escape velocity at base
69 self.v_esc = np.sqrt(2 * self.G * M_sun / R_sun)
70
71 def parker_equation_rhs(self, v, r):
72 """
73 Right-hand side of Parker wind equation.
74
75 dv/dr = (2c_s²/r - GM/r²) / (v - c_s²/v)
76
77 Parameters:
78 v (float): Velocity (m/s)
79 r (float): Radius (m)
80
81 Returns:
82 float: dv/dr
83 """
84 G = self.G
85 M = self.M_sun
86 c_s = self.c_s
87
88 # Numerator
89 numerator = 2 * c_s**2 / r - G * M / r**2
90
91 # Denominator
92 denominator = v - c_s**2 / v
93
94 # Avoid division by zero near sonic point
95 if abs(denominator) < 1e-10:
96 return 0.0
97
98 dv_dr = numerator / denominator
99
100 return dv_dr
101
102 def solve_parker_wind(self, r_array, v_init, solution_type='wind'):
103 """
104 Solve Parker wind equation.
105
106 Parameters:
107 r_array (ndarray): Radial grid
108 v_init (float): Initial velocity at r_array[0]
109 solution_type (str): 'wind', 'breeze', or 'accretion'
110
111 Returns:
112 ndarray: Velocity solution
113 """
114 # For wind solution: start just above critical point
115 # For breeze: start below critical point and integrate inward
116
117 if solution_type == 'wind':
118 # Start slightly above sonic point, integrate outward
119 r_start = self.r_c * 1.01
120 v_start = self.c_s * 1.01
121
122 r_integrate = r_array[r_array >= r_start]
123 v_solution = odeint(lambda v, r: self.parker_equation_rhs(v[0], r),
124 [v_start], r_integrate)[:, 0]
125
126 # Pad with NaN for r < r_c
127 v_full = np.full_like(r_array, np.nan)
128 v_full[r_array >= r_start] = v_solution
129
130 return v_full
131
132 elif solution_type == 'breeze':
133 # Start at base, integrate outward (stays subsonic)
134 v_solution = odeint(lambda v, r: self.parker_equation_rhs(v[0], r),
135 [v_init], r_array)[:, 0]
136 return v_solution
137
138 elif solution_type == 'accretion':
139 # Integrate inward from infinity (supersonic infall)
140 r_start = r_array[-1]
141 v_start = -self.c_s * 2.0 # Negative = inward
142
143 r_integrate = r_array[::-1]
144 v_solution = odeint(lambda v, r: self.parker_equation_rhs(v[0], r),
145 [v_start], r_integrate)[:, 0]
146
147 return v_solution[::-1]
148
149 return np.zeros_like(r_array)
150
151 def density_from_continuity(self, r, v, n0, r0):
152 """
153 Compute density from continuity equation.
154
155 ρ(r) v(r) r² = ρ₀ v₀ r₀²
156
157 Parameters:
158 r (float or ndarray): Radius
159 v (float or ndarray): Velocity
160 n0 (float): Base density
161 r0 (float): Base radius
162
163 Returns:
164 float or ndarray: Density
165 """
166 # Find velocity at r0
167 idx = np.argmin(np.abs(r - r0))
168 v0 = v[idx]
169
170 # Continuity
171 n = n0 * (v0 / v) * (r0 / r)**2
172
173 return n
174
175 def temperature_profile(self, r, T0, r0):
176 """
177 Temperature profile for polytropic case.
178
179 T(r) = T₀ (ρ/ρ₀)^(γ-1)
180
181 Parameters:
182 r (ndarray): Radius
183 T0 (float): Base temperature
184 r0 (float): Base radius
185
186 Returns:
187 ndarray: Temperature
188 """
189 if self.gamma == 1.0:
190 # Isothermal
191 return np.full_like(r, T0)
192 else:
193 # Polytropic (simplified)
194 return T0 * (r0 / r)**(2 * (self.gamma - 1))
195
196 def plot_solutions(self):
197 """
198 Plot Parker wind solutions: wind, breeze, accretion.
199 """
200 # Radial grid from 1 R_sun to 100 R_sun
201 r_array = np.linspace(self.R_sun, 100 * self.R_sun, 1000)
202
203 # Solve for different solution types
204 v_wind = self.solve_parker_wind(r_array, 0, 'wind')
205 v_breeze = self.solve_parker_wind(r_array, 1e3, 'breeze')
206 v_accretion = self.solve_parker_wind(r_array, 0, 'accretion')
207
208 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
209
210 # Velocity solutions
211 ax1.plot(r_array / self.R_sun, v_wind / 1e3, 'b-', linewidth=2,
212 label='Wind (transonic)')
213 ax1.plot(r_array / self.R_sun, v_breeze / 1e3, 'g--', linewidth=2,
214 label='Breeze (subsonic)')
215 ax1.plot(r_array / self.R_sun, np.abs(v_accretion) / 1e3, 'r:',
216 linewidth=2, label='Accretion (infall)')
217
218 # Sound speed
219 ax1.axhline(y=self.c_s / 1e3, color='k', linestyle='--', alpha=0.5,
220 label=f'Sound speed ({self.c_s/1e3:.0f} km/s)')
221
222 # Critical point
223 ax1.axvline(x=self.r_c / self.R_sun, color='orange', linestyle='--',
224 alpha=0.5, label=f'Critical radius ({self.r_c/self.R_sun:.1f} R☉)')
225
226 ax1.set_xlabel(r'Radius ($R_\odot$)', fontsize=12)
227 ax1.set_ylabel('Velocity (km/s)', fontsize=12)
228 ax1.set_title('Parker Solar Wind Solutions', fontsize=14)
229 ax1.set_xlim([1, 100])
230 ax1.set_ylim([0, 1000])
231 ax1.grid(True, alpha=0.3)
232 ax1.legend(fontsize=10)
233
234 # Mach number
235 M_wind = v_wind / self.c_s
236 M_breeze = v_breeze / self.c_s
237
238 ax2.semilogy(r_array / self.R_sun, M_wind, 'b-', linewidth=2,
239 label='Wind')
240 ax2.semilogy(r_array / self.R_sun, M_breeze, 'g--', linewidth=2,
241 label='Breeze')
242 ax2.axhline(y=1, color='k', linestyle='--', alpha=0.5,
243 label='Sonic point (M=1)')
244 ax2.axvline(x=self.r_c / self.R_sun, color='orange', linestyle='--',
245 alpha=0.5)
246
247 ax2.set_xlabel(r'Radius ($R_\odot$)', fontsize=12)
248 ax2.set_ylabel('Mach Number', fontsize=12)
249 ax2.set_title('Mach Number vs Radius', fontsize=14)
250 ax2.set_xlim([1, 100])
251 ax2.grid(True, alpha=0.3)
252 ax2.legend(fontsize=10)
253
254 plt.tight_layout()
255 return fig
256
257 def plot_with_observations(self):
258 """
259 Plot Parker wind solution with observational data points.
260 """
261 r_array = np.linspace(self.R_sun, 100 * self.R_sun, 1000)
262 v_wind = self.solve_parker_wind(r_array, 0, 'wind')
263
264 # Compute density and temperature
265 n_wind = self.density_from_continuity(r_array, v_wind,
266 self.n_base, self.R_sun)
267 T_wind = self.temperature_profile(r_array, self.T_corona, self.R_sun)
268
269 fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(14, 10))
270
271 # Velocity
272 ax1.plot(r_array / self.R_sun, v_wind / 1e3, 'b-', linewidth=2,
273 label='Parker model')
274
275 # Observational points (typical values)
276 r_obs = np.array([1, 10, 50, 100])
277 v_obs = np.array([0, 200, 350, 400])
278 ax1.plot(r_obs, v_obs, 'ro', markersize=8, label='Observations')
279
280 ax1.set_xlabel(r'Radius ($R_\odot$)', fontsize=11)
281 ax1.set_ylabel('Velocity (km/s)', fontsize=11)
282 ax1.set_title('Solar Wind Velocity', fontsize=13)
283 ax1.grid(True, alpha=0.3)
284 ax1.legend(fontsize=10)
285
286 # Density
287 ax2.semilogy(r_array / self.R_sun, n_wind, 'b-', linewidth=2)
288 ax2.set_xlabel(r'Radius ($R_\odot$)', fontsize=11)
289 ax2.set_ylabel(r'Density (m$^{-3}$)', fontsize=11)
290 ax2.set_title('Number Density', fontsize=13)
291 ax2.grid(True, alpha=0.3)
292
293 # Temperature
294 ax3.plot(r_array / self.R_sun, T_wind / 1e6, 'b-', linewidth=2)
295 ax3.set_xlabel(r'Radius ($R_\odot$)', fontsize=11)
296 ax3.set_ylabel('Temperature (MK)', fontsize=11)
297 ax3.set_title('Temperature Profile', fontsize=13)
298 ax3.grid(True, alpha=0.3)
299
300 # Mass flux
301 mass_flux = n_wind * self.m_p * v_wind * 4 * np.pi * r_array**2
302 ax4.plot(r_array / self.R_sun, mass_flux / 1e9, 'b-', linewidth=2)
303 ax4.set_xlabel(r'Radius ($R_\odot$)', fontsize=11)
304 ax4.set_ylabel(r'Mass Flux ($10^9$ kg/s)', fontsize=11)
305 ax4.set_title('Mass Loss Rate (should be constant)', fontsize=13)
306 ax4.grid(True, alpha=0.3)
307
308 plt.tight_layout()
309 return fig
310
311
312def main():
313 """
314 Main function demonstrating Parker solar wind model.
315 """
316 print("=" * 60)
317 print("Parker Solar Wind Model")
318 print("=" * 60)
319
320 # Create Parker wind model
321 parker = ParkerSolarWind(
322 M_sun=1.989e30,
323 R_sun=6.96e8,
324 T_corona=1.5e6, # 1.5 MK
325 n_base=1e14, # 10^8 cm^-3
326 gamma=1.0 # Isothermal
327 )
328
329 print(f"\nParameters:")
330 print(f" Coronal temperature: {parker.T_corona/1e6:.1f} MK")
331 print(f" Sound speed: {parker.c_s/1e3:.0f} km/s")
332 print(f" Critical radius: {parker.r_c/parker.R_sun:.2f} R☉")
333 print(f" Escape velocity at R☉: {parker.v_esc/1e3:.0f} km/s")
334
335 # Plot solutions
336 parker.plot_solutions()
337 parker.plot_with_observations()
338
339 plt.savefig('/tmp/parker_solar_wind.png', dpi=150, bbox_inches='tight')
340 print("\nPlots saved to /tmp/parker_solar_wind.png")
341
342 plt.show()
343
344
345if __name__ == "__main__":
346 main()