Velocità della luce

2026-11-09 18:20:54 +01:00
parent 8425ff90da
commit 4d25330361
9 changed files with 423 additions and 76 deletions
@@ -0,0 +1,71 @@
+f_finale;f_iniziale;x_finale;x_iniziale
+-907;-50;5,78;5,55
+-1444;-29;5,78;5,23
+-920;-76;5,62;5,37
+-900;-62;5,62;5,38
+-910;-110;5,61;5,35
+-923;-84;5,57;5,38
+-683;-56;5,58;5,4
+-850;-107;5,57;5,35
+-791;-129;5,56;5,37
+-834;-75;5,57;5,32
+-826;-112;5,52;5,32
+-811;-130;5,52;5,32
+-915;-70;5,54;5,29
+-894;-31;5,54;5,3
+-850;-100;5,52;5,31
+-797;-126;5,52;5,32
+-908;-16;5,55;5,29
+-921;-50;5,53;5,27
+-831;-42;5,52;5,31
+-900;-136;5,5;5,28
+-914;-58;5,52;5,26
+-1495;-75;5,1;4,69
+893;39;6,1;5,85
+862;55;6,09;5,85
+922;137;6,1;5,88
+829;25;5,98;5,75
+911;60;6,1;5,87
+849;81;6,08;5,86
+910;65;6,1;5,86
+713;30;6,04;5,84
+887;43;6,09;5,85
+901;64;6,1;5,86
+836;72;6,07;5,86
+890;29;6,09;5,85
+918;36;6,09;5,86
+927;77;6,11;5,87
+822;87;6,06;5,86
+789;136;6,05;5,87
+884;128;6,07;5,87
+722;100;6,03;5,85
+926;112;6,08;5,85
+842;135;6,09;5,88
+-748;767;5,63;6,08
+-631;532;5,64;5,98
+-920;900;5,55;6,05
+-910;915;5,55;6,07
+-757;858;5,6;6,06
+-882;866;5,55;6,05
+-758;920;5,59;6,08
+-790;693;5,6;6,01
+-840;929;5,59;6,08
+-905;929;5,55;6,05
+-679;897;5,59;6,04
+-714;808;5,59;6,02
+-673;796;5,57;6
+-612;935;5,59;6,04
+-817;927;5,52;6,02
+-607;805;5,56;5,97
+-769;833;5,51;5,97
+-880;937;5,48;5,99
+-920;900;5,5;6,01
+-798;801;5,54;5,99
+-789;731;5,52;5,94
+-765;935;5,52;5,99
+-935;787;5,47;5,96
+-935;781;5,45;5,94
+-739;938;5,51;5,99
+-936;948;5,48;6,02
+-754;930;5,53;6,01
+-492;896;5,62;6,01
@@ -0,0 +1,2 @@
+D1;D2;D3;D4;D5;D6;D7
+6,2;6,19;6,19;6,19;6,2;6,19;6,2
@@ -0,0 +1,2 @@
+D1;D2;D3;D4;D5;D6;D7
+6,66;6,66;6,66;6,65;6,67;6,67;6,67
@@ -0,0 +1,2 @@
+D1;D2;D3;D4;D5;D6;D7
+0,67;0,65;0,65;0,66;0,66;0,65;0,66
@@ -20,7 +20,6 @@ def _():
    Data = defs["Data"]
    compat = defs["compat"]
    media_pesata = defs["media_pesata"]
-    media_pesata_dersimonian_laird = defs["media_pesata_dersimonian_laird"]
    rimuovi_outlier = defs["rimuovi_outlier"]
    plot_gauss = defs["plot_gauss"]
    cartella_csv = Path(__file__).resolve().parent / "csv" / "focaleLente"
@@ -8,6 +8,7 @@ dependencies = [
    "altair>=6.3.0",
    "marimo>=0.25.0",
    "matplotlib>=3.11.2",
+    "mcp>=2.3.0",
    "numpy>=2.5.3",
    "pandas>=3.0.6",
    "scipy>=1.18.1",
@@ -22,3 +23,12 @@ lofm = "lofm:main"
 [build-system]
 requires = ["uv_build>=0.12.19,<0.13.0"]
 build-backend = "uv_build"
+
+[dependency-groups]
+dev = [
+    "python-lsp-ruff>=2.3.4",
+    "python-lsp-server>=1.15.0",
+    "vegafusion>=2.0.3",
+    "vl-convert-python>=1.9.0.post1",
+    "websockets>=17.3",
+]
@@ -185,9 +185,15 @@ def _(np, pd, stats):
                valori_puliti, rimossi = rimuovi_outlier(valori, err_strum)

                n = len(valori_puliti)
-                s = valori_puliti.std(ddof=1)
-                if err_della_media: s /= np.sqrt(n)
-                u = np.sqrt(s**2 + err_strum**2)
+                if n > 2:
+                    s = valori_puliti.std(ddof=1)
+                    if err_della_media:
+                        s /= np.sqrt(n)
+                    u = np.sqrt(s**2 + err_strum**2)
+                else:
+                    s = 0.0
+                    u = float(err_strum)
+            

                self.analisi_stat.loc[i, prefix]             = valori_puliti.mean()
                self.analisi_stat.loc[i, f"u{prefix}_strum"] = err_strum
@@ -199,7 +205,7 @@ def _(np, pd, stats):
        return self

    Data.stat_chauv = stat_chauv
-    return
+    return (trova_outlier,)


@app.cell(hide_code=True)
@@ -207,90 +213,83 @@ def _(mo):
    mo.md(r"""
    # Media pesata

-    input:
-        - array x: valori
-        - array ux: uvalori
-    output:
-        - xm: x medio
-        - uxm: sigma x medio
-        - chi2: $\chi^2$
-        - dof: gradi di libertà
-        - Pchi2: $P(\chi^2)$
+    **media_pesata( x, ux, chauvenet=True, info_outlier=False )**
+
+    *input:*
+      - array x: valori
+      - array ux: incertezze (tutte > 0)
+      - chauvenet: se True applica il criterio di Chauvenet gia' definito
+        (trova_outlier / p_t_student), con incertezza di_resolutione data
+        dall'RMS delle ux del campione
+      - info_outlier: se True aggiunge la lista dei punti scartati
+
+    *output:*
+      - xm: x medio
+      - uxm: sigma x medio
+      - chi2: chi quadro ridotto $\tilde\chi^2$
+      - dof: gradi di liberta
+      - Pchi2: $P(\chi^2)$
+      - out: [(valore, incertezza)] rimossi, solo se info_outlier=True

    data:
-        - prefix: prefisso del nome della variabile
+      - prefix: prefisso del nome della variabile
    """)
    return


@app.cell
-def _(np, stats):
-    def media_pesata(x, ux):
-        x, ux = np.asarray(x, float), np.asarray(ux, float)
-        w   = 1 / ux**2
-        xm  = np.sum(w * x) / np.sum(w)
-        uxm = 1 / np.sqrt(np.sum(w))
-        chi2 = np.sum(w * (x - xm) ** 2)
-        dof  = len(x) - 1
-        return xm, uxm, chi2/dof, dof, stats.chi2.sf(chi2, dof)
+def _(np, stats, trova_outlier):
+    def media_pesata(x, ux, chauvenet=True, info_outlier=False):
+        x, ux = np.asarray(x, float).ravel(), np.asarray(ux, float).ravel()
+
+        if x.shape != ux.shape:
+            raise ValueError("media_pesata: x e ux devono avere la stessa lunghezza")
+
+        validi = np.isfinite(x) & np.isfinite(ux)
+        x, ux  = x[validi], ux[validi]
+
+        if len(x) and np.any(ux <= 0):
+            raise ValueError("media_pesata: tutte le incertezze devono essere > 0")
+
+        out = []
+        if chauvenet:
+            # Chauvenet gia' presente nel notebook: trova_outlier -> p_t_student
+            # qui l'incertezza di risoluzione e' l'RMS delle ux del campione corrente
+            while len(x) > 2:
+                idx = trova_outlier(x, np.sqrt(np.mean(ux**2)))
+                if idx is None:
+                    break
+                out.append((x[idx], ux[idx]))
+                x  = np.delete(x, idx)
+                ux = np.delete(ux, idx)
+
+        if len(x) == 0:
+            return np.nan, np.nan, np.nan, 0, np.nan, []
+
+        w        = 1 / ux**2
+        xm       = np.sum(w * x) / np.sum(w)
+        uxm      = 1 / np.sqrt(np.sum(w))
+        chi2     = np.sum(w * (x - xm) ** 2)
+        dof      = len(x) - 1
+        chi2_red = chi2 / dof if dof > 0 else np.nan
+        P        = stats.chi2.sf(chi2, dof) if dof > 0 else np.nan
+
+        if info_outlier:
+            return xm, uxm, chi2_red, dof, P, out
+        return xm, uxm, chi2_red, dof, P
+
+
+    def media_pesata_data(self, prefix, chauvenet=True, info_outlier=False):
+        return media_pesata(self.analisi_stat[prefix],
+                            self.analisi_stat[f"u{prefix}"],
+                            chauvenet=chauvenet,
+                            info_outlier=info_outlier)

-    def media_pesata_data(self, prefix):
-        return media_pesata(self.analisi_stat[prefix], self.analisi_stat[f"u{prefix}"])

    Data.media_pesata = media_pesata_data
    return


-@app.cell
-def _(np, stats):
-    def media_pesata_dersimonian_laird(x, ux):
-        x, ux = np.asarray(x, float), np.asarray(ux, float)
-        n = len(x)
-        dof = n - 1
-
-        if n <= 1:
-            raise ValueError("Il dataset deve contenere almeno due misure.")
-        if np.any(ux <= 0):
-            raise ValueError("Le incertezze (ux) devono essere strettamente maggiori di zero.")
-
-        # --- STEP 1: Calcolo standard (Effetti Fissi) per stimare il Chi2 originale ---
-        w_fixed = 1 / ux**2
-        sum_w_fixed = np.sum(w_fixed)
-        xm_fixed = np.sum(w_fixed * x) / sum_w_fixed
-        chi2_fixed = np.sum(w_fixed * (x - xm_fixed) ** 2)
-
-        # --- STEP 2: Stima della varianza extra (Tau^2) ---
-        # Formula di DerSimonian-Laird
-        sum_w2_fixed = np.sum(w_fixed**2)
-        numerator = chi2_fixed - dof
-        denominator = sum_w_fixed - (sum_w2_fixed / sum_w_fixed)
-
-        # Tau^2 non può essere negativo. Se il Chi2 è minore del dof, Tau^2 = 0
-        tau2 = max(0.0, numerator / denominator)
-
-        # --- STEP 3: Ricalcolo con Effetti Casuali (Random Effects) ---
-        w_random = 1 / (ux**2 + tau2)
-        sum_w_random = np.sum(w_random)
-
-        xm_random = np.sum(w_random * x) / sum_w_random
-        uxm_random = 1 / np.sqrt(sum_w_random)
-
-        # Nota: Il p-value finale del Chi2 si riferisce alla consistenza iniziale dei dati
-        p_value = stats.chi2.sf(chi2_fixed, dof)
-
-        return {
-            "media": xm_random,
-            "errore_media": uxm_random,
-            "tau2": tau2,
-            "tau": np.sqrt(tau2),
-            "chi2_originale": chi2_fixed,
-            "dof": dof,
-            "p_value": p_value
-        }
-
-    return
-
-
@app.cell(hide_code=True)
 def _(mo):
    mo.md(r"""
@@ -0,0 +1,262 @@
+import marimo
+
+__generated_with = "0.25.2"
+app = marimo.App(width="medium")
+
+
+@app.cell
+def _():
+    import marimo as mo
+    import numpy as np
+    import pandas as pd
+    import sympy as sp
+    from scipy import stats
+
+    from pathlib import Path
+
+    import statistica
+
+    _outputs, defs = statistica.app.run()
+    Data = defs["Data"]
+    compat = defs["compat"]
+    media_pesata = defs["media_pesata"]
+    rimuovi_outlier = defs["rimuovi_outlier"]
+    plot_gauss = defs["plot_gauss"]
+    cartella_csv = Path(__file__).resolve().parent / "csv" / "velocita-C"
+    # cartella_immagini = Path(__file__).resolve().parent / "LaTeX" / "immagini"
+
+    def leggi_csv(nome):
+        dati = pd.read_csv(cartella_csv / nome, sep=";", decimal=",")
+        dati.columns = dati.columns.astype(str).str.strip()
+        dati = dati.dropna(axis=1, how="all")
+        return dati.apply(pd.to_numeric, errors="coerce")
+
+    return Data, compat, leggi_csv, mo, pd, sp
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    # Velocità della luce — metodo di Foucault
+
+    $$c = \frac{4 f_2 D^2 \left|\omega_{fin}-\omega_{in}\right|}{(D+a-f_2)\left|x_{fin}-x_{in}\right|},
+    \qquad \omega = 2\pi f_r$$
+
+    dove $\omega_{in},\ \omega_{fin}$ sono le pulsazioni di rotazione e
+    $x_{in},\ x_{fin}$ le posizioni del fascio (zero arbitrario) relative alle
+    due letture dello stesso punto sperimentale.
+
+    Le stime di $c$ dei singoli punti vengono combinate con una media pesata.
+    """)
+    return
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Costanti geometriche (parte sistematica)
+    """)
+    return
+
+
+@app.cell
+def _(Data, leggi_csv, pd, sp):
+    _D1 = Data(leggi_csv("D1.csv"))
+    _D1.stat_chauv("D", 0.01, False)
+
+    _D2 = Data(leggi_csv("D2.csv"))
+    _D2.stat_chauv("D", 0.01, False)
+
+    _D3 = Data(leggi_csv("D3.csv"))
+    _D3.stat_chauv("D", 0.01, False)
+
+    _D_tot = Data(None, analisi_stat = pd.DataFrame({
+        "L1": [_D1.analisi_stat.D[0]], "uL1": [_D1.analisi_stat.uD[0]],
+        "L2": [_D2.analisi_stat.D[0]], "uL2": [_D2.analisi_stat.uD[0]],
+        "L3": [_D3.analisi_stat.D[0]], "uL3": [_D3.analisi_stat.uD[0]]
+    }))
+
+    L1, L2, L3 = sp.symbols("L1 L2 L3", positive=True)
+    D = L1 + L2 + L3
+    _D_tot.calc_var(D, "D")
+
+    D, uD = _D_tot.analisi_stat.D[0], _D_tot.analisi_stat.uD[0]  # Distanza specchio rotante - specchio concavo [m]
+    a, ua = 0.376, 0.002     # Distanza lente - specchio rotante [m]
+    f2, uf2 = 0.2536, 0.0005 # Lunghezza focale lente L2 [m]
+    return D, a, f2, uD, ua, uf2
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Dati grezzi
+
+    Dati simulati, usati solo per verificare il procedimento.
+
+    - `C.csv`: una riga per punto sperimentale, con le quattro letture
+      `f_iniziale`, `f_finale` in Hz e `x_iniziale`, `x_finale` in mm;
+    - le frequenze sono negative per il verso di rotazione assunto;
+    - `D1.csv`, `D2.csv`, `D3.csv`: le tre tratte che compongono
+      $D=L_1+L_2+L_3$.
+
+    Dopo l'analisi statistica le posizioni vengono convertite da mm a m.
+    """)
+    return
+
+
+@app.cell
+def _(Data, leggi_csv):
+    err_strum = {
+        "f_finale"  : 1, 
+        "f_iniziale": 1,
+        "x_finale"  : 0.01,
+        "x_iniziale": 0.01
+    }  # risoluzione strumentale: Hz, mm
+
+    dati = Data(leggi_csv("C.csv"))
+    dati.stat_chauv(list(err_strum.keys()), list(err_strum.values()))
+    return (dati,)
+
+
+@app.cell
+def _(D, a, dati, f2, uD, ua, uf2):
+    dati.analisi_stat["D"]   = D
+    dati.analisi_stat["uD"]  = uD
+    dati.analisi_stat["a"]   = a
+    dati.analisi_stat["ua"]  = ua
+    dati.analisi_stat["f2"]  = f2
+    dati.analisi_stat["uf2"] = uf2
+
+    dati.analisi_stat
+    return
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Analisi con criterio di Chauvenet
+
+    Stessa analisi (`stat_chauv`) applicata sia ai punti con rotazione
+    sia alla baseline. Poi si converte `d` da mm a m e si riportano
+    $f_{r0}, d_0$ (baseline) come colonne costanti su tutte le righe dei
+    punti principali, dato che è la stessa misura per tutti.
+    """)
+    return
+
+
+@app.cell
+def _(dati):
+    # conversione mm -> m per le colonne di lunghezza
+    for _col in ["x_finale", "ux_finale", "x_iniziale", "ux_iniziale"]:
+        dati.analisi_stat[_col]  = dati.analisi_stat[_col]  * 1e-3
+
+
+    dati.analisi_stat
+    return
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Calcolo delle pulsazioni $\omega_{in}=2\pi f_{in}$ e $\omega_{fin}=2\pi f_{fin}$
+    """)
+    return
+
+
+@app.cell
+def _(dati, sp):
+    _f_iniziale, _f_finale = sp.symbols("f_iniziale f_finale")
+    dati.calc_var(2 * sp.pi * _f_iniziale, "omega_iniziale")
+    dati.calc_var(2 * sp.pi * _f_finale,   "omega_finale")
+
+    dati.analisi_stat
+    return
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Calcolo di $c$ punto per punto
+
+    $D$, $a$ e $f_2$ sono riportate come colonne costanti **con le loro
+    incertezze** (`uD`, `ua`, `uf2`), quindi `calc_var` propaga la parte
+    geometrica insieme a quella sperimentale. Le colonne
+    `%uc_f_iniziale`, `%uc_f_finale`, `uc_x_...` danno il peso percentuale di
+    ciascun termine e mostrano quale incertezza domina il risultato.
+    """)
+    return
+
+
+@app.cell(hide_code=True)
+def _(formula_c, mo, sp):
+    mo.md(f"""
+    $c={sp.latex(formula_c)}$
+    """)
+    return
+
+
+@app.cell
+def _(dati, sp):
+    _D, _a, _f2, _omega_iniziale, _omega_finale, _x_iniziale, _x_finale= sp.symbols("D a f2 omega_iniziale omega_finale x_iniziale x_finale", real=True)
+
+    formula_c = 4 * _f2 * _D**2 * sp.Abs(_omega_finale - _omega_iniziale) / ((_D + _a - _f2) * sp.Abs(_x_finale - _x_iniziale))
+    dati.calc_var(formula_c, "c")
+
+    dati.analisi_stat
+    return (formula_c,)
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Media pesata delle stime di $c$
+
+    Le stime indipendenti di $c$ — una per riga di `C.csv` — vengono combinate con
+    pesi $w_i = 1/u_{c,i}^2$. Prima di combinarle, `media_pesata` applica il
+    criterio di Chauvenet fra le stime; i punti rimossi sono riportati in `out`.
+
+    Si riportano $\tilde\chi^2$ (chi quadro ridotto), i gradi di libertà e
+    $P(\chi^2)$ riferiti ai soli punti mantenuti.
+    """)
+    return
+
+
+@app.cell
+def _(dati):
+    c_stat, uc_stat, chi2, dof, P, out= dati.media_pesata("c", info_outlier=True)
+
+    print(f"c (media pesata) = {c_stat:.6e} +- {uc_stat:.2e}")
+    print(f"chi2 = {chi2:.3f}, dof = {dof}, P(chi2) = {P:.3f}")
+    print(out)
+    return c_stat, uc_stat
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Risultato finale e confronto col valore vero
+    """)
+    return
+
+
+@app.cell
+def _(c_stat, compat, mo, uc_stat):
+    c_vero = 299792458.0
+    _k, _diff = compat(c_stat, c_vero, uc_stat, 0.0)
+
+    mo.md(
+        f"""
+        $$c = \\left({c_stat/1e8:.3f} \\pm {uc_stat/1e8:.3f}\\right)\\times 10^8\\ \\mathrm{{m/s}}$$
+
+        Confronto col valore tabulato
+        $c_{{tab}} = {c_vero/1e8:.3f}\\times 10^8\\ \\mathrm{{m/s}}$:
+        differenza ${_diff:.2E}$ m/s, fattore di compatibilita **k = {_k:.2f}**,
+        calcolato con $u_{{tab}}=0$, cioe $k=|c-c_{{tab}}|/u_c$.
+        {"✅ compatibile ($k<2$)" if _k < 2 else "⚠️ non compatibile ($k\\ge 2$)"}
+        """
+    )
+    return
+
+
+if __name__ == "__main__":
+    app.run()