[education/labplot] lib/examples: [scripting] added four new extensive examples to show-case the available functionality.
Alexander Semke <[email protected]>
| Newsgroups | gmane.comp.kde.cvs |
|---|---|
| Message-ID | <[email protected]> |
Git commit 60e027dd56725a586537f6ec5e45a66bae877f9a by Alexander Semke.
Committed on 07/08/2026 at 16:04.
Pushed by asemke into branch 'master'.
[scripting] added four new extensive examples to show-case the available functionality.
M +4 -0 lib/examples/CMakeLists.txt
A +3 -0 lib/examples/Peak Detection/CMakeLists.txt
A +325 -0 lib/examples/Peak Detection/README.md
A +10 -0 lib/examples/Peak Detection/main.py
A +298 -0 lib/examples/Peak Detection/script.py
A +3 -0 lib/examples/Publication Multi-Panel Figure/CMakeLists.txt
A +149 -0 lib/examples/Publication Multi-Panel Figure/README.md
A +10 -0 lib/examples/Publication Multi-Panel Figure/main.py
A +322 -0 lib/examples/Publication Multi-Panel Figure/script.py
A +3 -0 lib/examples/Quantum Wave Packet/CMakeLists.txt
A +109 -0 lib/examples/Quantum Wave Packet/README.md
A +203 -0 lib/examples/Quantum Wave Packet/main.cpp *
A +23 -0 lib/examples/Quantum Wave Packet/main.py
A +235 -0 lib/examples/Quantum Wave Packet/script.py
A +293 -0 lib/examples/Quantum Wave Packet/script_tunneling.py
A +3 -0 lib/examples/Statistical Comparison/CMakeLists.txt
A +220 -0 lib/examples/Statistical Comparison/README.md
A +10 -0 lib/examples/Statistical Comparison/main.py
A +285 -0 lib/examples/Statistical Comparison/script.py
The files marked with a * at the end have a non valid license. Please read: https://community.kde.org/Policies/Licensing_Policy and use the headers which are listed at that page.
https://invent.kde.org/education/labplot/-/commit/60e027dd56725a586537f6ec5e45a66bae877f9a
diff --git a/lib/examples/CMakeLists.txt b/lib/examples/CMakeLists.txt
index 0b4c5bd8df..97fe167f65 100644
--- a/lib/examples/CMakeLists.txt
+++ b/lib/examples/CMakeLists.txt
@@ -30,4 +30,8 @@ add_subdirectory("Data = Smooth + Rough")
add_subdirectory("SOS Morse Signal")
add_subdirectory("Space Debris")
add_subdirectory("Same Data Different Boxplots")
+add_subdirectory("Quantum Wave Packet")
+add_subdirectory("Publication Multi-Panel Figure")
+add_subdirectory("Statistical Comparison")
+add_subdirectory("Peak Detection")
add_subdirectory("Demo")
diff --git a/lib/examples/Peak Detection/CMakeLists.txt b/lib/examples/Peak Detection/CMakeLists.txt
new file mode 100644
index 0000000000..bb26378dff
--- /dev/null
+++ b/lib/examples/Peak Detection/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_executable(PeakDetection main.cpp)
+
+target_link_libraries(PeakDetection LabPlot::SDK)
diff --git a/lib/examples/Peak Detection/README.md b/lib/examples/Peak Detection/README.md
new file mode 100644
index 0000000000..a48c2f7cf3
--- /dev/null
+++ b/lib/examples/Peak Detection/README.md
@@ -0,0 +1,325 @@
+# Automated Peak Detection and Annotation
+
+This tutorial demonstrates automated peak finding and labeling in spectroscopic or chromatographic data - one of the most requested features in analytical chemistry software.
+
+## What You'll Learn
+
+- Automatic peak detection using scipy.signal
+- Peak annotation with TextLabels
+- Reference lines at peak positions
+- Customizable detection parameters
+- Exporting peak lists for further analysis
+
+## The Problem
+
+**Manual peak identification is tedious and error-prone!**
+
+In analytical chemistry and materials science, you often need to:
+1. Identify all peaks in a spectrum or chromatogram
+2. Measure peak positions, heights, and widths
+3. Label peaks for presentation/publication
+4. Export peak data for database comparison
+
+Doing this manually for dozens of spectra is time-consuming. **This script automates the entire workflow.**
+
+## What It Does
+
+### 1. Automatic Peak Detection
+Uses `scipy.signal.find_peaks()` with configurable parameters:
+- **Prominence**: Minimum height above baseline (filters noise)
+- **Distance**: Minimum separation between peaks (merges shoulders)
+- **Height**: Absolute minimum peak intensity
+
+### 2. Peak Information Extraction
+For each detected peak:
+- Position (x-coordinate)
+- Intensity (y-coordinate)
+- Prominence (height above surrounding baseline)
+- Width (FWHM - full width at half maximum)
+
+### 3. Automatic Annotation
+- Marks peaks with symbols
+- Adds vertical reference lines
+- Labels top N peaks with positions
+- Creates exportable peak table
+
+## Example: XRD Pattern Analysis
+
+The demo analyzes an X-ray diffraction pattern showing:
+- **Input**: 2θ angle vs intensity data
+- **Output**: Annotated pattern with labeled peaks
+- **Use**: Phase identification, crystallite size, lattice parameters
+
+## Use Cases by Field
+
+### Spectroscopy
+```python
+# Raman Spectroscopy
+x = wavenumber # cm⁻¹
+y = intensity
+# Find characteristic peaks for compound identification
+
+# IR Spectroscopy
+x = wavenumber # cm⁻¹
+y = transmittance # %
+# Identify functional groups
+
+# UV-Vis
+x = wavelength # nm
+y = absorbance
+# Find absorption maxima
+```
+
+### Chromatography
+```python
+# HPLC/GC
+x = retention_time # minutes
+y = absorbance or intensity
+# Quantify compounds, calculate areas
+
+# Mass Spectrometry
+x = mz_ratio # m/z
+y = abundance # %
+# Identify molecular fragments
+```
+
+### Materials Science
+```python
+# XRD (this example)
+x = two_theta # degrees
+y = intensity # counts
+# Phase identification, crystallinity
+
+# XPS
+x = binding_energy # eV
+y = counts
+# Element identification
+```
+
+### Signal Processing
+```python
+# ECG/EEG
+x = time # seconds
+y = voltage # mV
+# Detect R-peaks, spikes
+
+# Audio
+x = frequency # Hz
+y = amplitude
+# Harmonic analysis
+```
+
+## Customizing Peak Detection
+
+### Finding More Peaks (Lower Threshold)
+```python
+prominence_threshold = 50 # Lower value = more peaks
+```
+
+### Finding Fewer Peaks (Higher Threshold)
+```python
+prominence_threshold = 500 # Higher value = only major peaks
+```
+
+### Separate Close Peaks
+```python
+distance_threshold = 10 # Smaller = distinguish close peaks
+```
+
+### Merge Shoulders
+```python
+distance_threshold = 100 # Larger = treat shoulders as one peak
+```
+
+### Set Absolute Minimum Height
+```python
+peaks = find_peaks_with_info(
+ x, y,
+ prominence=200,
+ height=100 # Ignore peaks below this intensity
+)
+```
+
+## Advanced Techniques
+
+### Baseline Correction First
+```python
+from scipy.signal import savgol_filter
+
+# Smooth data
+y_smooth = savgol_filter(y, window_length=51, polyorder=3)
+
+# Estimate baseline (rolling minimum)
+from scipy.ndimage import minimum_filter
+baseline = minimum_filter(y_smooth, size=100)
+
+# Subtract baseline
+y_corrected = y - baseline
+
+# Then find peaks in corrected data
+peaks = find_peaks_with_info(x, y_corrected, ...)
+```
+
+### Peak Area Integration
+```python
+from scipy.integrate import simpson
+
+for peak in peaks:
+ # Define integration region (±3σ around peak)
+ left_idx = peak['index'] - int(3 * peak['width'] / (x[1] - x[0]))
+ right_idx = peak['index'] + int(3 * peak['width'] / (x[1] - x[0]))
+
+ # Integrate
+ x_region = x[left_idx:right_idx]
+ y_region = y[left_idx:right_idx]
+ area = simpson(y_region, x_region)
+
+ peak['area'] = area
+```
+
+### Peak Fitting (Gaussian, Lorentzian)
+```python
+from scipy.optimize import curve_fit
+
+def gaussian(x, amp, center, width):
+ return amp * np.exp(-(x - center)**2 / (2 * width**2))
+
+for peak in peaks:
+ # Extract region around peak
+ window = 50
+ idx = peak['index']
+ x_fit = x[idx-window:idx+window]
+ y_fit = y[idx-window:idx+window]
+
+ # Fit Gaussian
+ params, _ = curve_fit(
+ gaussian, x_fit, y_fit,
+ p0=[peak['height'], peak['x_pos'], peak['width']]
+ )
+
+ peak['fitted_center'] = params[1]
+ peak['fitted_width'] = params[2]
+```
+
+### Multi-Component Deconvolution
+```python
+# For overlapping peaks, fit sum of Gaussians
+def multi_gaussian(x, *params):
+ # params = [amp1, center1, width1, amp2, center2, width2, ...]
+ y = np.zeros_like(x)
+ for i in range(0, len(params), 3):
+ y += gaussian(x, params[i], params[i+1], params[i+2])
+ return y
+
+# Initial guess from detected peaks
+p0 = []
+for peak in peaks:
+ p0.extend([peak['height'], peak['x_pos'], peak['width']])
+
+# Fit all peaks simultaneously
+params_fitted, _ = curve_fit(multi_gaussian, x, y, p0=p0)
+```
+
+## Exporting Peak Data
+
+### To Spreadsheet/CSV
+```python
+# The peaks_spreadsheet contains all detected peaks
+# Export via LabPlot's UI or:
+# Save as CSV for database comparison, literature matching
+```
+
+### To Text Report
+```python
+with open('peak_report.txt', 'w') as f:
+ f.write("Peak Detection Report\n")
+ f.write("=" * 50 + "\n\n")
+ for i, peak in enumerate(peaks, 1):
+ f.write(f"Peak {i}:\n")
+ f.write(f" Position: {peak['x_pos']:.2f}\n")
+ f.write(f" Intensity: {peak['y_pos']:.1f}\n")
+ f.write(f" Prominence: {peak['prominence']:.1f}\n")
+ f.write(f" Width: {peak['width']:.3f}\n\n")
+```
+
+## Comparison with Competitors
+
+**OriginLab** - Peak Analyzer tool (commercial, $1000+)
+**PeakFit** - Dedicated peak fitting software (commercial, $500+)
+**Igor Pro** - Peak finding procedures (commercial, $700+)
+
+**LabPlot + This Script** - FREE and open source! 🎉
+
+## Validation and Quality Control
+
+### Check for False Positives
+```python
+# Visual inspection: do marked peaks look real?
+# Compare with known reference patterns
+# Check peak widths (too narrow = noise spike)
+```
+
+### Check for Missed Peaks
+```python
+# Lower prominence threshold gradually
+# Visually inspect spectrum for unmarked features
+```
+
+### Reproducibility
+```python
+# Run detection multiple times with same parameters
+# Should give identical results (deterministic algorithm)
+```
+
+## Best Practices
+
+### Data Preprocessing
+1. **Remove spikes** (cosmic rays in Raman, electrical noise)
+2. **Smooth if noisy** (Savitzky-Golay filter)
+3. **Correct baseline** (polynomial, rolling minimum)
+4. **Normalize** (to internal standard or max peak)
+
+### Parameter Selection
+1. **Start conservative** (high prominence, large distance)
+2. **Lower thresholds gradually** until false positives appear
+3. **Document final parameters** for reproducibility
+4. **Use same parameters** across sample series
+
+### Annotation Strategy
+1. **Label major peaks** only (top 5-10 by prominence)
+2. **Use consistent formatting** (font, color, position)
+3. **Include units** (degrees, cm⁻¹, ppm)
+4. **Reference literature** values if identifying phases
+
+## Troubleshooting
+
+### "Too many peaks detected (noise)"
+→ Increase `prominence_threshold`
+→ Smooth data first (`savgol_filter`)
+→ Apply baseline correction
+
+### "Missing real peaks"
+→ Decrease `prominence_threshold`
+→ Decrease `distance_threshold`
+→ Check for baseline drift
+
+### "Peaks too close together"
+→ Increase `distance_threshold`
+→ Or use peak deconvolution for genuine overlap
+
+### "Labels overlap"
+→ Reduce `n_labels` (only top N peaks)
+→ Adjust label positions manually
+→ Use leader lines/arrows
+
+## Related Examples
+
+- See `../NIST - Linear Regression/script.py` for calibration curves from peak areas
+- See `../Publication Multi-Panel Figure/script.py` for including annotated spectra in figures
+
+## Further Reading
+
+- scipy.signal.find_peaks documentation
+- "Automated Spectral Peak Finding" - Applied Spectroscopy reviews
+- NIST Chemistry WebBook for reference peak positions
+- International Centre for Diffraction Data (ICDD) for XRD patterns
diff --git a/lib/examples/Peak Detection/main.py b/lib/examples/Peak Detection/main.py
new file mode 100644
index 0000000000..a78c46afdb
--- /dev/null
+++ b/lib/examples/Peak Detection/main.py
@@ -0,0 +1,10 @@
+#!/usr/bin/env python3
+"""Runner for peak detection demo"""
+import sys
+import os
+
+script_dir = os.path.dirname(os.path.abspath(__file__))
+script_path = os.path.join(script_dir, 'script.py')
+
+with open(script_path, 'r') as f:
+ exec(f.read())
diff --git a/lib/examples/Peak Detection/script.py b/lib/examples/Peak Detection/script.py
new file mode 100644
index 0000000000..e4df8bbe94
--- /dev/null
+++ b/lib/examples/Peak Detection/script.py
@@ -0,0 +1,298 @@
+"""
+Automated Peak Detection and Annotation
+
+This script demonstrates automated peak finding and labeling in spectroscopic
+or chromatographic data - a common task in analytical chemistry, materials
+science, and signal processing.
+
+Features:
+- Automatic peak detection with prominence threshold
+- Peak annotation with TextLabels and arrows
+- Customizable detection parameters
+- Works with any signal: XRD, Raman, IR, HPLC, mass spec, etc.
+
+Use cases:
+- X-ray diffraction (XRD) pattern analysis
+- Raman/IR spectroscopy peak identification
+- Chromatography (HPLC, GC) peak labeling
+- Mass spectrometry peak annotation
+- Any 1D signal with peaks!
+"""
+
+import numpy as np
+from scipy import signal
+from PySide6.QtCore import QRectF, QPointF, Qt
+from PySide6.QtGui import QFont
+from PySide6.QtWidgets import QTextEdit
+from pylabplot import *
+
+# === Peak Detection Function ===
+def find_peaks_with_info(x, y, prominence=None, distance=None, height=None):
+ """
+ Find peaks in signal using scipy.signal.find_peaks
+
+ Parameters:
+ -----------
+ x : array
+ X-axis values (e.g., wavelength, retention time, 2-theta)
+ y : array
+ Y-axis values (intensity, absorbance, counts)
+ prominence : float
+ Minimum prominence (height above surrounding baseline)
+ distance : int
+ Minimum distance between peaks (in samples)
+ height : float
+ Minimum peak height
+
+ Returns:
+ --------
+ peak_info : list of dict
+ Each dict contains: index, x_pos, y_pos, prominence, width
+ """
+ # Find peaks
+ if prominence is None:
+ prominence = 0.1 * (np.max(y) - np.min(y))
+
+ peaks, properties = signal.find_peaks(
+ y,
+ prominence=prominence,
+ distance=distance,
+ height=height,
+ width=1 # Also calculate peak widths
+ )
+
+ # Extract peak information
+ peak_info = []
+ for i, peak_idx in enumerate(peaks):
+ info = {
+ 'index': peak_idx,
+ 'x_pos': x[peak_idx],
+ 'y_pos': y[peak_idx],
+ 'prominence': properties['prominences'][i],
+ 'width': properties['widths'][i] * (x[1] - x[0]) if len(x) > 1 else 0,
+ 'height': y[peak_idx]
+ }
+ peak_info.append(info)
+
+ return peak_info
+
+# === Generate Sample XRD Pattern ===
+# Simulated X-ray diffraction pattern with multiple peaks
+np.random.seed(42)
+
+# X-axis: 2-theta angle (degrees)
+two_theta = np.linspace(10, 80, 2000)
+
+# Background (polynomial baseline)
+baseline = 50 + 0.1 * two_theta + 0.001 * two_theta**2
+
+# Add multiple Gaussian peaks at different positions
+peak_positions = [22.5, 31.8, 38.4, 45.6, 48.2, 56.8, 62.3, 68.5, 74.2]
+peak_heights = [800, 1200, 600, 950, 400, 700, 500, 350, 450]
+peak_widths = [0.8, 0.6, 0.7, 0.5, 0.6, 0.7, 0.8, 0.6, 0.7]
+
+intensity = baseline.copy()
+for pos, height, width in zip(peak_positions, peak_heights, peak_widths):
+ intensity += height * np.exp(-((two_theta - pos) / width)**2)
+
+# Add noise
+intensity += np.random.normal(0, 20, len(intensity))
+
+# Ensure non-negative
+intensity = np.maximum(intensity, 0)
+
+# === Detect Peaks ===
+# Adjust prominence threshold to find significant peaks only
+prominence_threshold = 200 # Minimum peak prominence
+distance_threshold = 50 # Minimum distance between peaks (in samples)
+
+peaks = find_peaks_with_info(
+ two_theta,
+ intensity,
+ prominence=prominence_threshold,
+ distance=distance_threshold
+)
+
+print("=" * 70)
+print("AUTOMATED PEAK DETECTION")
+print("=" * 70)
+print(f"\nFound {len(peaks)} peaks in the spectrum\n")
+print(f"{'Peak #':<8} {'Position':<12} {'Intensity':<12} {'Prominence':<12} {'Width':<10}")
+print("-" * 70)
+for i, peak in enumerate(peaks, 1):
+ print(f"{i:<8} {peak['x_pos']:>10.2f}° {peak['y_pos']:>11.1f} {peak['prominence']:>11.1f} {peak['width']:>9.3f}°")
+print()
+
+# === Setup Project ===
+proj = project()
+
+# Create spreadsheet
+spreadsheet = Spreadsheet("XRD Pattern Data")
+spreadsheet.setColumnCount(2) # x and y data
+proj.addChild(spreadsheet)
+
+# Fill data
+spreadsheet.column(0).setName("2θ (degrees)")
+spreadsheet.column(1).setName("Intensity (counts)")
+spreadsheet.column(0).replaceValues(0, [float(x) for x in two_theta])
+spreadsheet.column(1).replaceValues(0, [float(x) for x in intensity])
+
+# Create peak positions spreadsheet
+peaks_spreadsheet = Spreadsheet("Detected Peaks")
+peaks_spreadsheet.setColumnCount(3) # position, intensity, prominence
+proj.addChild(peaks_spreadsheet)
+
+peaks_spreadsheet.column(0).setName("2θ (degrees)")
+peaks_spreadsheet.column(1).setName("Intensity")
+peaks_spreadsheet.column(2).setName("Prominence")
+
+peak_x = [peak['x_pos'] for peak in peaks]
+peak_y = [peak['y_pos'] for peak in peaks]
+peak_prom = [peak['prominence'] for peak in peaks]
+
+peaks_spreadsheet.column(0).replaceValues(0, [float(x) for x in peak_x])
+peaks_spreadsheet.column(1).replaceValues(0, [float(x) for x in peak_y])
+peaks_spreadsheet.column(2).replaceValues(0, [float(x) for x in peak_prom])
+
+# === Create Worksheet ===
+worksheet = Worksheet("XRD Pattern with Peak Labels")
+proj.addChild(worksheet)
+
+worksheet.setUseViewSize(False)
+w = Worksheet.convertToSceneUnits(22, Worksheet.Unit.Centimeter)
+h = Worksheet.convertToSceneUnits(14, Worksheet.Unit.Centimeter)
+worksheet.setPageRect(QRectF(0, 0, w, h))
+
+worksheet.setTheme("Bright")
+
+margin = Worksheet.convertToSceneUnits(0.5, Worksheet.Unit.Centimeter)
+worksheet.setLayoutTopMargin(margin)
+worksheet.setLayoutBottomMargin(margin)
+worksheet.setLayoutLeftMargin(margin)
+worksheet.setLayoutRightMargin(margin)
+
+# === Create Plot ===
+plot = CartesianPlot("XRD Pattern")
+plot.setType(CartesianPlot.Type.FourAxes)
+plot.title().setText("X-Ray Diffraction Pattern with Automatic Peak Detection")
+plot.setNiceExtend(True) # Auto-extend ranges for cleaner appearance
+
+worksheet.addChild(plot)
+
+# Configure axes - do this AFTER adding to worksheet
+x_axis = plot.horizontalAxis()
+x_axis.title().setText("2θ (degrees)")
+x_axis.majorGridLine().setStyle(Qt.PenStyle.NoPen)
+x_axis.minorGridLine().setStyle(Qt.PenStyle.NoPen)
+
+y_axis = plot.verticalAxis()
+y_axis.title().setText("Intensity (counts)")
+y_axis.majorGridLine().setStyle(Qt.PenStyle.NoPen)
+y_axis.minorGridLine().setStyle(Qt.PenStyle.NoPen)
+
+# === Plot Spectrum ===
+curve = XYCurve("XRD Pattern")
+curve.setXColumn(spreadsheet.column(0))
+curve.setYColumn(spreadsheet.column(1))
+curve.setLineType(XYCurve.LineType.Line)
+curve.line().setWidth(Worksheet.convertToSceneUnits(1.5, Worksheet.Unit.Point))
+curve.symbol().setStyle(Symbol.Style.NoSymbols)
+plot.addChild(curve)
+
+# === Mark Peaks with Symbols ===
+peaks_curve = XYCurve("Detected Peaks")
+peaks_curve.setXColumn(peaks_spreadsheet.column(0))
+peaks_curve.setYColumn(peaks_spreadsheet.column(1))
+peaks_curve.setLineType(XYCurve.LineType.NoLine)
+symbol = peaks_curve.symbol()
+symbol.setStyle(Symbol.Style.Circle)
+symbol.setSize(Worksheet.convertToSceneUnits(8, Worksheet.Unit.Point))
+pen = symbol.pen()
+pen.setWidth(Worksheet.convertToSceneUnits(2, Worksheet.Unit.Point))
+symbol.setPen(pen)
+plot.addChild(peaks_curve)
+
+# === Add Peak Labels ===
+# Create TextLabels for each peak
+label_font = QFont()
+label_font.setPointSize(9)
+
+te = QTextEdit()
+te.setFont(label_font)
+
+# Add reference lines and labels for major peaks (top N by prominence)
+# Sort peaks by prominence
+peaks_sorted = sorted(peaks, key=lambda p: p['prominence'], reverse=True)
+n_labels = min(10, len(peaks_sorted)) # Label top 10 peaks
+
+for i, peak in enumerate(peaks_sorted[:n_labels]):
+ # Add vertical reference line at peak position
+ ref_line = ReferenceLine(plot, f"Peak {i+1}")
+ plot.addChild(ref_line)
+ ref_line.setOrientation(ReferenceLine.Orientation.Vertical)
+ ref_line.setPositionLogical(QPointF(peak['x_pos'], 0))
+ ref_line.line().setStyle(Qt.PenStyle.DashLine)
+ ref_line.line().setWidth(Worksheet.convertToSceneUnits(0, Worksheet.Unit.Point))
+ ref_line.line().setOpacity(0.5)
+ ref_line.retransform() # Make line visible
+
+ # Add text label above peak - as child of plot, using logical coordinates
+ label = TextLabel(f"Peak Label {i+1}")
+ plot.addChild(label)
+
+ te.clear()
+ te.setPlainText(f"{peak['x_pos']:.1f}°")
+ label.setText(te.toHtml())
+
+ # Enable coordinate binding and position label above the peak in logical (data) coordinates
+ label.setCoordinateBindingEnabled(True)
+ label.setPositionLogical(QPointF(peak['x_pos'], peak['y_pos'] * 1.05)) # 5% above peak
+
+# === Add Summary Information ===
+info_label = TextLabel("Peak Detection Info")
+worksheet.addChild(info_label)
+
+te.clear()
+te.setFont(label_font)
+te.setPlainText(f"Detected {len(peaks)} peaks\nProminence threshold: {prominence_threshold}\nTop {n_labels} peaks labeled")
+info_label.setText(te.toHtml())
+info_label.setPositionScene(QPointF(
+ Worksheet.convertToSceneUnits(16, Worksheet.Unit.Centimeter),
+ Worksheet.convertToSceneUnits(11, Worksheet.Unit.Centimeter)
+))
+
+# Add legend
+plot.addLegend()
+
+# === Print Summary ===
+print("=" * 70)
+print("PEAK ANNOTATION COMPLETE")
+print("=" * 70)
+print()
+print(f"✓ Plotted spectrum with {len(two_theta)} data points")
+print(f"✓ Detected {len(peaks)} peaks automatically")
+print(f"✓ Labeled top {n_labels} peaks by prominence")
+print(f"✓ Added reference lines for major peaks")
+print()
+print("Customization options:")
+print(" • Adjust 'prominence_threshold' to find more/fewer peaks")
+print(" • Adjust 'distance_threshold' to merge/separate nearby peaks")
+print(" • Change 'n_labels' to annotate more/fewer peaks")
+print()
+print("Export options:")
+print(" • worksheet.exportToFile('xrd_pattern.pdf', Worksheet.ExportFormat.PDF)")
+print(" • peaks_spreadsheet can be exported to CSV for further analysis")
+print()
+print("=" * 70)
+print("TIP: This workflow works for ANY 1D signal!")
+print("=" * 70)
+print("Just replace the data with:")
+print(" • Raman/IR spectrum (wavenumber vs intensity)")
+print(" • HPLC chromatogram (time vs absorbance)")
+print(" • Mass spectrum (m/z vs abundance)")
+print(" • EEG/ECG signal (time vs voltage)")
+print()
+
+# Optional: Export
+# worksheet.exportToFile("peak_detection.pdf", Worksheet.ExportFormat.PDF)
+# worksheet.exportToFile("peak_detection.png", Worksheet.ExportFormat.PNG, 300)
diff --git a/lib/examples/Publication Multi-Panel Figure/CMakeLists.txt b/lib/examples/Publication Multi-Panel Figure/CMakeLists.txt
new file mode 100644
index 0000000000..1fb1d09a99
--- /dev/null
+++ b/lib/examples/Publication Multi-Panel Figure/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_executable(PublicationFigure main.cpp)
+
+target_link_libraries(PublicationFigure LabPlot::SDK)
diff --git a/lib/examples/Publication Multi-Panel Figure/README.md b/lib/examples/Publication Multi-Panel Figure/README.md
new file mode 100644
index 0000000000..bc2cf62e4d
--- /dev/null
+++ b/lib/examples/Publication Multi-Panel Figure/README.md
@@ -0,0 +1,149 @@
+# Publication-Ready Multi-Panel Figures
+
+This tutorial demonstrates how to create professional, publication-quality figures with multiple panels suitable for scientific journals.
+
+## What You'll Learn
+
+- Creating multi-panel layouts (2×2, 2×3, etc.)
+- Consistent styling across all panels
+- Adding panel labels (a), (b), (c), (d)
+- Different plot types in one figure
+- High-DPI export for journals (300+ DPI)
+- Proper sizing for journal requirements
+
+## The Result
+
+A complete "Figure 1" with four panels showing:
+- **(a)** Time series data with scatter plot
+- **(b)** Dose-response curve with sigmoid fit (log scale X-axis)
+- **(c)** Bar chart with error bars comparing treatment groups
+- **(d)** Scatter plot with linear regression and correlation
+
+## Why This Matters
+
+**Every scientific paper needs multi-panel figures!** Journals require:
+- Multiple related plots in one figure
+- Consistent styling (fonts, colors, sizes)
+- Proper panel labels for reference in text
+- Specific dimensions (single column ≈ 8.5 cm, double column ≈ 17 cm)
+- High resolution (≥300 DPI for print)
+
+This script automates the entire process, ensuring consistency and saving hours of manual work.
+
+## Key Features Demonstrated
+
+### Layout Control
+```python
+worksheet.setLayout(Worksheet.Layout.GridLayout)
+worksheet.setLayoutRowCount(2)
+worksheet.setLayoutColumnCount(2)
+```
+
+### Panel Labels
+Automatically positioned TextLabels with bold font:
+```python
+label_texts = ['(a)', '(b)', '(c)', '(d)']
+# Positioned at top-left of each panel
+```
+
+### Multiple Plot Types
+- Line plots with symbols (time series)
+- Scatter plots with log scale (dose-response)
+- Bar charts with error bars (grouped data)
+- Scatter with regression line (correlation)
+
+### High-Quality Export
+```python
+# PDF for print journals
+worksheet.exportToFile("Figure1.pdf", Worksheet.ExportFormat.PDF)
+
+# PNG at 300 DPI for online submissions
+worksheet.exportToFile("Figure1.png", Worksheet.ExportFormat.PNG, 300)
+
+# SVG for infinite scalability
+worksheet.exportToFile("Figure1.svg", Worksheet.ExportFormat.SVG)
+```
+
+## Usage
+
+1. **Run the script** - It generates synthetic data for demonstration
+2. **Customize for your data** - Replace the data generation section with your actual data import
+3. **Adjust layout** - Change row/column counts for different arrangements
+4. **Export** - Uncomment the export lines at the end
+
+## Customization Tips
+
+### Change Layout to 1×3 (horizontal strip)
+```python
+worksheet.setLayoutRowCount(1)
+worksheet.setLayoutColumnCount(3)
+```
+
+### Adjust Figure Size for Single Column
+```python
+w = Worksheet.convertToSceneUnits(8.5, Worksheet.Unit.Centimeter)
+h = Worksheet.convertToSceneUnits(12, Worksheet.Unit.Centimeter)
+```
+
+### Add Statistical Annotations
+Use `TextLabel` to add p-values, R², or other stats directly on plots:
+```python
+stat_label = TextLabel("R² = 0.92, p < 0.001")
+plot_d.addChild(stat_label)
+```
+
+### Consistent Colors Across Panels
+Define colors once and reuse:
+```python
+color1 = "#1f77b4" # Blue
+color2 = "#ff7f0e" # Orange
+curve_a.line().setColor(color1)
+bar_c.borderLineAt(0).setColor(color1)
+```
+
+## Journal Requirements
+
+Common specifications:
+- **Nature**: 89 mm (single) or 183 mm (double) width, 300 DPI minimum
+- **Science**: 5.5 cm (single) or 12 cm (double) width, CMYK or RGB
+- **PLOS**: 8.30 cm or 17.35 cm width, 300-600 DPI
+- **Cell**: 8.5 cm or 17.4 cm width, RGB, 300 DPI minimum
+
+This script uses standard double-column width (17 cm) that works for most journals.
+
+## Real-World Applications
+
+Replace the synthetic data with your actual experiments:
+
+### Example: Drug Study
+- Panel A: Cell viability over time
+- Panel B: Dose-response relationship
+- Panel C: Comparison across cell lines
+- Panel D: Drug A vs Drug B correlation
+
+### Example: Materials Science
+- Panel A: XRD pattern
+- Panel B: TEM image
+- Panel C: Particle size distribution
+- Panel D: Mechanical properties
+
+### Example: Climate Data
+- Panel A: Temperature time series
+- Panel B: Precipitation histogram
+- Panel C: Regional comparison
+- Panel D: Temperature vs CO₂ correlation
+
+## Tips for Publication Success
+
+1. **Use vector formats** (PDF, SVG) when possible - they scale perfectly
+2. **Check journal guidelines** before finalizing dimensions
+3. **Maintain consistent fonts** across all panels (usually Arial or Helvetica)
+4. **Keep panel labels bold** and positioned consistently
+5. **Test print** your figure at actual size before submission
+6. **Use colorblind-friendly palettes** for accessibility
+
+## Related Examples
+
+- See `../Demo/script.py` for data import and fitting basics
+- See `../Basic Plots/script.py` for more plot type examples
+- See `../Tufte's Minimal Ink Design/script.py` for minimalist publication style
diff --git a/lib/examples/Publication Multi-Panel Figure/main.py b/lib/examples/Publication Multi-Panel Figure/main.py
new file mode 100644
index 0000000000..cca0e19a3b
--- /dev/null
+++ b/lib/examples/Publication Multi-Panel Figure/main.py
@@ -0,0 +1,10 @@
+#!/usr/bin/env python3
+"""Runner for publication figure demo"""
+import sys
+import os
+
+script_dir = os.path.dirname(os.path.abspath(__file__))
+script_path = os.path.join(script_dir, 'script.py')
+
+with open(script_path, 'r') as f:
+ exec(f.read())
diff --git a/lib/examples/Publication Multi-Panel Figure/script.py b/lib/examples/Publication Multi-Panel Figure/script.py
new file mode 100644
index 0000000000..6927d819d8
--- /dev/null
+++ b/lib/examples/Publication Multi-Panel Figure/script.py
@@ -0,0 +1,322 @@
+"""
+Publication-Ready Multi-Panel Figure Generator
+
+This script demonstrates how to create a publication-quality figure with
+multiple panels (A, B, C, D) suitable for scientific journals.
+
+Features:
+- Automatic panel layout (2×2 grid)
+- Consistent styling across all panels
+- Panel labels (a), (b), (c), (d) in proper position
+- High-DPI export (300 DPI) for publication
+- Professional appearance following journal standards
+
+Use case: Creating "Figure 1" for your paper with multiple related plots
+"""
+
+import numpy as np
+from PySide6.QtCore import QRectF, QPointF, Qt
+from PySide6.QtGui import QFont
+from PySide6.QtWidgets import QTextEdit
+from pylabplot import *
+
+# === Generate Sample Data ===
+# In real use, you would import your actual data files
+np.random.seed(42)
+
+# Panel A: Time series data
+t = np.linspace(0, 10, 100)
+signal_a = np.sin(2 * np.pi * 0.5 * t) * np.exp(-t / 10) + 0.1 * np.random.randn(100)
+
+# Panel B: Dose-response curve
+dose = np.logspace(-2, 2, 20)
+response = 100 / (1 + (dose / 1.0)**(-2)) + 5 * np.random.randn(20)
+
+# Panel C: Grouped bar chart data
+groups = ['Control', 'Treatment A', 'Treatment B', 'Treatment C']
+values_c = [45, 67, 58, 72]
+errors_c = [5, 7, 6, 8]
+
+# Panel D: Scatter plot with correlation
+x_d = np.random.randn(50) * 2 + 5
+y_d = 1.5 * x_d + 3 + np.random.randn(50) * 2
+
+# === Setup Project ===
+proj = project()
+
+# Create spreadsheets for each panel
+spreadsheet_a = Spreadsheet("Panel A Data")
+spreadsheet_a.setColumnCount(2) # time, signal
+spreadsheet_b = Spreadsheet("Panel B Data")
+spreadsheet_b.setColumnCount(2) # dose, response
+spreadsheet_c = Spreadsheet("Panel C Data")
+spreadsheet_c.setColumnCount(3) # group, value, error
+spreadsheet_d = Spreadsheet("Panel D Data")
+spreadsheet_d.setColumnCount(2) # x, y
+
+proj.addChild(spreadsheet_a)
+proj.addChild(spreadsheet_b)
+proj.addChild(spreadsheet_c)
+proj.addChild(spreadsheet_d)
+
+# === Fill Data ===
+# Panel A
+spreadsheet_a.column(0).setName("Time (s)")
+spreadsheet_a.column(1).setName("Signal (mV)")
+spreadsheet_a.column(0).replaceValues(0, [float(x) for x in t])
+spreadsheet_a.column(1).replaceValues(0, [float(x) for x in signal_a])
+
+# Panel B
+spreadsheet_b.column(0).setName("Dose (μM)")
+spreadsheet_b.column(1).setName("Response (%)")
+spreadsheet_b.column(0).replaceValues(0, [float(x) for x in dose])
+spreadsheet_b.column(1).replaceValues(0, [float(x) for x in response])
+
+# Panel C
+spreadsheet_c.column(0).setName("Group")
+spreadsheet_c.column(1).setName("Value")
+spreadsheet_c.column(2).setName("Error")
+spreadsheet_c.column(0).replaceValues(0, [float(i+1) for i in range(len(groups))])
+spreadsheet_c.column(1).replaceValues(0, [float(x) for x in values_c])
+spreadsheet_c.column(2).replaceValues(0, [float(x) for x in errors_c])
+
+# Panel D
+spreadsheet_d.column(0).setName("X Variable")
+spreadsheet_d.column(1).setName("Y Variable")
+spreadsheet_d.column(0).replaceValues(0, [float(x) for x in x_d])
+spreadsheet_d.column(1).replaceValues(0, [float(x) for x in y_d])
+
+# === Create Worksheet ===
+worksheet = Worksheet("Figure 1")
+proj.addChild(worksheet)
+
+# Set size appropriate for publication (single column = 8.5 cm, double column = 17.4 cm)
+# Using double column width
+worksheet.setUseViewSize(False)
+w = Worksheet.convertToSceneUnits(17, Worksheet.Unit.Centimeter)
+h = Worksheet.convertToSceneUnits(17, Worksheet.Unit.Centimeter)
+worksheet.setPageRect(QRectF(0, 0, w, h))
+
+# Use publication-appropriate theme (clean, high contrast)
+worksheet.setTheme("Tufte")
+
+# Setup 2×2 grid layout
+worksheet.setLayout(Worksheet.Layout.GridLayout)
+worksheet.setLayoutRowCount(2)
+worksheet.setLayoutColumnCount(2)
+
+# Tight margins for publication (journals prefer compact figures)
+margin = Worksheet.convertToSceneUnits(0.3, Worksheet.Unit.Centimeter)
+worksheet.setLayoutTopMargin(margin)
+worksheet.setLayoutBottomMargin(margin)
+worksheet.setLayoutLeftMargin(margin)
+worksheet.setLayoutRightMargin(margin)
+worksheet.setLayoutHorizontalSpacing(Worksheet.convertToSceneUnits(0.8, Worksheet.Unit.Centimeter))
+worksheet.setLayoutVerticalSpacing(Worksheet.convertToSceneUnits(0.8, Worksheet.Unit.Centimeter))
+
+# === Panel A: Time Series ===
+plot_a = CartesianPlot("Panel A")
+plot_a.setType(CartesianPlot.Type.FourAxes)
+plot_a.setNiceExtend(True)
+
+worksheet.addChild(plot_a)
+
+# Configure axes
+x_axis = plot_a.horizontalAxis()
+x_axis.title().setText("Time (s)")
+x_axis.majorGridLine().setStyle(Qt.PenStyle.NoPen)
+x_axis.minorGridLine().setStyle(Qt.PenStyle.NoPen)
+
+y_axis = plot_a.verticalAxis()
+y_axis.title().setText("Signal (mV)")
+y_axis.majorGridLine().setStyle(Qt.PenStyle.NoPen)
+y_axis.minorGridLine().setStyle(Qt.PenStyle.NoPen)
+
+curve_a = XYCurve("Time Series")
+curve_a.setXColumn(spreadsheet_a.column(0))
+curve_a.setYColumn(spreadsheet_a.column(1))
+curve_a.setLineType(XYCurve.LineType.Line)
+curve_a.symbol().setStyle(Symbol.Style.Circle)
+curve_a.symbol().setSize(Worksheet.convertToSceneUnits(3, Worksheet.Unit.Point))
+plot_a.addChild(curve_a)
+
+# === Panel B: Dose-Response (log scale) ===
+plot_b = CartesianPlot("Panel B")
+plot_b.setType(CartesianPlot.Type.FourAxes)
+plot_b.setNiceExtend(True)
+
+worksheet.addChild(plot_b)
+
+# Configure axis titles
+x_axis = plot_b.horizontalAxis()
+x_axis.title().setText("Dose (μM)")
+x_axis.majorGridLine().setStyle(Qt.PenStyle.NoPen)
+x_axis.minorGridLine().setStyle(Qt.PenStyle.NoPen)
+
+y_axis = plot_b.verticalAxis()
+y_axis.title().setText("Response (%)")
+y_axis.majorGridLine().setStyle(Qt.PenStyle.NoPen)
+y_axis.minorGridLine().setStyle(Qt.PenStyle.NoPen)
+
+# Set X axis to logarithmic scale (via Range object)
+plot_b.enableAutoScale(CartesianCoordinateSystem.Dimension.X, 0, False)
+rangeX_b = plot_b.range(CartesianCoordinateSystem.Dimension.X, 0)
+rangeX_b.setScale(RangeT.Scale.Log10)
+rangeX_b.setRange(0.01, 100) # Set appropriate range for log scale
+plot_b.setRange(CartesianCoordinateSystem.Dimension.X, 0, rangeX_b)
+
+curve_b = XYCurve("Dose Response")
+curve_b.setXColumn(spreadsheet_b.column(0))
+curve_b.setYColumn(spreadsheet_b.column(1))
+curve_b.setLineType(XYCurve.LineType.NoLine)
+curve_b.symbol().setStyle(Symbol.Style.Square)
+curve_b.symbol().setSize(Worksheet.convertToSceneUnits(5, Worksheet.Unit.Point))
+plot_b.addChild(curve_b)
+
+# Add sigmoid fit
+fit_b = XYFitCurve("Sigmoid Fit")
+fit_b.setXDataColumn(spreadsheet_b.column(0))
+fit_b.setYDataColumn(spreadsheet_b.column(1))
+
+fitData_b = fit_b.fitData()
+fitData_b.modelCategory = nsl_fit_model_category.nsl_fit_model_growth
+fitData_b.modelType = nsl_fit_model_type_growth.nsl_fit_model_sigmoid
+fitData_b = XYFitCurve.initFitData(fitData_b)
+fitData_b = fit_b.initStartValues(fitData_b)
+fit_b.setFitData(fitData_b)
+fit_b.recalculate()
+plot_b.addChild(fit_b)
+
+# === Panel C: Grouped Bar Chart with Error Bars ===
+plot_c = CartesianPlot("Panel C")
+plot_c.setType(CartesianPlot.Type.FourAxes)
+plot_c.setNiceExtend(True)
+
+worksheet.addChild(plot_c)
+
+# Configure axes
+x_axis = plot_c.horizontalAxis()
+x_axis.title().setText("Treatment Group")
+x_axis.majorGridLine().setStyle(Qt.PenStyle.NoPen)
+x_axis.minorGridLine().setStyle(Qt.PenStyle.NoPen)
+# Set custom labels
+x_axis.setMajorTicksType(Axis.TicksType.CustomColumn)
+x_axis.setMajorTicksColumn(spreadsheet_c.column(0))
+
+y_axis = plot_c.verticalAxis()
+y_axis.title().setText("Effect Size (%)")
+y_axis.majorGridLine().setStyle(Qt.PenStyle.NoPen)
+y_axis.minorGridLine().setStyle(Qt.PenStyle.NoPen)
+
+bar_c = BarPlot("Treatment Effects")
+bar_c.setDataColumns([spreadsheet_c.column(1)])
+bar_c.setOrientation(BarPlot.Orientation.Vertical)
+bar_c.setType(BarPlot.Type.Grouped)
+
+plot_c.addChild(bar_c)
+
+# Add error bars via ErrorBar object
+error_bar = bar_c.errorBarAt(0) # Get error bar for first (only) data series
+error_bar.setYErrorType(ErrorBar.ErrorType.Symmetric)
+error_bar.setYPlusColumn(spreadsheet_c.column(2))
+
+# === Panel D: Scatter with Linear Regression ===
+plot_d = CartesianPlot("Panel D")
+plot_d.setType(CartesianPlot.Type.FourAxes)
+plot_d.setNiceExtend(True)
+
+worksheet.addChild(plot_d)
+
+x_axis = plot_d.horizontalAxis()
+x_axis.title().setText("X Variable (a.u.)")
+x_axis.majorGridLine().setStyle(Qt.PenStyle.NoPen)
+x_axis.minorGridLine().setStyle(Qt.PenStyle.NoPen)
+
+y_axis = plot_d.verticalAxis()
+y_axis.title().setText("Y Variable (a.u.)")
+y_axis.majorGridLine().setStyle(Qt.PenStyle.NoPen)
+y_axis.minorGridLine().setStyle(Qt.PenStyle.NoPen)
+
+scatter_d = XYCurve("Data Points")
+scatter_d.setXColumn(spreadsheet_d.column(0))
+scatter_d.setYColumn(spreadsheet_d.column(1))
+scatter_d.setLineType(XYCurve.LineType.NoLine)
+scatter_d.symbol().setStyle(Symbol.Style.Circle)
+scatter_d.symbol().setSize(Worksheet.convertToSceneUnits(4, Worksheet.Unit.Point))
+scatter_d.symbol().setOpacity(0.6)
+plot_d.addChild(scatter_d)
+
+# Linear regression
+fit_d = XYFitCurve("Linear Fit")
+fit_d.setXDataColumn(spreadsheet_d.column(0))
+fit_d.setYDataColumn(spreadsheet_d.column(1))
+
+fitData_d = fit_d.fitData()
+fitData_d.modelCategory = nsl_fit_model_category.nsl_fit_model_basic
+fitData_d.modelType = nsl_fit_model_type_basic.nsl_fit_model_polynomial
+fitData_d.degree = 1 # Linear
+fitData_d = XYFitCurve.initFitData(fitData_d)
+fitData_d = fit_d.initStartValues(fitData_d)
+fit_d.setFitData(fitData_d)
+fit_d.recalculate()
+plot_d.addChild(fit_d)
+
+# === Add Panel Labels (a), (b), (c), (d) ===
+# These are positioned at top-left of each panel
+label_font = QFont()
+label_font.setPointSize(12)
+label_font.setBold(True)
+
+label_positions = [
+ (-0.5, 8.0), # Panel (a) - top left
+ (8.0, 8.0), # Panel (b) - top right
+ (-0.5, -0.5), # Panel (c) - bottom left
+ (8.0, -0.5) # Panel (d) - bottom right
+]
+
+label_texts = ['(a)', '(b)', '(c)', '(d)']
+
+te = QTextEdit()
+for i, (label_text, (x, y)) in enumerate(zip(label_texts, label_positions)):
+ label = TextLabel(f"Label {label_text}")
+ worksheet.addChild(label)
+
+ te.clear()
+ te.setFont(label_font)
+ te.setPlainText(label_text)
+ label.setText(te.toHtml())
+
+ x_pos = Worksheet.convertToSceneUnits(x, Worksheet.Unit.Centimeter)
+ y_pos = Worksheet.convertToSceneUnits(y, Worksheet.Unit.Centimeter)
+ label.setPositionScene(QPointF(x_pos, y_pos))
+
+print("=" * 60)
+print("Publication Figure Created Successfully!")
+print("=" * 60)
+print()
+print("Figure Layout: 2×2 grid with panels (a), (b), (c), (d)")
+print()
+print("Panel (a): Time series with noise")
+print("Panel (b): Dose-response curve with sigmoid fit")
+print("Panel (c): Grouped bar chart with error bars")
+print("Panel (d): Scatter plot with linear regression")
+print()
+print("=" * 60)
+print("Export Options for Publication:")
+print("=" * 60)
+print()
+print("Use worksheet.exportToFile() for high-quality export:")
+print(' worksheet.exportToFile("Figure1.pdf", Worksheet.ExportFormat.PDF)')
+print(' worksheet.exportToFile("Figure1.png", Worksheet.ExportFormat.PNG, 300) # 300 DPI')
+print(' worksheet.exportToFile("Figure1.svg", Worksheet.ExportFormat.SVG)')
+print()
+print("Most journals prefer:")
+print(" - PDF or TIFF for print")
+print(" - 300 DPI minimum resolution")
+print(" - RGB color mode")
+print(" - Single or double column width (8.5 or 17.4 cm)")
+print()
+
+# Uncomment to auto-export:
+# worksheet.exportToFile("Figure1.pdf", Worksheet.ExportFormat.PDF)
+# worksheet.exportToFile("Figure1_300dpi.png", Worksheet.ExportFormat.PNG, 300)
diff --git a/lib/examples/Quantum Wave Packet/CMakeLists.txt b/lib/examples/Quantum Wave Packet/CMakeLists.txt
new file mode 100644
index 0000000000..c1a5e889c3
--- /dev/null
+++ b/lib/examples/Quantum Wave Packet/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_executable(QuantumWavePacket main.cpp)
+
+target_link_libraries(QuantumWavePacket LabPlot::SDK)
diff --git a/lib/examples/Quantum Wave Packet/README.md b/lib/examples/Quantum Wave Packet/README.md
new file mode 100644
index 0000000000..ce4455e3db
--- /dev/null
+++ b/lib/examples/Quantum Wave Packet/README.md
@@ -0,0 +1,109 @@
+# Quantum Wave Packet Evolution
+
+This example demonstrates the dynamic capabilities of LabPlot's Python scripting by simulating the time evolution of quantum mechanical wave packets.
+
+**Two demos are included:**
+1. **script.py** - Free particle wave packet dispersion
+2. **script_tunneling.py** - Quantum tunneling through a potential barrier
+
+## What it demonstrates
+
+**Physics:**
+- Gaussian wave packet propagating as a free particle
+- Time evolution according to the Schrödinger equation
+- Wave packet dispersion (spreading over time)
+
+**LabPlot Scripting Features:**
+- Creating project structure (Worksheet, CartesianPlot, Spreadsheet) via Python
+- Dynamic data generation with NumPy
+- Real-time column data updates for animation
+- Multiple synchronized curves on a single plot
+- Theme application and layout control
+
+## The Physics
+
+A free particle described by a Gaussian wave packet:
+
+```
+ψ(x,0) = (1/√(2πσ₀²))^(1/2) exp(ik₀x) exp(-(x-x₀)²/(4σ₀²))
+```
+
+evolves according to the Schrödinger equation. The analytical solution shows:
+- The wave packet center moves with group velocity v = ℏk₀/m
+- The width increases as σ(t) = σ₀√(1 + (ℏt/mσ₀²)²)
+- The probability density |ψ(x,t)|² spreads out over time
+
+## What you'll see
+
+### Demo 1: Free Particle (script.py)
+
+The animation shows three curves:
+- **Re(ψ)** - Real part of the wave function (blue)
+- **Im(ψ)** - Imaginary part of the wave function (orange)
+- **|ψ|²** - Probability density (green, bold)
+
+As time progresses, you'll observe:
+1. The wave packet moves to the right (positive momentum k₀ > 0)
+2. The oscillations maintain constant phase velocity
+3. The envelope spreads out (quantum dispersion)
+4. The peak of |ψ|² decreases (probability conservation with spreading)
+
+### Demo 2: Quantum Tunneling (script_tunneling.py)
+
+The animation shows:
+- **|ψ|²** - Probability density (approaching the barrier)
+- **Re(ψ)** - Real part showing interference
+- **Barrier (V/E)** - Potential energy normalized to particle energy (shaded region)
+
+You'll observe the fascinating quantum phenomenon:
+1. Wave packet approaches barrier from the left
+2. Partial reflection at the barrier entrance
+3. **Tunneling** - part of the wave appears beyond the barrier despite E < V₀!
+4. Transmitted and reflected components separate
+5. Final transmission/reflection probabilities printed at end
+
+## Parameters
+
+**Demo 1 (Free particle)** - Edit script.py to explore:
+- `k0`: Initial momentum (try negative for leftward motion)
+- `sigma0`: Initial width (smaller = faster dispersion)
+- `N_steps`: Number of frames (more = smoother animation)
+- `t_max`: Total simulation time
+
+**Demo 2 (Tunneling)** - Edit script_tunneling.py to explore:
+- `V0`: Barrier height relative to energy (try `V0 = 0.5*E` or `V0 = 1.2*E`)
+- `barrier_x2 - barrier_x1`: Barrier width (wider = less tunneling)
+- `k0`: Particle momentum (higher energy = more tunneling)
+- Watch the transmission/reflection probabilities!
+
+## Running the scripts
+
+**Method 1: Within LabPlot (Recommended)**
+1. Open LabPlot
+2. Create a new Script (File → New → Script)
+3. Select Python as the language
+4. Copy the contents of `script.py` or `script_tunneling.py`
+5. Click "Run" and watch the evolution of the quantum wave packet!
+
+**Method 2: Standalone C++ executable**
+1. Build with CMake (if examples are enabled)
+2. Run `./QuantumWavePacket` from the build directory
+
+**Method 3: Python script (if pylabplot installed)**
+```bash
+python3 main.py # Runs script.py
+python3 script_tunneling.py # Runs tunneling demo
+```
+
+## Educational value
+
+This demo is perfect for:
+- Quantum mechanics courses (visualizing abstract wave functions)
+- Demonstrating LabPlot's scientific computing capabilities
+- Teaching students about wave-particle duality
+- Showing how to create interactive scientific visualizations
+
+## References
+
+- Griffiths, "Introduction to Quantum Mechanics", Chapter 2
+- The blog post that inspired this: https://ben.land/post/2022/03/09/quantum-mechanics-simulation/
diff --git a/lib/examples/Quantum Wave Packet/main.cpp b/lib/examples/Quantum Wave Packet/main.cpp
new file mode 100644
index 0000000000..2259ad4a46
--- /dev/null
+++ b/lib/examples/Quantum Wave Packet/main.cpp
@@ -0,0 +1,203 @@
+#include <QApplication>
+#include <QTimer>
+#include <QElapsedTimer>
+#include <cmath>
+#include <complex>
+#include <vector>
+
+#include <labplot.h>
+
+// Physical constants (natural units: ℏ = m = 1)
+constexpr double HBAR = 1.0;
+constexpr double MASS = 1.0;
+
+// Wave packet parameters
+constexpr double X0 = 0.0; // Initial position
+constexpr double K0 = 5.0; // Initial wave number
+constexpr double SIGMA0 = 1.0; // Initial width
+
+// Spatial grid
+constexpr int N_POINTS = 512;
+constexpr double X_MIN = -10.0;
+constexpr double X_MAX = 10.0;
+
+// Time parameters
+constexpr int N_STEPS = 100;
+constexpr double T_MAX = 3.0;
+constexpr double DT = T_MAX / N_STEPS;
+
+// Animation delay (ms)
+constexpr int FRAME_DELAY = 50;
+
+/**
+ * Compute the wave function ψ(x,t) for a free Gaussian wave packet
+ * Returns complex-valued wave function
+ */
+std::vector<std::complex<double>> computeWaveFunction(const std::vector<double>& x, double t) {
+ std::vector<std::complex<double>> psi(x.size());
+
+ // Time-dependent width
+ double sigma_t = SIGMA0 * std::sqrt(1.0 + std::pow(HBAR * t / (MASS * SIGMA0 * SIGMA0), 2));
+
+ // Normalization factor
+ double norm = 1.0 / std::sqrt(sigma_t * std::sqrt(M_PI));
+
+ for (size_t i = 0; i < x.size(); ++i) {
+ double dx = x[i] - X0;
+ double dx_shifted = dx - HBAR * K0 * t / MASS;
+
+ // Phase factors
+ std::complex<double> phase1(0.0, K0 * dx);
+ std::complex<double> phase2(-dx_shifted * dx_shifted / (2.0 * sigma_t * sigma_t), 0.0);
+ std::complex<double> phase3(0.0, HBAR * K0 * K0 * t / (2.0 * MASS));
+
+ psi[i] = norm * std::exp(phase1 + phase2 + phase3);
+ }
+
+ return psi;
+}
+
+int main(int argc, char** argv) {
+ QApplication app(argc, argv);
+
+ // Setup spatial grid
+ std::vector<double> x(N_POINTS);
+ double dx = (X_MAX - X_MIN) / (N_POINTS - 1);
+ for (int i = 0; i < N_POINTS; ++i) {
+ x[i] = X_MIN + i * dx;
+ }
+
+ // Create project structure
+ auto* project = new Project();
+
+ // Create spreadsheet with data columns
+ auto* spreadsheet = new Spreadsheet(QStringLiteral("Wave Data"));
+ project->addChild(spreadsheet);
+
+ // Get columns (automatically created)
+ auto* colX = spreadsheet->column(0);
+ auto* colReal = spreadsheet->column(1);
+ auto* colImag = spreadsheet->column(2);
+ auto* colProb = spreadsheet->column(3);
+
+ colX->setName(QStringLiteral("x"));
+ colReal->setName(QStringLiteral("Re(ψ)"));
+ colImag->setName(QStringLiteral("Im(ψ)"));
+ colProb->setName(QStringLiteral("|ψ|²"));
+
+ // Initialize with t=0 data
+ QVector<double> x_data, real_data, imag_data, prob_data;
+ x_data.reserve(N_POINTS);
+ real_data.reserve(N_POINTS);
+ imag_data.reserve(N_POINTS);
+ prob_data.reserve(N_POINTS);
+
+ auto psi0 = computeWaveFunction(x, 0.0);
+ for (int i = 0; i < N_POINTS; ++i) {
+ x_data.append(x[i]);
+ real_data.append(psi0[i].real());
+ imag_data.append(psi0[i].imag());
+ prob_data.append(std::norm(psi0[i]));
+ }
+
+ colX->replaceValues(0, x_data);
+ colReal->replaceValues(0, real_data);
+ colImag->replaceValues(0, imag_data);
+ colProb->replaceValues(0, prob_data);
+
+ // Create worksheet
+ auto* worksheet = new Worksheet(QStringLiteral("Quantum Wave Packet"));
+ project->addChild(worksheet);
+
+ worksheet->setUseViewSize(false);
+ double w = Worksheet::convertToSceneUnits(20, Worksheet::Unit::Centimeter);
+ double h = Worksheet::convertToSceneUnits(15, Worksheet::Unit::Centimeter);
+ worksheet->setPageRect(QRectF(0, 0, w, h));
+ worksheet->setTheme(QStringLiteral("Tufte"));
+
+ // Create plot area
+ auto* plotArea = new CartesianPlot(QStringLiteral("Wave Function Plot"));
+ plotArea->setType(CartesianPlot::Type::FourAxes);
+ plotArea->title()->setText(QStringLiteral("Quantum Wave Packet Evolution"));
+ worksheet->addChild(plotArea);
+
+ // Configure axes
+ for (auto* axis : plotArea->children<Axis>()) {
+ if (axis->orientation() == WorksheetElement::Orientation::Horizontal
+ && axis->position() == Axis::Position::Bottom) {
+ axis->title()->setText(QStringLiteral("Position x"));
+ } else if (axis->orientation() == WorksheetElement::Orientation::Vertical
+ && axis->position() == Axis::Position::Left) {
+ axis->title()->setText(QStringLiteral("ψ(x,t)"));
+ }
+ }
+
+ // Create curves
+ auto* curveReal = new XYCurve(QStringLiteral("Re(ψ)"));
+ curveReal->setXColumn(colX);
+ curveReal->setYColumn(colReal);
+ curveReal->setLineType(XYCurve::LineType::Line);
+ curveReal->symbol()->setStyle(Symbol::Style::NoSymbols);
+ plotArea->addChild(curveReal);
+
+ auto* curveImag = new XYCurve(QStringLiteral("Im(ψ)"));
+ curveImag->setXColumn(colX);
+ curveImag->setYColumn(colImag);
+ curveImag->setLineType(XYCurve::LineType::Line);
+ curveImag->symbol()->setStyle(Symbol::Style::NoSymbols);
+ plotArea->addChild(curveImag);
+
+ auto* curveProb = new XYCurve(QStringLiteral("|ψ|²"));
+ curveProb->setXColumn(colX);
+ curveProb->setYColumn(colProb);
+ curveProb->setLineType(XYCurve::LineType::Line);
+ curveProb->symbol()->setStyle(Symbol::Style::NoSymbols);
+ curveProb->line()->setWidth(Worksheet::convertToSceneUnits(3, Worksheet::Unit::Point));
+ plotArea->addChild(curveProb);
+
+ plotArea->addLegend();
+
+ // Show the worksheet
+ worksheet->view()->show();
+
+ // Animation using QTimer
+ int currentStep = 0;
+ QTimer* timer = new QTimer();
+
+ QObject::connect(timer, &QTimer::timeout, [&]() {
+ if (currentStep >= N_STEPS) {
+ timer->stop();
+ qDebug() << "Animation complete!";
+ return;
+ }
+
+ double t = currentStep * DT;
+ auto psi_t = computeWaveFunction(x, t);
+
+ // Update column data
+ real_data.clear();
+ imag_data.clear();
+ prob_data.clear();
+
+ for (const auto& val : psi_t) {
+ real_data.append(val.real());
+ imag_data.append(val.imag());
+ prob_data.append(std::norm(val));
+ }
+
+ colReal->replaceValues(0, real_data);
+ colImag->replaceValues(0, imag_data);
+ colProb->replaceValues(0, prob_data);
+
+ if ((currentStep + 1) % 10 == 0) {
+ qDebug() << "Frame" << (currentStep + 1) << "/" << N_STEPS << "(t =" << t << ")";
+ }
+
+ currentStep++;
+ });
+
+ // Start animation
+ timer->start(FRAME_DELAY);
+
+ return app.exec();
+}
diff --git a/lib/examples/Quantum Wave Packet/main.py b/lib/examples/Quantum Wave Packet/main.py
new file mode 100644
index 0000000000..4bcd7b67e9
--- /dev/null
+++ b/lib/examples/Quantum Wave Packet/main.py
@@ -0,0 +1,23 @@
+#!/usr/bin/env python3
+"""
+Standalone runner for the Quantum Wave Packet demo.
+
+This script can be executed directly with Python if LabPlot's Python
+bindings are installed in your Python environment.
+"""
+
+import sys
+import os
+
+# Add the path to pylabplot if needed
+# sys.path.insert(0, '/path/to/pylabplot')
+
+# Import the script
+script_dir = os.path.dirname(os.path.abspath(__file__))
+script_path = os.path.join(script_dir, 'script.py')
+
+with open(script_path, 'r') as f:
+ script_code = f.read()
+
+# Execute the script
+exec(script_code)
diff --git a/lib/examples/Quantum Wave Packet/script.py b/lib/examples/Quantum Wave Packet/script.py
new file mode 100644
index 0000000000..f28e0d44b5
--- /dev/null
+++ b/lib/examples/Quantum Wave Packet/script.py
@@ -0,0 +1,235 @@
+"""
+Quantum Wave Packet Evolution Demo
+
+This script demonstrates the time evolution of a Gaussian wave packet
+in quantum mechanics, showing:
+- Real and imaginary parts of the wave function ψ(x,t)
+- Probability density |ψ(x,t)|²
+- Wave packet propagation and dispersion
+
+Physics:
+A free particle Gaussian wave packet disperses as it propagates.
+The wave function evolves according to the Schrödinger equation.
+"""
+
+import sys
+import time
+import numpy as np
+from PySide6.QtCore import QCoreApplication
+from pylabplot import *
+
+# === Physical Parameters ===
+hbar = 1.0 # Reduced Planck constant (natural units)
+m = 1.0 # Particle mass (natural units)
+
+# Wave packet initial parameters
+x0 = 0.0 # Initial center position
+k0 = 5.0 # Initial wave number (momentum p = hbar*k)
+sigma0 = 1.0 # Initial width (standard deviation)
+
+# Spatial grid
+N_points = 512
+x_min, x_max = -10.0, 10.0
+x = np.linspace(x_min, x_max, N_points)
+
+# Time parameters
+N_steps = 100 # Number of animation frames
+t_max = 3.0 # Maximum time
+dt = t_max / N_steps
+
+# === Wave Function Evolution ===
+def psi(x, t):
+ """
+ Analytical solution for free particle Gaussian wave packet.
+ Returns complex wave function ψ(x,t).
+ """
+ # Time-dependent width
+ sigma_t = sigma0 * np.sqrt(1 + (hbar * t / (m * sigma0**2))**2)
+
+ # Normalization
+ norm = 1.0 / np.sqrt(sigma_t * np.sqrt(np.pi))
+
+ # Phase factors
+ phase1 = 1j * k0 * (x - x0)
+ phase2 = -(x - x0 - hbar * k0 * t / m)**2 / (2 * sigma_t**2)
+ phase3 = 1j * (hbar * k0**2 * t) / (2 * m)
+
+ return norm * np.exp(phase1 + phase2 + phase3)
+
+# === Setup LabPlot Project ===
+proj = project()
+
+# Check if objects already exist (script was run before)
+existing_spreadsheet = None
+existing_worksheet = None
+
+for child in proj.children(AspectType.Spreadsheet):
+ if child.name() == "Wave Data":
+ existing_spreadsheet = child
+ break
+
+for child in proj.children(AspectType.Worksheet):
+ if child.name() == "Quantum Wave Packet":
+ existing_worksheet = child
+ break
+
+# Determine if we need to create objects or just update data
+create_objects = (existing_spreadsheet is None or existing_worksheet is None)
+
+if create_objects:
+ print("First run - creating project structure...")
+
+ # Create spreadsheet with columns
+ spreadsheet = Spreadsheet("Wave Data")
+ spreadsheet.setColumnCount(4) # We need 4 columns: x, real(ψ), imag(ψ), |ψ|²
+ proj.addChild(spreadsheet)
+
+ # Get column references
+ col_x = spreadsheet.column(0)
+ col_real = spreadsheet.column(1)
+ col_imag = spreadsheet.column(2)
+ col_prob = spreadsheet.column(3)
+
+ # Set column names
+ col_x.setName("x")
+ col_real.setName("Re(ψ)")
+ col_imag.setName("Im(ψ)")
+ col_prob.setName("|ψ|²")
+
+ # Initialize with t=0 data
+ psi_t0 = psi(x, 0)
+ col_x.replaceValues(0, [float(xi) for xi in x])
+ col_real.replaceValues(0, [float(val.real) for val in psi_t0])
+ col_imag.replaceValues(0, [float(val.imag) for val in psi_t0])
+ col_prob.replaceValues(0, [float(abs(val)**2) for val in psi_t0])
+
+ # === Create Worksheet ===
+ worksheet = Worksheet("Quantum Wave Packet")
+ proj.addChild(worksheet)
+
+ # Set worksheet size
+ worksheet.setUseViewSize(False)
+ w = Worksheet.convertToSceneUnits(20, Worksheet.Unit.Centimeter)
+ h = Worksheet.convertToSceneUnits(15, Worksheet.Unit.Centimeter)
+ from PySide6.QtCore import QRectF
+ worksheet.setPageRect(QRectF(0, 0, w, h))
+
+ # Apply a nice theme
+ worksheet.setTheme("Tufte")
+
+ # Setup layout
+ worksheet.setLayout(Worksheet.Layout.VerticalLayout)
+ ms = Worksheet.convertToSceneUnits(0.5, Worksheet.Unit.Centimeter)
+ worksheet.setLayoutTopMargin(ms)
+ worksheet.setLayoutBottomMargin(ms)
+ worksheet.setLayoutLeftMargin(ms)
+ worksheet.setLayoutRightMargin(ms)
+ worksheet.setLayoutVerticalSpacing(ms)
+
+ # === Create Plot Area ===
+ plotArea = CartesianPlot("Wave Function Plot")
+ plotArea.setType(CartesianPlot.Type.FourAxes)
+ plotArea.title().setText("Quantum Wave Packet Evolution")
+
+ # Configure axes
+ for axis in plotArea.children(AspectType.Axis):
+ if axis.orientation() == WorksheetElement.Orientation.Horizontal and axis.position() == Axis.Position.Bottom:
+ axis.title().setText("Position x")
+ elif axis.orientation() == WorksheetElement.Orientation.Vertical and axis.position() == Axis.Position.Left:
+ axis.title().setText("ψ(x,t)")
+
+ worksheet.addChild(plotArea)
+
+ # Set fixed Y range to reduce flickering during animation
+ plotArea.enableAutoScale(CartesianCoordinateSystem.Dimension.Y, 0, False)
+ rangeY = plotArea.range(CartesianCoordinateSystem.Dimension.Y, 0)
+ rangeY.setRange(-1.0, 1.0)
+ plotArea.setRange(CartesianCoordinateSystem.Dimension.Y, 0, rangeY)
+
+ # === Create Curves ===
+ # Real part
+ curve_real = XYCurve("Re(ψ)")
+ curve_real.setXColumn(col_x)
+ curve_real.setYColumn(col_real)
+ curve_real.setLineType(XYCurve.LineType.Line)
+ curve_real.symbol().setStyle(Symbol.Style.NoSymbols)
+ curve_real.line().setWidth(Worksheet.convertToSceneUnits(2, Worksheet.Unit.Point))
+ plotArea.addChild(curve_real)
+
+ # Imaginary part
+ curve_imag = XYCurve("Im(ψ)")
+ curve_imag.setXColumn(col_x)
+ curve_imag.setYColumn(col_imag)
+ curve_imag.setLineType(XYCurve.LineType.Line)
+ curve_imag.symbol().setStyle(Symbol.Style.NoSymbols)
+ curve_imag.line().setWidth(Worksheet.convertToSceneUnits(2, Worksheet.Unit.Point))
+ plotArea.addChild(curve_imag)
+
+ # Probability density
+ curve_prob = XYCurve("|ψ|²")
+ curve_prob.setXColumn(col_x)
+ curve_prob.setYColumn(col_prob)
+ curve_prob.setLineType(XYCurve.LineType.Line)
+ curve_prob.symbol().setStyle(Symbol.Style.NoSymbols)
+ curve_prob.line().setWidth(Worksheet.convertToSceneUnits(3, Worksheet.Unit.Point))
+ plotArea.addChild(curve_prob)
+
+ # Add legend
+ plotArea.addLegend()
+else:
+ print("Objects already exist - reusing and resetting animation...")
+
+ # Reuse existing objects
+ spreadsheet = existing_spreadsheet
+ worksheet = existing_worksheet
+
+ # Get columns
+ col_x = spreadsheet.column(0)
+ col_real = spreadsheet.column(1)
+ col_imag = spreadsheet.column(2)
+ col_prob = spreadsheet.column(3)
+
+ # Get plot area (first CartesianPlot child of worksheet)
+ plotArea = None
+ for child in worksheet.children(AspectType.CartesianPlot):
+ if child.name() == "Wave Function Plot":
+ plotArea = child
+ break
+
+ if plotArea is None:
+ print("Warning: Could not find plot area!")
+
+# Reset to initial state (t=0) for animation
+psi_t0 = psi(x, 0)
+col_real.replaceValues(0, [float(val.real) for val in psi_t0])
+col_imag.replaceValues(0, [float(val.imag) for val in psi_t0])
+col_prob.replaceValues(0, [float(abs(val)**2) for val in psi_t0])
+
+# === Animation Loop ===
+print("Starting wave packet evolution animation...")
+print(f"Time step: {dt:.3f}, Total frames: {N_steps}")
+
+for i in range(N_steps):
+ t = i * dt
+
+ # Compute wave function at time t
+ psi_t = psi(x, t)
+
+ # Update column data
+ col_real.replaceValues(0, [float(val.real) for val in psi_t])
+ col_imag.replaceValues(0, [float(val.imag) for val in psi_t])
+ col_prob.replaceValues(0, [float(abs(val)**2) for val in psi_t])
+
+ # Process Qt events to update the UI
+ QCoreApplication.processEvents()
+
+ # Small delay for animation effect
+ time.sleep(0.05)
+
+ # Progress indicator
+ if (i + 1) % 10 == 0:
+ print(f"Frame {i+1}/{N_steps} (t = {t:.2f})")
+
+print("Animation complete!")
+print(f"Final time: t = {t_max:.2f}")
+print(f"Wave packet has dispersed by factor: {np.sqrt(1 + (hbar * t_max / (m * sigma0**2))**2):.2f}x")
diff --git a/lib/examples/Quantum Wave Packet/script_tunneling.py b/lib/examples/Quantum Wave Packet/script_tunneling.py
new file mode 100644
index 0000000000..fe5f5d704f
--- /dev/null
+++ b/lib/examples/Quantum Wave Packet/script_tunneling.py
@@ -0,0 +1,293 @@
+"""
+Quantum Tunneling Demo - Wave Packet Encounters a Potential Barrier
+
+This script demonstrates quantum tunneling by simulating a Gaussian wave packet
+encountering a rectangular potential barrier.
+
+Physics:
+- A wave packet with definite momentum approaches a barrier
+- Classical: particle reflects if E < V₀
+- Quantum: finite probability of tunneling through the barrier
+- Shows transmission and reflection coefficients visually
+
+This is a more advanced demo showing:
+- Potential energy visualization
+- Split-operator method for time evolution
+- Complex quantum phenomena (tunneling)
+"""
+
+import sys
+import time
+import numpy as np
+from PySide6.QtCore import QCoreApplication
+from pylabplot import *
+
+# === Physical Parameters ===
+hbar = 1.0 # Reduced Planck constant
+m = 1.0 # Particle mass
+
+# Wave packet initial parameters
+x0 = -5.0 # Initial center position (left side)
+k0 = 8.0 # Initial wave number (momentum)
+sigma0 = 0.8 # Initial width
+E = (hbar * k0)**2 / (2 * m) # Energy
+
+# Potential barrier
+V0 = 0.7 * E # Barrier height (< E for partial tunneling)
+barrier_x1 = 0.0 # Barrier start
+barrier_x2 = 2.0 # Barrier end
+
+# Spatial grid
+N_points = 1024
+x_min, x_max = -15.0, 15.0
+x = np.linspace(x_min, x_max, N_points)
+dx = x[1] - x[0]
+
+# Potential energy function
+def V(x_val):
+ """Rectangular potential barrier"""
+ return V0 if (barrier_x1 <= x_val <= barrier_x2) else 0.0
+
+V_array = np.array([V(xi) for xi in x])
+
+# Time parameters
+N_steps = 200
+dt = 0.02
+t_max = N_steps * dt
+
+# === Split-Operator Method for Time Evolution ===
+# This method is more accurate than analytical free-particle solution
+# when a potential is present
+
+def evolve_split_operator(psi, V_array, dt):
+ """
+ Evolve wave function by one time step using split-operator method.
+ U(dt) ≈ exp(-iV dt/2ℏ) · exp(-iT dt/ℏ) · exp(-iV dt/2ℏ)
+ where T is kinetic energy operator.
+ """
+ # Momentum space grid
+ dk = 2 * np.pi / (N_points * dx)
+ k_vals = np.fft.fftfreq(N_points, dx) * 2 * np.pi
+
+ # Half step in position space (potential energy)
+ psi = psi * np.exp(-1j * V_array * dt / (2 * hbar))
+
+ # Full step in momentum space (kinetic energy)
+ psi_k = np.fft.fft(psi)
+ psi_k = psi_k * np.exp(-1j * hbar * k_vals**2 * dt / (2 * m))
+ psi = np.fft.ifft(psi_k)
+
+ # Half step in position space (potential energy)
+ psi = psi * np.exp(-1j * V_array * dt / (2 * hbar))
+
+ return psi
+
+# === Initial Wave Function ===
+# Gaussian wave packet with momentum k0
+psi = np.exp(1j * k0 * (x - x0)) * np.exp(-(x - x0)**2 / (2 * sigma0**2))
+psi = psi / np.sqrt(np.sum(np.abs(psi)**2) * dx) # Normalize
+
+# === Setup LabPlot Project ===
+proj = project()
+
+# Check if objects already exist (script was run before)
+existing_spreadsheet = None
+existing_worksheet = None
+
+for child in proj.children(AspectType.Spreadsheet):
+ if child.name() == "Tunneling Data":
+ existing_spreadsheet = child
+ break
+
+for child in proj.children(AspectType.Worksheet):
+ if child.name() == "Quantum Tunneling":
+ existing_worksheet = child
+ break
+
+# Determine if we need to create objects or just update data
+create_objects = (existing_spreadsheet is None or existing_worksheet is None)
+
+if create_objects:
+ print("First run - creating project structure...")
+
+ # Create spreadsheet
+ spreadsheet = Spreadsheet("Tunneling Data")
+ spreadsheet.setColumnCount(5)
+ proj.addChild(spreadsheet)
+
+ # Columns: x, Re(ψ), Im(ψ), |ψ|², V(x)
+ col_x = spreadsheet.column(0)
+ col_real = spreadsheet.column(1)
+ col_imag = spreadsheet.column(2)
+ col_prob = spreadsheet.column(3)
+ col_potential = spreadsheet.column(4)
+
+ col_x.setName("x")
+ col_real.setName("Re(ψ)")
+ col_imag.setName("Im(ψ)")
+ col_prob.setName("|ψ|²")
+ col_potential.setName("V(x)/E")
+
+ # Initialize x and potential (these don't change)
+ col_x.replaceValues(0, [float(xi) for xi in x])
+ col_potential.replaceValues(0, [float(V(xi) / E) for xi in x])
+
+ # === Create Worksheet ===
+ worksheet = Worksheet("Quantum Tunneling")
+ proj.addChild(worksheet)
+
+ worksheet.setUseViewSize(False)
+ w = Worksheet.convertToSceneUnits(22, Worksheet.Unit.Centimeter)
+ h = Worksheet.convertToSceneUnits(15, Worksheet.Unit.Centimeter)
+ from PySide6.QtCore import QRectF
+ worksheet.setPageRect(QRectF(0, 0, w, h))
+
+ worksheet.setTheme("Bright")
+
+ ms = Worksheet.convertToSceneUnits(0.5, Worksheet.Unit.Centimeter)
+ worksheet.setLayoutTopMargin(ms)
+ worksheet.setLayoutBottomMargin(ms)
+ worksheet.setLayoutLeftMargin(ms)
+ worksheet.setLayoutRightMargin(ms)
+
+ # === Create Plot Area ===
+ plotArea = CartesianPlot("Tunneling Simulation")
+ plotArea.setType(CartesianPlot.Type.FourAxes)
+ plotArea.title().setText(f"Quantum Tunneling (E = {E:.2f}, V₀ = {V0:.2f})")
+ plotArea.setNiceExtend(True) # Auto-extend ranges for cleaner appearance
+
+ # Set fixed Y range to reduce flickering during animation
+ plotArea.enableAutoScale(CartesianCoordinateSystem.Dimension.Y, 0, False)
+ rangeY = plotArea.range(CartesianCoordinateSystem.Dimension.Y, 0)
+ rangeY.setRange(-1.0, 1.0)
+ plotArea.setRange(CartesianCoordinateSystem.Dimension.Y, 0, rangeY)
+
+ # Configure axes
+ for axis in plotArea.children(AspectType.Axis):
+ if axis.orientation() == WorksheetElement.Orientation.Horizontal and axis.position() == Axis.Position.Bottom:
+ axis.title().setText("Position x")
+ elif axis.orientation() == WorksheetElement.Orientation.Vertical and axis.position() == Axis.Position.Left:
+ axis.title().setText("Wave Function / Potential")
+
+ worksheet.addChild(plotArea)
+
+ # === Create Potential Barrier Visualization ===
+ # Use a filled area to show the barrier
+ curve_potential = XYCurve("Barrier (V/E)")
+ curve_potential.setXColumn(col_x)
+ curve_potential.setYColumn(col_potential)
+ curve_potential.setLineType(XYCurve.LineType.Line)
+ curve_potential.symbol().setStyle(Symbol.Style.NoSymbols)
+ curve_potential.line().setWidth(Worksheet.convertToSceneUnits(2, Worksheet.Unit.Point))
+
+ # Fill under the potential curve
+ from PySide6.QtCore import Qt
+ curve_potential.background().setEnabled(True)
+ curve_potential.background().setType(Background.Type.Color)
+ curve_potential.background().setPosition(Background.Position.Below)
+ curve_potential.background().setOpacity(0.3)
+
+ plotArea.addChild(curve_potential)
+
+ # === Create Wave Function Curves ===
+ curve_prob = XYCurve("|ψ|²")
+ curve_prob.setXColumn(col_x)
+ curve_prob.setYColumn(col_prob)
+ curve_prob.setLineType(XYCurve.LineType.Line)
+ curve_prob.symbol().setStyle(Symbol.Style.NoSymbols)
+ curve_prob.line().setWidth(Worksheet.convertToSceneUnits(3, Worksheet.Unit.Point))
+ plotArea.addChild(curve_prob)
+
+ curve_real = XYCurve("Re(ψ)")
+ curve_real.setXColumn(col_x)
+ curve_real.setYColumn(col_real)
+ curve_real.setLineType(XYCurve.LineType.Line)
+ curve_real.symbol().setStyle(Symbol.Style.NoSymbols)
+ curve_real.line().setWidth(Worksheet.convertToSceneUnits(1.5, Worksheet.Unit.Point))
+ curve_real.line().setOpacity(0.6)
+ plotArea.addChild(curve_real)
+
+ plotArea.addLegend()
+
+else:
+ print("Objects already exist - reusing and resetting animation...")
+
+ # Reuse existing objects
+ spreadsheet = existing_spreadsheet
+ worksheet = existing_worksheet
+
+ # Get columns
+ col_x = spreadsheet.column(0)
+ col_real = spreadsheet.column(1)
+ col_imag = spreadsheet.column(2)
+ col_prob = spreadsheet.column(3)
+ col_potential = spreadsheet.column(4)
+
+ # Get plot area (first CartesianPlot child of worksheet)
+ plotArea = None
+ for child in worksheet.children(AspectType.CartesianPlot):
+ if child.name() == "Tunneling Simulation":
+ plotArea = child
+ break
+
+ if plotArea is None:
+ print("Warning: Could not find plot area!")
+
+# Reset wave function to initial state for animation
+psi = np.exp(1j * k0 * (x - x0)) * np.exp(-(x - x0)**2 / (2 * sigma0**2))
+psi = psi / np.sqrt(np.sum(np.abs(psi)**2) * dx) # Normalize
+
+# Initialize wave function data
+col_real.replaceValues(0, [float(val.real) for val in psi])
+col_imag.replaceValues(0, [float(val.imag) for val in psi])
+col_prob.replaceValues(0, [float(abs(val)**2) for val in psi])
+
+# === Animation Loop ===
+print("Starting quantum tunneling simulation...")
+print(f"Energy: E = {E:.3f}")
+print(f"Barrier: V₀ = {V0:.3f} ({V0/E*100:.1f}% of E)")
+print(f"Barrier width: {barrier_x2 - barrier_x1:.2f}")
+print(f"Time step: {dt:.4f}, Total frames: {N_steps}")
+print()
+
+for i in range(N_steps):
+ t = i * dt
+
+ # Evolve wave function
+ psi = evolve_split_operator(psi, V_array, dt)
+
+ # Update column data
+ col_real.replaceValues(0, [float(val.real) for val in psi])
+ col_imag.replaceValues(0, [float(val.imag) for val in psi])
+ col_prob.replaceValues(0, [float(abs(val)**2) for val in psi])
+
+ # Process Qt events
+ QCoreApplication.processEvents()
+
+ # Animation delay
+ time.sleep(0.03)
+
+ # Progress and analysis
+ if (i + 1) % 20 == 0:
+ # Compute transmission and reflection
+ prob_transmitted = np.sum(np.abs(psi[x > barrier_x2])**2) * dx
+ prob_reflected = np.sum(np.abs(psi[x < barrier_x1])**2) * dx
+ prob_barrier = np.sum(np.abs(psi[(x >= barrier_x1) & (x <= barrier_x2)])**2) * dx
+
+ print(f"Frame {i+1}/{N_steps} (t={t:.2f}) | "
+ f"R={prob_reflected:.3f}, T={prob_transmitted:.3f}, "
+ f"B={prob_barrier:.3f}")
+
+print()
+print("Simulation complete!")
+print()
+print("=== Final Analysis ===")
+prob_transmitted = np.sum(np.abs(psi[x > barrier_x2])**2) * dx
+prob_reflected = np.sum(np.abs(psi[x < barrier_x1])**2) * dx
+total_prob = np.sum(np.abs(psi)**2) * dx
+print(f"Transmitted: {prob_transmitted:.4f} ({prob_transmitted/total_prob*100:.1f}%)")
+print(f"Reflected: {prob_reflected:.4f} ({prob_reflected/total_prob*100:.1f}%)")
+print(f"Total probability: {total_prob:.4f} (should be ≈1.0)")
+print()
+print("Quantum tunneling observed! The particle has a finite probability")
+print("of appearing beyond the barrier despite having insufficient classical energy.")
diff --git a/lib/examples/Statistical Comparison/CMakeLists.txt b/lib/examples/Statistical Comparison/CMakeLists.txt
new file mode 100644
index 0000000000..ea08c32b89
--- /dev/null
+++ b/lib/examples/Statistical Comparison/CMakeLists.txt
@@ -0,0 +1,3 @@
+add_executable(StatisticalComparison main.cpp)
+
+target_link_libraries(StatisticalComparison LabPlot::SDK)
diff --git a/lib/examples/Statistical Comparison/README.md b/lib/examples/Statistical Comparison/README.md
new file mode 100644
index 0000000000..c41313b03f
--- /dev/null
+++ b/lib/examples/Statistical Comparison/README.md
@@ -0,0 +1,220 @@
+# Statistical Comparison with Box Plots
+
+This tutorial demonstrates a complete statistical comparison workflow - one of the most common tasks in experimental science, clinical research, and data analysis.
+
+## What You'll Learn
+
+- Importing and organizing multiple datasets
+- Computing descriptive statistics (mean, median, std, SEM, etc.)
+- Creating publication-quality box plots
+- Visualizing individual data points with jittering
+- Generating summary statistics tables
+- Interpreting box plot elements
+
+## The Complete Workflow
+
+```
+Data Import → Statistical Analysis → Visualization → Export
+```
+
+### Step 1: Data Import
+- Load data from multiple groups/conditions
+- Can use CSV files, Excel, or generate synthetic data
+- Each group gets its own spreadsheet
+
+### Step 2: Statistical Analysis
+Computes for each group:
+- **N** - Sample size
+- **Mean ± SEM** - Average and standard error
+- **Median** - 50th percentile
+- **Std Dev** - Standard deviation (variability)
+- **IQR** - Interquartile range (25th to 75th percentile)
+- **Min/Max** - Range of values
+
+### Step 3: Visualization
+Creates a box plot showing:
+- **Box** - Interquartile range (middle 50% of data)
+- **Line in box** - Median value
+- **Diamond** - Mean value
+- **Notches** - 95% confidence interval around median
+- **Whiskers** - Extend to 1.5 × IQR
+- **Individual points** - All raw data (jittered for visibility)
+- **Outliers** - Points beyond whiskers
+- **Rug marks** - Data density at bottom
+
+### Step 4: Export
+- High-resolution figure for publication
+- Statistics table as spreadsheet
+- Summary statistics in console
+
+## Why Box Plots?
+
+Box plots are **ideal for comparing groups** because they show:
+
+1. **Central tendency** (median and mean)
+2. **Spread** (IQR, range)
+3. **Skewness** (asymmetric box)
+4. **Outliers** (unusual values)
+5. **Statistical significance** (via notches)
+
+They're much more informative than just bar charts with error bars!
+
+## Use Cases
+
+### Biological Sciences
+```python
+# Compare drug treatments
+groups = ['Vehicle', 'Drug 10μM', 'Drug 50μM', 'Drug 100μM']
+# Response: cell viability, gene expression, etc.
+```
+
+### Clinical Research
+```python
+# Compare patient groups
+groups = ['Healthy', 'Stage I', 'Stage II', 'Stage III']
+# Response: biomarker levels, symptoms scores, etc.
+```
+
+### Quality Control
+```python
+# Compare production batches
+groups = ['Batch A', 'Batch B', 'Batch C', 'Batch D']
+# Response: purity, yield, defect rate, etc.
+```
+
+### A/B Testing
+```python
+# Compare website versions
+groups = ['Control', 'Variant A', 'Variant B']
+# Response: conversion rate, time on page, etc.
+```
+
+## Interpreting the Results
+
+### Visual Assessment
+
+**Non-overlapping boxes** → Groups likely different
+```
+ Control: [----]
+Treatment A: [-----]
+```
+
+**Overlapping notches** → No significant difference
+```
+ Group 1: [---(===)---]
+ Group 2: [---(===)---]
+ Notches overlap
+```
+
+**Many outliers** → Check data quality or use robust methods
+```
+ Group X: [-----] o o oo (many outliers)
+```
+
+### Statistical Testing
+
+The box plot is **exploratory** - follow up with formal tests:
+
+**For normally distributed data:**
+- **ANOVA** (3+ groups) or **t-test** (2 groups)
+- Post-hoc: Tukey HSD, Bonferroni
+
+**For non-normal data:**
+- **Kruskal-Wallis** (3+ groups) or **Mann-Whitney** (2 groups)
+- Post-hoc: Dunn's test with correction
+
+**Check assumptions:**
+```python
+# Normality: Shapiro-Wilk test
+# Equal variances: Levene's test
+```
+
+## Customization
+
+### Import Real Data from CSV
+Replace the data generation section:
+```python
+import pandas as pd
+
+data_dict = {}
+for group in ['control', 'treatment_a', 'treatment_b']:
+ df = pd.read_csv(f"{group}.csv")
+ data_dict[group] = df['value'].values
+```
+
+### Change Box Plot Style
+
+**Show only medians (no individual points):**
+```python
+boxplot.symbolData().setStyle(Symbol.Style.NoSymbols)
+boxplot.setJitteringEnabled(False)
+```
+
+**Horizontal orientation:**
+```python
+boxplot.setOrientation(BoxPlot.Orientation.Horizontal)
+```
+
+**Colored boxes by group:**
+```python
+colors = ['#1f77b4', '#ff7f0e', '#2ca02c', '#d62728']
+for i, color in enumerate(colors):
+ boxplot.backgroundAt(i).setFirstColor(color)
+ boxplot.backgroundAt(i).setOpacity(0.3)
+```
+
+### Add Statistical Annotations
+
+```python
+# Add significance markers (*, **, ***) between groups
+# Requires implementing statistical tests first
+sig_label = TextLabel("*** p < 0.001")
+plot.addChild(sig_label)
+# Position above compared groups
+```
+
+## Best Practices
+
+### Sample Size
+- **N < 10** per group: Show all individual points
+- **N = 10-50**: Box plot with jittered points (this example)
+- **N > 50**: Box plot without individual points (too crowded)
+
+### Reporting
+Always report:
+- Sample sizes (N)
+- Central tendency (mean or median)
+- Variability (SEM or SD)
+- Statistical test used and p-values
+
+### Publication
+- Use **median ± IQR** for skewed data
+- Use **mean ± SEM** for symmetric data
+- Report **outliers** but don't automatically exclude them
+- Mention **normality tests** if using parametric statistics
+
+## Common Mistakes to Avoid
+
+❌ **Using only mean ± SD without showing distribution**
+✅ Show box plots to reveal skewness, outliers, distribution shape
+
+❌ **Excluding outliers without justification**
+✅ Report outliers, investigate causes, use robust statistics
+
+❌ **Claiming significance from overlapping notches alone**
+✅ Perform proper statistical tests (ANOVA, t-test)
+
+❌ **Bar charts for non-normal data**
+✅ Box plots show distribution shape, not just mean
+
+## Related Examples
+
+- See `../Basic Plots/script.py` for basic box plot creation
+- See `../Same Data Different Boxplots/script.py` for style variations
+- See `../Publication Multi-Panel Figure/script.py` for combining with other plots
+
+## Further Reading
+
+- Tukey, J.W. (1977). "Exploratory Data Analysis"
+- McGill et al. (1978). "Variations of Box Plots" (notched boxes)
+- Wickham, H. (2009). "ggplot2: Elegant Graphics for Data Analysis"
diff --git a/lib/examples/Statistical Comparison/main.py b/lib/examples/Statistical Comparison/main.py
new file mode 100644
index 0000000000..bc99c36caf
--- /dev/null
+++ b/lib/examples/Statistical Comparison/main.py
@@ -0,0 +1,10 @@
+#!/usr/bin/env python3
+"""Runner for statistical comparison demo"""
+import sys
+import os
+
+script_dir = os.path.dirname(os.path.abspath(__file__))
+script_path = os.path.join(script_dir, 'script.py')
+
+with open(script_path, 'r') as f:
+ exec(f.read())
diff --git a/lib/examples/Statistical Comparison/script.py b/lib/examples/Statistical Comparison/script.py
new file mode 100644
index 0000000000..230aac8eb7
--- /dev/null
+++ b/lib/examples/Statistical Comparison/script.py
@@ -0,0 +1,285 @@
+"""
+Statistical Comparison with Box Plots
+
+This script demonstrates a complete workflow for comparing multiple datasets
+with statistical analysis and visualization - a common task in experimental
+science, clinical trials, and data analysis.
+
+Workflow:
+1. Import data from multiple CSV files (or generate synthetic data)
+2. Compute descriptive statistics for each group
+3. Create box plots with individual data points
+4. Create summary statistics table
+5. Export publication-ready figure and statistics
+
+Use cases:
+- Comparing treatment groups in experiments
+- A/B testing analysis
+- Quality control across batches
+- Regional or temporal comparisons
+"""
+
+import numpy as np
+import sys
+from PySide6.QtCore import QRectF, QPointF
+from PySide6.QtGui import QFont
+from PySide6.QtWidgets import QTextEdit
+from pylabplot import *
+
+# === Generate Sample Data ===
+# In real use, import from CSV files with:
+# filter = AsciiFilter()
+# filter.readDataFromFile("group1.csv", spreadsheet)
+
+np.random.seed(42)
+
+# Simulate experimental data: 4 groups with different distributions
+groups = ['Control', 'Treatment A', 'Treatment B', 'Treatment C']
+
+# Generate realistic experimental data (e.g., cell viability %)
+data_dict = {
+ 'Control': np.random.normal(50, 8, 25),
+ 'Treatment A': np.random.normal(65, 10, 28),
+ 'Treatment B': np.random.normal(75, 9, 26),
+ 'Treatment C': np.random.normal(82, 7, 30)
+}
+
+# === Compute Statistics ===
+def compute_stats(data):
+ """Calculate descriptive statistics for a dataset"""
+ return {
+ 'n': len(data),
+ 'mean': np.mean(data),
+ 'median': np.median(data),
+ 'std': np.std(data, ddof=1), # Sample std
+ 'sem': np.std(data, ddof=1) / np.sqrt(len(data)), # Standard error
+ 'min': np.min(data),
+ 'max': np.max(data),
+ 'q25': np.percentile(data, 25),
+ 'q75': np.percentile(data, 75),
+ 'iqr': np.percentile(data, 75) - np.percentile(data, 25)
+ }
+
+stats_dict = {group: compute_stats(data) for group, data in data_dict.items()}
+
+# === Setup Project ===
+proj = project()
+
+# === Create Spreadsheets for Each Group ===
+spreadsheets = {}
+for group_name, data in data_dict.items():
+ ss = Spreadsheet(f"{group_name} Data")
+ ss.setColumnCount(1) # Single value column per group
+ proj.addChild(ss)
+
+ # Fill data column
+ ss.column(0).setName("Value")
+ ss.column(0).replaceValues(0, [float(x) for x in data])
+
+ spreadsheets[group_name] = ss
+
+# === Create Statistics Summary Spreadsheet ===
+stats_spreadsheet = Spreadsheet("Summary Statistics")
+stats_spreadsheet.setColumnCount(11) # 11 columns for all statistics
+proj.addChild(stats_spreadsheet)
+
+# Setup columns
+col_names = ['Group', 'N', 'Mean', 'Median', 'Std Dev', 'SEM', 'Min', 'Max', 'Q25', 'Q75', 'IQR']
+for i, name in enumerate(col_names):
+ stats_spreadsheet.column(i).setName(name)
+
+# Fill statistics table
+for row_idx, group_name in enumerate(groups):
+ stats = stats_dict[group_name]
+
+ # Group name (as index)
+ stats_spreadsheet.column(0).replaceValues(row_idx, [float(row_idx + 1)])
+
+ # Statistics
+ stats_values = [
+ stats['n'],
+ stats['mean'],
+ stats['median'],
+ stats['std'],
+ stats['sem'],
+ stats['min'],
+ stats['max'],
+ stats['q25'],
+ stats['q75'],
+ stats['iqr']
+ ]
+
+ for col_idx, value in enumerate(stats_values, start=1):
+ stats_spreadsheet.column(col_idx).replaceValues(row_idx, [float(value)])
+
+# === Create Worksheet ===
+worksheet = Worksheet("Statistical Comparison")
+proj.addChild(worksheet)
+
+worksheet.setUseViewSize(False)
+w = Worksheet.convertToSceneUnits(30, Worksheet.Unit.Centimeter)
+h = Worksheet.convertToSceneUnits(20, Worksheet.Unit.Centimeter)
+worksheet.setPageRect(QRectF(0, 0, w, h))
+
+worksheet.setTheme("Bright")
+
+# Vertical layout: box plot on top, statistics table below
+worksheet.setLayout(Worksheet.Layout.VerticalLayout)
+margin = Worksheet.convertToSceneUnits(0.5, Worksheet.Unit.Centimeter)
+worksheet.setLayoutTopMargin(margin)
+worksheet.setLayoutBottomMargin(margin)
+worksheet.setLayoutLeftMargin(margin)
+worksheet.setLayoutRightMargin(margin)
+worksheet.setLayoutVerticalSpacing(Worksheet.convertToSceneUnits(1.0, Worksheet.Unit.Centimeter))
+
+# === Create Box Plot ===
+plot = CartesianPlot("Box Plot Comparison")
+plot.setType(CartesianPlot.Type.FourAxes)
+plot.title().setText("Treatment Group Comparison")
+plot.setNiceExtend(True) # Auto-extend ranges for cleaner appearance
+
+worksheet.addChild(plot)
+
+# Configure axes - do this AFTER adding to worksheet
+x_axis = plot.horizontalAxis()
+x_axis.title().setText("Treatment Group")
+x_axis.majorGridLine().setStyle(Qt.PenStyle.NoPen)
+x_axis.minorGridLine().setStyle(Qt.PenStyle.NoPen)
+
+y_axis = plot.verticalAxis()
+y_axis.title().setText("Response (%)")
+y_axis.majorGridLine().setStyle(Qt.PenStyle.NoPen)
+y_axis.minorGridLine().setStyle(Qt.PenStyle.NoPen)
+
+# Create box plot with all groups
+boxplot = BoxPlot("Group Comparison")
+
+# Set data columns (one per group)
+data_columns = [spreadsheets[group].column(0) for group in groups]
+boxplot.setDataColumns(data_columns)
+
+# Customize box plot appearance
+boxplot.setOrientation(BoxPlot.Orientation.Vertical)
+boxplot.setWhiskersType(BoxPlot.WhiskersType.IQR)
+boxplot.setWhiskersRangeParameter(1.5) # Standard 1.5×IQR rule
+boxplot.setNotchesEnabled(True) # Show confidence interval around median
+boxplot.setVariableWidth(False) # Equal width boxes
+
+# Show individual data points with jittering
+boxplot.symbolData().setStyle(Symbol.Style.Circle)
+boxplot.symbolData().setSize(Worksheet.convertToSceneUnits(3, Worksheet.Unit.Point))
+boxplot.symbolData().setOpacity(0.4)
+boxplot.setJitteringEnabled(True)
+
+# Show mean as well as median
+boxplot.symbolMean().setStyle(Symbol.Style.Diamond)
+boxplot.symbolMean().setSize(Worksheet.convertToSceneUnits(6, Worksheet.Unit.Point))
+
+# Highlight outliers
+boxplot.symbolOutlier().setStyle(Symbol.Style.Circle)
+boxplot.symbolOutlier().setSize(Worksheet.convertToSceneUnits(4, Worksheet.Unit.Point))
+
+# Enable rug plot for better data distribution visualization
+boxplot.setRugEnabled(True)
+boxplot.setRugLength(Worksheet.convertToSceneUnits(8, Worksheet.Unit.Point))
+
+plot.addChild(boxplot)
+boxplot.recalc()
+
+# Add legend
+plot.addLegend()
+
+# === Create Statistics Table as Text ===
+# Create formatted statistics table as child of plot (bottom right corner)
+table_label = TextLabel("Statistics Table")
+plot.addChild(table_label)
+
+# Build HTML table
+# font-family: Arial; font-size: 6pt;
+html_table = "<html><body><table border='1' cellpadding='3' cellspacing='0' style='border-collapse: collapse;'>"
+html_table += "<tr style='background-color: #e0e0e0; font-weight: bold;'>"
+html_table += "<th>Group</th><th>N</th><th>Mean ± SEM</th><th>Median</th><th>Std Dev</th><th>Range</th></tr>"
+
+for group_name in groups:
+ stats = stats_dict[group_name]
+ html_table += "<tr>"
+ html_table += f"<td>{group_name}</td>"
+ html_table += f"<td>{stats['n']}</td>"
+ html_table += f"<td>{stats['mean']:.2f} ± {stats['sem']:.2f}</td>"
+ html_table += f"<td>{stats['median']:.2f}</td>"
+ html_table += f"<td>{stats['std']:.2f}</td>"
+ html_table += f"<td>{stats['min']:.1f} - {stats['max']:.1f}</td>"
+ html_table += "</tr>"
+
+html_table += "</table></body></html>"
+
+table_label.setText(html_table)
+
+# Position in bottom right corner of plot using logical coordinates
+table_label.setCoordinateBindingEnabled(True)
+# Get plot data range to position in bottom right
+rangeX = plot.range(CartesianCoordinateSystem.Dimension.X, 0)
+rangeY = plot.range(CartesianCoordinateSystem.Dimension.Y, 0)
+# Position at 95% of X range and 5% of Y range (bottom right)
+table_label.setPositionLogical(QPointF(rangeX.end() * 0.95, rangeY.start() + (rangeY.end() - rangeY.start()) * 0.05))
+table_label.setHorizontalAlignment(WorksheetElement.HorizontalAlignment.Right)
+table_label.setVerticalAlignment(WorksheetElement.VerticalAlignment.Bottom)
+
+# === Print Statistics Summary ===
+print("=" * 80)
+print("STATISTICAL COMPARISON ANALYSIS")
+print("=" * 80)
+print()
+
+for group_name in groups:
+ stats = stats_dict[group_name]
+ print(f"{group_name}:")
+ print(f" N = {stats['n']}")
+ print(f" Mean ± SEM = {stats['mean']:.2f} ± {stats['sem']:.2f}")
+ print(f" Median (IQR) = {stats['median']:.2f} ({stats['q25']:.2f} - {stats['q75']:.2f})")
+ print(f" Std Dev = {stats['std']:.2f}")
+ print(f" Range = {stats['min']:.2f} - {stats['max']:.2f}")
+ print()
+
+print("=" * 80)
+print("INTERPRETATION GUIDE")
+print("=" * 80)
+print()
+print("Box Plot Elements:")
+print(" • Box = Interquartile range (IQR, 25th to 75th percentile)")
+print(" • Line in box = Median")
+print(" • Diamond symbol = Mean")
+print(" • Notches = 95% confidence interval around median")
+print(" • Whiskers = 1.5 × IQR (standard Tukey method)")
+print(" • Individual points = Raw data with jitter for visibility")
+print(" • Circles beyond whiskers = Outliers")
+print(" • Rug marks at bottom = Data density")
+print()
+print("Statistical Notes:")
+print(" • Use SEM (Standard Error of Mean) for comparing means")
+print(" • Use Std Dev for describing variability within groups")
+print(" • Notched boxes: non-overlapping notches suggest significant difference")
+print(" • For formal testing, use ANOVA or Kruskal-Wallis test")
+print()
+print("=" * 80)
+print("NEXT STEPS")
+print("=" * 80)
+print()
+print("1. Visual Assessment:")
+print(" - Do the boxes overlap? (suggests similarity)")
+print(" - Do the notches overlap? (suggests no significant difference)")
+print(" - Are outliers present? (check data quality)")
+print()
+print("2. Statistical Testing:")
+print(" - Run ANOVA for normally distributed data")
+print(" - Use Kruskal-Wallis for non-normal data")
+print(" - Perform post-hoc tests (Tukey HSD, Dunn's test)")
+print()
+print("3. Export:")
+print(" - Save figure: worksheet.exportToFile('comparison.pdf', ...)")
+print(" - Save statistics: export stats_spreadsheet to CSV")
+print()
+
+# Optional: Export
+# worksheet.exportToFile("statistical_comparison.pdf", Worksheet.ExportFormat.PDF)
+# worksheet.exportToFile("statistical_comparison.png", Worksheet.ExportFormat.PNG, 300)