1#!/usr/bin/env python3
2"""
3Simplified 2D Magnetosphere Model
4
5This module creates a simplified 2D magnetosphere model showing the interaction
6between Earth's dipole magnetic field and the solar wind.
7
8Key features:
9- Dipole field inside magnetopause
10- Uniform solar wind flow (potential flow approximation)
11- Chapman-Ferraro magnetopause boundary (pressure balance)
12- Visualization of field lines, stagnation point, and cusps
13
14The magnetopause is where magnetic pressure balances solar wind ram pressure:
15 B²/(2μ₀) = ½ρ_sw v_sw²
16
17Author: MHD Course Examples
18License: MIT
19"""
20
21import numpy as np
22import matplotlib.pyplot as plt
23
24
25class Magnetosphere2D:
26 """
27 2D magnetosphere model with dipole field and solar wind.
28
29 Attributes:
30 B_dipole (float): Dipole moment strength (T·m³)
31 R_planet (float): Planet radius (m)
32 v_sw (float): Solar wind velocity (m/s)
33 n_sw (float): Solar wind density (m^-3)
34 rho_sw (float): Solar wind mass density (kg/m³)
35 """
36
37 def __init__(self, B_dipole=3.1e-5, R_planet=6.4e6,
38 v_sw=4e5, n_sw=5e6):
39 """
40 Initialize magnetosphere model.
41
42 Parameters:
43 B_dipole (float): Dipole field at equator at R_planet (T)
44 R_planet (float): Planet radius (m)
45 v_sw (float): Solar wind velocity (m/s)
46 n_sw (float): Solar wind number density (m^-3)
47 """
48 self.B_dipole = B_dipole
49 self.R_planet = R_planet
50 self.v_sw = v_sw
51 self.n_sw = n_sw
52
53 # Constants
54 self.mu_0 = 4 * np.pi * 1e-7
55 self.m_p = 1.67e-27
56
57 # Solar wind mass density
58 self.rho_sw = n_sw * self.m_p
59
60 # Dipole moment: M = B_eq * R³
61 self.M_dipole = B_dipole * R_planet**3
62
63 # Compute magnetopause standoff distance
64 self.R_mp = self.calculate_magnetopause_distance()
65
66 def calculate_magnetopause_distance(self):
67 """
68 Calculate magnetopause standoff distance from pressure balance.
69
70 B²/(2μ₀) = ½ρ_sw v_sw²
71
72 At subsolar point: B = M/(R_mp)³
73
74 Returns:
75 float: Magnetopause distance (m)
76 """
77 # Dynamic pressure
78 P_ram = 0.5 * self.rho_sw * self.v_sw**2
79
80 # Magnetic pressure at distance R: P_mag = B²/(2μ₀) = M²/(2μ₀ R⁶)
81 # Solve: M²/(2μ₀ R_mp⁶) = P_ram
82
83 R_mp = (self.M_dipole**2 / (2 * self.mu_0 * P_ram))**(1/6)
84
85 return R_mp
86
87 def dipole_field(self, x, y):
88 """
89 Compute dipole magnetic field components.
90
91 B = (μ₀/4π) × (3(m·r)r/r⁵ - m/r³)
92
93 For dipole in z-direction in x-y plane:
94 B_x = 3μ₀M xy / (4π r⁵)
95 B_y = μ₀M(3y² - r²) / (4π r⁵)
96
97 Parameters:
98 x, y (ndarray): Coordinates (m)
99
100 Returns:
101 tuple: (B_x, B_y) field components (T)
102 """
103 r = np.sqrt(x**2 + y**2)
104
105 # Avoid singularity at origin
106 r = np.maximum(r, 0.1 * self.R_planet)
107
108 # Dipole field (simplified 2D version)
109 # In x-y plane with dipole in z-direction
110 B_x = (3 * self.M_dipole * x * y) / (r**5)
111 B_y = self.M_dipole * (3 * y**2 - r**2) / (r**5)
112
113 return B_x, B_y
114
115 def is_inside_magnetopause(self, x, y):
116 """
117 Check if point is inside magnetopause.
118
119 Magnetopause shape: Chapman-Ferraro approximation
120 r(θ) ≈ R_mp (2/(1+cos(θ)))^(1/6)
121
122 Parameters:
123 x, y (ndarray): Coordinates
124
125 Returns:
126 ndarray: Boolean mask
127 """
128 r = np.sqrt(x**2 + y**2)
129 theta = np.arctan2(y, x)
130
131 # Chapman-Ferraro shape (approximate)
132 r_mp_theta = self.R_mp * (2 / (1 + np.cos(theta)))**(1/6)
133
134 # Points inside magnetopause
135 inside = r < r_mp_theta
136
137 return inside
138
139 def magnetopause_boundary(self, theta_array):
140 """
141 Compute magnetopause boundary shape.
142
143 Parameters:
144 theta_array (ndarray): Angles from 0 to 2π
145
146 Returns:
147 tuple: (x, y) coordinates of boundary
148 """
149 # Chapman-Ferraro model
150 r_mp = self.R_mp * (2 / (1 + np.cos(theta_array)))**(1/6)
151
152 x = r_mp * np.cos(theta_array)
153 y = r_mp * np.sin(theta_array)
154
155 return x, y
156
157 def plot_magnetosphere(self):
158 """
159 Plot 2D magnetosphere structure.
160 """
161 fig, ax = plt.subplots(figsize=(12, 10))
162
163 # Grid
164 x_max = 15 * self.R_planet
165 y_max = 15 * self.R_planet
166 nx, ny = 150, 150
167
168 x = np.linspace(-5 * self.R_planet, x_max, nx)
169 y = np.linspace(-y_max, y_max, ny)
170 X, Y = np.meshgrid(x, y)
171
172 # Compute magnetic field
173 B_x, B_y = self.dipole_field(X, Y)
174 B_mag = np.sqrt(B_x**2 + B_y**2)
175
176 # Mask outside magnetopause
177 inside = self.is_inside_magnetopause(X, Y)
178
179 # Plot field magnitude (inside magnetopause only)
180 B_plot = np.where(inside, B_mag, np.nan)
181 im = ax.contourf(X / self.R_planet, Y / self.R_planet,
182 np.log10(B_plot * 1e9), levels=20,
183 cmap='plasma', alpha=0.6)
184 plt.colorbar(im, ax=ax, label=r'$\log_{10}(B)$ [nT]')
185
186 # Plot field lines (streamplot)
187 # Sample on coarser grid for clarity
188 skip = 5
189 ax.streamplot(X[::skip, ::skip] / self.R_planet,
190 Y[::skip, ::skip] / self.R_planet,
191 B_x[::skip, ::skip],
192 B_y[::skip, ::skip],
193 color='white', linewidth=0.5, density=1.0,
194 arrowsize=0.8)
195
196 # Plot magnetopause
197 theta = np.linspace(0, 2*np.pi, 200)
198 x_mp, y_mp = self.magnetopause_boundary(theta)
199 ax.plot(x_mp / self.R_planet, y_mp / self.R_planet,
200 'r-', linewidth=2.5, label='Magnetopause')
201
202 # Plot Earth
203 earth = plt.Circle((0, 0), 1, color='blue', alpha=0.8,
204 label='Planet')
205 ax.add_patch(earth)
206
207 # Mark stagnation point (subsolar point)
208 ax.plot(self.R_mp / self.R_planet, 0, 'yo', markersize=12,
209 label=f'Stagnation point ({self.R_mp/self.R_planet:.1f} R)')
210
211 # Mark magnetic cusps (approximate)
212 cusp_angle = np.pi / 3 # ~60 degrees
213 theta_cusp = np.array([cusp_angle, -cusp_angle])
214 x_cusp, y_cusp = self.magnetopause_boundary(theta_cusp)
215 ax.plot(x_cusp / self.R_planet, y_cusp / self.R_planet,
216 'go', markersize=10, label='Cusps')
217
218 # Solar wind direction
219 ax.arrow(-4, 12, 3, 0, head_width=1, head_length=0.5,
220 fc='orange', ec='orange', linewidth=2)
221 ax.text(-2.5, 13, 'Solar Wind', fontsize=12, color='orange',
222 weight='bold')
223
224 ax.set_xlabel(r'X ($R_{\rm planet}$)', fontsize=12)
225 ax.set_ylabel(r'Y ($R_{\rm planet}$)', fontsize=12)
226 ax.set_title('2D Magnetosphere Model\n(Dipole + Solar Wind)',
227 fontsize=14, weight='bold')
228 ax.set_aspect('equal')
229 ax.grid(True, alpha=0.3)
230 ax.legend(loc='upper right', fontsize=10)
231 ax.set_xlim([-5, 15])
232 ax.set_ylim([-15, 15])
233
234 plt.tight_layout()
235 return fig
236
237 def plot_pressure_profiles(self):
238 """
239 Plot magnetic and ram pressure along x-axis.
240 """
241 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(10, 8))
242
243 # Array along x-axis (y=0)
244 x_array = np.linspace(0.5 * self.R_planet, 20 * self.R_planet, 200)
245 y_array = np.zeros_like(x_array)
246
247 # Magnetic field and pressure
248 B_x, B_y = self.dipole_field(x_array, y_array)
249 B_mag = np.sqrt(B_x**2 + B_y**2)
250 P_mag = B_mag**2 / (2 * self.mu_0)
251
252 # Solar wind ram pressure
253 P_ram = 0.5 * self.rho_sw * self.v_sw**2
254
255 # Plot pressures
256 ax1.loglog(x_array / self.R_planet, P_mag * 1e9, 'b-',
257 linewidth=2, label='Magnetic pressure')
258 ax1.axhline(y=P_ram * 1e9, color='r', linestyle='--',
259 linewidth=2, label='Solar wind ram pressure')
260 ax1.axvline(x=self.R_mp / self.R_planet, color='g',
261 linestyle='--', alpha=0.7, label='Magnetopause')
262
263 ax1.set_xlabel(r'Distance ($R_{\rm planet}$)', fontsize=11)
264 ax1.set_ylabel('Pressure (nPa)', fontsize=11)
265 ax1.set_title('Pressure Balance along X-axis', fontsize=13)
266 ax1.grid(True, alpha=0.3)
267 ax1.legend(fontsize=10)
268
269 # Plot beta = P_gas / P_mag
270 # Assuming thermal pressure ~ ram pressure
271 beta = P_ram / P_mag
272
273 ax2.loglog(x_array / self.R_planet, beta, 'purple',
274 linewidth=2)
275 ax2.axhline(y=1, color='k', linestyle='--', alpha=0.5,
276 label='β = 1')
277 ax2.axvline(x=self.R_mp / self.R_planet, color='g',
278 linestyle='--', alpha=0.7, label='Magnetopause')
279
280 ax2.set_xlabel(r'Distance ($R_{\rm planet}$)', fontsize=11)
281 ax2.set_ylabel(r'Plasma Beta $\beta$', fontsize=11)
282 ax2.set_title(r'Plasma Beta ($\beta = P_{\rm ram}/P_{\rm mag}$)',
283 fontsize=13)
284 ax2.grid(True, alpha=0.3)
285 ax2.legend(fontsize=10)
286
287 plt.tight_layout()
288 return fig
289
290
291def main():
292 """
293 Main function demonstrating 2D magnetosphere model.
294 """
295 print("=" * 60)
296 print("2D Magnetosphere Model")
297 print("=" * 60)
298
299 # Create magnetosphere (Earth parameters)
300 mag = Magnetosphere2D(
301 B_dipole=3.1e-5, # 31,000 nT at equator
302 R_planet=6.4e6, # Earth radius
303 v_sw=4e5, # 400 km/s
304 n_sw=5e6 # 5 cm^-3
305 )
306
307 print(f"\nParameters:")
308 print(f" Dipole field at surface: {mag.B_dipole*1e9:.0f} nT")
309 print(f" Solar wind velocity: {mag.v_sw/1e3:.0f} km/s")
310 print(f" Solar wind density: {mag.n_sw/1e6:.1f} cm^-3")
311 print(f" Ram pressure: {0.5*mag.rho_sw*mag.v_sw**2*1e9:.2f} nPa")
312 print(f"\nMagnetopause standoff distance: {mag.R_mp/mag.R_planet:.1f} R_planet")
313
314 # Plot magnetosphere
315 mag.plot_magnetosphere()
316
317 # Plot pressure profiles
318 mag.plot_pressure_profiles()
319
320 plt.savefig('/tmp/magnetosphere_2d.png', dpi=150, bbox_inches='tight')
321 print("\nPlots saved to /tmp/magnetosphere_2d.png")
322
323 plt.show()
324
325
326if __name__ == "__main__":
327 main()