import netCDF4 as nc
import matplotlib.pyplot as plt
import numpy as np

file_path = '04008.atmos_daily.nc'
TAG = "sim00"
output_dir = "/cantiere/opendap_data/emm_deliverable/"

try:
    dataset = nc.Dataset(file_path, mode='r')
    time_raw = dataset.variables['time'][:]
    areo_raw = dataset.variables['areo'][:, 0] 
    dataset.close()

    # 1. IDENTIFICAZIONE PUNTI BUONI
    valore_finale = areo_raw[-1]
    valore_iniziale = areo_raw[0]
    mask_valore = (areo_raw <= valore_finale) & (areo_raw >= valore_iniziale)
    
    time_clean = time_raw[mask_valore]
    areo_clean = areo_raw[mask_valore]

    # 2. CALCOLO DERIVATA (Velocità Orbitale)
    delta_areo = np.diff(areo_clean)
    delta_time = np.diff(time_clean)
    pendenza = delta_areo / delta_time
    
    # Filtro pendenza per eliminare residui di spike
    mask_pendenza = np.abs(pendenza) < 1.0 
    
    t_plot = time_clean[:-1][mask_pendenza]
    p_plot = pendenza[mask_pendenza]
    # Prendiamo gli Ls corrispondenti ai punti della pendenza
    ls_plot = areo_clean[:-1][mask_pendenza]

    # --- PLOTTING (3 Figure sovrapposte) ---
    fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(12, 15), sharex=True)

    # Subplot 1: Areo cumulativo vs Time (Allineato al tempo!)
    ax1.plot(time_clean, areo_clean, color='blue', linewidth=1)
    ax1.set_title(f'1. Ls Cumulativo vs Tempo - {TAG}')
    ax1.set_ylabel('Ls Totale (degrees)')
    ax1.grid(True)

    # Subplot 2: Velocità Orbitale dLs/dt vs Time
    ax2.plot(t_plot, p_plot, color='red', marker='o', markersize=2, linestyle='-')
    ax2.set_title('2. Velocità Orbitale dLs/dt (PULITA)')
    ax2.set_ylabel('Degrees / Day')
    ax2.grid(True)

    # Subplot 3: Ls Stagionale (Modulo 360) vs Time
    # Questo ti serve per leggere l'Ls del perielio (0-360)
    ax3.plot(time_clean, areo_clean % 360, color='green', linewidth=1)
    ax3.set_title('3. Ls Stagionale (Modulo 360) vs Tempo')
    ax3.set_xlabel('Time (days)')
    ax3.set_ylabel('Ls (0-360°)')
    ax3.grid(True)

    plt.tight_layout()
    
    fn = f"check_time_areo_3PANE_{TAG}.png"
    plt.savefig(output_dir + fn)
    
    print(f"Dataset ripulito. Generata figura a 3 pannelli.")
    print(f"Lancio calcolo perielio teorico...")
    
    # Identifica il punto di max velocità nel grafico
    idx_max = np.argmax(p_plot)
    time_peri = t_plot[idx_max]
    ls_peri = ls_plot[idx_max] % 360
    
    print(f"\n--- RISULTATO ANALISI ---")
    print(f"Picco velocità trovato al giorno: {time_peri:.2f}")
    print(f"Corrispondente a Ls: {ls_peri:.2f}°")

except Exception as e:
    print(f"Errore: {e}")

print(f"\nLink: https://opendap.le.isac.cnr.it/opendap/emm_deliverable/{fn}")