import matplotlib.pyplot as plt
import seaborn as sns

def plot_inertia_density(points, tuning_c=3.0):
    # Ponownie obliczamy odchylenia (Z-scores) jak w g07dbc
    theta = np.median(points, axis=0)
    diff = np.abs(points - theta)
    mad = np.median(diff, axis=0) / 0.67448
    mad[mad == 0] = 1e-6
    modified_z_score = 0.6745 * diff / mad
    deviation_magnitude = np.linalg.norm(modified_z_score, axis=1)

    # Tworzenie wykresu
    plt.figure(figsize=(12, 6))
    
    # Histogram gęstości
    sns.histplot(deviation_magnitude, kde=True, color="blue", bins=50, label="Gęstość inercji")
    
    # Linia progu relatywizmu (tuning_c)
    plt.axvline(tuning_c, color='red', linestyle='--', linewidth=2, label=f'Próg relatywizmu (c={tuning_c})')
    
    # Opisy zgodne z Twoją fizyką
    plt.fill_betweenx([0, plt.gca().get_ylim()[1]], 0, tuning_c, color='blue', alpha=0.1, label='Strefa Klasycznej Inercji')
    plt.fill_betweenx([0, plt.gca().get_ylim()[1]], tuning_c, max(deviation_magnitude), color='red', alpha=0.1, label='Strefa Relatywistyczna (Singularity)')
    
    plt.title('Rozkład Gęstości Inercji w Układzie', fontsize=15)
    plt.xlabel('Magnituda odchylenia (Z-score / Bezwładność)', fontsize=12)
    plt.ylabel('Liczba cząstek danych', fontsize=12)
    plt.legend()
    plt.grid(alpha=0.3)
    plt.show()

# Uruchomienie analizy dla punktów ze sceny Nerfies
plot_inertia_density(new_scene_manager.points, tuning_c=3.0)