1#!/usr/bin/env python3
2"""
3Magnetic Field Polarity Reversal
4
5This module demonstrates magnetic field polarity reversals similar to those
6observed in Earth's paleomagnetic record and the Sun's 22-year magnetic cycle.
7
8The model uses an α-Ω dynamo with stochastic noise to trigger reversals.
9Key features:
10- Bistable dipole states (normal/reversed polarity)
11- Stochastic transitions between states
12- Waiting time distribution analysis
13- Comparison with paleomagnetic observations
14
15Author: MHD Course Examples
16License: MIT
17"""
18
19import numpy as np
20import matplotlib.pyplot as plt
21from scipy.integrate import odeint
22
23
24class ReversalDynamo:
25 """
26 Stochastic dynamo model with field reversals.
27
28 This implements a low-order model of the geodynamo/solar dynamo
29 that exhibits polarity reversals due to stochastic fluctuations.
30
31 Attributes:
32 alpha (float): Alpha effect strength
33 omega (float): Omega effect (differential rotation)
34 eta (float): Magnetic diffusivity
35 noise_strength (float): Amplitude of stochastic forcing
36 """
37
38 def __init__(self, alpha=1.0, omega=10.0, eta=0.1, noise_strength=0.5):
39 """
40 Initialize reversal dynamo.
41
42 Parameters:
43 alpha (float): Alpha effect strength
44 omega (float): Omega effect
45 eta (float): Diffusivity
46 noise_strength (float): Noise amplitude
47 """
48 self.alpha = alpha
49 self.omega = omega
50 self.eta = eta
51 self.noise_strength = noise_strength
52
53 # Dynamo number
54 self.D = alpha * omega / eta**2
55
56 # State variables (dipole + quadrupole components)
57 self.dipole = 0.1
58 self.quadrupole = 0.1
59
60 # History
61 self.time_history = []
62 self.dipole_history = []
63 self.quadrupole_history = []
64 self.reversal_times = []
65
66 def rhs(self, state, t, noise):
67 """
68 Right-hand side of dynamo equations.
69
70 Simplified model based on Rikitake dynamo / Lorenz-like system.
71
72 Parameters:
73 state (list): [dipole, quadrupole]
74 t (float): Time
75 noise (float): Noise realization
76
77 Returns:
78 list: Time derivatives
79 """
80 d, q = state
81
82 # Coupled equations with cubic nonlinearity (saturation)
83 # and stochastic forcing
84 dd_dt = self.alpha * q - self.eta * d - d * (d**2 + q**2) + noise
85 dq_dt = self.omega * d - self.eta * q - q * (d**2 + q**2)
86
87 return [dd_dt, dq_dt]
88
89 def detect_reversal(self, dipole_current, dipole_previous):
90 """
91 Detect if a reversal has occurred.
92
93 Parameters:
94 dipole_current (float): Current dipole value
95 dipole_previous (float): Previous dipole value
96
97 Returns:
98 bool: True if reversal detected
99 """
100 # Reversal = sign change of dipole
101 return dipole_current * dipole_previous < 0
102
103 def run_simulation(self, t_end=1000.0, dt=0.01):
104 """
105 Run stochastic dynamo simulation.
106
107 Parameters:
108 t_end (float): End time
109 dt (float): Time step
110 """
111 n_steps = int(t_end / dt)
112 t_array = np.linspace(0, t_end, n_steps)
113
114 print("Running magnetic reversal simulation...")
115 print(f"Dynamo number D = {self.D:.2f}")
116 print(f"Noise strength = {self.noise_strength:.2f}")
117
118 # Initialize
119 state = np.array([self.dipole, self.quadrupole])
120 self.time_history = []
121 self.dipole_history = []
122 self.quadrupole_history = []
123 self.reversal_times = []
124
125 previous_dipole = state[0]
126
127 for i, t in enumerate(t_array):
128 # Generate noise
129 noise = self.noise_strength * np.random.randn()
130
131 # Integrate one step
132 state_new = odeint(self.rhs, state, [t, t + dt], args=(noise,))[-1]
133 state = state_new
134
135 # Store
136 self.time_history.append(t)
137 self.dipole_history.append(state[0])
138 self.quadrupole_history.append(state[1])
139
140 # Detect reversal
141 if self.detect_reversal(state[0], previous_dipole):
142 self.reversal_times.append(t)
143 print(f"Reversal at t = {t:.2f}")
144
145 previous_dipole = state[0]
146
147 if i % 10000 == 0:
148 print(f"Step {i}/{n_steps}, t={t:.1f}")
149
150 print(f"\nSimulation complete!")
151 print(f"Total reversals: {len(self.reversal_times)}")
152
153 if len(self.reversal_times) > 1:
154 waiting_times = np.diff(self.reversal_times)
155 print(f"Mean waiting time: {np.mean(waiting_times):.1f}")
156 print(f"Std waiting time: {np.std(waiting_times):.1f}")
157
158 def plot_field_evolution(self):
159 """
160 Plot magnetic field evolution and reversals.
161 """
162 fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(12, 10))
163
164 # Dipole evolution
165 ax1.plot(self.time_history, self.dipole_history, 'b-', linewidth=1.5,
166 label='Dipole component')
167 ax1.axhline(y=0, color='k', linestyle='--', alpha=0.5)
168
169 # Mark reversals
170 for t_rev in self.reversal_times:
171 ax1.axvline(x=t_rev, color='r', alpha=0.3, linewidth=0.5)
172
173 ax1.set_ylabel('Dipole Field', fontsize=11)
174 ax1.set_title('Magnetic Field Reversals (Dipole Component)', fontsize=13)
175 ax1.grid(True, alpha=0.3)
176 ax1.legend(fontsize=10)
177
178 # Quadrupole evolution
179 ax2.plot(self.time_history, self.quadrupole_history, 'g-', linewidth=1.5,
180 label='Quadrupole component')
181 ax2.axhline(y=0, color='k', linestyle='--', alpha=0.5)
182 ax2.set_ylabel('Quadrupole Field', fontsize=11)
183 ax2.grid(True, alpha=0.3)
184 ax2.legend(fontsize=10)
185
186 # Phase space
187 ax3.plot(self.dipole_history, self.quadrupole_history, 'b-',
188 alpha=0.5, linewidth=0.5)
189 ax3.plot(self.dipole_history[0], self.quadrupole_history[0], 'go',
190 markersize=8, label='Start')
191 ax3.plot(self.dipole_history[-1], self.quadrupole_history[-1], 'ro',
192 markersize=8, label='End')
193 ax3.set_xlabel('Dipole', fontsize=11)
194 ax3.set_ylabel('Quadrupole', fontsize=11)
195 ax3.set_title('Phase Space Trajectory', fontsize=13)
196 ax3.grid(True, alpha=0.3)
197 ax3.legend(fontsize=10)
198 ax3.axhline(y=0, color='k', linestyle='--', alpha=0.3)
199 ax3.axvline(x=0, color='k', linestyle='--', alpha=0.3)
200
201 plt.tight_layout()
202 return fig
203
204 def plot_reversal_statistics(self):
205 """
206 Plot reversal waiting time statistics.
207 """
208 if len(self.reversal_times) < 2:
209 print("Not enough reversals for statistics")
210 return None
211
212 waiting_times = np.diff(self.reversal_times)
213
214 fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(12, 5))
215
216 # Histogram of waiting times
217 ax1.hist(waiting_times, bins=20, density=True, alpha=0.7,
218 color='blue', edgecolor='black')
219
220 # Theoretical exponential (Poisson process)
221 mean_wait = np.mean(waiting_times)
222 x_theory = np.linspace(0, np.max(waiting_times), 100)
223 y_theory = (1 / mean_wait) * np.exp(-x_theory / mean_wait)
224 ax1.plot(x_theory, y_theory, 'r-', linewidth=2,
225 label=f'Exponential (λ={1/mean_wait:.3f})')
226
227 ax1.set_xlabel('Waiting Time', fontsize=11)
228 ax1.set_ylabel('Probability Density', fontsize=11)
229 ax1.set_title('Distribution of Reversal Waiting Times', fontsize=13)
230 ax1.grid(True, alpha=0.3)
231 ax1.legend(fontsize=10)
232
233 # Cumulative distribution
234 sorted_waits = np.sort(waiting_times)
235 cumulative = np.arange(1, len(sorted_waits) + 1) / len(sorted_waits)
236
237 ax2.plot(sorted_waits, cumulative, 'b-', linewidth=2,
238 label='Empirical CDF')
239
240 # Theoretical CDF
241 cdf_theory = 1 - np.exp(-x_theory / mean_wait)
242 ax2.plot(x_theory, cdf_theory, 'r--', linewidth=2,
243 label='Exponential CDF')
244
245 ax2.set_xlabel('Waiting Time', fontsize=11)
246 ax2.set_ylabel('Cumulative Probability', fontsize=11)
247 ax2.set_title('Cumulative Distribution Function', fontsize=13)
248 ax2.grid(True, alpha=0.3)
249 ax2.legend(fontsize=10)
250
251 plt.tight_layout()
252 return fig
253
254 def plot_paleomagnetic_comparison(self):
255 """
256 Plot comparison with schematic paleomagnetic data.
257 """
258 fig, (ax1, ax2) = plt.subplots(2, 1, figsize=(12, 8))
259
260 # Our simulation (zoom to show detail)
261 t_max = min(200, self.time_history[-1])
262 idx = np.searchsorted(self.time_history, t_max)
263
264 dipole_sign = np.sign(self.dipole_history[:idx])
265 ax1.fill_between(self.time_history[:idx], -1, 1,
266 where=(dipole_sign > 0),
267 alpha=0.7, color='black', label='Normal polarity')
268 ax1.fill_between(self.time_history[:idx], -1, 1,
269 where=(dipole_sign < 0),
270 alpha=0.7, color='white', edgecolor='black',
271 label='Reversed polarity')
272 ax1.set_ylabel('Polarity', fontsize=11)
273 ax1.set_title('Simulated Magnetic Polarity History', fontsize=13)
274 ax1.set_ylim([-1.5, 1.5])
275 ax1.set_yticks([-1, 0, 1])
276 ax1.set_yticklabels(['Reversed', 'Transition', 'Normal'])
277 ax1.legend(loc='upper right', fontsize=10)
278 ax1.grid(True, alpha=0.3, axis='x')
279
280 # Schematic Earth paleomagnetic record (simplified)
281 # Based on geomagnetic polarity time scale
282 earth_times = np.array([0, 30, 60, 80, 100, 125, 150, 170, 200])
283 earth_polarity = np.array([1, 1, -1, -1, 1, -1, 1, 1, -1])
284
285 for i in range(len(earth_times) - 1):
286 color = 'black' if earth_polarity[i] > 0 else 'white'
287 ax2.fill_between([earth_times[i], earth_times[i+1]], -1, 1,
288 alpha=0.7, color=color, edgecolor='black')
289
290 ax2.set_xlabel('Time (relative units)', fontsize=11)
291 ax2.set_ylabel('Polarity', fontsize=11)
292 ax2.set_title('Schematic Earth Paleomagnetic Record', fontsize=13)
293 ax2.set_ylim([-1.5, 1.5])
294 ax2.set_yticks([-1, 0, 1])
295 ax2.set_yticklabels(['Reversed', 'Transition', 'Normal'])
296 ax2.grid(True, alpha=0.3, axis='x')
297
298 plt.tight_layout()
299 return fig
300
301
302def main():
303 """
304 Main function demonstrating magnetic field reversals.
305 """
306 print("=" * 60)
307 print("Magnetic Field Polarity Reversal Simulation")
308 print("=" * 60)
309
310 # Create dynamo
311 dynamo = ReversalDynamo(
312 alpha=1.0,
313 omega=10.0,
314 eta=0.1,
315 noise_strength=0.8 # Higher noise → more reversals
316 )
317
318 # Run simulation
319 dynamo.run_simulation(t_end=1000.0, dt=0.01)
320
321 # Plot results
322 dynamo.plot_field_evolution()
323 if len(dynamo.reversal_times) >= 2:
324 dynamo.plot_reversal_statistics()
325 dynamo.plot_paleomagnetic_comparison()
326
327 plt.savefig('/tmp/field_reversal.png', dpi=150, bbox_inches='tight')
328 print("\nPlots saved to /tmp/field_reversal.png")
329
330 plt.show()
331
332
333if __name__ == "__main__":
334 main()