-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathmain.py
More file actions
655 lines (552 loc) · 27.9 KB
/
main.py
File metadata and controls
655 lines (552 loc) · 27.9 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
import argparse
import json
import matplotlib.pyplot as plt
import matplotlib.animation as animation
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import os
import time
from typing import Tuple, List, Dict, Any, Optional
from matplotlib.figure import Figure
from matplotlib.axes import Axes
from matplotlib.patches import Rectangle
def calculate_derivatives(t: float, x: np.ndarray, G: float, masses: np.ndarray, epsilon: float, dimensions: int = 2) -> np.ndarray:
"""
Calcola la derivata dello stato per un sistema a N corpi in 2D o 3D.
Args:
t: Tempo corrente (secondi)
x: Vettore di stato appiattito contenente posizioni e velocità
G: Costante gravitazionale universale (m³/kg·s²)
masses: Array delle masse dei corpi (kg)
epsilon: Parametro di smoothing per evitare singolarità (m)
dimensions: Numero di dimensioni spaziali (2 o 3)
Returns:
Array appiattito contenente velocità e accelerazioni
"""
num_bodies = len(masses)
# Reshape dinamico: 2D = (N, 4), 3D = (N, 6)
state_size = dimensions * 2
state_matrix = x.reshape(num_bodies, state_size)
positions = state_matrix[:, :dimensions]
velocities = state_matrix[:, dimensions:]
accelerations = np.zeros_like(positions)
for i in range(num_bodies):
for j in range(num_bodies):
if i == j:
continue
r_vec = positions[j] - positions[i]
dist = np.sqrt(np.sum(r_vec**2) + epsilon**2)
accelerations[i] += G * masses[j] * r_vec / dist**3
derivatives_matrix = np.hstack((velocities, accelerations))
return derivatives_matrix.flatten()
def RK4_step(func, t: float, x: np.ndarray, dt: float, *args) -> np.ndarray:
"""
Esegue un singolo passo di integrazione con il metodo Runge-Kutta di 4° ordine.
Args:
func: Funzione che calcola le derivate
t: Tempo corrente (secondi)
x: Vettore di stato corrente
dt: Passo temporale (secondi)
*args: Argomenti aggiuntivi da passare a func
Returns:
Nuovo vettore di stato dopo un passo RK4
"""
k1 = dt * func(t, x, *args)
k2 = dt * func(t + 0.5 * dt, x + 0.5 * k1, *args)
k3 = dt * func(t + 0.5 * dt, x + 0.5 * k2, *args)
k4 = dt * func(t + dt, x + k3, *args)
return x + (k1 + 2*k2 + 2*k3 + k4) / 6.0
def run_simulation(config_filename: str) -> None:
"""
Carica la configurazione ed esegue la simulazione, salvandola come GIF se richiesto.
Args:
config_filename: Percorso del file di configurazione JSON
"""
# --- Caricamento e parsing della configurazione ---
with open(config_filename, 'r') as f:
config = json.load(f)
sim_params = config['simulation']
output_params = config['output']
G = sim_params['G']
epsilon = sim_params['epsilon']
# Quale oggetto seguire? (es. "origin", "com", "Earth")
camera_target_name = sim_params.get('camera_target', 'origin')
# Zoom fisso (usato se il target è 'origin' o un oggetto)
view_radius = sim_params.get('view_radius', 1.5e11)
# Fattore di zoom dinamico (usato solo se il target è 'com')
zoom_factor = sim_params.get('zoom_factor', 1.5)
step_days = sim_params['step_days']
total_days = sim_params['total_days']
dt_seconds = step_days * 24 * 3600
# Non abbiamo più bisogno di total_time_seconds, possiamo calcolare num_steps direttamente
num_steps = int(total_days / step_days)
object_names = list(config['objects'].keys())
objects_data = list(config['objects'].values())
num_bodies = len(objects_data)
masses = np.array([obj['mass'] for obj in objects_data])
initial_positions = np.array([obj['position'] for obj in objects_data])
initial_velocities = np.array([obj['velocity'] for obj in objects_data])
colors = [obj.get('color', 'white') for obj in objects_data]
# Rileva automaticamente se la simulazione è 2D o 3D
dimensions = initial_positions.shape[1]
if dimensions not in [2, 3]:
raise ValueError(f"Le posizioni devono essere 2D o 3D, trovate {dimensions} dimensioni")
is_3d = (dimensions == 3)
print(f"Modalità: {'3D' if is_3d else '2D'}")
# Parametri grafici configurabili
display_params = config.get('display', {})
show_labels = display_params.get('show_labels', True)
show_legend = display_params.get('show_legend', True)
show_com = display_params.get('show_com', False)
use_au = display_params.get('use_au', False)
max_trail_length = display_params.get('max_trail_length', 1000)
proportional_sizes = display_params.get('proportional_sizes', True)
show_info_panel = display_params.get('show_info_panel', True)
show_tooltip = display_params.get('show_tooltip', True)
# Calcola dimensioni proporzionali alle masse (scala logaritmica)
if proportional_sizes and len(masses) > 1:
min_mass = np.min(masses)
max_mass = np.max(masses)
if max_mass / min_mass > 10: # Solo se c'è differenza significativa
marker_sizes = [5 + 20 * (np.log10(m) - np.log10(min_mass)) /
(np.log10(max_mass) - np.log10(min_mass)) for m in masses]
else:
marker_sizes = [10] * num_bodies
else:
marker_sizes = [10] * num_bodies
# Costante per conversione metri -> AU
AU = 1.496e11 # 1 Unità Astronomica in metri
x = np.hstack((initial_positions, initial_velocities)).flatten()
target_index = -1 # Default: 'origin'
if camera_target_name in object_names:
target_index = object_names.index(camera_target_name)
elif camera_target_name == 'com':
target_index = -2 # Codice speciale per Center-of-Mass (COM)
# --- Inizializzazione della figura ---
plt.style.use('dark_background')
fig = plt.figure(figsize=(12, 10))
if is_3d:
ax = fig.add_subplot(111, projection='3d')
ax.grid(True, linestyle=':', alpha=0.3, linewidth=0.5)
# Imposta etichette assi con unità appropriate per 3D
if use_au:
ax.set_xlabel('x (AU)', fontsize=11)
ax.set_ylabel('y (AU)', fontsize=11)
ax.set_zlabel('z (AU)', fontsize=11)
else:
ax.set_xlabel('x (m)', fontsize=11)
ax.set_ylabel('y (m)', fontsize=11)
ax.set_zlabel('z (m)', fontsize=11)
else:
ax = fig.add_subplot(111)
ax.set_aspect('equal')
ax.grid(True, linestyle=':', alpha=0.3, linewidth=0.5)
# Imposta etichette assi con unità appropriate per 2D
if use_au:
ax.set_xlabel('x (AU)', fontsize=11)
ax.set_ylabel('y (AU)', fontsize=11)
else:
ax.set_xlabel('x (m)', fontsize=11)
ax.set_ylabel('y (m)', fontsize=11)
# Crea plot per i corpi con dimensioni proporzionali
if is_3d:
body_plots = [ax.plot([], [], [], 'o', color=c, markersize=marker_sizes[i],
markeredgewidth=0.5, markeredgecolor='white')[0]
for i, c in enumerate(colors)]
trajectory_plots = [ax.plot([], [], [], '-', color=c, linewidth=1, alpha=0.6)[0]
for c in colors]
else:
body_plots = [ax.plot([], [], 'o', color=c, markersize=marker_sizes[i],
markeredgewidth=0.5, markeredgecolor='white')[0]
for i, c in enumerate(colors)]
trajectory_plots = [ax.plot([], [], '-', color=c, linewidth=1, alpha=0.6)[0]
for c in colors]
trajectory_data = [[pos] for pos in initial_positions]
# Etichette per i nomi degli oggetti
if show_labels:
if is_3d:
label_texts = [ax.text(0, 0, 0, name, color='white', fontsize=8,
ha='left', va='bottom', alpha=0.9)
for name in object_names]
else:
label_texts = [ax.text(0, 0, name, color='white', fontsize=8,
ha='left', va='bottom', alpha=0.9)
for name in object_names]
else:
label_texts = []
# Indicatore del centro di massa
if show_com:
if is_3d:
com_plot = ax.plot([], [], [], '+', color='red', markersize=12,
markeredgewidth=2, label='Centro di Massa')[0]
else:
com_plot = ax.plot([], [], '+', color='red', markersize=12,
markeredgewidth=2, label='Centro di Massa')[0]
else:
com_plot = None
# Info panel
if is_3d:
# In 3D, ax.text richiede coordinate 2D (usa text2D)
if show_info_panel:
info_text = ax.text2D(0.02, 0.98, '', transform=ax.transAxes,
color='white', fontsize=9, verticalalignment='top',
family='monospace',
bbox=dict(boxstyle='round', facecolor='black', alpha=0.6))
else:
info_text = ax.text2D(0.02, 0.95, '', transform=ax.transAxes, color='white')
else:
if show_info_panel:
info_text = ax.text(0.02, 0.98, '', transform=ax.transAxes,
color='white', fontsize=9, verticalalignment='top',
family='monospace',
bbox=dict(boxstyle='round', facecolor='black', alpha=0.6))
else:
info_text = ax.text(0.02, 0.95, '', transform=ax.transAxes, color='white')
# Legenda
if show_legend:
legend_elements = [plt.Line2D([0], [0], marker='o', color='w',
markerfacecolor=colors[i], markersize=8,
label=object_names[i], linestyle='None')
for i in range(num_bodies)]
if show_com and com_plot:
legend_elements.append(plt.Line2D([0], [0], marker='+', color='red',
markersize=10, label='Centro di Massa',
linestyle='None', markeredgewidth=2))
ax.legend(handles=legend_elements, loc='upper right', fontsize=8,
framealpha=0.7)
# Tooltip interattivo (disabilitato in 3D per complessità)
if show_tooltip and not is_3d:
tooltip_box = ax.text(0, 0, '', bbox=dict(boxstyle='round,pad=0.5',
facecolor='black', alpha=0.85,
edgecolor='white', linewidth=1.5),
fontsize=9, visible=False, zorder=1000,
family='monospace', color='white')
# Variabili per tracciare lo stato del tooltip
tooltip_state = {'current_positions': initial_positions,
'current_velocities': initial_velocities,
'current_time': 0.0}
def format_number(num: float, use_scientific: bool = False) -> str:
"""Formatta un numero in modo leggibile"""
if use_scientific or abs(num) >= 1e6 or (abs(num) < 0.01 and num != 0):
return f'{num:.2e}'
return f'{num:.2f}'
def on_hover(event) -> None:
"""Gestisce l'evento hover del mouse"""
if event.inaxes != ax:
tooltip_box.set_visible(False)
fig.canvas.draw_idle()
return
# Ottieni posizioni correnti
positions = tooltip_state['current_positions']
velocities = tooltip_state['current_velocities']
t = tooltip_state['current_time']
# Converti coordinate mouse
mouse_x, mouse_y = event.xdata, event.ydata
if mouse_x is None or mouse_y is None:
tooltip_box.set_visible(False)
fig.canvas.draw_idle()
return
# Easter egg: angolo basso a destra
xlim = ax.get_xlim()
ylim = ax.get_ylim()
x_range = xlim[1] - xlim[0]
y_range = ylim[1] - ylim[0]
# Area easter egg: 15% dall'angolo basso a destra, più in alto per non essere tagliato
easter_egg_x_min = xlim[1] - 0.35 * x_range
easter_egg_y_min = ylim[0] + 0.05 * y_range
easter_egg_y_max = ylim[0] + 0.20 * y_range
if mouse_x >= easter_egg_x_min and easter_egg_y_min <= mouse_y <= easter_egg_y_max:
tooltip_box.set_text('Tommaso sei gay')
tooltip_box.set_position((mouse_x, mouse_y))
tooltip_box.set_visible(True)
fig.canvas.draw_idle()
return
# Trova il corpo più vicino al mouse
min_dist = float('inf')
closest_idx = -1
for i in range(num_bodies):
pos_x = positions[i, 0] / AU if use_au else positions[i, 0]
pos_y = positions[i, 1] / AU if use_au else positions[i, 1]
# Calcola distanza in pixel (approssimativa)
dx = pos_x - mouse_x
dy = pos_y - mouse_y
dist = np.sqrt(dx**2 + dy**2)
# Considera anche la dimensione del marker
threshold = marker_sizes[i] * 0.01 * (ax.get_xlim()[1] - ax.get_xlim()[0])
if dist < threshold and dist < min_dist:
min_dist = dist
closest_idx = i
if closest_idx >= 0:
# Calcola statistiche del corpo
vel_magnitude = np.linalg.norm(velocities[closest_idx])
pos_magnitude = np.linalg.norm(positions[closest_idx])
# Calcola energia cinetica del corpo
kinetic_e = 0.5 * masses[closest_idx] * vel_magnitude**2
# Crea testo tooltip
unit_dist = 'AU' if use_au else 'm'
pos_display = positions[closest_idx] / AU if use_au else positions[closest_idx]
tooltip_text = f'{object_names[closest_idx]}\n'
tooltip_text += f'─────────────────\n'
tooltip_text += f'Massa: {format_number(masses[closest_idx], True)} kg\n'
tooltip_text += f'Pos X: {format_number(pos_display[0])} {unit_dist}\n'
tooltip_text += f'Pos Y: {format_number(pos_display[1])} {unit_dist}\n'
if is_3d:
tooltip_text += f'Pos Z: {format_number(pos_display[2])} {unit_dist}\n'
tooltip_text += f'Vel: {format_number(vel_magnitude)} m/s\n'
tooltip_text += f'Vel X: {format_number(velocities[closest_idx, 0])} m/s\n'
tooltip_text += f'Vel Y: {format_number(velocities[closest_idx, 1])} m/s\n'
if is_3d:
tooltip_text += f'Vel Z: {format_number(velocities[closest_idx, 2])} m/s\n'
tooltip_text += f'E_cin: {format_number(kinetic_e, True)} J\n'
tooltip_text += f'Distanza origine: {format_number(pos_magnitude / AU if use_au else pos_magnitude)} {unit_dist}'
tooltip_box.set_text(tooltip_text)
tooltip_box.set_position((mouse_x, mouse_y))
tooltip_box.set_visible(True)
else:
tooltip_box.set_visible(False)
fig.canvas.draw_idle()
# Connetti l'evento hover
fig.canvas.mpl_connect('motion_notify_event', on_hover)
else:
tooltip_box = None
tooltip_state = None
# Funzione di inizializzazione per l'animazione
def init() -> List:
# Imposta la vista iniziale
if is_3d:
center_x, center_y, center_z = 0.0, 0.0, 0.0
if target_index >= 0: # Se il target è un oggetto
center_x = initial_positions[target_index, 0]
center_y = initial_positions[target_index, 1]
center_z = initial_positions[target_index, 2]
# Converti in AU se richiesto
if use_au:
ax.set_xlim((center_x - view_radius) / AU, (center_x + view_radius) / AU)
ax.set_ylim((center_y - view_radius) / AU, (center_y + view_radius) / AU)
ax.set_zlim((center_z - view_radius) / AU, (center_z + view_radius) / AU)
else:
ax.set_xlim(center_x - view_radius, center_x + view_radius)
ax.set_ylim(center_y - view_radius, center_y + view_radius)
ax.set_zlim(center_z - view_radius, center_z + view_radius)
for j in range(num_bodies):
body_plots[j].set_data([], [])
body_plots[j].set_3d_properties([])
trajectory_plots[j].set_data([], [])
trajectory_plots[j].set_3d_properties([])
if show_com and com_plot:
com_plot.set_data([], [])
com_plot.set_3d_properties([])
else:
center_x, center_y = 0.0, 0.0
if target_index >= 0: # Se il target è un oggetto
center_x = initial_positions[target_index, 0]
center_y = initial_positions[target_index, 1]
# Converti in AU se richiesto
if use_au:
ax.set_xlim((center_x - view_radius) / AU, (center_x + view_radius) / AU)
ax.set_ylim((center_y - view_radius) / AU, (center_y + view_radius) / AU)
else:
ax.set_xlim(center_x - view_radius, center_x + view_radius)
ax.set_ylim(center_y - view_radius, center_y + view_radius)
for j in range(num_bodies):
body_plots[j].set_data([], [])
trajectory_plots[j].set_data([], [])
if show_com and com_plot:
com_plot.set_data([], [])
if show_labels:
for label in label_texts:
label.set_text('')
if show_tooltip and tooltip_box:
tooltip_box.set_visible(False)
info_text.set_text('')
artists = body_plots + trajectory_plots + [info_text]
if show_labels:
artists += label_texts
if show_com and com_plot:
artists.append(com_plot)
if show_tooltip and tooltip_box:
artists.append(tooltip_box)
return artists
# Funzione di aggiornamento per ogni frame dell'animazione
def update(frame: int) -> List:
nonlocal x
t = frame * dt_seconds
x = RK4_step(calculate_derivatives, t, x, dt_seconds, G, masses, epsilon, dimensions)
# Estrai posizioni e velocità in base alla dimensionalità
if is_3d:
current_positions = x.reshape(num_bodies, 6)[:, :3]
current_velocities = x.reshape(num_bodies, 6)[:, 3:]
else:
current_positions = x.reshape(num_bodies, 4)[:, :2]
current_velocities = x.reshape(num_bodies, 4)[:, 2:]
# --- Logica di centratura telecamera ---
if is_3d:
center_point = np.array([0.0, 0.0, 0.0]) # Default per 'origin'
else:
center_point = np.array([0.0, 0.0]) # Default per 'origin'
radius = view_radius # Default zoom
com_position = None
if target_index >= 0:
# Target è un oggetto (es. "Earth")
center_point = current_positions[target_index]
# radius rimane view_radius
elif target_index == -2:
# Target è il Centro di Massa (COM) con zoom dinamico
total_mass = np.sum(masses)
center_point = np.sum(masses[:, np.newaxis] * current_positions, axis=0) / total_mass
com_position = center_point
max_dist_from_com = np.max(np.sqrt(np.sum((current_positions - center_point)**2, axis=1)))
radius = max_dist_from_com * zoom_factor
if radius == 0: radius = 1e11 # Failsafe
# Calcola COM anche se non è il target (per show_com)
if show_com and com_position is None:
total_mass = np.sum(masses)
com_position = np.sum(masses[:, np.newaxis] * current_positions, axis=0) / total_mass
# Imposta limiti assi (con conversione AU se richiesto)
if is_3d:
if use_au:
ax.set_xlim((center_point[0] - radius) / AU, (center_point[0] + radius) / AU)
ax.set_ylim((center_point[1] - radius) / AU, (center_point[1] + radius) / AU)
ax.set_zlim((center_point[2] - radius) / AU, (center_point[2] + radius) / AU)
else:
ax.set_xlim(center_point[0] - radius, center_point[0] + radius)
ax.set_ylim(center_point[1] - radius, center_point[1] + radius)
ax.set_zlim(center_point[2] - radius, center_point[2] + radius)
else:
if use_au:
ax.set_xlim((center_point[0] - radius) / AU, (center_point[0] + radius) / AU)
ax.set_ylim((center_point[1] - radius) / AU, (center_point[1] + radius) / AU)
else:
ax.set_xlim(center_point[0] - radius, center_point[0] + radius)
ax.set_ylim(center_point[1] - radius, center_point[1] + radius)
# Aggiorna posizioni corpi e traiettorie
for j in range(num_bodies):
# Converti posizioni in AU se richiesto
pos_x = current_positions[j, 0] / AU if use_au else current_positions[j, 0]
pos_y = current_positions[j, 1] / AU if use_au else current_positions[j, 1]
if is_3d:
pos_z = current_positions[j, 2] / AU if use_au else current_positions[j, 2]
body_plots[j].set_data([pos_x], [pos_y])
body_plots[j].set_3d_properties([pos_z])
# Gestione traiettorie con lunghezza massima
trajectory_data[j].append(current_positions[j])
if len(trajectory_data[j]) > max_trail_length:
trajectory_data[j] = trajectory_data[j][-max_trail_length:]
# Se la telecamera si muove, la scia deve essere ricalcolata relativa alla telecamera
if target_index >= 0:
# Converte i dati storici della traiettoria nel frame di riferimento attuale
relative_trajectory = np.array(trajectory_data[j]) - center_point
if use_au:
relative_trajectory = relative_trajectory / AU
trajectory_plots[j].set_data(relative_trajectory[:, 0], relative_trajectory[:, 1])
trajectory_plots[j].set_3d_properties(relative_trajectory[:, 2])
else:
# Comportamento normale (frame fisso o COM)
traj_array = np.array(trajectory_data[j])
if use_au:
traj_array = traj_array / AU
trajectory_plots[j].set_data(traj_array[:, 0], traj_array[:, 1])
trajectory_plots[j].set_3d_properties(traj_array[:, 2])
else:
body_plots[j].set_data([pos_x], [pos_y])
# Gestione traiettorie con lunghezza massima
trajectory_data[j].append(current_positions[j])
if len(trajectory_data[j]) > max_trail_length:
trajectory_data[j] = trajectory_data[j][-max_trail_length:]
# Se la telecamera si muove, la scia deve essere ricalcolata relativa alla telecamera
if target_index >= 0:
# Converte i dati storici della traiettoria nel frame di riferimento attuale
relative_trajectory = np.array(trajectory_data[j]) - center_point
if use_au:
relative_trajectory = relative_trajectory / AU
trajectory_plots[j].set_data(relative_trajectory.T)
else:
# Comportamento normale (frame fisso o COM)
traj_array = np.array(trajectory_data[j])
if use_au:
traj_array = traj_array / AU
trajectory_plots[j].set_data(traj_array.T)
# Aggiorna etichette
if show_labels:
for j in range(num_bodies):
pos_x = current_positions[j, 0] / AU if use_au else current_positions[j, 0]
pos_y = current_positions[j, 1] / AU if use_au else current_positions[j, 1]
if is_3d:
pos_z = current_positions[j, 2] / AU if use_au else current_positions[j, 2]
label_texts[j].set_position((pos_x, pos_y))
label_texts[j].set_3d_properties(pos_z, 'z')
else:
label_texts[j].set_position((pos_x, pos_y))
label_texts[j].set_text(object_names[j])
# Aggiorna indicatore centro di massa
if show_com and com_plot and com_position is not None:
com_x = com_position[0] / AU if use_au else com_position[0]
com_y = com_position[1] / AU if use_au else com_position[1]
if is_3d:
com_z = com_position[2] / AU if use_au else com_position[2]
com_plot.set_data([com_x], [com_y])
com_plot.set_3d_properties([com_z])
else:
com_plot.set_data([com_x], [com_y])
# Calcola statistiche per info panel
years_passed = t / (365.25 * 24 * 3600)
# current_velocities già estratto sopra
# Calcola energia totale (cinetica + potenziale)
kinetic_energy = 0.5 * np.sum(masses[:, np.newaxis] * current_velocities**2)
potential_energy = 0.0
for i in range(num_bodies):
for j in range(i + 1, num_bodies):
r = np.linalg.norm(current_positions[j] - current_positions[i])
potential_energy -= G * masses[i] * masses[j] / (r + epsilon)
total_energy = kinetic_energy + potential_energy
# Info panel migliorato
if show_info_panel:
info_lines = [
f'Tempo: {years_passed:.2f} anni',
f'Frame: {frame}/{num_steps}',
f'Energia: {total_energy:.2e} J',
f'E_cin: {kinetic_energy:.2e} J',
f'E_pot: {potential_energy:.2e} J'
]
info_text.set_text('\n'.join(info_lines))
else:
info_text.set_text(f'Anno: {years_passed:.2f}')
# Aggiorna stato tooltip
if show_tooltip and tooltip_state is not None:
tooltip_state['current_positions'] = current_positions
tooltip_state['current_velocities'] = current_velocities
tooltip_state['current_time'] = t
if frame % 10 == 0:
print(f"Generando frame: {frame}/{num_steps} ({frame/num_steps*100:.1f}%)", end='\r')
artists = body_plots + trajectory_plots + [info_text]
if show_labels:
artists += label_texts
if show_com and com_plot:
artists.append(com_plot)
if show_tooltip and tooltip_box:
artists.append(tooltip_box)
return artists
# --- Creazione e Salvataggio dell'Animazione ---
ani = animation.FuncAnimation(fig, update, frames=num_steps,
init_func=init, blit=True, interval=20)
# 1. Salva il GIF se richiesto
if output_params.get("save_gif", False):
filename = output_params.get("filename", "simulation.gif")
fps = output_params.get("fps", 30)
print(f"\nSalvataggio dell'animazione come '{filename}' a {fps} FPS. Potrebbe richiedere tempo...")
start_time = time.time()
ani.save(filename, writer='pillow', fps=fps)
end_time = time.time()
print(f"\nSalvataggio completato in {end_time - start_time:.2f} secondi.")
# 2. Mostra l'animazione
print("Mostrando l'animazione in tempo reale...")
plt.show()
if __name__ == '__main__':
parser = argparse.ArgumentParser(description='Avvia una simulazione N-corpi.')
default_config_path = os.path.join('configs', 'config.json')
parser.add_argument('config_file',
nargs='?',
default='configs/config.json',
help='Percorso del file di configurazione JSON (default: config.json)')
args = parser.parse_args()
run_simulation(args.config_file)