diff --git a/focaleLente.py b/focaleLente.py
new file mode 100644
index 0000000..77f0167
--- /dev/null
+++ b/focaleLente.py
@@ -0,0 +1,404 @@
+import marimo
+
+__generated_with = "0.25.0"
+app = marimo.App(width="medium")
+
+
+@app.cell
+def _():
+    import marimo as mo
+
+    return (mo,)
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    # Misura della distanza focale di una lente convergente
+
+    Quattro stime, tutte basate sull'equazione della lente sottile $\frac{1}{f}=\frac{1}{p}+\frac{1}{q}$ e sull'ingrandimento $m=\frac{q}{p}=\frac{l_{Freccia}}{L}$ (con $L$ = lunghezza dell'oggetto):
+
+    1. **Punti coniugati**: $f=\dfrac{p\,q}{p+q}$
+    2. **a) $p$ e ingrandimento**: $f=\dfrac{p\,m}{1+m}$
+    2. **b) $q$ e ingrandimento**: $f=\dfrac{q}{1+m}$
+    3. **Bessel**: $f=\dfrac{D^2-d^2}{4D}$ (vedi sotto)
+
+    Le incertezze sono propagate con `Data.calc_var` di `statistica.py`.
+    """)
+    return
+
+
+@app.cell
+def _():
+    import numpy as np
+    import pandas as pd
+    import sympy as sp
+    from scipy import stats
+
+    import statistica
+
+    # Esegue statistica.py
+    _outputs, defs = statistica.app.run()
+
+    Data = defs["Data"]
+    compat = defs["compat"]
+    plot_gauss = defs["plot_gauss"]
+    return Data, compat, np, pd, plot_gauss, sp, stats
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Dati
+
+    **Le righe sono a coppie (0-1, 2-3, 4-5, ...)** per il metodo di Bessel: le due righe di una coppia hanno la **stessa distanza sorgente-schermo** $D=p+q$ e corrispondono alle due posizioni della lente con immagine nitida.
+
+    In una posizione la lente è vicina alla sorgente (immagine ingrandita, $p$ piccolo), nell'altra è vicina allo schermo (immagine rimpicciolita, $p$ grande): le due posizioni sono **coniugate**, cioè $p$ e $q$ si scambiano ($p_B=q_A$, $q_B=p_A$) e gli ingrandimenti sono uno il reciproco dell'altro.
+
+    Ordine delle righe dentro la coppia indifferente.
+    """)
+    return
+
+
+@app.cell
+def _(Data, pd):
+    dfFenditura = pd.DataFrame({
+        "l1": [10.05,],
+        "l2": [10.00,],
+        "l3": [10.03,],
+        "l4": [9.94, ],
+        "l5": [9.99, ],
+    })
+
+    datiFenditura = Data(dfFenditura)
+    datiFenditura.stat_chauv(["l"], [0.02])
+    datiFenditura
+
+    lunghezzaOggetto  = datiFenditura.analisi_stat.l[0]     # variabile che contiene la lunghezza dell'oggetto
+    ulunghezzaOggetto = datiFenditura.analisi_stat.ul[0]    # variabile che contiene l'errore sulla lunghezza dell'oggetto
+    return lunghezzaOggetto, ulunghezzaOggetto
+
+
+@app.cell
+def _(pd):
+    df = pd.DataFrame({
+        "p1": [23.82, 41.22, 20.58, 53.91, 19.38, 65.51],
+        "p2": [23.39, 41.62, 20.94, 54.35, 19.07, 65.62],
+        "p3": [23.49, 41.54, 20.40, 54.72, 19.43, 65.70],
+        "p4": [23.57, 41.57, 20.84, 54.32, 19.79, 65.51],
+        "p5": [23.33, 41.21, 20.32, 54.17, 19.37, 65.19],
+
+        "q1": [41.51, 23.98, 54.07, 21.30, 65.28, 19.65],
+        "q2": [41.51, 23.53, 53.91, 20.80, 65.19, 19.72],
+        "q3": [40.99, 23.37, 54.71, 20.76, 65.23, 19.58],
+        "q4": [41.82, 24.09, 54.80, 20.81, 65.46, 19.17],
+        "q5": [41.69, 23.47, 54.17, 20.69, 65.19, 19.39],
+
+        "lFreccia1": [17.64, 5.58, 26.22, 3.80, 33.77, 2.94],
+        "lFreccia2": [17.67, 5.64, 26.17, 3.75, 33.69, 2.95],
+        "lFreccia3": [17.70, 5.54, 26.21, 3.84, 33.70, 3.03],
+        "lFreccia4": [17.66, 5.71, 26.14, 3.81, 33.77, 3.04],
+        "lFreccia5": [17.66, 5.64, 26.09, 3.88, 33.77, 3.00],
+    })
+
+    # incertezze strumentali (risoluzione) di p, q, lFreccia
+    err_strum = {"p": 0.1, "q": 0.1, "lFreccia": 0.05}
+    df
+    return df, err_strum
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Analisi statistica (criterio di Chauvenet)
+
+    Per ogni configurazione: media, incertezza statistica + strumentale, ed
+    eventuali outlier scartati. Alla tabella si aggiunge la lunghezza
+    dell'oggetto, comune a tutte le righe.
+    """)
+    return
+
+
+@app.cell
+def _(Data, df, err_strum, lunghezzaOggetto, ulunghezzaOggetto):
+    dati = Data(df)
+    dati.stat_chauv(list(err_strum.keys()), list(err_strum.values()))
+
+    dati.analisi_stat["lunghezzaOggetto"]  = lunghezzaOggetto
+    dati.analisi_stat["ulunghezzaOggetto"] = ulunghezzaOggetto
+
+    dati.analisi_stat
+    return (dati,)
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Metodo 1 — punti coniugati
+    """)
+    return
+
+
+@app.cell(hide_code=True)
+def _(f1, mo, sp):
+    mo.md(f"""
+    $f = {sp.latex(f1)}$
+    """)
+    return
+
+
+@app.cell
+def _(dati, sp):
+    _p, _q = sp.symbols("p q", positive=True)
+
+    f1 = _p * _q / (_p + _q)
+
+    dati.calc_var(f1, "f1")
+
+    res_1 = dati.analisi_stat[["f1", "uf1", "%uf1_p", "%uf1_q"]]
+    res_1
+    return f1, res_1
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Metodo 2a — $p$ e ingrandimento
+    """)
+    return
+
+
+@app.cell(hide_code=True)
+def _(f2a, mo, sp):
+    mo.md(f"""
+    $f = {sp.latex(f2a)}$
+    """)
+    return
+
+
+@app.cell
+def _(dati, sp):
+    _p, _lF, _L = sp.symbols("p lFreccia lunghezzaOggetto", positive=True)
+    _m = _lF / _L
+
+    f2a = _p * _m / (1 + _m)
+
+    dati.calc_var(f2a, "f2a")
+
+    res_2a = dati.analisi_stat[
+        ["f2a", "uf2a", "%uf2a_p", "%uf2a_lFreccia", "%uf2a_lunghezzaOggetto"]
+    ]
+    res_2a
+    return f2a, res_2a
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Metodo 2b — $q$ e ingrandimento
+    """)
+    return
+
+
+@app.cell(hide_code=True)
+def _(f2b, mo, sp):
+    mo.md(f"""
+    $f = {sp.latex(f2b)}$
+    """)
+    return
+
+
+@app.cell
+def _(dati, sp):
+    _q, _lF, _L = sp.symbols("q lFreccia lunghezzaOggetto", positive=True)
+    _m = _lF / _L
+
+    f2b = _q / (1 + _m)
+
+    dati.calc_var(f2b, "f2b")
+
+    res_2b = dati.analisi_stat[
+        ["f2b", "uf2b", "%uf2b_q", "%uf2b_lFreccia", "%uf2b_lunghezzaOggetto"]
+    ]
+    res_2b
+    return f2b, res_2b
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Metodo 3 — Bessel
+
+    Sorgente e schermo sono a distanza fissa $D>4f$. Esistono due posizioni della lente, distanti $d$ tra loro, con immagine nitida.
+
+    Per la prima (A) $p_A=\frac{D-d}{2}$, $q_A=\frac{D+d}{2}$, quindi
+    $$f=\frac{p_Aq_A}{p_A+q_A}=\frac{D^2-d^2}{4D}$$
+
+    Ogni coppia di righe dà **una** stima di $f$.
+    Le quattro misure $p_A,q_A,p_B,q_B$ misurano tutte $D$ e $d$ (con A = posizione a $p$ minore):
+
+    $$D=\frac{p_A+q_A+p_B+q_B}{2},\qquad d=\frac{p_B-p_A+q_A-q_B}{2}$$
+
+    **Attenzione alla propagazione:** $D$ e $d$ sono combinazioni delle stesse quattro misure, quindi *non* sono indipendenti.
+
+    Propagare $f(D,d)$ trattandoli come indipendenti sbaglierebbe l'incertezza: $f$ viene perciò propagata direttamente da $p_A,q_A,p_B,q_B$ (che sono misure indipendenti). $D$ e $d$ (con le loro incertezze) sono calcolati solo per i controlli.
+    """)
+    return
+
+
+@app.cell
+def _(dati, pd):
+    _s = dati.analisi_stat
+    assert len(_s) % 2 == 0, "Bessel: servono righe a coppie (0-1, 2-3, ...)"
+
+    _righe, _etichette = [], []
+    for _k in range(0, len(_s), 2):
+        _a, _b = _s.iloc[_k], _s.iloc[_k + 1]
+        if float(_a["p"]) > float(_b["p"]):   # A = posizione con p minore
+            _a, _b = _b, _a
+        _righe.append({
+            "pA": float(_a["p"]), "upA": float(_a["up"]),
+            "qA": float(_a["q"]), "uqA": float(_a["uq"]),
+            "pB": float(_b["p"]), "upB": float(_b["up"]),
+            "qB": float(_b["q"]), "uqB": float(_b["uq"]),
+        })
+        _etichette.append(f"{_k}-{_k + 1}")
+
+    coppie = pd.DataFrame(_righe, index=_etichette)
+    coppie
+    return (coppie,)
+
+
+@app.cell
+def _(Data, coppie, dati, sp):
+    _pA, _qA, _pB, _qB = sp.symbols("pA qA pB qB", positive=True)
+    _D = (_pA + _qA + _pB + _qB) / 2
+    _d = (_pB - _pA + _qA - _qB) / 2
+
+    f3 = (_D**2 - _d**2) / (4 * _D)
+
+    datiBessel = Data(dati.campione, analisi_stat=coppie.copy())
+    datiBessel.calc_var(_D, "D") 
+    datiBessel.calc_var(_d, "d") 
+    datiBessel.calc_var(f3, "f3")
+
+    res_3 = datiBessel.analisi_stat[
+        ["D", "uD", "d", "ud", "f3", "uf3", "%uf3_pA", "%uf3_qA", "%uf3_pB", "%uf3_qB"]
+    ]
+    res_3.index = coppie.index
+    res_3
+    return (res_3,)
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ### Controlli di consistenza del metodo di Bessel
+
+    - $D$ deve essere uguale nelle due righe: $p_A+q_A$ vs $p_B+q_B$
+    - $d$ misurato con $p$ ($p_B-p_A$) e con $q$ ($q_A-q_B$) deve coincidere
+    - condizione di applicabilità $D>4f$
+
+    Compatibilità con `compat` (k < 2).
+    """)
+    return
+
+
+@app.cell
+def _(compat, coppie, np, pd, res_3):
+    _c = coppie
+    _DA, _uDA = _c.pA + _c.qA, np.sqrt(_c.upA**2 + _c.uqA**2)
+    _DB, _uDB = _c.pB + _c.qB, np.sqrt(_c.upB**2 + _c.uqB**2)
+    _dp, _udp = _c.pB - _c.pA, np.sqrt(_c.upA**2 + _c.upB**2)
+    _dq, _udq = _c.qA - _c.qB, np.sqrt(_c.uqA**2 + _c.uqB**2)
+
+    _kD, _ = compat(_DA, _DB, _uDA, _uDB)
+    _kd, _ = compat(_dp, _dq, _udp, _udq)
+
+    controlli_bessel = pd.DataFrame({
+        "D_A": _DA, "D_B": _DB, "k(D)": _kD, "D compatibili": _kD < 2,
+        "d da p": _dp, "d da q": _dq, "k(d)": _kd, "d compatibili": _kd < 2,
+        "D > 4f": res_3["D"] > 4 * res_3["f3"],
+    })
+    controlli_bessel
+    return
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Risultati finali
+
+    Per ogni metodo si combinano le stime (una per configurazione; una per
+    coppia nel metodo 3) con la **media pesata** ($w_i = 1/u_i^2$). Si
+    riportano anche $\chi^2$ e $P(\chi^2)$ di compatibilità delle stime con
+    un valore costante: se $\chi^2/\mathrm{dof}>1$ l'incertezza va ingrandita
+    di $\sqrt{\chi^2/\mathrm{dof}}$ (colonna `u_corr`).
+    """)
+    return
+
+
+@app.cell
+def _(np, pd, res_1, res_2a, res_2b, res_3, stats):
+    def media_pesata(f, uf):
+        f, uf = np.asarray(f, float), np.asarray(uf, float)
+        w   = 1 / uf**2
+        fm  = np.sum(w * f) / np.sum(w)
+        ufm = 1 / np.sqrt(np.sum(w))
+        chi2 = np.sum(w * (f - fm) ** 2)
+        dof  = len(f) - 1
+        return fm, ufm, chi2, dof, stats.chi2.sf(chi2, dof)
+
+    _righe = {}
+    for _nome, _res, _c in [
+        ("1  - punti coniugati", res_1, "f1"),
+        ("2a - p e ingrandimento", res_2a, "f2a"),
+        ("2b - q e ingrandimento", res_2b, "f2b"),
+        ("3  - Bessel", res_3, "f3"),
+    ]:
+        _fm, _u, _chi2, _dof, _P = media_pesata(_res[_c], _res[f"u{_c}"])
+        _righe[_nome] = {
+            "f [cm]": _fm, "u": _u,
+            "chi2": _chi2, "dof": _dof, "P": _P,
+            "u_corr": _u * np.sqrt(max(1.0, _chi2 / _dof)),
+        }
+
+    risultati = pd.DataFrame(_righe).T
+    risultati
+    return (risultati,)
+
+
+@app.cell(hide_code=True)
+def _(mo):
+    mo.md(r"""
+    ## Compatibilità tra i metodi e confronto grafico
+    """)
+    return
+
+
+@app.cell
+def _(compat, pd, risultati):
+    _nomi = list(risultati.index)
+    _righe = []
+    for _i in range(len(_nomi)):
+        for _j in range(_i + 1, len(_nomi)):
+            _a, _b = risultati.loc[_nomi[_i]], risultati.loc[_nomi[_j]]
+            _k, _diff = compat(_a["f [cm]"], _b["f [cm]"], _a["u_corr"], _b["u_corr"])
+            _righe.append({"metodi": f"{_nomi[_i][:2].strip()} vs {_nomi[_j][:2].strip()}",
+                           "differenza": _diff, "k": _k,
+                           "compatibili (k<2)": bool(_k < 2)})
+    pd.DataFrame(_righe).set_index("metodi")
+    return
+
+
+@app.cell
+def _(plot_gauss, risultati):
+    plot_gauss([
+        (r["f [cm]"], r["u_corr"], _i, _n.strip())
+        for _i, (_n, r) in enumerate(risultati.iterrows())
+    ])
+    return
+
+
+if __name__ == "__main__":
+    app.run()
diff --git a/src/lofm/__init__.py b/src/lofm/__init__.py
new file mode 100644
index 0000000..a4b66b4
--- /dev/null
+++ b/src/lofm/__init__.py
@@ -0,0 +1,2 @@
+def main() -> None:
+    print("Hello from lofm!")