diff --git a/docs/environment_docs.yml b/docs/environment_docs.yml
index e591e8d2..d32eeb2c 100644
--- a/docs/environment_docs.yml
+++ b/docs/environment_docs.yml
@@ -7,6 +7,7 @@ dependencies:
   - xarray
   - sphinx
   - pydata-sphinx-theme
+  - arm-test-data
   - pip:
       - sphinx_gallery
       - sphinx-copybutton
diff --git a/pysp2/testing/sample_files.py b/pysp2/testing/sample_files.py
index 5f5dec23..13b8e421 100644
--- a/pysp2/testing/sample_files.py
+++ b/pysp2/testing/sample_files.py
@@ -2,11 +2,28 @@
 These are the sample files used for testing PySP2.
 """
 import os
+from pathlib import Path
+from arm_test_data import DATASETS
 
-DATA_PATH = os.path.join(os.path.dirname(__file__), 'data')
-EXAMPLE_SP2B = os.path.join(
-    DATA_PATH, 'mosaossp2M1.00.20191216.130601.raw.20191216x193.sp2b')
-EXAMPLE_INI = os.path.join(
-    DATA_PATH, 'mosaossp2M1.00.20191216.000601.raw.20191216000000.ini')
-EXAMPLE_HK = os.path.join(
-    DATA_PATH, 'mosaossp2auxM1.00.20191217.010801.raw.20191216000000.hk')
+DATA_PATH = Path(__file__).resolve().parent / "data"
+
+
+def _sample_file(local_name: str) -> str:
+    local_path = DATA_PATH / local_name
+    if local_path.exists():
+        return str(local_path)
+    return DATASETS.fetch(local_name)
+
+#DATA_PATH = os.path.join(os.path.dirname(__file__), 'data')
+EXAMPLE_SP2B = _sample_file(
+    'mosaossp2M1.00.20191216.130601.raw.20191216x193.sp2b')
+EXAMPLE_INI = _sample_file(
+    'mosaossp2M1.00.20191216.000601.raw.20191216000000.ini')
+EXAMPLE_HK = _sample_file(
+    'mosaossp2auxM1.00.20191217.010801.raw.20191216000000.hk')
+
+print("Fetching files from ARM test data repository...")
+print("DATASETS:", DATASETS.registry_files)
+EXAMPLE_SP2B_PSL = _sample_file("20230721x002.sp2b")
+EXAMPLE_INI_PSL = _sample_file("20230721121710.ini")
+EXAMPLE_HK_PSL = _sample_file("20230721121711.hk")
\ No newline at end of file
diff --git a/pysp2/util/normalized_derivative_method.py b/pysp2/util/normalized_derivative_method.py
index 3cd22573..f539984d 100644
--- a/pysp2/util/normalized_derivative_method.py
+++ b/pysp2/util/normalized_derivative_method.py
@@ -2,6 +2,7 @@
 
 import numpy as np
 import xarray as xr
+import matplotlib.pyplot as plt
 
 from dataclasses import dataclass
 from typing import Optional, Union
@@ -54,7 +55,7 @@ def central_difference(S, num_records=None, normalize=True, baseline_to_zero=Tru
     if num_records is None:
         num_records = S.sizes['event_index']
 
-    dt = 200e-9
+    dt = 0.4
     channels = ['Data_ch0', 'Data_ch4']
     dSdt = {}
 
@@ -120,7 +121,6 @@ def plot_normalized_derivative(S, ds, record_no, chn=0, plot_scattering_signal=F
     ax: matplotlib Axes
         The axes object containing the plot.
     """
-    import matplotlib.pyplot as plt
 
     if chn not in [0, 4]:
         raise ValueError("Channel number must be 0 or 4.")
@@ -131,8 +131,8 @@ def plot_normalized_derivative(S, ds, record_no, chn=0, plot_scattering_signal=F
     inp_data['time'] = xr.DataArray(np.array(time[np.newaxis]),
                                     dims=['time'])
 
-    bins = np.arange(0, 0.00004-0.3e-6, 0.4e-6)  # 0 to 0.0004 microseconds in steps of 0.4e-6 seconds
-    bins = bins*1e6  # convert to microseconds for plotting
+    bins = np.arange(0, 40, 0.4)  # 0 to 39 in steps of 0.4
+    #bins = bins*1e6  # convert to microseconds for plotting
 
     ch_name = f'Data_ch{chn}'
     plt.figure(figsize=(10, 6))
@@ -145,8 +145,9 @@ def plot_normalized_derivative(S, ds, record_no, chn=0, plot_scattering_signal=F
     )
 
     ax.set_xlim([bins[0], bins[-1]])
-    ax.set_xlabel('Time ($\\mu$s)')
-    ax.set_ylabel('Normalized Derivative')
+    ax.set_xlabel(r'Time ($\mu$s)')
+    ax.set_ylim([-1.0,1.0])
+    ax.set_ylabel(r'Normalized Derivative ($\rm \mu s^{-1}$)')
 
     # --- Secondary axis: scattering signal ---
     if plot_scattering_signal:
@@ -189,15 +190,13 @@ def _resolve_peakfit_window(
     """
     Build the MLE fit window from PySP2 peak-fit outputs.
 
-    Left bound:
-        PkStart_ch{peak_ch}
-
-    Right bound:
-        PkStart_ch{peak_ch} + width
+    The window is restricted to the leading edge only:
+        - Left bound: max(min_start, PkStart_chX)
+        - Right bound: min(PkPos_chX, PkStart_chX + width)
 
     where width is either:
-        - PkFWHM_ch{peak_ch}, or
-        - PkFWHM_ch{peak_ch} converted to FWTM if requested
+        - PkFWHM_chX, or
+        - PkFWHM_chX converted to FWTM if requested
 
     Returns
     -------
@@ -205,7 +204,6 @@ def _resolve_peakfit_window(
     """
     # --- Normalize channel input ---
     if isinstance(ch, str):
-        # Expect something like "Data_ch0"
         if "ch" not in ch:
             raise ValueError(f"Cannot infer channel from {ch}")
         ch_num = int(ch.split("ch")[-1])
@@ -213,18 +211,24 @@ def _resolve_peakfit_window(
         ch_num = int(ch)
 
     start_var = f"PkStart_ch{ch_num}"
+    pos_var = f"PkPos_ch{ch_num}"
     width_var = f"PkFWHM_ch{ch_num}"
 
     if start_var not in peak_ds:
         raise ValueError(f"{start_var} not found in dataset")
+    if pos_var not in peak_ds:
+        raise ValueError(f"{pos_var} not found in dataset")
     if width_var not in peak_ds:
         raise ValueError(f"{width_var} not found in dataset")
 
     pk_start = float(peak_ds[start_var].isel({event_dim: event_index}).values)
+    pk_pos = float(peak_ds[pos_var].isel({event_dim: event_index}).values)
     pk_fwhm = float(peak_ds[width_var].isel({event_dim: event_index}).values)
 
     if not np.isfinite(pk_start):
         raise ValueError(f"Invalid peak start: {pk_start}")
+    if not np.isfinite(pk_pos):
+        raise ValueError(f"Invalid peak position: {pk_pos}")
     if not np.isfinite(pk_fwhm) or pk_fwhm <= 0:
         raise ValueError(f"Invalid peak width: {pk_fwhm}")
 
@@ -235,8 +239,17 @@ def _resolve_peakfit_window(
     else:
         raise ValueError("width_metric must be 'fwhm' or 'fwtm'")
 
+    # Left bound from the fitted peak start, but never before min_start.
     fit_start = max(min_start, int(np.floor(pk_start)))
-    fit_stop = int(np.ceil(fit_start + width))
+
+    # Width-based upper bound.
+    width_stop = int(np.ceil(fit_start + width))
+
+    # Leading-edge-only upper bound: stop before the peak maximum.
+    peak_stop = int(np.floor(pk_pos))
+
+    # Use the earlier of the two, so the window stays on the leading edge.
+    fit_stop = min(width_stop, peak_stop)
 
     if n_samples is not None:
         fit_stop = min(fit_stop, n_samples)
@@ -247,6 +260,7 @@ def _resolve_peakfit_window(
     return fit_start, fit_stop, {
         "ch": ch_num,
         "pk_start": pk_start,
+        "pk_pos": pk_pos,
         "pk_fwhm": pk_fwhm,
         "width_metric": width_metric,
     }
@@ -458,7 +472,7 @@ def mle_tau_moteki_kondo(
     # --- Peak-based fit window ---
     # This ensures the MLE only evaluates subsets within the Gaussian peak region.
     fit_start, fit_stop, _ = _resolve_peakfit_window(
-        S_original,
+        peak_ds=S_original,
         event_index=event_index,
         ch=ch,
         event_dim=event_dim,
@@ -667,7 +681,7 @@ def compute_d2_moteki_kondo(
 
    # --- Same peak-based window as MLE ---
     fit_start, fit_stop, _ = _resolve_peakfit_window(
-        S_original,
+        peak_ds=S_original,
         event_index=event_index,
         ch=ch,
         event_dim=event_dim,
@@ -721,4 +735,404 @@ def compute_d2_moteki_kondo(
         coords={"k": k_values},
         name="d2",
         attrs={"fit_start": fit_start, "fit_stop": fit_stop}
-    )
\ No newline at end of file
+    )
+
+def compute_sigma_moteki_kondo(
+    S: Union[xr.DataArray, xr.Dataset],
+    norm_deriv: Union[xr.DataArray, xr.Dataset],
+    tau_hat: Union[np.ndarray, xr.DataArray],
+    d2: Union[np.ndarray, xr.DataArray],
+    p: int = 11,
+    *,
+    ch: Optional[str] = None,
+    event_index: int,
+    event_dim: str = "event_index",
+    S_sample_dim: Optional[str] = None,
+    y_sample_dim: Optional[str] = None,
+    min_start: int = 15,
+    width_metric: str = "fwhm",
+    d2_threshold: float = 80000.0,
+    config: Optional[MLEConfig] = None,
+) -> xr.Dataset:
+    """
+    Estimate Gaussian width sigma using the Moteki & Kondo method.
+
+    1. Use the same peak-based window as the tau/d2 routines.
+       - Left bound: PkStart_chX
+       - Right bound: PkStart_chX + width derived from PkFWHM_chX
+    2. Determine kbest as the subset index that minimizes d²(k).
+    3. Use tau_hat[kbest] as tau_best.
+    4. Fit the linear relation y_i = a (t_i - tau_best) on the kbest sub-array
+       using weighted least squares with weights 1 / (δy_i)_ran².
+    5. Convert slope to width by sigma = sqrt(-1 / a).
+
+    Notes
+    -----
+    The paper applies the sigma estimate after requiring d²(kbest) < 200000.
+    For consistency, this function returns sigma_hat = NaN when the threshold
+    is not met, while still returning diagnostic fields.
+    """
+    if config is None:
+        raise ValueError("config must be provided.")
+
+    # Convert datasets to DataArrays if needed.
+    S_original = S
+    S = _to_dataarray(S, "S", ch=ch)
+    norm_deriv = _to_dataarray(norm_deriv, "norm_deriv", ch=ch)
+
+    # Require event dimension.
+    if event_dim not in S.dims:
+        raise ValueError(f"{event_dim!r} not found in S.dims={S.dims}")
+    if event_dim not in norm_deriv.dims:
+        raise ValueError(f"{event_dim!r} not found in norm_deriv.dims={norm_deriv.dims}")
+
+    # Infer sample dimensions if not provided.
+    if S_sample_dim is None:
+        s_non_event_dims = [d for d in S.dims if d != event_dim]
+        if len(s_non_event_dims) != 1:
+            raise ValueError(
+                f"Could not infer S sample dim. Non-event dims in S: {s_non_event_dims}"
+            )
+        S_sample_dim = s_non_event_dims[0]
+
+    if y_sample_dim is None:
+        y_non_event_dims = [d for d in norm_deriv.dims if d != event_dim]
+        if len(y_non_event_dims) != 1:
+            raise ValueError(
+                f"Could not infer norm_deriv sample dim. Non-event dims in norm_deriv: {y_non_event_dims}"
+            )
+        y_sample_dim = y_non_event_dims[0]
+
+    # Standardize sample dimension name.
+    S_std = S.rename({S_sample_dim: "sample"})
+    y_std = norm_deriv.rename({y_sample_dim: "sample"})
+    S_std, y_std = xr.align(S_std, y_std, join="inner")
+
+    if event_index < 0 or event_index >= S_std.sizes[event_dim]:
+        raise ValueError(
+            f"event_index must be in [0, {S_std.sizes[event_dim] - 1}], got {event_index}"
+        )
+
+    n_samples = S_std.sizes["sample"]
+
+    # Same internal peak-based window used by the tau and d2 routines.
+    fit_start, fit_stop, _ = _resolve_peakfit_window(
+        peak_ds=S_original,
+        event_index=event_index,
+        ch=ch if ch is not None else "Data_ch0",
+        event_dim=event_dim,
+        min_start=min_start,
+        width_metric=width_metric,
+        n_samples=n_samples,
+    )
+
+    k_values = np.arange(fit_start, fit_stop - p + 1)
+    if k_values.size == 0:
+        raise ValueError(
+            f"No valid k subsets for p={p} inside window [{fit_start}, {fit_stop})."
+        )
+
+    tau_hat_np = np.asarray(
+        tau_hat.data if isinstance(tau_hat, xr.DataArray) else tau_hat,
+        dtype=float,
+    )
+    d2_np = np.asarray(
+        d2.data if isinstance(d2, xr.DataArray) else d2,
+        dtype=float,
+    )
+
+    if tau_hat_np.ndim != 1:
+        raise ValueError("tau_hat must be 1D.")
+    if d2_np.ndim != 1:
+        raise ValueError("d2 must be 1D.")
+    if tau_hat_np.size != k_values.size:
+        raise ValueError(
+            f"tau_hat length mismatch. Expected {k_values.size}, got {tau_hat_np.size}."
+        )
+    if d2_np.size != k_values.size:
+        raise ValueError(
+            f"d2 length mismatch. Expected {k_values.size}, got {d2_np.size}."
+        )
+
+    finite = np.isfinite(d2_np) & np.isfinite(tau_hat_np)
+    if not np.any(finite):
+        raise ValueError("No finite d2/tau_hat values available to determine kbest.")
+
+    # Moteki & Kondo Appendix A parameters.
+    h = float(config.h)
+    sigma_bar = float(config.sigma_bar)
+    delta_sigma = float(config.delta_sigma)
+    A1, A2, A3 = float(config.A1), float(config.A2), float(config.A3)
+
+    if h <= 0:
+        raise ValueError("config.h must be positive.")
+    if sigma_bar <= 0:
+        raise ValueError("config.sigma_bar must be positive.")
+    if delta_sigma < 0:
+        raise ValueError("config.delta_sigma must be >= 0.")
+
+    # Time axis in the same units used by the MLE routines.
+    t = np.arange(n_samples, dtype=float) * h
+
+    # Select the requested event.
+    s_event = np.asarray(S_std.sel({event_dim: event_index}).values, dtype=float)
+    y_event = np.asarray(y_std.sel({event_dim: event_index}).values, dtype=float)
+
+    if not (np.all(np.isfinite(s_event)) and np.all(np.isfinite(y_event))):
+        raise ValueError(f"Selected event_index={event_index} contains non-finite values.")
+
+    # kbest is the subset index that minimizes d²(k) [Appendix A.5].
+    kbest_local = int(np.nanargmin(np.where(finite, d2_np, np.nan)))
+    kbest = int(k_values[kbest_local])
+    tau_best = float(tau_hat_np[kbest_local])
+    d2_best = float(d2_np[kbest_local])
+
+    # Paper-consistent acceptance test: only proceed when d²(kbest) is small enough.
+    accepted = bool(np.isfinite(d2_best) and (d2_best < d2_threshold))
+
+    # The kbest sub-array used in Appendix A.6.
+    yk = y_event[kbest : kbest + p]
+    sk = s_event[kbest : kbest + p]
+    tk = t[kbest : kbest + p]
+
+    if not (np.all(np.isfinite(yk)) and np.all(np.isfinite(sk))):
+        raise ValueError("kbest sub-array contains non-finite values.")
+
+    # Random-error model from Appendix A.3:
+    #   (δy_i)_ran = Af_d / h * (1 / S_i) * δS_i
+    # with δS_i = sqrt(A1^2 + A2^2 S_i + A3^2 S_i^2)
+    Af_d = np.sqrt(130.0) / 12.0
+    deltaS = np.sqrt(A1 * A1 + (A2 * A2) * sk + (A3 * A3) * (sk * sk))
+
+    with np.errstate(divide="ignore", invalid="ignore"):
+        var_rand = (Af_d * Af_d) / (h * h) * (deltaS * deltaS) / (sk * sk)
+
+    if not np.all(np.isfinite(var_rand)) or np.any(var_rand <= 0):
+        raise ValueError("Invalid random-variance weights in the kbest sub-array.")
+
+    # Appendix A.6:
+    # Fit y_i = a (t_i - tau_best) with weights 1 / (δy_i)_ran^2.
+    x = tk - tau_best
+    denom = np.sum((x * x) / var_rand)
+    if denom <= 0 or not np.isfinite(denom):
+        raise ValueError("Degenerate weighted least-squares system for sigma estimation.")
+
+    # Weighted slope estimate for the zero-intercept linear model.
+    a = np.sum((x * yk) / var_rand) / denom
+
+    # Appendix A.6: sigma = sqrt(-1 / a)
+    if accepted and np.isfinite(a) and a < 0:
+        sigma_hat = float(np.sqrt(-1.0 / a))
+    else:
+        sigma_hat = np.nan
+
+    out = xr.Dataset(
+        data_vars={
+            "sigma_hat": xr.DataArray(
+                sigma_hat,
+                attrs={
+                    "long_name": "Moteki & Kondo Gaussian width sigma",
+                    "units": "time units of t",
+                },
+            ),
+            "slope": xr.DataArray(
+                float(a),
+                attrs={"long_name": "Weighted least-squares slope a"},
+            ),
+            "tau_best": xr.DataArray(
+                tau_best,
+                attrs={"long_name": "tau associated with kbest"},
+            ),
+            "kbest": xr.DataArray(
+                kbest,
+                attrs={"long_name": "Subset index minimizing d2"},
+            ),
+            "d2_best": xr.DataArray(
+                d2_best,
+                attrs={"long_name": "Minimum d2 value at kbest"},
+            ),
+            "accepted": xr.DataArray(
+                accepted,
+                attrs={"long_name": f"True when d2_best < {d2_threshold}"},
+            ),
+            "fit_start": xr.DataArray(fit_start),
+            "fit_stop": xr.DataArray(fit_stop),
+        },
+        attrs={
+            "method": "Moteki & Kondo weighted least-squares sigma estimate",
+            "width_metric": width_metric,
+            "d2_threshold": d2_threshold,
+        },
+    )
+    return out
+
+def plot_incident_irradiance(
+    S: xr.Dataset,
+    ds: xr.Dataset,
+    record_no: int,
+    chn: int = 0,
+    plot_scattering_signal: bool = True,
+    sigma_ds: Optional[xr.Dataset] = None,
+    tau: Optional[float] = None,
+    sigma: Optional[float] = None,
+    h: float = 0.4,
+    time_units: str = "us",
+    show_fit_window: bool = True,
+):
+    """
+    Plot normalized derivative S'(t)/S(t), expected I'(t)/I(t), and optionally
+    the scattering signal, all against the same bins-based time axis.
+
+    Parameters
+    ----------
+    S : xr.Dataset
+        Original scattering signal dataset.
+    ds : xr.Dataset
+        Dataset containing the normalized derivative.
+    record_no : int
+        Event index to plot.
+    chn : int
+        Channel number (0 or 4).
+    plot_scattering_signal : bool
+        If True, overlay the scattering signal on a secondary y-axis.
+    sigma_ds : xr.Dataset, optional
+        Output of compute_sigma_moteki_kondo(). If provided, tau/sigma are
+        taken from sigma_ds["tau_best"] and sigma_ds["sigma_hat"].
+    tau : float, optional
+        Beam-center time in seconds.
+    sigma : float, optional
+        Gaussian width in seconds.
+    h : float
+        Sampling interval in seconds.
+    time_units : {"us", "s"}
+        Units for the x-axis.
+    show_fit_window : bool
+        If True, shade the fitted leading-edge window when available.
+
+    Returns
+    -------
+    ax : matplotlib Axes
+        Primary axes object.
+    """
+    if chn not in [0, 4]:
+        raise ValueError("Channel number must be 0 or 4.")
+
+    ch_name = f"Data_ch{chn}"
+
+    if ch_name not in S.data_vars:
+        raise ValueError(f"{ch_name!r} not found in S.data_vars.")
+    if ch_name not in ds.data_vars:
+        raise ValueError(f"{ch_name!r} not found in ds.data_vars.")
+
+    spectra = ds.isel(event_index=record_no)
+
+    y_norm = np.asarray(spectra[ch_name].values, dtype=float)
+    y_scatter = np.asarray(S[ch_name].isel(event_index=record_no).values, dtype=float)
+
+    n_samples = y_norm.shape[-1]
+    if y_scatter.shape[-1] != n_samples:
+        raise ValueError(
+            f"Normalized derivative and scattering signal have different lengths: "
+            f"{n_samples} vs {y_scatter.shape[-1]}"
+        )
+
+    # Use the same bins convention everywhere.
+    t = np.arange(n_samples, dtype=float) * h
+    if time_units == "us":
+        t_plot = t
+        tau_scale = 1.0
+        x_label = r"Time ($\rm \mu$s)"
+    elif time_units == "s":
+        t_plot = t * 1e-6
+        tau_scale = 1.0e-6
+        x_label = r"Time (s)"
+    else:
+        raise ValueError("time_units must be 'us' or 's'.")
+
+    # Pull tau and sigma from sigma_ds if supplied.
+    if sigma_ds is not None:
+        if tau is None:
+            tau = float(sigma_ds["tau_best"].item())
+        if sigma is None:
+            sigma = float(sigma_ds["sigma_hat"].item())
+
+    if tau is None or sigma is None:
+        raise ValueError("Provide either sigma_ds or both tau and sigma.")
+
+    if not np.isfinite(tau) or not np.isfinite(sigma) or sigma <= 0:
+        raise ValueError(f"Invalid tau/sigma values: tau={tau}, sigma={sigma}")
+
+    tau_plot = tau * tau_scale
+    sigma_plot = sigma * tau_scale
+
+    # Expected I'/I line from Moteki & Kondo.
+    i_ratio_expected = -(t_plot - tau_plot) / (sigma_plot ** 2)
+
+    fig, ax = plt.subplots(figsize=(10, 6))
+
+    # Normalized derivative.
+    line1, = ax.plot(
+        t_plot,
+        y_norm,  # Scale for visibility
+        label=f"{ch_name} (Normalized dS/dt)",
+        linewidth=1.2,
+    )
+
+    # Expected I'/I.
+    line2, = ax.plot(
+        t_plot,
+        i_ratio_expected,
+        linestyle="--",
+        linewidth=2.0,
+        label=r"Expected $I'(t)/I(t)$",
+    )
+
+    ax.set_xlabel(x_label)
+    ax.set_ylim(-1.0, 1.0)
+    ax.set_xlim(t_plot[10], t_plot[-30])
+    ax.set_ylabel(r"Normalized Derivative ($\rm \mu s^{-1}$)")
+    ax.grid(True, alpha=0.3)
+
+    # Optional fit window shading.
+    if show_fit_window and sigma_ds is not None:
+        fit_start = int(sigma_ds["fit_start"].item())
+        fit_stop = int(sigma_ds["fit_stop"].item())
+
+        fit_start = max(0, min(fit_start, n_samples - 1))
+        fit_stop = max(fit_start + 1, min(fit_stop, n_samples))
+
+        ax.axvspan(
+            t_plot[fit_start],
+            t_plot[fit_stop - 1],
+            color="gray",
+            alpha=0.12,
+            label="Fit window",
+        )
+
+    # Optional scattering signal overlay.
+    if plot_scattering_signal:
+        ax2 = ax.twinx()
+        y_scatter_shifted = y_scatter - np.nanmin(y_scatter)
+
+        line3, = ax2.plot(
+            t_plot,
+            y_scatter_shifted,
+            color="black",
+            linestyle="--",
+            linewidth=1.2,
+            label=f"{ch_name} (Scattering Signal)",
+        )
+        ax2.set_ylabel("Scattering Signal (baseline shifted)")
+
+        lines = [line1,line2, line3]
+        labels = [l.get_label() for l in lines]
+        ax.legend(lines, labels, loc="best")
+    else:
+        ax.legend(loc="best")
+
+    ax.set_title(
+        f"Normalized Derivative, Expected I'(t)/I(t), and Scattering Signal - "
+        f"Channel {chn} Record {record_no}"
+    )
+
+    return ax
\ No newline at end of file
diff --git a/tests/baseline/test_plot_incident_irradiance.png b/tests/baseline/test_plot_incident_irradiance.png
new file mode 100644
index 00000000..a07f445b
Binary files /dev/null and b/tests/baseline/test_plot_incident_irradiance.png differ
diff --git a/tests/baseline/test_plot_normalized_derivative.png b/tests/baseline/test_plot_normalized_derivative.png
index 9a8846c9..6faac7ef 100644
Binary files a/tests/baseline/test_plot_normalized_derivative.png and b/tests/baseline/test_plot_normalized_derivative.png differ
diff --git a/tests/baseline/test_plot_wave.png b/tests/baseline/test_plot_wave.png
index 941839ac..5eddb820 100644
Binary files a/tests/baseline/test_plot_wave.png and b/tests/baseline/test_plot_wave.png differ
diff --git a/tests/test_ndm.py b/tests/test_ndm.py
index 34014c09..00c830f3 100644
--- a/tests/test_ndm.py
+++ b/tests/test_ndm.py
@@ -3,41 +3,40 @@
 np.set_printoptions(threshold=np.inf)
 
 from pysp2.util.normalized_derivative_method import MLEConfig, mle_tau_moteki_kondo, compute_d2_moteki_kondo
+from pysp2.util.normalized_derivative_method import compute_sigma_moteki_kondo
+event=152
+my_sp2b = pysp2.io.read_sp2(pysp2.testing.EXAMPLE_SP2B_PSL, arm_convention=False)
+my_ini = pysp2.io.read_config(pysp2.testing.EXAMPLE_INI_PSL)
+my_binary = pysp2.util.gaussian_fit(my_sp2b, my_ini, parallel=False)
 
 def test_central_difference():
-    my_sp2b = pysp2.io.read_sp2(pysp2.testing.EXAMPLE_SP2B)
-    my_ini = pysp2.io.read_config(pysp2.testing.EXAMPLE_INI)
-    my_binary = pysp2.util.gaussian_fit(my_sp2b, my_ini, parallel=False)
     dSdt = pysp2.util.central_difference(my_binary, normalize=False, baseline_to_zero=False)
     
-    np.testing.assert_almost_equal(dSdt['Data_ch4'].isel(event_index=5876, time=0).item(), 
-                                   8.3333333333e6, decimal=2)
-    np.testing.assert_almost_equal(dSdt['Data_ch4'].isel(event_index=5876, time=99).item(), 
-                                   7.166666666e7, decimal=2)
-    np.testing.assert_almost_equal(dSdt['Data_ch4'].isel(event_index=5876, time=19).item(), 
-                                   1.5e7, decimal=2)
+    np.testing.assert_almost_equal(dSdt['Data_ch4'].isel(event_index=event, time=0).item(), 
+                                   2.50, decimal=2)
+    np.testing.assert_almost_equal(dSdt['Data_ch4'].isel(event_index=event, time=99).item(), 
+                                   13.33, decimal=2)
+    np.testing.assert_almost_equal(dSdt['Data_ch4'].isel(event_index=event, time=2).item(), 
+                                    0.83, decimal=2)
     assert np.isfinite(dSdt).all()
 
-    dSdt_norm = pysp2.util.central_difference(my_binary, normalize=True, baseline_to_zero=False)
-    np.testing.assert_almost_equal(dSdt_norm['Data_ch4'].isel(event_index=5876, time=0).item(), 
-                                   8.3333333333e6/-30168, decimal=2)
-    np.testing.assert_almost_equal(dSdt_norm['Data_ch4'].isel(event_index=5876, time=99).item(), 
-                                   7.166666666e7/-30152, decimal=2)
-    np.testing.assert_almost_equal(dSdt_norm['Data_ch4'].isel(event_index=5876, time=19).item(), 
-                                   1.5e7/-30132, decimal=2)
+    dSdt_norm = pysp2.util.central_difference(my_binary, normalize=True, baseline_to_zero=True)
+    np.testing.assert_almost_equal(dSdt_norm['Data_ch0'].isel(event_index=event, time=0).item(), 
+                                   -0.70, decimal=2)
+    np.testing.assert_almost_equal(dSdt_norm['Data_ch0'].isel(event_index=event, time=99).item(), 
+                                   0.799, decimal=2)
+    np.testing.assert_almost_equal(dSdt_norm['Data_ch0'].isel(event_index=event, time=2).item(), 
+                                    -0.078, decimal=2)
     
-def test_mle_estimate_tau():
-    my_sp2b = pysp2.io.read_sp2(pysp2.testing.EXAMPLE_SP2B)
-    my_ini = pysp2.io.read_config(pysp2.testing.EXAMPLE_INI)
-    my_binary = pysp2.util.gaussian_fit(my_sp2b, my_ini, parallel=False, baseline_to_zero=True)
-    dSdt = pysp2.util.central_difference(my_binary, normalize=True, baseline_to_zero=True)
+def test_ndm_moteki_kondo():
+    dSdt = pysp2.util.central_difference(my_binary, normalize=True, baseline_to_zero=True)  
 
     cfg = MLEConfig(
-    h=0.4e-6,           # example: 0.4 microseconds
-    sigma_bar= 16.6*0.4e-6,  # example; use your measured average width
-    delta_sigma=1.2*0.4e-6, # example; use your measured width std dev
-    A1=0.37,
-    A2=1.6e-2,
+    h=0.4,           # example: 0.4 microseconds
+    sigma_bar= 18.5*0.4,  # example; use your measured average width
+    delta_sigma=1.2*0.4, # example; use your measured width std dev
+    A1=0.37*2.44,
+    A2=(1.6e-2)*2.44**(1/2),
     A3=6.2e-4,
 )
     
@@ -45,28 +44,28 @@ def test_mle_estimate_tau():
     tau = mle_tau_moteki_kondo(
         S=my_binary,
         norm_deriv=dSdt,
-        p=10,
+        p=13,
         ch="Data_ch0",
-        event_index=499,
+        event_index=event,
         min_start=15,
-        width_metric="fwhm",
+        width_metric="fwtm",
         config=cfg,
     )
 
-    tau_val_true = my_binary['Data_ch0'].isel(event_index=499).argmax().item()*0.4e-6
+    tau_val_true = my_binary['Data_ch0'].isel(event_index=event).argmax().item()*0.4
     # Test that the estimated tau for a subset of results is close to the true value for the event
-    for i in range(13, 17):
-        np.testing.assert_almost_equal(tau[i], tau_val_true, decimal=6)
+    for i in range(10, 15):
+        np.testing.assert_allclose(tau[i], tau_val_true, atol=0.3)
     
     d2 = compute_d2_moteki_kondo(
         S=my_binary,
         norm_deriv=dSdt,
         tau_hat=tau,
-        p=10,
+        p=13,
         ch="Data_ch0",
-        event_index=499,
+        event_index=event,
         min_start=15,
-        width_metric="fwhm",
+        width_metric="fwtm",
         config=cfg,
     )
 
@@ -78,86 +77,28 @@ def test_mle_estimate_tau():
     np.testing.assert_allclose(
         tau_best,
         tau_val_true,
-        atol=0.01e-05,  # absolute tolerance = 1e-7
-    )
-
-    ## Test another event ##################################################
-    tau = mle_tau_moteki_kondo(
-        S=my_binary,
-        norm_deriv=dSdt,
-        p=6,
-        ch="Data_ch0",
-        event_index=1040,
-        min_start=15,
-        width_metric="fwtm",
-        config=cfg,
+        atol=0.3,  # absolute tolerance = 5e-7
     )
 
-    tau_val = my_binary['Data_ch0'].isel(event_index=1040).argmax().item()*0.4e-6
-    # Test that the estimated tau for a subset of results is close to the true value for the event
-    for i in range(18, 28):
-        np.testing.assert_allclose(tau[i], tau_val, atol=0.08e-05)
-
-    d2 = compute_d2_moteki_kondo(
+    sigma_ds = compute_sigma_moteki_kondo(
         S=my_binary,
         norm_deriv=dSdt,
         tau_hat=tau,
-        p=6,
+        d2=d2,
+        p=13,
         ch="Data_ch0",
-        event_index=1040,
+        event_index=event,
         min_start=15,
         width_metric="fwtm",
         config=cfg,
     )
 
-    k_min_local = int(d2.argmin(dim="k").item())
-    # Get corresponding tau value
-    tau_best = tau.isel(k=k_min_local).item()
+    # example: use the best sigma value from your analysis, divided by 4.29193 to convert FWTM value of
+    # 18.51*np.sqrt(np.log(10)/np.log(2)) to sigma where 18.51 is the average FWHM in us
+    sigma_best = (33.7366*0.4)/4.29193
 
-    # Assert closeness
     np.testing.assert_allclose(
-        tau_best,
-        tau_val,
-        atol=0.04e-05,  # absolute tolerance = 4e-7
-    )
-    
-    ## Test another event ##################################################
-    tau = mle_tau_moteki_kondo(
-        S=my_binary,
-        norm_deriv=dSdt,
-        p=8,
-        ch="Data_ch4",
-        event_index=2008,
-        min_start=15,
-        width_metric="fwtm",
-        config=cfg,
-    )
-
-    d2 = compute_d2_moteki_kondo(
-        S=my_binary,
-        norm_deriv=dSdt,
-        tau_hat=tau,
-        p=8,
-        ch="Data_ch4",
-        event_index=2008,
-        min_start=15,
-        width_metric="fwtm",
-        config=cfg,
-    )
-
-    # Test that the estimated tau for a subset of results is close to the true value for the event
-    for i in range(2, 4):
-        np.testing.assert_allclose(tau[i], tau_val_true, atol=0.08e-05)
-    for i in range(10,15):
-        np.testing.assert_allclose(tau[i], tau_val_true, atol=0.05e-05)
-
-    k_min_local = int(d2.argmin(dim="k").item())
-    # Get corresponding tau value
-    tau_best = tau.isel(k=k_min_local).item()
- 
-    # Assert closeness
-    np.testing.assert_allclose(
-        tau_best,
-        tau_val_true,
-        atol=0.01e-05,  # absolute tolerance = 2e-6 (larger tolerance for evaporation events)
-    )
\ No newline at end of file
+         sigma_ds['sigma_hat'].values,
+         sigma_best,
+         atol=0.12,  # absolute tolerance = 1.5 microseconds
+     )
\ No newline at end of file
diff --git a/tests/test_vis.py b/tests/test_vis.py
index a858ace1..9c1c1805 100644
--- a/tests/test_vis.py
+++ b/tests/test_vis.py
@@ -2,23 +2,28 @@
 import pytest
 import xarray as xr
 import numpy as np
+import matplotlib.pyplot as plt
 
 import pysp2
 from pysp2.util.normalized_derivative_method import plot_normalized_derivative
+from pysp2.util.normalized_derivative_method import MLEConfig, mle_tau_moteki_kondo, compute_d2_moteki_kondo
+from pysp2.util.normalized_derivative_method import compute_sigma_moteki_kondo
+from pysp2.util.normalized_derivative_method import plot_incident_irradiance
 from pysp2.vis.plot_wave import plot_wave
 
 matplotlib.use("Agg")
 
+my_sp2b = pysp2.io.read_sp2(pysp2.testing.EXAMPLE_SP2B_PSL, arm_convention=False)
+my_ini = pysp2.io.read_config(pysp2.testing.EXAMPLE_INI_PSL)
+event = 152
+
 @pytest.mark.mpl_image_compare(tolerance=10)
 def test_plot_normalized_derivative():
-
-    my_sp2b = pysp2.io.read_sp2(pysp2.testing.EXAMPLE_SP2B)
-    my_ini = pysp2.io.read_config(pysp2.testing.EXAMPLE_INI)
     my_binary = pysp2.util.gaussian_fit(my_sp2b, my_ini, parallel=False, baseline_to_zero=True)
     dSdt_norm = pysp2.util.central_difference(my_binary, normalize=True, baseline_to_zero=True)
 
-    # Test the plotting function for channel 0 and record number 499
-    ax = plot_normalized_derivative(my_sp2b, dSdt_norm, record_no=499, chn=0, plot_scattering_signal=True)
+    # Test the plotting function for channel 0 and record number 720
+    ax = plot_normalized_derivative(my_binary, dSdt_norm, record_no=event, chn=0, plot_scattering_signal=True)
     fig = ax.figure
     
     return fig
@@ -26,11 +31,75 @@ def test_plot_normalized_derivative():
 @pytest.mark.mpl_image_compare(tolerance=10)
 def test_plot_wave():
 
-    my_sp2b = pysp2.io.read_sp2(pysp2.testing.EXAMPLE_SP2B)
-    my_ini = pysp2.io.read_config(pysp2.testing.EXAMPLE_INI)
-    my_binary = pysp2.util.gaussian_fit(my_sp2b, my_ini, parallel=False, baseline_to_zero=True)
+    my_binary = pysp2.util.gaussian_fit(my_sp2b, my_ini, parallel=False, baseline_to_zero=False)
 
-    # Test the plotting function for channel 0 and record number 499
-    display = plot_wave(my_binary, record_no=499, chn=0)
+    # Test the plotting function for channel 0 and record number 720
+    display = plot_wave(my_binary, record_no=event, chn=0, plot_fit=True)
     fig = display.axes[0].figure
+    return fig
+
+@pytest.mark.mpl_image_compare(tolerance=10)
+def test_plot_incident_irradiance():
+
+    my_binary = pysp2.util.gaussian_fit(my_sp2b, my_ini, parallel=False, baseline_to_zero=True)
+    dSdt = pysp2.util.central_difference(my_binary, normalize=True, baseline_to_zero=True)
+    print("dSdt dimensions:", dSdt.dims)
+
+    cfg = MLEConfig(
+    h=0.4,           # example: 0.4 microseconds
+    sigma_bar= 18.5*0.4,  # example; use your measured average width
+    delta_sigma=1.2*0.4, # example; use your measured width std dev
+    A1=0.37*2.44,
+    A2=(1.6e-2)*2.44**(1/2),
+    A3=6.2e-4,
+)
+    
+    tau = mle_tau_moteki_kondo(
+        S=my_binary,
+        norm_deriv=dSdt,
+        p=13,
+        ch="Data_ch0",
+        event_index=event,
+        min_start=15,
+        width_metric="fwtm",
+        config=cfg,
+    )
+
+    d2 = compute_d2_moteki_kondo(
+        S=my_binary,
+        norm_deriv=dSdt,
+        tau_hat=tau,
+        p=13,
+        ch="Data_ch0",
+        event_index=event,
+        min_start=15,
+        width_metric="fwtm",
+        config=cfg,
+    )
+
+    sigma_ds = compute_sigma_moteki_kondo(
+        S=my_binary,
+        norm_deriv=dSdt,
+        tau_hat=tau,
+        d2=d2,
+        p=13,
+        ch="Data_ch0",
+        event_index=event,
+        min_start=15,
+        width_metric="fwtm",
+        config=cfg,
+    )
+
+    # Test the plotting function for channel 0 and record number 720
+    ax = plot_incident_irradiance(
+        S=my_binary,
+        ds=dSdt,
+        record_no=event,
+        chn=0,
+        plot_scattering_signal=True,
+        sigma_ds=sigma_ds,
+        time_units="us",
+    )
+    fig = ax.figure
+    
     return fig
\ No newline at end of file