diff --git a/doc/locale/fr/LC_MESSAGES/release_notes/release_1.01.po b/doc/locale/fr/LC_MESSAGES/release_notes/release_1.01.po
index 241c5a7..e9c9d04 100644
--- a/doc/locale/fr/LC_MESSAGES/release_notes/release_1.01.po
+++ b/doc/locale/fr/LC_MESSAGES/release_notes/release_1.01.po
@@ -30,8 +30,39 @@ msgstr "Les annotations PlotPy existantes restent lisibles sans modifier l'objet
 msgid "Canonical annotation identifiers, metadata, extensions, and persistent lock state are preserved across PlotPy editing round trips. Application-specific opaque annotation entries continue to coexist with graphical annotations."
 msgstr "Les identifiants d'annotation canoniques, les métadonnées, les extensions et l'état de verrouillage persistant sont conservés lors des allers-retours d'édition PlotPy. Les entrées d'annotation opaques spécifiques à l'application continuent de coexister avec les annotations graphiques."
 
+msgid "Main window"
+msgstr ""
+
+msgid "With Qt 5, the dock widget sizes saved on exit are now restored at startup: the window is sized and laid out before its state is restored, so tabified docks no longer shrink to the default window size and leave the extra width to the central widget."
+msgstr ""
+
+msgid "Qt helpers"
+msgstr ""
+
+msgid "`block_signals` now restores the previous signal-blocking state of the widget and of its children instead of unblocking them unconditionally: nested blocking contexts and widgets already blocked by the caller keep their state."
+msgstr ""
+
+msgid "`sigimax_app_context` no longer schedules the unattended/screenshot close timer when the Qt event loop is not executed (`exec_loop=False`): the pending timer previously closed unrelated top-level widgets created later, such as the window of a subsequent test."
+msgstr ""
+
+msgid "Curve fitting dialogs"
+msgstr "Fenêtres d'ajustement de courbes"
+
+msgid "Interactive fitting dialogs now show the formula of the fitted model below the fit parameters, and parameter labels show the symbols used in the formula (e.g. \"Scale factor (x₀)\"), so that the parameters stored in the fitted signal metadata can be related to the model. `guifit` accepts a new `formula` argument for this purpose. This closes [Issue #9](https://github.com/DataLab-Platform/SigimaX/issues/9)."
+msgstr "Les fenêtres d'ajustement interactif affichent désormais la formule du modèle ajusté sous les paramètres d'ajustement, et les libellés des paramètres indiquent les symboles utilisés dans la formule (par exemple « Facteur d'échelle (x₀) »), afin de pouvoir relier au modèle les paramètres enregistrés dans les métadonnées du signal ajusté. `guifit` accepte pour cela un nouvel argument `formula`. Cela corrige [l'issue #9](https://github.com/DataLab-Platform/SigimaX/issues/9)."
+
+msgid "Polynomial and linear fit coefficients are now labelled `a`, `b`, `c`... as in the model formula and in the fitted signal metadata, instead of `c2`, `c1`, `c0`."
+msgstr "Les coefficients des ajustements polynomial et linéaire sont désormais nommés `a`, `b`, `c`... comme dans la formule du modèle et dans les métadonnées du signal ajusté, au lieu de `c2`, `c1`, `c0`."
+
+msgid "The Planckian fit dialog starts from a curve matching the data peak (its initial amplitude was about 147 times too low on blackbody data in physical units), and its \"width factor\" is now labelled \"temperature factor\": for wavelengths in meters, the temperature is `h c σ / (5 k_B x₀)`."
+msgstr "La fenêtre d'ajustement planckien part d'une courbe correspondant au pic des données (son amplitude initiale était environ 147 fois trop faible sur des données de corps noir en unités physiques), et son « facteur de largeur » est désormais nommé « facteur de température » : pour des longueurs d'onde en mètres, la température vaut `h c σ / (5 k_B x₀)`."
+
+msgid "The exponential fit dialog now evaluates the same model as Sigima, and its rate slider covers short time constants whatever the X unit (e.g. 1 ms with X in seconds)."
+msgstr "La fenêtre d'ajustement exponentiel évalue désormais le même modèle que Sigima, et son curseur de taux couvre les constantes de temps courtes quelle que soit l'unité de X (par exemple 1 ms avec X en secondes)."
+
 msgid "Requirements"
 msgstr "Exigences"
 
-msgid "Sigima 1.3.0 or later is required for the portable annotation model and PlotPy conversion helpers."
-msgstr "Sigima 1.3.0 ou une version ultérieure est requise pour le modèle d'annotation portable et les helpers de conversion PlotPy."
+msgid "Sigima 1.3.0 or later is required for the portable annotation model, PlotPy conversion helpers and curve fitting model formulas."
+msgstr "Sigima 1.3.0 ou une version ultérieure est requise pour le modèle d'annotation portable, les helpers de conversion PlotPy et les formules des modèles d'ajustement de courbes."
+
diff --git a/doc/release_notes/release_1.01.md b/doc/release_notes/release_1.01.md
index dd76ddc..a8b1e0a 100644
--- a/doc/release_notes/release_1.01.md
+++ b/doc/release_notes/release_1.01.md
@@ -24,7 +24,13 @@
 * `block_signals` now restores the previous signal-blocking state of the widget and of its children instead of unblocking them unconditionally: nested blocking contexts and widgets already blocked by the caller keep their state.
 * `sigimax_app_context` no longer schedules the unattended/screenshot close timer when the Qt event loop is not executed (`exec_loop=False`): the pending timer previously closed unrelated top-level widgets created later, such as the window of a subsequent test.
 
+### Curve fitting dialogs ###
+
+* Interactive fitting dialogs now show the formula of the fitted model below the fit parameters, and parameter labels show the symbols used in the formula (e.g. "Scale factor (x₀)"), so that the parameters stored in the fitted signal metadata can be related to the model. `guifit` accepts a new `formula` argument for this purpose. This closes [Issue #9](https://github.com/DataLab-Platform/SigimaX/issues/9).
+* Polynomial and linear fit coefficients are now labelled `a`, `b`, `c`... as in the model formula and in the fitted signal metadata, instead of `c2`, `c1`, `c0`.
+* The Planckian fit dialog starts from a curve matching the data peak (its initial amplitude was about 147 times too low on blackbody data in physical units), and its "width factor" is now labelled "temperature factor": for wavelengths in meters, the temperature is `h c σ / (5 k_B x₀)`.
+* The exponential fit dialog now evaluates the same model as Sigima, and its rate slider covers short time constants whatever the X unit (e.g. 1 ms with X in seconds).
+
 ### Requirements ###
 
-* Sigima 1.3.0 or later is required for the portable annotation model and
-  PlotPy conversion helpers.
\ No newline at end of file
+* Sigima 1.3.0 or later is required for the portable annotation model, PlotPy conversion helpers and curve fitting model formulas.
\ No newline at end of file
diff --git a/sigimax/locale/fr/LC_MESSAGES/sigimax.po b/sigimax/locale/fr/LC_MESSAGES/sigimax.po
index df17ee0..810d59e 100644
--- a/sigimax/locale/fr/LC_MESSAGES/sigimax.po
+++ b/sigimax/locale/fr/LC_MESSAGES/sigimax.po
@@ -332,6 +332,10 @@ msgstr "Connexion échouée : %s"
 msgid "Contents of file"
 msgstr "Contenu du fichier"
 
+#, python-format
+msgid "Model: %s"
+msgstr "Modèle : %s"
+
 msgid "Polymomial fit"
 msgstr "Ajustement polynomial"
 
@@ -395,8 +399,8 @@ msgstr "Ajustement CDF"
 msgid "Scale factor"
 msgstr "Facteur d'échelle"
 
-msgid "Width factor"
-msgstr "Facteur de largeur"
+msgid "Temperature factor"
+msgstr "Facteur de température"
 
 msgid "Planckian fit"
 msgstr "Ajustement planckien"
diff --git a/sigimax/tests/widgets/test_fitdialog.py b/sigimax/tests/widgets/test_fitdialog.py
index 1d88319..ca48b04 100644
--- a/sigimax/tests/widgets/test_fitdialog.py
+++ b/sigimax/tests/widgets/test_fitdialog.py
@@ -11,7 +11,13 @@
 import numpy as np
 import pytest
 from guidata.qthelpers import qt_app_context
-from sigima.objects import NormalDistribution1DParam
+from qtpy import QtCore as QC
+from qtpy import QtWidgets as QW
+from sigima.objects import (
+    NormalDistribution1DParam,
+    PlanckParam,
+    create_signal_from_param,
+)
 from sigima.tests.data import create_noisy_signal, get_test_signal
 from sigima.tools.signal import fitting, pulse
 from sigima.tools.signal.peakdetection import peak_indices
@@ -258,5 +264,104 @@ def test_fit_dialog_bounds_contain_the_optimum(
         )
 
 
+FORMULA_CASES = (
+    (fdlg.linear_fit, None, fitting.LinearFitComputer.FORMULA),
+    (
+        lambda x, y: fdlg.polynomial_fit(x, y, 3),
+        None,
+        "y = a x<sup>3</sup> + b x<sup>2</sup> + c x + d",
+    ),
+    (fdlg.gaussian_fit, None, fitting.GaussianFitComputer.FORMULA),
+    (fdlg.lorentzian_fit, None, fitting.LorentzianFitComputer.FORMULA),
+    (fdlg.voigt_fit, None, fitting.VoigtFitComputer.FORMULA),
+    (fdlg.multigaussian_fit, "peaks", fitting.MultiGaussianFitComputer.FORMULA),
+    (fdlg.multilorentzian_fit, "peaks", fitting.MultiLorentzianFitComputer.FORMULA),
+    (fdlg.exponential_fit, None, fitting.ExponentialFitComputer.FORMULA),
+    (fdlg.sinusoidal_fit, None, fitting.SinusoidalFitComputer.FORMULA),
+    (fdlg.cdf_fit, None, fitting.CDFFitComputer.FORMULA),
+    (fdlg.planckian_fit, None, fitting.PlanckianFitComputer.FORMULA),
+    (fdlg.twohalfgaussian_fit, None, fitting.TwoHalfGaussianFitComputer.FORMULA),
+    (
+        fdlg.piecewiseexponential_fit,
+        None,
+        fitting.DoubleExponentialFitComputer.FORMULA,
+    ),
+)
+
+
+@pytest.mark.parametrize(("dialog", "extra", "formula"), FORMULA_CASES)
+def test_fit_dialogs_show_model_formula(monkeypatch, dialog, extra, formula):
+    """Every interactive fit dialog displays the formula of its Sigima model."""
+    captured = {}
+
+    def accept_initial_values(_x, _y, _fitfunc, fitparams, **kwargs):
+        captured.update(kwargs)
+        return [param.value for param in fitparams]
+
+    monkeypatch.setattr(fdlg, "guifit", accept_initial_values)
+    signal = get_test_signal("paracetamol.txt")
+    if extra == "peaks":
+        assert dialog(signal.x, signal.y, peak_indices(signal.y)) is not None
+    else:
+        assert dialog(signal.x, signal.y) is not None
+    assert captured["formula"] == formula
+
+
+def test_guifit_displays_formula(monkeypatch):
+    """`guifit` adds a selectable label showing the model formula."""
+    dialogs = []
+
+    def reject_dialog(dialog):
+        dialogs.append(dialog)
+        return False
+
+    monkeypatch.setattr(fdlg, "exec_dialog", reject_dialog)
+    x = np.linspace(0.0, 1.0, 50)
+    y = 2.0 * x + 1.0
+    with qt_app_context():
+        params = [fdlg.FitParam("a", 1.0, -5.0, 5.0), fdlg.FitParam("b", 0.0, -5, 5)]
+        formula = fitting.LinearFitComputer.FORMULA
+        result = fdlg.guifit(
+            x, y, lambda x, p: p[0] * x + p[1], params, formula=formula
+        )
+        assert result is None
+        label = dialogs[0].findChild(QW.QLabel, "fit_formula")
+        assert label is not None
+        assert formula in label.text()
+        assert label.textInteractionFlags() & QC.Qt.TextSelectableByMouse
+
+
+def test_exponential_dialog_reaches_short_time_constants(monkeypatch):
+    """The rate slider reaches b = -1000 for a 1 ms time constant (X in s).
+
+    The dialog model is Sigima's, so the stored parameters reproduce the curve.
+    """
+    captured = __capture_fit_params(monkeypatch)
+    x = np.linspace(0.0, 5e-3, 500)
+    y = 2.0 * np.exp(-x / 1e-3) + 0.1
+
+    y_fitted, _params, fit_params = fdlg.exponential_fit(x, y)
+
+    rate = captured[1]
+    assert rate.min <= -1000.0 <= rate.max
+    np.testing.assert_allclose(fitting.evaluate_fit(x, **fit_params), y_fitted)
+
+
+def test_planckian_dialog_starts_on_blackbody_peak(monkeypatch):
+    """On blackbody data, the initial Planckian curve matches the data peak.
+
+    The initial amplitude used to be about 147 times too low, and the initial
+    scale factor could lie outside its slider range.
+    """
+    captured = __capture_fit_params(monkeypatch)
+    signal = create_signal_from_param(PlanckParam.create(T=296.0))
+
+    y_initial, _params, _fit_params = fdlg.planckian_fit(signal.x, signal.y)
+
+    assert y_initial.max() == pytest.approx(signal.y.max(), rel=0.01)
+    for param in captured:
+        assert param.min <= param.value <= param.max, param.name
+
+
 if __name__ == "__main__":
     test_fit_dialog()
diff --git a/sigimax/widgets/fitdialog.py b/sigimax/widgets/fitdialog.py
index 6799d9b..142b718 100644
--- a/sigimax/widgets/fitdialog.py
+++ b/sigimax/widgets/fitdialog.py
@@ -3,6 +3,7 @@
 """Curve fitting dialog widgets
 
 .. autofunction:: guifit
+.. autofunction:: create_formula_label
 .. autofunction:: linear_fit
 .. autofunction:: polynomial_fit
 .. autofunction:: gaussian_fit
@@ -20,11 +21,15 @@
 
 # pylint: disable=invalid-name  # Allows short reference names like x, y, ...
 
+from __future__ import annotations
+
 import numpy as np
 from guidata.configtools import get_icon
 from guidata.qthelpers import exec_dialog
 from plotpy.plot import PlotOptions
 from plotpy.widgets.fit import FitDialog, FitParam
+from qtpy import QtCore as QC
+from qtpy import QtWidgets as QW
 from scipy.special import erf  # pylint: disable=no-name-in-module
 from sigima.tools.checks import check_1d_arrays
 from sigima.tools.signal import fitting, fourier, pulse
@@ -33,6 +38,7 @@
 
 __all__ = [
     "cdf_fit",
+    "create_formula_label",
     "exponential_fit",
     "gaussian_fit",
     "guifit",
@@ -59,6 +65,25 @@ def create_interactive_fit_params(fit_type, values, y, y_fitted):
     )
 
 
+def create_formula_label(formula: str, parent: QW.QWidget | None = None) -> QW.QLabel:
+    """Create a selectable rich-text label showing a fit model formula.
+
+    Args:
+        formula: Model formula, as Qt rich text (see
+         :attr:`sigima.tools.signal.fitting.FitComputer.FORMULA`).
+        parent: Parent widget.
+
+    Returns:
+        Formula label.
+    """
+    label = QW.QLabel(_("Model: %s") % formula, parent)
+    label.setObjectName("fit_formula")
+    label.setTextFormat(QC.Qt.RichText)
+    label.setTextInteractionFlags(QC.Qt.TextSelectableByMouse)
+    label.setWordWrap(True)
+    return label
+
+
 def guifit(
     x,
     y,
@@ -76,8 +101,12 @@ def guifit(
     winpos=None,
     parent=None,
     name=None,
+    formula=None,
 ):  # pylint: disable=too-many-positional-arguments
-    """GUI-based curve fitting tool"""
+    """GUI-based curve fitting tool
+
+    The optional ``formula`` (Qt rich text) is displayed below the fit
+    parameters, so that the user knows which model is fitted."""
     win = FitDialog(
         edit=True,
         title=wintitle,
@@ -96,6 +125,8 @@ def guifit(
         auto_fit=auto_fit,
     )
     win.setObjectName(name)
+    if formula:
+        win.add_widget(create_formula_label(formula, win))
     win.set_data(x, y, fitfunc, fitparams, fitargs, fitkwargs)
     try:
         win.autofit()  # TODO: [P3] make this optional
@@ -122,24 +153,26 @@ def polynomial_fit(x, y, degree, parent=None, name=None, fit_type="polynomial"):
     computer = fitting.PolynomialFitComputer(x, y, degree)
     ivals = np.polyfit(x, y, degree)
 
+    # Parameters are named after the model formula (and the stored metadata):
+    # `a` for the highest degree coefficient, then `b`, `c`, ...
     params = []
-    for index in range(degree + 1):
-        val = ivals[index]
+    for pname, val in zip(computer.get_params_names(), ivals):
         vmax = max(1.0, np.abs(val))
-        param = FitParam(
-            f"c{(len(ivals) - index - 1):d}",
-            val,
-            -2 * vmax,
-            2 * vmax,
-            format=DEFAULT_FORMAT,
-        )
+        param = FitParam(pname, val, -2 * vmax, 2 * vmax, format=DEFAULT_FORMAT)
         params.append(param)
 
     def fitfunc(x, params):
         return np.polyval(params, x)
 
     values = guifit(
-        x, y, fitfunc, params, parent=parent, wintitle=_("Polymomial fit"), name=name
+        x,
+        y,
+        fitfunc,
+        params,
+        parent=parent,
+        wintitle=_("Polymomial fit"),
+        name=name,
+        formula=computer.get_formula(),
     )
     if values:
         y_fitted = fitfunc(x, values)
@@ -178,14 +211,18 @@ def gaussian_fit(x, y, parent=None, name=None):
     dy = np.max(y) - np.min(y)
     max_amplitude = max(2.0 * dy, 2.0 * abs(amplitude_guess))
     amplitude = FitParam(
-        _("Amplitude"),
+        _("Amplitude") + " (A)",
         amplitude_guess,
         -max_amplitude,
         max_amplitude,
         format=DEFAULT_FORMAT,
     )
     b = FitParam(
-        _("Base line"), b_guess, np.min(y) - 0.1 * dy, np.max(y), format=DEFAULT_FORMAT
+        _("Base line") + " (y₀)",
+        b_guess,
+        np.min(y) - 0.1 * dy,
+        np.max(y),
+        format=DEFAULT_FORMAT,
     )
     sigma = FitParam(
         _("Std-dev") + " (σ)",
@@ -195,7 +232,7 @@ def gaussian_fit(x, y, parent=None, name=None):
         format=DEFAULT_FORMAT,
     )
     mu = FitParam(
-        _("Mean") + " (μ)", mu_guess, np.min(x), np.max(x), format=DEFAULT_FORMAT
+        _("Mean") + " (x₀)", mu_guess, np.min(x), np.max(x), format=DEFAULT_FORMAT
     )
 
     params = [amplitude, sigma, mu, b]
@@ -204,7 +241,14 @@ def fitfunc(x, params):
         return pulse.GaussianModel.evaluate(x, *params)
 
     values = guifit(
-        x, y, fitfunc, params, parent=parent, wintitle=_("Gaussian fit"), name=name
+        x,
+        y,
+        fitfunc,
+        params,
+        parent=parent,
+        wintitle=_("Gaussian fit"),
+        name=name,
+        formula=computer.get_formula(),
     )
     if values:
         y_fitted = fitfunc(x, values)
@@ -233,14 +277,18 @@ def lorentzian_fit(x, y, parent=None, name=None):
 
     max_amplitude = max(2.0 * dy, 2.0 * abs(amplitude_guess))
     amplitude = FitParam(
-        _("Amplitude"),
+        _("Amplitude") + " (A)",
         amplitude_guess,
         -max_amplitude,
         max_amplitude,
         format=DEFAULT_FORMAT,
     )
     b = FitParam(
-        _("Base line"), b_guess, np.min(y) - 0.1 * dy, np.max(y), format=DEFAULT_FORMAT
+        _("Base line") + " (y₀)",
+        b_guess,
+        np.min(y) - 0.1 * dy,
+        np.max(y),
+        format=DEFAULT_FORMAT,
     )
     sigma = FitParam(
         _("Std-dev") + " (σ)",
@@ -250,7 +298,7 @@ def lorentzian_fit(x, y, parent=None, name=None):
         format=DEFAULT_FORMAT,
     )
     mu = FitParam(
-        _("Mean") + " (μ)", mu_guess, np.min(x), np.max(x), format=DEFAULT_FORMAT
+        _("Mean") + " (x₀)", mu_guess, np.min(x), np.max(x), format=DEFAULT_FORMAT
     )
 
     params = [amplitude, sigma, mu, b]
@@ -259,7 +307,14 @@ def fitfunc(x, params):
         return pulse.LorentzianModel.evaluate(x, *params)
 
     values = guifit(
-        x, y, fitfunc, params, parent=parent, wintitle=_("Lorentzian fit"), name=name
+        x,
+        y,
+        fitfunc,
+        params,
+        parent=parent,
+        wintitle=_("Lorentzian fit"),
+        name=name,
+        formula=computer.get_formula(),
     )
     if values:
         y_fitted = fitfunc(x, values)
@@ -288,14 +343,18 @@ def voigt_fit(x, y, parent=None, name=None):
 
     max_amplitude = max(2.0 * dy, 2.0 * abs(amplitude_guess))
     amplitude = FitParam(
-        _("Amplitude"),
+        _("Amplitude") + " (A)",
         amplitude_guess,
         -max_amplitude,
         max_amplitude,
         format=DEFAULT_FORMAT,
     )
     b = FitParam(
-        _("Base line"), b_guess, np.min(y) - 0.1 * dy, np.max(y), format=DEFAULT_FORMAT
+        _("Base line") + " (y₀)",
+        b_guess,
+        np.min(y) - 0.1 * dy,
+        np.max(y),
+        format=DEFAULT_FORMAT,
     )
     sigma = FitParam(
         _("Std-dev") + " (σ)",
@@ -305,7 +364,7 @@ def voigt_fit(x, y, parent=None, name=None):
         format=DEFAULT_FORMAT,
     )
     mu = FitParam(
-        _("Mean") + " (μ)", mu_guess, np.min(x), np.max(x), format=DEFAULT_FORMAT
+        _("Mean") + " (x₀)", mu_guess, np.min(x), np.max(x), format=DEFAULT_FORMAT
     )
 
     params = [amplitude, sigma, mu, b]
@@ -314,7 +373,14 @@ def fitfunc(x, params):
         return pulse.VoigtModel.evaluate(x, *params)
 
     values = guifit(
-        x, y, fitfunc, params, parent=parent, wintitle=_("Voigt fit"), name=name
+        x,
+        y,
+        fitfunc,
+        params,
+        parent=parent,
+        wintitle=_("Voigt fit"),
+        name=name,
+        formula=computer.get_formula(),
     )
     if values:
         y_fitted = fitfunc(x, values)
@@ -409,6 +475,7 @@ def fitfunc(xi, params):
         parent=parent,
         name=name,
         wintitle=_("Multi-Gaussian fit"),
+        formula=computer.get_formula(),
     )
     if values:
         y_fitted = fitfunc(x, values)
@@ -508,6 +575,7 @@ def fitfunc(xi, params):
         parent=parent,
         name=name,
         wintitle=_("Multi-Lorentzian fit"),
+        formula=computer.get_formula(),
     )
     if values:
         y_fitted = fitfunc(x, values)
@@ -532,7 +600,7 @@ def exponential_fit(x: np.ndarray, y: np.ndarray, parent=None, name=None):
 
     Returns (yfit, params, fit_params), where yfit is the fitted curve, params are
     the fitting parameters and fit_params is the canonical metadata dictionary"""
-    # Get initial parameter estimates from Sigima ExponentialFitComputer
+    # Get initial parameter estimates and bounds from Sigima ExponentialFitComputer
     computer = fitting.ExponentialFitComputer(x, y)
     initial_params = computer.compute_initial_params()
     oa = initial_params["a"]
@@ -540,26 +608,31 @@ def exponential_fit(x: np.ndarray, y: np.ndarray, parent=None, name=None):
     oc = initial_params["y0"]
 
     # Create parameter bounds
-    moa, mob, moc = np.maximum(1, [abs(oa), abs(ob), abs(oc)])
+    moa, _mob, moc = np.maximum(1, [abs(oa), abs(ob), abs(oc)])
     a_p = FitParam(
         _("A coefficient"), oa, -2 * moa, 2 * moa, logscale=True, format=DEFAULT_FORMAT
     )
     # B must be free to change sign: a positive-only range makes every decaying
-    # exponential unreachable. Sigima uses (-10, 10) for the same parameter.
-    mob = max(10.0, 2 * mob)
-    b_p = FitParam(_("B coefficient"), ob, -mob, mob, format=DEFAULT_FORMAT)
+    # exponential unreachable. Its range scales with the X range, as in Sigima,
+    # so that short time constants remain reachable whatever the X unit.
+    b_min, b_max = computer.compute_bounds(**initial_params)[1]
+    b_p = FitParam(_("B coefficient"), ob, b_min, b_max, format=DEFAULT_FORMAT)
     c_p = FitParam(_("y0 constant"), oc, -2 * moc, 2 * moc, format=DEFAULT_FORMAT)
 
     params = [a_p, b_p, c_p]
 
-    def modelfunc(x, a, b, c):
-        return a * np.exp(b * x) + c
-
     def fitfunc(x, params):
-        return modelfunc(x, *params)
+        return fitting.ExponentialFitComputer.evaluate(x, *params)
 
     values = guifit(
-        x, y, fitfunc, params, parent=parent, wintitle=_("Exponential fit"), name=name
+        x,
+        y,
+        fitfunc,
+        params,
+        parent=parent,
+        wintitle=_("Exponential fit"),
+        name=name,
+        formula=computer.get_formula(),
     )
     if values:
         y_fitted = fitfunc(x, values)
@@ -603,11 +676,19 @@ def sinusoidal_fit(x: np.ndarray, y: np.ndarray, parent=None, name=None):
     # Create parameter bounds
     abs_values = [abs(guess_a), abs(guess_f), abs(guess_ph), abs(guess_c)]
     moa, mof, _mop, moc = np.maximum(1, abs_values)
-    a_p = FitParam(_("Amplitude"), guess_a, -2 * moa, 2 * moa, format=DEFAULT_FORMAT)
-    f_p = FitParam(_("Frequency"), guess_f, 0, 2 * mof, format=DEFAULT_FORMAT)
-    p_p = FitParam(_("Phase"), guess_ph, -360, 360, format=DEFAULT_FORMAT)
+    a_p = FitParam(
+        _("Amplitude") + " (A)", guess_a, -2 * moa, 2 * moa, format=DEFAULT_FORMAT
+    )
+    f_p = FitParam(_("Frequency") + " (f)", guess_f, 0, 2 * mof, format=DEFAULT_FORMAT)
+    p_p = FitParam(
+        _("Phase") + " (φ)", guess_ph, -360, 360, format=DEFAULT_FORMAT, unit="°"
+    )
     c_p = FitParam(
-        _("Continuous component"), guess_c, -2 * moc, 2 * moc, format=DEFAULT_FORMAT
+        _("Continuous component") + " (y₀)",
+        guess_c,
+        -2 * moc,
+        2 * moc,
+        format=DEFAULT_FORMAT,
     )
 
     params = [a_p, f_p, p_p, c_p]
@@ -619,7 +700,14 @@ def fitfunc(x, params):
         return modelfunc(x, *params)
 
     values = guifit(
-        x, y, fitfunc, params, parent=parent, wintitle=_("Sinusoidal fit"), name=name
+        x,
+        y,
+        fitfunc,
+        params,
+        parent=parent,
+        wintitle=_("Sinusoidal fit"),
+        name=name,
+        formula=computer.get_formula(),
     )
     if values:
         y_fitted = fitfunc(x, values)
@@ -662,10 +750,18 @@ def cdf_fit(x: np.ndarray, y: np.ndarray, parent=None, name=None):
     # Amplitude must be free to change sign, otherwise a descending transition
     # cannot be fitted at all.
     a = FitParam(
-        _("Amplitude"), a_guess, -iamp * 2.0, iamp * 2.0, format=DEFAULT_FORMAT
+        _("Amplitude") + " (A)",
+        a_guess,
+        -iamp * 2.0,
+        iamp * 2.0,
+        format=DEFAULT_FORMAT,
     )
     b = FitParam(
-        _("Base line"), b_guess, np.min(y) - 0.1 * dy, np.max(y), format=DEFAULT_FORMAT
+        _("Base line") + " (y₀)",
+        b_guess,
+        np.min(y) - 0.1 * dy,
+        np.max(y),
+        format=DEFAULT_FORMAT,
     )
     # Bound sigma and mu to the abscissa range rather than to the magnitude of
     # the initial guess, which excluded negative and near-zero means.
@@ -694,6 +790,7 @@ def fitfunc(x, params):
         parent=parent,
         wintitle=_("CDF fit"),
         name=name,
+        formula=computer.get_formula(),
     )
     if values:
         y_fitted = fitfunc(x, values)
@@ -709,34 +806,36 @@ def planckian_fit(x: np.ndarray, y: np.ndarray, parent=None, name=None):
 
     Returns (yfit, params, fit_params), where yfit is the fitted curve, params are
     the fitting parameters and fit_params is the canonical metadata dictionary"""
-    # Get initial parameter estimates from Sigima PlanckianFitComputer
+    # Get initial parameter estimates and bounds from Sigima PlanckianFitComputer
     computer = fitting.PlanckianFitComputer(x, y)
     initial_params = computer.compute_initial_params()
-    amp_guess = initial_params["amp"]
-    x0_guess = initial_params["x0"]
-    sigma_guess = initial_params["sigma"]
-    y0_guess = initial_params["y0"]
-
-    # Create parameter bounds
-    dy = np.max(y) - np.min(y)
+    bounds = dict(
+        zip(computer.get_params_names(), computer.compute_bounds(**initial_params))
+    )
 
-    # Parameter bounds with appropriate ranges for Planckian fitting
+    # `sigma / x0` is proportional to the temperature (see PlanckianFitComputer)
     amp = FitParam(
-        _("Amplitude"),
-        amp_guess,
-        amp_guess * 0.01,
-        amp_guess * 100,
+        _("Amplitude") + " (A)",
+        initial_params["amp"],
+        *bounds["amp"],
         format=DEFAULT_FORMAT,
     )
     x0 = FitParam(
-        _("Scale factor"), x0_guess, np.min(x), np.max(x), format=DEFAULT_FORMAT
+        _("Scale factor") + " (x₀)",
+        initial_params["x0"],
+        *bounds["x0"],
+        format=DEFAULT_FORMAT,
+    )
+    sigma = FitParam(
+        _("Temperature factor") + " (σ)",
+        initial_params["sigma"],
+        *bounds["sigma"],
+        format=DEFAULT_FORMAT,
     )
-    sigma = FitParam(_("Width factor"), sigma_guess, 0.1, 5.0, format=DEFAULT_FORMAT)
     y0 = FitParam(
-        _("Base line"),
-        y0_guess,
-        y0_guess - 0.2 * dy,
-        y0_guess + 0.2 * dy,
+        _("Base line") + " (y₀)",
+        initial_params["y0"],
+        *bounds["y0"],
         format=DEFAULT_FORMAT,
     )
 
@@ -747,7 +846,14 @@ def fitfunc(x, params: list[float]) -> np.ndarray:
         return fitting.PlanckianFitComputer.evaluate(x, *params)
 
     values = guifit(
-        x, y, fitfunc, params, parent=parent, wintitle=_("Planckian fit"), name=name
+        x,
+        y,
+        fitfunc,
+        params,
+        parent=parent,
+        wintitle=_("Planckian fit"),
+        name=name,
+        formula=computer.get_formula(),
     )
     if values:
         y_fitted = fitfunc(x, values)
@@ -782,20 +888,28 @@ def twohalfgaussian_fit(x: np.ndarray, y: np.ndarray, parent=None, name=None):
     # New model signature: func(x, amp_left, amp_right, sigma_left,
     #                            sigma_right, x0, y0_left, y0_right)
     amp_left = FitParam(
-        _("Left amplitude"), amp_left_guess, dy * 0.1, dy * 3, format=DEFAULT_FORMAT
+        _("Left amplitude") + " (A<sub>L</sub>)",
+        amp_left_guess,
+        dy * 0.1,
+        dy * 3,
+        format=DEFAULT_FORMAT,
     )
     amp_right = FitParam(
-        _("Right amplitude"), amp_right_guess, dy * 0.1, dy * 3, format=DEFAULT_FORMAT
+        _("Right amplitude") + " (A<sub>R</sub>)",
+        amp_right_guess,
+        dy * 0.1,
+        dy * 3,
+        format=DEFAULT_FORMAT,
     )
     sigma_left = FitParam(
-        _("Left width") + " (σL)",
+        _("Left width") + " (σ<sub>L</sub>)",
         sigma_left_guess,
         dx * 0.001,  # Very small minimum
         dx * 0.5,  # Reasonable maximum
         format=DEFAULT_FORMAT,
     )
     sigma_right = FitParam(
-        _("Right width") + " (σR)",
+        _("Right width") + " (σ<sub>R</sub>)",
         sigma_right_guess,
         dx * 0.001,  # Very small minimum
         dx * 0.5,  # Reasonable maximum
@@ -805,14 +919,14 @@ def twohalfgaussian_fit(x: np.ndarray, y: np.ndarray, parent=None, name=None):
         _("Center") + " (x₀)", x0_guess, np.min(x), np.max(x), format=DEFAULT_FORMAT
     )
     y0_left = FitParam(
-        _("Left baseline"),
+        _("Left baseline") + " (y<sub>0,L</sub>)",
         y0_left_guess,
         y0_left_guess - 0.2 * dy,
         y0_left_guess + 0.2 * dy,
         format=DEFAULT_FORMAT,
     )
     y0_right = FitParam(
-        _("Right baseline"),
+        _("Right baseline") + " (y<sub>0,R</sub>)",
         y0_right_guess,
         y0_right_guess - 0.2 * dy,
         y0_right_guess + 0.2 * dy,
@@ -832,6 +946,7 @@ def fitfunc(x, params):
         parent=parent,
         wintitle=_("Two half-Gaussian fit"),
         name=name,
+        formula=computer.get_formula(),
     )
     if values:
         y_fitted = fitfunc(x, values)
@@ -875,38 +990,42 @@ def piecewiseexponential_fit(x: np.ndarray, y: np.ndarray, parent=None, name=Non
     # assumes a rise-then-decay shape and makes the opposite shape unreachable.
     rate_bound = 100.0 / x_range
     x_center = FitParam(
-        _("Center position"), x_center_guess, x_min, x_max, format=DEFAULT_FORMAT
+        _("Center position") + " (x<sub>c</sub>)",
+        x_center_guess,
+        x_min,
+        x_max,
+        format=DEFAULT_FORMAT,
     )
     a_left = FitParam(
-        _("Left amplitude"),
+        _("Left amplitude") + " (a<sub>L</sub>)",
         a_left_guess,
         -amp_bound,
         amp_bound,
         format=DEFAULT_FORMAT,
     )
     b_left = FitParam(
-        _("Left rate") + " (bL)",
+        _("Left rate") + " (b<sub>L</sub>)",
         b_left_guess,  # Already in coefficient form
         -rate_bound,
         rate_bound,
         format=DEFAULT_FORMAT,
     )
     a_right = FitParam(
-        _("Right amplitude"),
+        _("Right amplitude") + " (a<sub>R</sub>)",
         a_right_guess,
         -amp_bound,
         amp_bound,
         format=DEFAULT_FORMAT,
     )
     b_right = FitParam(
-        _("Right rate") + " (bR)",
+        _("Right rate") + " (b<sub>R</sub>)",
         b_right_guess,  # Already in coefficient form
         -rate_bound,
         rate_bound,
         format=DEFAULT_FORMAT,
     )
     y0 = FitParam(
-        _("Base line"),
+        _("Base line") + " (y₀)",
         y0_guess,
         y0_guess - 0.2 * y_range,
         y0_guess + 0.2 * y_range,
@@ -926,6 +1045,7 @@ def fitfunc(x, params):
         parent=parent,
         wintitle=_("Piecewise exponential (raise-decay) fit"),
         name=name,
+        formula=computer.get_formula(),
     )
     if values:
         y_fitted = fitfunc(x, values)