Metadata-Version: 2.4
Name: ramansep
Version: 0.13.1
Summary: Separate strain from carrier density in 2D-material Raman maps using a two-mode linear inversion
Author-email: "Tanvir M. Mahim" <tanvir.mahim@bracu.ac.bd>
License-Expression: Apache-2.0
Project-URL: Homepage, https://tanvir-mahmud-mahim.github.io/software/
Project-URL: Repository, https://github.com/TaN-MM-Org/ramansep
Project-URL: Issues, https://github.com/TaN-MM-Org/ramansep/issues
Project-URL: Changelog, https://github.com/TaN-MM-Org/ramansep/releases
Keywords: Raman spectroscopy,2D materials,MoS2,strain,carrier density,TMD
Classifier: Development Status :: 3 - Alpha
Classifier: Intended Audience :: Science/Research
Classifier: Programming Language :: Python :: 3
Classifier: Topic :: Scientific/Engineering :: Physics
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.22
Requires-Dist: scipy>=1.8
Provides-Extra: plot
Requires-Dist: matplotlib; extra == "plot"
Provides-Extra: test
Requires-Dist: pytest; extra == "test"
Dynamic: license-file

# ramansep

[![PyPI](https://img.shields.io/pypi/v/ramansep)](https://pypi.org/project/ramansep/) [![DOI](https://img.shields.io/badge/DOI-10.5281%2Fzenodo.22014913-blue)](https://doi.org/10.5281/zenodo.22014913) [![tests](https://github.com/TaN-MM-Org/ramansep/actions/workflows/ci.yml/badge.svg)](https://github.com/TaN-MM-Org/ramansep/actions)

`ramansep` is a Python package that tells **strain** apart from
**charge** in Raman maps of 2D materials (sheets only one or a few atoms
thick, such as graphene or monolayer MoS2).

A Raman peak moves when the material is stretched, and it also moves
when extra electrons or holes are added. One peak alone cannot say
which of the two happened. Two peaks that respond differently can: if
you know how much each peak moves per unit of strain and per unit of
charge, the two measured shifts give two equations with two unknowns.
`ramansep` solves those equations at every pixel of a map and gives a
strain map and a charge-density map, each with an error bar. It also
covers the steps around that:

- How do I read the text files my instrument exports, and get the
  peak shifts out of the measured spectra? Which pixels failed or
  carry no peak at all?
- What if two peaks overlap?
- Which pair of peaks separates the two causes best?
- With three or more peaks, does "strain plus charge" actually explain
  what I measured, pixel by pixel?
- Can I trade some resolution for smaller error bars if the fields are
  smooth, and fill the pixels whose fits failed?
- How do I measure the response numbers on my own instrument, even
  when the starting strain and charge of my sample are unknown, and
  how should I plan that measurement?
- How much of my final error bar comes from that calibration, with two
  peaks or more?
- What if laser heating also shifts the peaks?

The results are checked by automated tests against independent
calculations and exact limits (see
[How the results are checked](#how-the-results-are-checked)). When an
input cannot give a trustworthy answer, the package stops with an error
that says why, rather than returning a number that looks fine.

## Contents

- [A short guide to the words used here](#a-short-guide-to-the-words-used-here)
- [Install, requirements and units](#install-requirements-and-units)
- [Examples](#examples) (each with the output it prints)
- [What is in the package](#what-is-in-the-package)
- [Cited coefficient sets](#cited-coefficient-sets)
- [When it refuses, and why](#when-it-refuses-and-why)
- [How the results are checked](#how-the-results-are-checked)
- [Corrections in earlier versions](#corrections-in-earlier-versions)
- [Limits, and what is deliberately not included](#limits-and-what-is-deliberately-not-included)
- [Where it comes from](#where-it-comes-from)
- [Citing, support and license](#citing-support-and-license)

## A short guide to the words used here

- **Raman spectrum** -- the light scattered by a material, plotted
  against its energy loss in **cm^-1** (wavenumbers). Each lattice
  vibration (**phonon**, or **mode**) gives a peak.
- **Mode names** -- labels such as A'1, E', 2LA(M) (MoS2) or G, 2D
  (graphene). They only name which peak is meant.
- **Shift** (`dw`) -- how far a peak has moved from its position in
  the undisturbed ("pristine") material, in cm^-1. Negative means the
  peak moved to lower wavenumber (the mode "softened").
- **Strain** -- stretching of the sheet, in percent. **Biaxial** means
  stretched equally in both directions; **uniaxial** means along one
  direction. Positive is tension.
- **Carrier density** (or **doping**, **charge**) -- extra electrons
  or holes (missing electrons) per area, in the MoS2 sets in units of
  1e13 cm^-2.
- **Gate, gated** -- an electrode that adds a known amount of charge to
  the sheet when a voltage is applied; a **strain stage** stretches the
  sample by a known amount. Both are used to measure how far a peak
  moves per unit of charge or strain.
- **Lever arm** -- how much one peak shifts per unit of one cause, for
  example -5.1 cm^-1 per percent strain. The lever arms of all the
  peaks form the **lever-arm matrix** `K` (one row per peak, one column
  per cause). The package assumes the shifts are linear in the causes:
  `shifts = K @ (strain, density)`.
- **Condition number** -- a single number that says how close the
  matrix is to being unsolvable. Large values mean the peaks respond
  too similarly, so small shift errors become large strain or density
  errors.
- **Sigma / error bar** -- one standard deviation of a value.
  **Correlation** between two estimates near +1 or -1 means they are
  hard to tell apart. The **covariance** matrix collects the squared
  error bars (on its diagonal) and how each pair of estimates varies
  together (off the diagonal).
- **Weighted least squares** -- the standard way to solve more
  equations than unknowns: find the values that make the weighted sum
  of squared mismatches smallest, each measurement weighted by
  1/sigma^2 so that more precise peaks count more.
- **Chi-square and p-value** -- with more peaks than causes, the
  equations cannot all be met exactly when something else is going on.
  Chi-square measures the mismatch, weighted by the error bars; the
  p-value is the probability of a mismatch at least that large if
  strain and charge alone were at work. A tiny p-value flags a pixel.
  The **degrees of freedom** (`dof`) of this check are the number of
  peaks minus the number of causes solved for.
- **Hyperspectral map (cube)** -- one spectrum per pixel, stored as an
  array of shape (rows H, columns W, spectral points L).
- **Lorentzian / Voigt** -- peak shapes. A Lorentzian is the natural
  shape of a phonon peak; a Voigt is a Lorentzian blurred by a
  Gaussian, for peaks broadened by the instrument.
- **FWHM** -- full width at half maximum, the width of a peak measured
  halfway up its height.
- **Baseline** -- the background under a peak. A sloped background
  (for example from fluorescence) pulls a fitted peak centre sideways
  unless the slope is fitted too.
- **Stokes / anti-Stokes** -- Raman peaks on the energy-loss side and
  the energy-gain side. The ratio of their heights depends on
  temperature, so it works as a thermometer.
- **DFT (density functional theory)** -- a standard quantum-mechanical
  computer calculation of a material. "Frozen-phonon" DFT computes a
  vibration's frequency by displacing the atoms along that vibration.
  The source paper used it to compute the MoS2 strain lever arms.
- **Grueneisen parameter** (gamma) -- a number that says how strongly
  a vibration's frequency w changes when the sheet is stretched; for
  biaxial strain the package uses gamma = -(1/(2w)) dw/d(strain), with
  strain as a fraction (not percent).

## Install, requirements and units

```
pip install ramansep
```

For development: clone the repository and `pip install -e .[test]`.

It needs Python 3.9 or newer, NumPy 1.22 or newer and SciPy 1.8 or
newer. `matplotlib` is optional (`pip install ramansep[plot]`) and is
used only by the demo script `examples/synthetic_map.py`.

Units and signs:

- Shifts are in cm^-1, relative to the pristine reference frequency of
  each mode (you supply those references, e.g. `ref1`, `ref2`).
- The inversion returns strain and density in whatever units your
  lever arms use. The shipped MoS2 sets use percent biaxial strain
  (positive = tension) and electron density in 1e13 cm^-2. The
  graphene set uses percent uniaxial strain and, on purpose, a doping
  coordinate in cm^-1 (see [Cited coefficient sets](#cited-coefficient-sets)).
- Every `ModeCoefficients` must carry a `reference` string that says
  where its numbers came from; it is a required field.

## Examples

Each example below runs as written, and the output shown is what it
printed with ramansep 0.13.0. Unless a value is said to come from a
cited source, it is an illustrative value chosen for the example. The
repository also has a longer demo, `examples/synthetic_map.py`, that
inverts a noisy synthetic 128 x 128 map and plots it with matplotlib.

### 1. Strain and charge at two spots of a MoS2 ribbon

The shifts here are the ones quoted by the source paper (Mahim and
Rahman, under review) from tip-enhanced Raman maps (maps taken with a
sharp metal tip that gives nanometre resolution; Krayev et al., Appl.
Phys. Lett. 128, 203102 (2026)); the 0.1 cm^-1 shift error is
illustrative.

```python
import numpy as np
import ramansep as rs

model = rs.SeparationModel(rs.mos2_a1_2la())    # monolayer MoS2, A'1 + 2LA(M)

# Peak shifts (cm^-1) quoted by the source paper at two spots:
# the ribbon edge (A'1 -0.5, 2LA(M) 0.0) and an interior spot (-0.6, -2.8).
dw_a1 = np.array([-0.5, -0.6])
dw_2la = np.array([0.0, -2.8])
res = model.invert(dw_a1, dw_2la, sigma1=0.1, sigma2=0.1)   # 0.1 cm^-1: illustrative

for name, i in [("edge", 0), ("interior", 1)]:
    print(f"{name:8s} strain {res.strain[i]:+.3f} +/- {res.strain_sigma[i]:.3f} %   "
          f"electrons {res.density[i] * 1e13:+.2e} +/- {res.density_sigma[i] * 1e13:.1e} cm^-2")
print(f"condition number of the matrix: {res.condition_number:.1f}")
```

```
edge     strain +0.000 +/- 0.005 %   electrons +2.27e+12 +/- 4.7e+11 cm^-2
interior strain +0.134 +/- 0.005 %   electrons -3.78e+11 +/- 4.7e+11 cm^-2
condition number of the matrix: 10.1
```

The edge carries about 2.3e12 electrons per cm^2 and no strain; the
interior spot is under 0.134 % tension, and its density is within one
error bar of zero. These are the paper's published separations; the
test suite checks them (see below). `res.correlation` gives the
per-pixel correlation of the two estimates.

### 2. Choosing a mode pair

```python
import numpy as np
import ramansep as rs

sig = 0.1                                        # illustrative shift error, cm^-1
paper = rs.SeparationModel(rs.mos2_a1_2la()).invert(0.0, 0.0, sig, sig)
optical = rs.SeparationModel(rs.mos2_eprime_a1()).invert(0.0, 0.0, sig, sig)
print(f"A'1 + 2LA(M): strain error {float(paper.strain_sigma):.4f} %")
print(f"E'  + A'1   : strain error {float(optical.strain_sigma):.4f} %")
print(f"ratio: {float(optical.strain_sigma / paper.strain_sigma):.1f}")
```

```
A'1 + 2LA(M): strain error 0.0048 %
E'  + A'1   : strain error 0.0257 %
ratio: 5.4
```

Both optical modes E' and A'1 move only a little with strain, so with
the same shift errors that pair gives a strain error about five times
larger than the A'1 + 2LA(M) pair. The error bars depend only on the
lever arms and the shift errors, so this comparison can be made before
measuring anything.

### 3. From a map file to strain and charge maps

```python
import os, tempfile
import numpy as np
import ramansep as rs
from ramansep import lorentzian

coeffs = rs.mos2_a1_2la()
model = rs.SeparationModel(coeffs)

# An illustrative 4 x 5 map: known strain (%) and electron density (1e13 cm^-2).
rng = np.random.default_rng(3)
strain = rng.uniform(-0.1, 0.1, (4, 5))
density = rng.uniform(0.0, 0.2, (4, 5))
dw1, dw2 = model.forward(strain, density)

# Make a spectrum for every pixel: two Lorentzian peaks plus noise.
x = np.linspace(380.0, 480.0, 401)
cube = np.empty((4, 5, x.size))
for i in range(4):
    for j in range(5):
        cube[i, j] = (lorentzian(x, 404.7 + dw1[i, j], 3.0, 800.0, 50.0)
                      + lorentzian(x, 452.0 + dw2[i, j], 4.0, 500.0, 0.0)
                      + rng.normal(0.0, 3.0, x.size))
cube[3, 4, 100] = np.nan                         # one broken pixel

path = os.path.join(tempfile.mkdtemp(), "map.csv")
rs.save_map_csv(path, x, cube)                   # the documented file format

x2, cube2 = rs.load_map_csv(path)
fits = rs.fit_map(x2, cube2, window1=(395, 415), window2=(440, 465),
                  ref1=404.7, ref2=452.0)
print("masked pixels:", fits.n_masked, "at", np.argwhere(~fits.ok).tolist())

inv = model.invert(fits.dw1, fits.dw2, fits.sigma1, fits.sigma2)
ok = fits.ok
z_s = np.abs(inv.strain - strain)[ok] / inv.strain_sigma[ok]
z_d = np.abs(inv.density - density)[ok] / inv.density_sigma[ok]
print(f"typical error bars: strain {np.median(inv.strain_sigma[ok]):.4f} %, "
      f"density {np.median(inv.density_sigma[ok]) * 1e13:.1e} cm^-2")
print(f"pixels whose error is within 3 error bars: strain {(z_s < 3).sum()} of {ok.sum()}, "
      f"density {(z_d < 3).sum()} of {ok.sum()}")
print("strain at the broken pixel:", inv.strain[3, 4])
```

```
masked pixels: 1 at [[3, 4]]
typical error bars: strain 0.0002 %, density 1.3e+10 cm^-2
pixels whose error is within 3 error bars: strain 19 of 19, density 19 of 19
strain at the broken pixel: nan
```

The map file is plain CSV with the header `row,col,wavenumber_cm1,counts`
(0-based pixel indices, one shared wavenumber axis). A single spectrum
uses `wavenumber_cm1,counts` (`save_spectrum_csv`, `load_spectrum_csv`).
To read the plain text tables that instruments export, see example 16.
`fit_map` fits both peaks at every pixel. A pixel is masked (NaN in all
four output maps, `False` in `fits.ok`) when its counts are not finite,
a fit stops with an error or does not converge, a fitted centre ends
up outside its window, or a fitted peak height is zero or negative (a
dip is not a Raman peak; this check is new in 0.12.0). `fits.reason`
records, per pixel, why it was masked (`rs.MASK_REASONS` spells out
the codes). Maps with pixels that carry no peak at all need one more
check, `min_snr` (example 12). Mistakes in the settings themselves
(window limits, `baseline` name, a wavenumber axis containing NaN) are
refused with an error before any pixel is fitted.
The per-pixel inversion keeps masked pixels as NaN and does not touch
their neighbours. For spectra on a sloped background, pass
`baseline="linear"`; for peaks close enough to overlap, pass
`joint=True` (example 13).

### 4. Three peaks, and a check that the model holds

```python
import numpy as np
import ramansep as rs

# Lever arms (cm^-1 per unit strain, per unit density) of three modes.
# Rows 1-2 are the MoS2 A'1 and 2LA(M) values; row 3 is an invented,
# illustrative third mode.
K = np.array([[-5.1, -2.2],
              [-20.9, 0.0],
              [-8.0, -1.1]])
sig = [0.15, 0.10, 0.12]                         # illustrative shift errors, cm^-1
mm = rs.MultiModeModel(K, mode_names=["A'1", "2LA(M)", "X"])

rng = np.random.default_rng(3)
strain = rng.normal(0, 0.2, (30, 30))
density = rng.normal(0, 1.0, (30, 30))
shifts = [s + rng.normal(0, e, s.shape) for s, e in zip(mm.forward(strain, density), sig)]
shifts[2][10:15, :] += 2.0                       # rows 10-14: a cause the model does not know

res = mm.invert(shifts, sigmas=sig)
print("degrees of freedom of the check:", res.dof)
print(f"median p-value, rows 10-14 : {np.median(res.p_value[10:15]):.1e}")
print(f"median p-value, other rows : {np.median(res.p_value[:10]):.2f}")

for r in rs.compare_mode_sets(K, sig, mode_names=["A'1", "2LA(M)", "X"]):
    print(r["names"], f"var_strain {r['var_strain']:.2e}  var_density {r['var_density']:.2e}")
print(f"all three     var_strain {res.strain_sigma[0, 0]**2:.2e}  var_density {res.density_sigma[0, 0]**2:.2e}")
```

```
degrees of freedom of the check: 1
median p-value, rows 10-14 : 1.0e-43
median p-value, other rows : 0.50
["A'1", '2LA(M)'] var_strain 2.29e-05  var_density 4.77e-03
['2LA(M)', 'X'] var_strain 2.29e-05  var_density 1.31e-02
["A'1", 'X'] var_strain 6.74e-04  var_density 1.26e-02
all three     var_strain 2.21e-05  var_density 3.65e-03
```

With more peaks than causes, `MultiModeModel` uses weighted least
squares (each peak weighted by 1/sigma^2). The extra peak shrinks the
error bars and makes the chi-square check possible: the rows with the
hidden extra shift stand out with p-values near zero, while elsewhere
the p-values scatter around 0.5 as they should. With exactly two peaks
there is no check (`dof` is 0, `chi2_map` and `p_value` are `None`).
`compare_mode_sets` ranks every subset of the modes (pairs by default)
that can separate the two causes by the sum of the two variances
(squared error bars), smallest first.

### 5. Smoothing a whole map

```python
import numpy as np
import ramansep as rs

K = np.array([[-2.3, 1.1], [-0.9, -1.7], [0.4, 2.2]])   # illustrative, non-physical
sig = np.array([0.03, 0.05, 0.02])                      # cm^-1

# A smooth illustrative truth on a 12 x 12 map, measured with noise.
yy, xx = np.mgrid[0:12, 0:12]
strain = 0.002 * xx                                     # gentle gradient
density = 0.01 + 0.001 * yy
rng = np.random.default_rng(7)
shifts = np.einsum("mj,jhw->mhw", K, np.stack([strain, density])) \
    + rng.normal(size=(3, 12, 12)) * sig[:, None, None]

for lam in (0.0, 1000.0, 10000.0):
    r = rs.bayesian_map_inversion(K, shifts, sig, lam_strain=lam, posterior_sigma=True)
    err = np.sqrt(np.mean((r.strain - strain) ** 2))
    print(f"lam = {lam:6.0f}: mean strain sigma {r.strain_sigma.mean():.4f}, "
          f"RMS strain error {err:.4f}")

shifts[0, 2, 3] = np.nan                                # a masked pixel
try:
    rs.bayesian_map_inversion(K, shifts, sig, lam_strain=100.0)
except ValueError as err:
    print("refused:", str(err)[:60], "...")
```

```
lam =      0: mean strain sigma 0.0123, RMS strain error 0.0115
lam =   1000: mean strain sigma 0.0101, RMS strain error 0.0078
lam =  10000: mean strain sigma 0.0057, RMS strain error 0.0039
refused: shifts contain non-finite values; the smoothness prior coupl ...
```

`bayesian_map_inversion` solves the whole map at once, adding a
penalty `lam` times the squared differences between neighbouring
pixels (a "smoothness prior"; `lam_strain` and `lam_density` can
differ). It is one sparse linear system, solved exactly. At `lam = 0`
it gives the same answer as the pixel-by-pixel inversion, and a larger
`lam` never increases the error bars. Smoothing helps only when the
real fields are smooth, as they are here by construction; you choose
`lam`, and the package does not estimate it for you. Masked (NaN)
pixels are refused by default because the prior links every pixel to
its neighbours; `fill_masked=True` fills them from their neighbours
instead, and `sigmas` may be given per pixel (example 14). Exact
per-pixel error bars (`posterior_sigma=True`) need a dense matrix
inverse and are refused above `max_dense` = 4096 unknowns (two per
pixel) unless you raise that limit.

### 6. Measuring the lever arms on your own instrument

```python
import numpy as np
import ramansep as rs

# Reference states you prepared: strain (%) from a strain stage,
# electron density (1e13 cm^-2) from a gate. Illustrative values.
strain = np.array([0.0, 0.4, 0.8, 0.0, 0.3, 0.6, 0.2, 0.0])
density = np.array([0.0, 0.0, 0.0, 0.9, 0.6, 0.3, 1.2, 0.5])

# Pretend the true lever arms are the MoS2 A'1 / 2LA(M) values and
# measure the shifts with 0.1 cm^-1 noise.
K_true = np.array([[-5.1, -2.2], [-20.9, 0.0]])
rng = np.random.default_rng(1)
shifts = np.column_stack([strain, density]) @ K_true.T + rng.normal(0, 0.1, (8, 2))

cal = rs.calibrate_lever_arms(strain, density, shifts, sigmas=0.1,
                              mode_names=["A'1", "2LA(M)"])
np.set_printoptions(precision=2, suppress=True)
print("fitted lever arms:\n", cal.K)
print("their error bars:\n", cal.K_sigma)
print("chi2 / dof per mode:", cal.chi2 / cal.chi2_dof)

model = rs.SeparationModel(cal.coefficients(reference="stage + gate calibration, illustrative"))
print(model.coefficients.reference)

try:                                             # a strain-only sweep
    rs.calibrate_lever_arms(strain[:3], density[:3], shifts[:3], sigmas=0.1)
except ValueError as err:
    print("refused:", err)
```

```
fitted lever arms:
 [[ -5.    -2.27]
 [-20.9    0.03]]
their error bars:
 [[0.09 0.06]
 [0.09 0.06]]
chi2 / dof per mode: [0.09 0.56]
stage + gate calibration, illustrative
refused: reference states are collinear in the (strain, density) plane, so the strain and density lever arms are not separately identifiable; add a reference point off that line (e.g. a gated point to a strain-only sweep)
```

`calibrate_lever_arms` fits each mode's two lever arms by weighted
least squares from reference states of known strain and density. With
`sigmas` given, the error bars are the exact result for that noise
level, and `chi2 / dof` near 1 says the straight-line model and your
stated errors agree. Without `sigmas`, the error bars are estimated
from the scatter of the fit, which needs at least 3 reference points.
A sweep that changes strain only (or strain and charge in fixed
proportion) cannot tell the two lever arms apart, so it is refused.
This fit assumes the shifts are measured from the true pristine
frequencies and that strain and density are known in absolute terms.
If they are known only relative to the starting state of the sample,
use `fit_offset=True` (example 17).
`cal.K` of a calibration with three or more modes goes straight into
`MultiModeModel`.

### 7. Planning the calibration before measuring it

```python
import numpy as np
from ramansep import plan_calibration, design_references, repeats_for_sigma

# Four planned reference states, 0.15 cm^-1 expected shift error (illustrative).
plan = plan_calibration(strain=[0.0, 0.3, 0.6, 0.0],
                        density=[0.0, 0.0, 0.1, 0.8], sigmas=0.15)
print("identifiable:", plan["identifiable"])
print("expected error bars (strain arm, density arm):", plan["K_sigma"].round(3))

# States the stage and gate can reach; which 4 are worth measuring?
eps = np.array([0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.0, 0.0, 0.25])
rho = np.array([0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.4, 0.8, 0.8])
pick = design_references(eps, rho, n_pick=4, sigmas=0.15)
print("measure candidates:", pick["indices"])

r, plan_r = repeats_for_sigma(0.1, eps[pick["indices"]], rho[pick["indices"]], sigmas=0.15)
print(f"repeat that set {r} times -> worst error bar {plan_r['K_sigma'].max():.3f}")

line = plan_calibration([0.1, 0.2, 0.4], [0.0, 0.0, 0.0], sigmas=0.15)
print("strain-only sweep identifiable:", line["identifiable"])
```

```
identifiable: True
expected error bars (strain arm, density arm): [0.225 0.187]
measure candidates: [8, 5, 7, 4]
repeat that set 6 times -> worst error bar 0.092
strain-only sweep identifiable: False
```

Because the calibration is a linear fit, the error bars it will
report depend only on which reference states you prepare and on the
shift errors, not on the measured shifts. `plan_calibration` computes
them in advance. `design_references` picks states one at a time, each
time taking the one that shrinks the joint uncertainty most (the
"D-optimal" rule; F. Pukelsheim, Optimal Design of Experiments, SIAM
(2006)); this is a good heuristic, not a guarantee of the best
possible subset. `repeats_for_sigma` uses the fact that repeating the
whole set `r` times divides the covariance by `r`. `design_references`
and `repeats_for_sigma` refuse a candidate set that lies on one line.

### 8. How much of the error bar comes from the calibration

```python
import numpy as np
from ramansep import calibrate_lever_arms, separation_with_calibration

# Five reference states (illustrative) and noise-free shifts from the
# MoS2 A'1 / 2LA(M) lever arms; the calibration is told they carry 0.2 cm^-1 errors.
K_true = np.array([[-5.1, -2.2], [-20.9, 0.0]])
eps = np.array([0.0, 0.4, 0.8, 0.0, 0.3])
rho = np.array([0.0, 0.0, 0.2, 0.9, 0.6])
cal = calibrate_lever_arms(eps, rho, np.column_stack([eps, rho]) @ K_true.T, sigmas=0.2)

out = separation_with_calibration(cal, dw1=-2.0, dw2=-8.0, sigma1=0.1, sigma2=0.1)
print(f"strain  {float(out['strain']):.4f} %")
print(f"  error from the spectra     {float(out['strain_sigma_shifts']):.4f}")
print(f"  error from the calibration {float(out['strain_sigma_calibration']):.4f}")
print(f"  total                      {float(out['strain_sigma']):.4f}")
print(f"density {float(out['density']):.4f} x 1e13 cm^-2")
print(f"  error from the spectra     {float(out['density_sigma_shifts']):.4f}")
print(f"  error from the calibration {float(out['density_sigma_calibration']):.4f}")
print(f"  total                      {float(out['density_sigma']):.4f}")
```

```
strain  0.3828 %
  error from the spectra     0.0048
  error from the calibration 0.0040
  total                      0.0063
density 0.0217 x 1e13 cm^-2
  error from the spectra     0.0468
  error from the calibration 0.0396
  total                      0.0613
```

Other inversions treat the lever arms as exact.
`separation_with_calibration` also carries the calibration's own
uncertainty into the result, and reports the two parts separately; the
total is their quadrature sum (square root of the sum of squares). It
uses the first-order (linearised) error propagation, based on the
identity dx/dK_mj = -K^-1 E_mj x for the two-mode solution x = K^-1 y
(E_mj is the matrix with a single 1 at row m, column j). It takes
two-mode calibrations only; for three or more modes use
`multimode_with_calibration` (example 15). The calibration part comes
from one set of lever arms shared by every pixel, so it is an error
common to the whole map: unlike the spectral noise, it does not shrink
when you average many pixels.

### 9. Temperature from the anti-Stokes / Stokes ratio

```python
from ramansep import (HC_OVER_KB_CM_K, anti_stokes_ratio,
                      calibrate_anti_stokes, temperature_from_anti_stokes)

omega = 385.0                  # the mode's Raman shift, cm^-1
print(f"hc/kB = {HC_OVER_KB_CM_K:.7f} cm K")

# Calibrate once, on a spot you know is at 295 K (illustrative ratio).
C = calibrate_anti_stokes(ratio_measured=0.13, omega_cm1=omega, temperature_K=295.0)
print(f"calibration constant C = {C:.4f}")

# Later, under stronger laser power, the ratio has risen.
out = temperature_from_anti_stokes(0.17, omega, C, ratio_sigma=0.005)
print(f"T = {out['temperature']:.1f} +/- {out['temperature_sigma']:.1f} K")
print(f"with C wrongly set to 1: T = {temperature_from_anti_stokes(0.17, omega, 1.0):.1f} K")

try:
    temperature_from_anti_stokes(1.2 * C, omega, C)
except ValueError as err:
    print("refused:", str(err)[:60], "...")
```

```
hc/kB = 1.4387769 cm K
calibration constant C = 0.8500
T = 344.2 +/- 6.3 K
with C wrongly set to 1: T = 312.6 K
refused: ratio >= calibration constant: the Bose factor is strictly < ...
```

The model is `I_AS / I_S = C exp(-hc omega / (kB T))`, where `omega`
is the mode's Raman shift. The constant `C` absorbs the instrument's
response and the scattering factors; it is close to, but not exactly,
1, and assuming 1 gives a wrong temperature (as the last printed line
shows). So `C` is measured once at a known temperature. The constant
`hc/kB` is computed from the exact SI values of h, c and kB, not typed
in by hand. A ratio at or above `C` has no finite positive temperature
and is refused; it points to a wrong calibration or background
subtraction. `anti_stokes_ratio` is the forward model.

### 10. Strain, charge and temperature together

```python
import numpy as np
from ramansep import ThreeCauseModel

# Four modes; columns: shift per % strain, per 1e13 cm^-2, per kelvin.
# Invented, illustrative values -- no temperature coefficients ship with the package.
K = np.array([[-2.1, -0.8, -0.011],
              [-1.0, -2.3, -0.016],
              [-3.2, -0.4, -0.007],
              [-0.6, -1.1, -0.021]])
model = ThreeCauseModel(K)
shifts = model.forward(strain=0.2, density=0.5, temperature=40.0)   # 40 K of heating
res = model.invert(shifts, sigmas=[0.05, 0.05, 0.05, 0.05])
print(f"strain {float(res.strain):.3f} +/- {float(res.strain_sigma):.3f} %")
print(f"density {float(res.density):.3f} +/- {float(res.density_sigma):.3f} x 1e13 cm^-2")
print(f"temperature {float(res.temperature):.1f} +/- {float(res.temperature_sigma):.1f} K")
```

```
strain 0.200 +/- 0.016 %
density 0.500 +/- 0.042 x 1e13 cm^-2
temperature 40.0 +/- 4.2 K
```

Heating shifts peaks too, and a two-cause analysis would count that
shift as strain or charge. `ThreeCauseModel` adds a third column of
lever arms, the temperature coefficients (cm^-1 per kelvin) that you
have calibrated for your material, and needs at least three modes.
Because shifts are measured from the reference frequencies, the
"temperature" it returns is the change from the temperature at which
those references hold. With four or more modes it also returns the
chi-square check (`chi2_map`, `p_value`) and the full 3 x 3 covariance
per pixel. A pixel with a NaN shift or error bar (a pixel masked by
`fit_map`) comes out as NaN without affecting the others. Measuring the temperature with example 9 and comparing it
with this estimate is a direct cross-check.

### 11. Graphene: the G + 2D decomposition

```python
import warnings
import ramansep as rs

coeffs = rs.graphene_g_2d_lee2012()
with warnings.catch_warnings(record=True) as caught:
    warnings.simplefilter("always")
    rs.SeparationModel(coeffs)                   # default condition_warn = 30
print("warning:", str(caught[0].message)[:62], "...")

model = rs.SeparationModel(coeffs, condition_warn=100)
res = model.invert(dw1=0.3, dw2=-6.84)          # G and 2D shifts, cm^-1 (illustrative)
print(f"strain {float(res.strain):.3f} %, hole-doping G shift {float(res.density):.2f} cm^-1")
```

```
warning: coefficient matrix condition number is 91.5; the mode pair is  ...
strain 0.200 %, hole-doping G shift 5.00 cm^-1
```

The G peak has barely moved here, but only because strain and doping
pushed it in opposite directions; the 2D peak separates the two. The
condition number of this set (91.5) is high mainly because its two
columns are in different units (percent and cm^-1), so the docstring
advises `condition_warn=100`. The noise amplification is moderate: a
1 cm^-1 error on both shifts costs about 0.035 % of strain and
1.6 cm^-1 of doping coordinate (asserted in the tests). The "density"
here is a G-band shift, not a carrier density; see below.

### 12. Pixels with no peak at all

```python
import warnings
import numpy as np
import ramansep as rs
from ramansep import lorentzian

# An illustrative 2 x 20 map: row 0 is off the flake (background and
# noise only, no Raman peak), row 1 carries both MoS2 peaks.
x = np.linspace(380.0, 480.0, 401)
rng = np.random.default_rng(6)
cube = np.empty((2, 20, x.size))
cube[0] = 50.0 + rng.normal(0.0, 3.0, (20, x.size))
cube[1] = (lorentzian(x, 404.7, 3.0, 40.0, 50.0) + lorentzian(x, 452.0, 4.0, 25.0, 0.0)
           + rng.normal(0.0, 3.0, (20, x.size)))

with warnings.catch_warnings():
    warnings.simplefilter("ignore", RuntimeWarning)   # overflow in runaway noise fits
    plain = rs.fit_map(x, cube, (395, 415), (440, 465), 404.7, 452.0)
    checked = rs.fit_map(x, cube, (395, 415), (440, 465), 404.7, 452.0, min_snr=3.0)

for name, f in (("without min_snr", plain), ("with min_snr=3 ", checked)):
    print(f"{name}: kept off-flake {f.ok[0].sum()} of 20, on-flake {f.ok[1].sum()} of 20")
for j in np.flatnonzero(plain.ok[0]):
    print(f"  off-flake pixel {j:2d} kept without min_snr: "
          f"dw1 {plain.dw1[0, j]:+.2f} +/- {plain.sigma1[0, j]:.2f} cm^-1")
for code in np.unique(checked.reason[0]):
    print(f"  off-flake reason {code}: {(checked.reason[0] == code).sum():2d} x {rs.MASK_REASONS[code]}")
```

```
without min_snr: kept off-flake 3 of 20, on-flake 20 of 20
with min_snr=3 : kept off-flake 0 of 20, on-flake 20 of 20
  off-flake pixel 11 kept without min_snr: dw1 +9.92 +/- 0.05 cm^-1
  off-flake pixel 15 kept without min_snr: dw1 -7.72 +/- 0.11 cm^-1
  off-flake pixel 17 kept without min_snr: dw1 -0.32 +/- 1.14 cm^-1
  off-flake reason 3: 13 x fit did not converge
  off-flake reason 4:  1 x fitted centre outside its window
  off-flake reason 5:  3 x fitted peak height not positive (a dip or no peak)
  off-flake reason 6:  3 x peak height below min_snr times its error bar
```

A fit to a spectrum that holds no peak (off the flake, over a hole)
sometimes converges on a noise spike and returns a shift that looks
like a measurement: pixel 11 above reports +9.92 +/- 0.05 cm^-1 from
pure noise. `min_snr=3` masks a pixel unless each fitted peak height is
at least 3 times its own error bar (a signal-to-noise test), and
removes all such pixels here while keeping every real one; the right
threshold for your data is your choice. `fwhm_range=(lo, hi)` adds a
check on the fitted widths. Both are off by default. Before 0.12.0,
peaks with a negative fitted height were kept too: on this map 6
off-flake pixels were kept instead of 3.

### 13. Two peaks that overlap

```python
import numpy as np
import ramansep as rs
from ramansep import lorentzian

# Two broad, overlapping peaks at known positions 385 and 404 cm^-1
# (illustrative widths), no noise: both shifts should come out as zero.
x = np.linspace(360.0, 430.0, 561)
y = lorentzian(x, 385.0, 5.0, 600.0, 20.0) + lorentzian(x, 404.0, 6.0, 1000.0, 0.0)

for joint in (False, True):
    dw1, dw2, s1, s2, f1, f2 = rs.fit_two_modes(x, y, (375, 394), (395, 415),
                                                ref1=385.0, ref2=404.0, joint=joint)
    print(f"joint={joint!s:5}: dw1 {dw1:+.4f}  dw2 {dw2:+.4f} cm^-1")
print(f"correlation of the two fitted centres: {f1.center_correlation:.3f}")
```

```
joint=False: dw1 +0.0795  dw2 -0.0279 cm^-1
joint=True : dw1 +0.0000  dw2 +0.0000 cm^-1
correlation of the two fitted centres: 0.014
```

Fitting each peak in its own window leaves the tail of the other peak
in that window, and the tail pulls the fitted centre: here by 0.08
cm^-1, comparable to the 0.1 cm^-1 shift error used in the examples
above, and as a systematic offset rather than random scatter. `joint=True` fits both peaks together, as two
Lorentzians on one shared baseline, over the whole range from the lower
edge of the first window to the upper edge of the second; that range
must not contain a third peak. The two fitted centres are then slightly
correlated (`center_correlation`); the inversion treats the two shifts
as independent, which is a good approximation when this number is
small. `fit_map(..., joint=True)` does the same at every pixel.

### 14. Smoothing a map that has gaps

```python
import numpy as np
import ramansep as rs

K = np.array([[-2.3, 1.1], [-0.9, -1.7], [0.4, 2.2]])   # illustrative, non-physical
rng = np.random.default_rng(7)
yy, xx = np.mgrid[0:12, 0:12]
strain = 0.002 * xx
density = 0.01 + 0.001 * yy

# One error bar per mode AND pixel (as fit_map returns them), and a
# 3 x 3 patch of pixels whose fits failed.
sig = np.array([0.03, 0.05, 0.02])[:, None, None] * rng.uniform(0.5, 2.0, (3, 12, 12))
shifts = np.einsum("mj,jhw->mhw", K, np.stack([strain, density])) \
    + rng.normal(size=(3, 12, 12)) * sig
shifts[:, 4:7, 4:7] = np.nan

r = rs.bayesian_map_inversion(K, shifts, sig, lam_strain=1000.0,
                              fill_masked=True, posterior_sigma=True)
hole = np.zeros((12, 12), bool)
hole[4:7, 4:7] = True
print(f"RMS strain error, measured pixels: {np.sqrt(np.mean((r.strain - strain)[~hole] ** 2)):.4f}")
print(f"RMS strain error, filled pixels  : {np.sqrt(np.mean((r.strain - strain)[hole] ** 2)):.4f}")
print(f"strain sigma: hole centre {r.strain_sigma[5, 5]:.4f}, "
      f"measured pixels (median) {np.median(r.strain_sigma[~hole]):.4f}")
print("filled centre = mean of its 4 neighbours:",
      np.isclose(r.strain[5, 5], (r.strain[4, 5] + r.strain[6, 5] + r.strain[5, 4] + r.strain[5, 6]) / 4))
```

```
RMS strain error, measured pixels: 0.0077
RMS strain error, filled pixels  : 0.0064
strain sigma: hole centre 0.0197, measured pixels (median) 0.0116
filled centre = mean of its 4 neighbours: True
```

`bayesian_map_inversion` accepts (since 0.12.0) one error bar per mode
and per pixel (an array of the same shape as `shifts`), such as the
`sigma1`, `sigma2` maps of `fit_map`. With `fill_masked=True`, a NaN
shift (or NaN error bar) counts as a missing measurement: it gets zero
weight, and the smoothness prior fills the gap from the neighbouring
pixels. A pixel with no data at all ends up exactly at the average of
its neighbours, and its error bar is larger than that of a measured
pixel. The filled values are interpolations, not measurements; they
were accurate here because the true fields are smooth by construction.
Filling needs both smoothness weights above zero.

### 15. Calibration uncertainty with three or more modes

```python
import numpy as np
import ramansep as rs

# Three modes calibrated together on six reference states (illustrative);
# the shifts are generated from K_true with 0.1 cm^-1 noise.
K_true = np.array([[-5.1, -2.2], [-20.9, 0.0], [-8.0, -1.1]])
eps = np.array([0.0, 0.4, 0.8, 0.0, 0.3, 0.6])
rho = np.array([0.0, 0.0, 0.2, 0.9, 0.6, 0.3])
rng = np.random.default_rng(4)
ref_shifts = np.column_stack([eps, rho]) @ K_true.T + rng.normal(0, 0.1, (6, 3))
cal = rs.calibrate_lever_arms(eps, rho, ref_shifts, sigmas=0.1)

out = rs.multimode_with_calibration(cal, shifts=[-2.0, -8.0, -3.4],
                                    sigmas=[0.15, 0.10, 0.12])
print(f"strain  {float(out['strain']):.4f} %")
print(f"  error from the spectra     {float(out['strain_sigma_shifts']):.4f}")
print(f"  error from the calibration {float(out['strain_sigma_calibration']):.4f}")
print(f"  total                      {float(out['strain_sigma']):.4f}")
print(f"density {float(out['density']):.4f} x 1e13 cm^-2 +/- {float(out['density_sigma']):.4f}")
print(f"p-value of the model check: {float(out['p_value']):.2f}")
```

```
strain  0.3821 %
  error from the spectra     0.0046
  error from the calibration 0.0016
  total                      0.0049
density 0.1036 x 1e13 cm^-2 +/- 0.0642
p-value of the model check: 0.06
```

`multimode_with_calibration` is the counterpart of example 8 for any
number of modes: the weighted least-squares separation of example 4,
with the calibration's own uncertainty carried into the error bars, to
first order. With more modes than causes the derivative of the answer
with respect to a lever arm also depends on the leftover mismatch
(residual) of the fit; the docstring gives the formula, and the tests
check it against a direct numerical derivative. The chi-square check
(`chi2_map`, `p_value`) is the one of `MultiModeModel`, which treats the
lever arms as exact.

### 16. Reading a map exported by the instrument

```python
import os, tempfile
import numpy as np
import ramansep as rs
from ramansep import lorentzian

# Write an illustrative 3 x 4 map the way mapping software often exports
# it: tab-separated "x y wavenumber counts", stage positions in
# micrometers, a serpentine scan, and the wavenumber running high to low.
rng = np.random.default_rng(2)
x_um = np.array([0.0, 0.5, 1.0, 1.5])
y_um = np.array([10.0, 10.5, 11.0])
wn = np.linspace(480.0, 380.0, 401)
path = os.path.join(tempfile.mkdtemp(), "map_export.txt")
with open(path, "w") as fh:
    fh.write("#X\t#Y\t#Wave\t#Intensity\n")
    for i, y in enumerate(y_um):
        for x in (x_um if i % 2 == 0 else x_um[::-1]):
            counts = (lorentzian(wn, 404.7 + 0.4 * x, 3.0, 800.0, 50.0)
                      + lorentzian(wn, 452.0, 4.0, 500.0, 0.0)
                      + rng.normal(0.0, 3.0, wn.size))
            for w, c in zip(wn, counts):
                fh.write(f"{x}\t{y}\t{w:.2f}\t{c:.1f}\n")

wavenumber, cube, xs, ys = rs.load_map_table(path)
print("cube:", cube.shape, " axis:", wavenumber[0], "to", wavenumber[-1], "cm^-1")
print("x (um):", xs, " y (um):", ys)
fits = rs.fit_map(wavenumber, cube, (395, 415), (440, 465), ref1=404.7, ref2=452.0)
print("A'1 shift at y = 10 um (cm^-1):", " ".join(f"{v:+.2f}" for v in fits.dw1[0]))

with open(path) as fh:                           # a scan that stopped early
    lines = fh.readlines()
with open(path, "w") as fh:
    fh.writelines(lines[:-401])
try:
    rs.load_map_table(path)
except ValueError as err:
    print("refused:", str(err).split(";")[0])
```

```
cube: (3, 4, 401)  axis: 380.0 to 480.0 cm^-1
x (um): [0.  0.5 1.  1.5]  y (um): [10.  10.5 11. ]
A'1 shift at y = 10 um (cm^-1): -0.00 +0.20 +0.40 +0.60
refused: the stage positions do not form a complete rectangular grid: 1 of the 3 x 4 combinations of distinct y and x values are missing (first: x = 1.5, y = 11.0)
```

Mapping software rarely writes the CSV layout of example 3. A common
plain-text export has four columns, the stage position x, the stage
position y, the wavenumber and the counts, with one line per pixel and
spectral point. `load_map_table` reads such a file in any line order
(raster or serpentine scans) and with the wavenumber running in either
direction. It rebuilds the pixel grid from the stage positions:
`cube[i, j]` is the spectrum at `x[j]`, `y[i]`, and both position
arrays increase. Columns may be separated by tabs, spaces or commas.
Lines starting with `#` and one text header line are skipped. If your
export has its columns in another order, say which with `usecols`.
A scan that is not a complete rectangle (here, one stopped before the
last pixel), a pixel with its own wavenumber axis, or a repeated point
is refused rather than padded or interpolated. Stage positions that
differ only in the last digits can be merged with `decimals=`.
`load_spectrum_table` reads a single two-column spectrum the same
way. Check that your export matches this layout before relying on it;
the package does not read vendor binary formats.

### 17. Calibrating when the starting state is unknown

```python
import numpy as np
import ramansep as rs

# A stage + gate calibration (illustrative). Strain (%) and electron
# density (1e13 cm^-2) are known only relative to the starting state, and
# the peak positions are absolute frequencies (cm^-1).
eps = np.array([0.0, 0.2, 0.4, 0.6, 0.0, 0.3, 0.1])
rho = np.array([0.0, 0.0, 0.0, 0.1, 0.8, 0.5, 1.1])
K_true = np.array([[-5.1, -2.2], [-20.9, 0.0]])
start = np.array([403.6, 451.1])        # the starting state: unknown in practice
rng = np.random.default_rng(2)
freq = start + np.column_stack([eps, rho]) @ K_true.T + rng.normal(0, 0.1, (7, 2))

np.set_printoptions(precision=2, suppress=True)
old = rs.calibrate_lever_arms(eps, rho, freq - [403.0, 452.0], sigmas=0.1)
print("through zero, shifts from 403.0 and 452.0 cm^-1:\n", old.K)
print("  chi2 / dof:", old.chi2 / old.chi2_dof)
new = rs.calibrate_lever_arms(eps, rho, freq, sigmas=0.1, fit_offset=True)
print("with a fitted offset per mode:\n", new.K)
print("  error bars:\n", new.K_sigma)
print("  chi2 / dof:", new.chi2 / new.chi2_dof)
print("  starting frequencies:", new.offset, "+/-", new.offset_sigma)
print("planned error bars:", rs.plan_calibration(eps, rho, 0.1, fit_offset=True)["K_sigma"])
```

```
through zero, shifts from 403.0 and 452.0 cm^-1:
 [[ -3.8   -1.67]
 [-22.62  -0.87]]
  chi2 / dof: [13.5  34.69]
with a fitted offset per mode:
 [[ -5.05  -2.23]
 [-20.57   0.03]]
  error bars:
 [[0.2 0.1]
 [0.2 0.1]]
  chi2 / dof: [0.95 1.35]
  starting frequencies: [403.63 450.98] +/- [0.08 0.08]
planned error bars: [0.2 0.1]
```

In a real calibration, the stage reports strain relative to where it
started and the gate reports charge relative to the starting voltage.
The starting strain and doping of the flake, and so the pristine
frequencies, are unknown. The fit of example 6 goes through zero and
needs them. Here the shifts were measured from guessed references
(403.0 and 452.0 cm^-1, both off), and the lever arms come out wrong
by up to 1.7 cm^-1 per percent. The fit itself says so: chi2 / dof is
far above 1. With `fit_offset=True` each mode gets one more free
parameter, a constant offset (the shift at the calibration's own
zero). The lever arms then depend only on how the shifts change
between the reference states, and they come out within two error
bars of the true values, with chi2 / dof near 1. The offset is reported with its error
bar, but it is not the pristine frequency, and the inversions do not
use it. The error bars of the lever arms already include the
uncertainty of the offset, so `separation_with_calibration` and
`multimode_with_calibration` take the result as it is. The price of
the extra parameter: the reference states must not lie on any one
straight line (not only a line through zero), and the fit needs one
more reference point. `plan_calibration` and `repeats_for_sigma` take
`fit_offset=True` too, so the error bars can be planned in advance.

## What is in the package

Every name below is exported from `ramansep` (its `__all__`). Each
docstring (`help(ramansep.fit_map)`, for example) gives inputs, units
and conventions.

**Two-mode inversion**

- `SeparationModel(coefficients, condition_warn=30.0)` -- inverts two
  shift maps (`invert(dw1, dw2, sigma1, sigma2)`) into a
  `SeparationResult` with `strain`, `density`, `strain_sigma`,
  `density_sigma`, `correlation` and `condition_number`; `forward`
  predicts shifts. Error bars are computed only when both `sigma1` and
  `sigma2` are given.
- `ModeCoefficients` -- the four lever arms of two modes plus the
  required `reference` string.

**Coefficient sets**

- `mos2_a1_2la()`, `mos2_eprime_a1()`, `graphene_g_2d_lee2012()` --
  the cited sets described in the next section.
- `synthetic_demo()` -- arbitrary, non-physical values for tests and
  demos only.

**More modes, and whole maps**

- `MultiModeModel(K, mode_names)`, `MultiModeResult` -- weighted
  least squares for any number (two or more) of modes, with `chi2_map`,
  `p_value` and `dof`.
- `compare_mode_sets(K, sigmas, mode_names, subset_size=2)` -- ranks
  mode subsets by delivered variance.
- `bayesian_map_inversion(K, shifts, sigmas, lam_strain, lam_density,
  posterior_sigma, max_dense, fill_masked)`, `BayesianMapResult` --
  joint inversion of a whole map with a smoothness prior (a Gaussian Markov random
  field, Rue and Held 2005: a standard statistical way of saying that
  each pixel is probably close to its four neighbours).

**From spectra to shifts**

- `lorentzian`, `fit_lorentzian`, `PeakFit` -- Lorentzian peak shape
  and a least-squares fit by the Levenberg-Marquardt method (a standard
  step-by-step curve-fitting method), using exact formulas for the
  derivatives; `baseline="constant"` or `"linear"`. Parameter errors
  are the standard linearised estimates, which are reliable when the
  residuals are dominated by uncorrelated noise.
- `voigt`, `fit_voigt`, `VoigtFit` -- the Voigt shape, computed
  exactly with the Faddeeva function (a standard special function,
  `scipy.special.wofz`), normalised to peak height.
- `fit_two_modes(x, y, window1, window2, ref1, ref2, baseline,
  joint=False)` -- fits both peaks, each in its own non-overlapping
  window, or with `joint=True` both together (example 13), and returns
  `(dw1, dw2, sigma1, sigma2, fit1, fit2)`.
- `fit_map(..., joint=False, min_snr=None, fwhm_range=None)`,
  `MapFitResult`, `MASK_REASONS` -- the same at every pixel of a cube,
  with masking of failed pixels and a per-pixel reason code
  (example 12).
- `load_spectrum_csv`, `save_spectrum_csv`, `load_map_csv`,
  `save_map_csv` -- the documented CSV formats of example 3.
- `load_spectrum_table(path, usecols, delimiter)`,
  `load_map_table(path, usecols, delimiter, decimals)`: plain text
  tables as instruments export them, with stage positions for maps
  (example 16, new in 0.13.0).

**Calibration**

- `calibrate_lever_arms(strain, density, shifts, sigmas, mode_names,
  fit_offset=False)`, `CalibrationResult` (`K`, `K_sigma`, `cov`,
  `chi2`, `chi2_dof`, `n_points`, `condition_number`, `mode_names`,
  `coefficients(reference)` for two modes, and, with
  `fit_offset=True`, `offset`, `offset_sigma` and `full_cov`; example
  17, new in 0.13.0).
- `plan_calibration`, `design_references`, `repeats_for_sigma` --
  planning, example 7 (`plan_calibration` and `repeats_for_sigma` also
  take `fit_offset=True`).
- `separation_with_calibration` -- example 8;
  `multimode_with_calibration` -- the same for three or more modes,
  example 15.

**Temperature**

- `ThreeCauseModel`, `ThreeCauseResult` -- example 10.
- `anti_stokes_ratio`, `calibrate_anti_stokes`,
  `temperature_from_anti_stokes`, `HC_OVER_KB_CM_K` (hc/kB in cm K) --
  example 9.

## Cited coefficient sets

Three sets ship with full provenance in their `reference` field and
docstring. Each docstring states which number comes from which source
and the conditions it applies under. Check that your sample matches
before use.

- `mos2_a1_2la()`: the A'1 + 2LA(M) pair for monolayer MoS2. Strain
  arms -5.1 and -20.9 cm^-1 per percent biaxial strain, computed by
  frozen-phonon DFT in the source paper (Mahim and Rahman, under
  review; data doi:10.5281/zenodo.21778170). A'1 doping response
  -2.2 cm^-1 per 1e13 cm^-2 of electrons, measured with a gate by
  Chakraborty et al., Phys. Rev. B 85, 161403(R) (2012). The 2LA(M)
  doping response has never been measured and is set to zero; the
  source paper bounds the consequence at 16% of a recovered edge
  charge. Stated conditions: 1H monolayer (the common semiconducting
  crystal form), 532 nm laser, on an SiO2 substrate, 300 K. The tests
  reproduce the source paper's published separations from this set
  (example 1).
- `mos2_eprime_a1()`: the historically used all-optical E' + A'1 pair
  (Michail et al., Appl. Phys. Lett. 108, 173102 (2016)). Strain arms
  -4.31 (E') and -2.50 (A'1) cm^-1 per percent, from the measured
  biaxial Grueneisen parameters 0.56 and 0.31 of Michail et al., ACS
  Appl. Mater. Interfaces 16, 49602 (2024), converted with the mode
  frequencies 385.0 and 403.0 cm^-1. These are rounded monolayer
  values (see Mignuzzi et al., Phys. Rev. B 91, 195411 (2015)); Michail
  et al. 2016 place their strain-doping origin at 384.6 and
  402.7 cm^-1, which gives the same arms to the digits shown. Doping responses -0.33 (E') and
  -2.2 (A'1) cm^-1 per 1e13 cm^-2, measured by Chakraborty et al.
  (2012). Stated conditions: 1H monolayer, on SiO2, 300 K. Both modes
  respond weakly to strain, so it carries about five times the strain
  uncertainty (example 2) -- included for comparison.
- `graphene_g_2d_lee2012()`: the G + 2D decomposition of Lee et al.,
  Nat. Commun. 3, 1024 (2012). Its "doping" output is deliberately a
  G-band shift, not a density: graphene's doping response is
  nonlinear, so no universal linear rate exists to ship, and the
  docstring points to the gated calibrations needed for conversion.
  Stated for hole doping, uniaxial strain of random orientation.

**Tungsten-based materials (2024-2026).** For W-based monolayers the
literature now provides verified strain and temperature responses --
monolayer WSe2 biaxial strain rates (Michail et al., ACS Appl. Mater.
Interfaces 16, 49602 (2024)), monolayer WS2 strain rates (Roy, Yang
and Gao, Sci. Rep. 14, 3860 (2024)) and temperature coefficients
(Huang et al., Sci. Rep. 6, 32236 (2016)) -- while gate-calibrated
doping arms remain without a clean linear coefficient (the measured
behavior is an electron-only, threshold-like softening of the
out-of-plane modes, the vibrations in which atoms move perpendicular
to the sheet; Sohier et al., PRX 9, 031019 (2019)). That is why
no W-material set ships as constants here: the strain column is
citable, the doping column is not yet, and this package does not ship
half a lever-arm matrix. Calibrate your own with
`calibrate_lever_arms` -- the doping arm from your own gated reference
points -- and the provenance travels with the analysis.

For monolayer MoS2 temperature coefficients, the `thermal` module
docstring points to S. Sahoo et al., J. Phys. Chem. C 117, 9042
(2013); no temperature coefficient is shipped.

## When it refuses, and why

`ramansep` raises `ValueError` instead of guessing when:

- a lever-arm matrix contains NaN or infinity; a two-mode lever-arm
  matrix is not 2 x 2 or is exactly singular (the two modes respond
  identically); a merely poorly conditioned one
  (condition number above `condition_warn`, default 30) gives a
  warning, not an error;
- a multimode or map lever-arm matrix has rank below 2 (below 3 for
  `ThreeCauseModel`, which also needs at least three modes), or the
  weighted problem is singular at some pixel;
- the shift maps do not match in shape or number;
- a shift error (`sigmas`) is zero or negative in `MultiModeModel`,
  `ThreeCauseModel`, `bayesian_map_inversion`, the calibration and
  planning tools, is zero, negative, NaN or infinite in
  `compare_mode_sets` (which also refuses a `K` that is not (m, 2)),
  or negative in
  `SeparationModel` (zero is allowed there; NaN marks a masked pixel
  and is passed through as NaN);
- `bayesian_map_inversion` gets NaN or infinite shifts or error bars
  without `fill_masked=True`, missing values with a smoothness weight
  of zero, so little data that the unmasked shifts cannot fix both
  fields, per-pixel error bars of the wrong shape, a negative, NaN or
  infinite smoothness weight (NaN and infinity are refused since
  0.13.0), or asks for exact error bars above `max_dense` unknowns;
- `MultiModeModel` or `ThreeCauseModel` gets a number of error bars
  other than one per mode, or shift maps of different shapes (since
  0.13.0 the message says which);
- a peak fit gets too few points (5 for a Lorentzian with a constant
  baseline, 6 with a linear baseline or for a Voigt), non-finite data,
  an unknown `baseline` name, or a non-positive starting width;
- the two fit windows are reversed, overlap, or hold too few points,
  the `baseline` name is unknown, or a reference frequency `ref1`,
  `ref2` is NaN or infinite (`fit_two_modes`, and `fit_map` before
  fitting any pixel; the reference check is new in 0.13.0); `fit_map` also refuses a wavenumber axis
  that contains NaN or infinity, a `min_snr` that is not positive, or
  a `fwhm_range` that is not (lo, hi) with 0 <= lo < hi;
- a CSV file has the wrong header or field count, a wavenumber axis
  that is not strictly increasing or contains NaN or infinity, a pixel
  grid with holes, or a pixel with its own wavenumber axis;
- a text table (`load_spectrum_table`, `load_map_table`) has lines
  without the requested numeric columns, a NaN or infinite wavenumber
  or stage position, a wavenumber axis that repeats a value or is not
  monotonic, stage positions that do not form a complete rectangle, or
  a pixel with its own wavenumber axis or number of points;
- calibration reference states lie on one line through the origin of
  the (strain, density) plane (the two lever arms cannot be told
  apart), there are fewer than 2 of them, or fewer than 3 when no
  `sigmas` are given; with `fit_offset=True`, the states lie on any
  one straight line, or there are fewer than 4 of them without
  `sigmas`; `sigmas` whose shape fits neither the points nor the
  modes;
- `CalibrationResult.coefficients` is asked to package other than two
  modes, or `separation_with_calibration` gets such a calibration
  (`multimode_with_calibration` takes any number);
- the anti-Stokes ratio is at or above the calibration constant, or a
  temperature, frequency, ratio or constant is not positive.

## How the results are checked

114 automated tests run on every push and pull request, on Python 3.9,
3.10, 3.11, 3.12, 3.13 and 3.14, and once more on Python 3.9 with the
oldest NumPy (1.22.0) and SciPy (1.8.0) the package allows. The
numerical checks compare the package with an independent calculation,
a closed form, a published number, or a simulation with fixed random
seeds; the others check that refusals fire. The main checks, with the
tolerance each test really uses. "To 1e-12" means the difference may be
at most 1e-12; "plus a relative 1e-5" means that 1e-5 times the size
of the expected value is also allowed. Where a test does not set both
parts, NumPy's own defaults apply, and these are listed too:

**Two-mode inversion and the cited sets**

- Forward-then-invert round trips recover strain and density to
  1e-12 plus a relative 1e-7 (synthetic set, both MoS2 sets) and to
  1e-12 plus a relative 1e-5 (graphene set).
- Reported error bars match the scatter of 200 000 seeded noisy
  inversions to within 2 %.
- MoS2 A'1 + 2LA(M), paper edge spot: zero strain to 1e-12, density
  0.5/2.2 x 1e13 cm^-2 to a relative 1e-12, within 0.5e11 of the
  published 2.3e12 cm^-2. Interior spot: strain 0.134 % within 0.0005,
  density below 7.6e11 cm^-2 in size.
- The E' + A'1 strain arms equal -2 gamma omega / 100 from the cited
  Grueneisen parameters within 0.005; its strain error bar is more
  than 4 times that of the A'1 + 2LA(M) pair.
- Graphene: pure-strain and pure-doping trajectories decompose to zero
  doping and zero strain to 1e-12; the noise amplification equals its
  closed form to a relative 1e-12 plus 1e-8.

**Peak fitting and maps**

- Analytic derivatives match finite differences (Lorentzian: 1e-7
  plus a relative 1e-5; Voigt: 1e-6).
- A noiseless Lorentzian is recovered from automatic starting values
  (centre and width to 1e-8, height and offset to 1e-6); a noiseless
  Voigt to 1e-8.
- On 40 seeded noisy spectra, (centre error / reported sigma) has a
  standard deviation between 0.5 and 2.0 and a mean below 0.6 in size.
- The Voigt with zero Lorentzian width equals the Gaussian to 1e-14;
  as the Gaussian width shrinks tenfold (1e-3 to 1e-4), the gap to
  the Lorentzian shrinks by a factor of more than 9 (and is below 1e-5
  at width 1e-3).
- A sloped background biases a constant-baseline centre by more than
  3 of its sigmas; the linear baseline recovers the centre to 1e-6.
- `fit_map` gives exactly (bit for bit) the `fit_two_modes` result at
  the pixel checked; on a noiseless synthetic cube the maps invert back
  to the truth (strain to 1e-4, density to 1e-3); a dead pixel is
  masked and the other pixels' shifts are bit-identical to the clean
  run; with noise, every strain error is below 5 error bars + 1e-3.
- The CSV save/load round trips are exact (bit-identical arrays).
- Configuration errors in `fit_map` (unknown `baseline` name,
  reversed or overlapping windows, too few points in a window, a NaN
  in the wavenumber axis) are refused, and valid settings still give
  the `fit_two_modes` result bit for bit (new in 0.11.1).

**Multimode and map inversion**

- With two modes, `MultiModeModel` equals `SeparationModel` (maps to
  1e-12 plus a relative 1e-5, error bars to a relative 1e-12 plus
  1e-8).
- On 1600 seeded pixels the mean chi-square is within 0.1 of its
  degrees of freedom and the mean p-value within 0.05 of 0.5; a hidden
  extra shift gives median p-values below 1e-6, against above 0.1
  elsewhere.
- `bayesian_map_inversion` at `lam = 0` equals the per-pixel result
  (maps to 1e-12, error bars to 1e-13); a constant noiseless truth is
  recovered to 1e-12 at every `lam` tried (0, 1, 250); at `lam = 1e9`
  it matches an independently computed pooled estimate to 1e-5; error
  bars never grow as `lam` increases (checked at `lam` = 0, 0.5, 2
  and 10). A rank-deficient matrix is refused at `lam` = 0 and 1 (new
  in 0.11.1).

**Calibration and planning**

- Noise-free calibration recovers the lever arms to 1e-12, and agrees
  with an independent `numpy.linalg.lstsq` solve to 1e-10.
- Over 400 seeded calibrations the scatter of the fitted arms matches
  the reported error bars within 15 %.
- The planned covariance equals the one the calibration reports to
  1e-15; repeating a design 7 times divides it by 7 to a relative
  1e-12; on a test set of candidates where only one is off the line,
  the greedy design includes it, and it is never worse than any of 30
  random subsets of the same size.
- `separation_with_calibration`: with zero calibration covariance it
  equals the plain inversion to 1e-12 plus a relative 1e-5; its
  calibration part matches
  400 seeded lever-arm draws within 20 %; the derivative identity it
  uses matches finite differences to a relative 1e-5 plus 1e-10.

**Temperature**

- `ThreeCauseModel` recovers noise-free strain and density to 1e-9
  and temperature to 1e-8, agrees with an independent per-pixel
  least-squares solve to the same tolerances, and with a known
  temperature subtracted matches the two-cause solver to 1e-10.
- hc/kB agrees with SciPy's CODATA table to a relative 1e-9.
- The anti-Stokes thermometer round trip holds to a relative 1e-9,
  the calibration constant to 1e-12, and the error propagation matches
  a finite difference to a relative 1e-5.
- With a NaN shift at one pixel and a NaN error bar at another,
  `ThreeCauseModel` returns NaN at exactly those two pixels and, at
  every other pixel, bit-identical results to the run without them
  (new in 0.12.0).

**New in 0.12.0** (`tests/test_v012.py`, 19 tests)

- Joint two-peak fit: its exact derivatives match central finite
  differences to 1e-6 of each derivative's largest value (plus 1e-9),
  with both baselines. On a noiseless pair 19 cm^-1 apart (widths 5
  and 6 cm^-1) the separate-window fits are off by more than 0.05 and
  0.02 cm^-1; the joint fit recovers both centres and widths to 1e-8
  cm^-1 and heights and offset to 1e-6, also on a sloped background
  (slope to 1e-8) and with the windows given in reverse order. On 600
  seeded noisy spectra of a pair 12 cm^-1 apart, the scatter of each
  fitted centre is within 12 % of the reported error bar, the mean
  error is below 4 of its standard errors, and the measured
  correlation of the two centres is within 0.12 of the reported one
  (which is above 0.05). `fit_map(..., joint=True)` equals
  `fit_two_modes(..., joint=True)` bit for bit at every pixel checked.
- `fit_map` checks: on 30 seeded peak-free spectra, without `min_snr`
  at least one pixel is kept and at least one is masked for a negative
  height; with `min_snr=3` none is kept. All 10 seeded spectra with
  real (weak) peaks are kept, with shifts and error bars bit-identical
  to the run without the checks. Constructed spectra produce the
  reason codes 0, 1, 5 and 7.
- `bayesian_map_inversion`: with per-pixel error bars and `lam = 0`
  it equals `MultiModeModel` given the same error bars (maps to 1e-12,
  error bars to 1e-13); per-pixel error bars that are constant per
  mode equal the scalar path (maps to 1e-12, error bars to 1e-13).
  With a pixel without data, a missing mode at another pixel and a
  missing error bar at a third, `fill_masked=True` agrees with an
  independent dense least-squares solve of the same objective
  (`numpy.linalg.lstsq`) to 1e-10. A filled pixel equals the average
  of its neighbours to 1e-12 (interior and corner pixel), with a larger
  error bar than its measured neighbours; a constant truth is recovered
  through a 2 x 3 hole to 1e-12 at `lam` = 0.5 and 250.
- `multimode_with_calibration`: with a two-mode calibration it equals
  `separation_with_calibration` (all values and error bars, relative
  1e-10 plus 1e-13); its derivative formula matches central finite
  differences of the `MultiModeModel` re-solve to a relative 1e-6 plus
  1e-10 at a point with a non-zero residual; its calibration part
  matches 400 seeded lever-arm draws within 20 %; with zero
  calibration covariance it gives the `MultiModeModel` error bar to a
  relative 1e-12; masked pixels stay NaN.
- New refusals fire: non-finite lever arms (`SeparationModel`,
  `MultiModeModel`, `compare_mode_sets`), a negative shift error in
  `SeparationModel` (NaN still passes as a masked pixel), a zero,
  negative, NaN or infinite sigma and a `K` that is not (m, 2) in
  `compare_mode_sets`.

**New in 0.13.0** (`tests/test_v013.py`, 18 tests)

- Each fix below has a test that fails on 0.12.0 and passes now.
- Text tables: files are written line by line from arrays built
  independently in the test (a serpentine scan with the wavenumber
  running high to low; a comma-separated file with a text header and
  the columns in another order). `load_map_table` returns those arrays
  exactly (bit for bit), and `fit_map` on the loaded cube equals
  `fit_map` on the original cube bit for bit. Each refusal fires.
  `load_spectrum_table` reverses a descending export exactly.
- `calibrate_lever_arms(..., fit_offset=True)`: noise-free arms and
  offsets are recovered to 1e-12. With per-point error bars, the arms
  and offsets agree with an independent `numpy.linalg.lstsq` solve to
  1e-10, and the covariance agrees to a relative 1e-10 with the closed
  form of the Frisch-Waugh-Lovell theorem (arms from the
  weighted-centered design; offset variance 1/sum(w) plus the
  propagated mean state). Without error bars the covariance equals the
  residual-variance formula with n - 3 degrees of freedom (relative
  1e-9). The through-origin fit of offset data is off by exactly the
  closed-form projection of the offsets (to 1e-12), more than 1 cm^-1
  per percent here; the offset fit removes it. Over 400 seeded
  calibrations the scatter of the arms and offsets matches the
  reported error bars within 15 %, and the mean chi-square is within
  0.4 of its 4 degrees of freedom. The calibration part of
  `separation_with_calibration` matches 400 seeded offset calibrations
  within 20 %. `plan_calibration(..., fit_offset=True)` equals the
  reported covariance to 1e-15, and repeating the design r times
  divides it by r (relative 1e-12).

## Corrections in earlier versions

**0.13.0 fixed three silent failures** (before and after values from
running 0.12.0 and 0.13.0 on the same input).

- `load_spectrum_csv` returned a wavenumber axis containing NaN or
  infinity, although it promises to refuse an axis that is not
  strictly increasing: a file with the axis 400, nan, 402, 403 loaded
  as `[400, nan, 402, 403]`. It is now refused. (`load_map_csv`
  refused such an axis already, but with the message that pixel
  (0, 0) "carries a different wavenumber axis than pixel (0, 0)"; the
  message now says what is wrong.)
- `fit_map` with a NaN reference frequency reported every pixel as
  fitted (`ok` True at all pixels, `n_masked` 0) while every shift was
  NaN, and `fit_two_modes` returned a NaN shift. Both now refuse a
  reference that is not a finite number.
- `bayesian_map_inversion` with a NaN or infinite smoothness weight
  returned all-NaN maps with only SciPy's "Matrix is exactly singular"
  warning. It now refuses such a weight and suggests a large finite
  one for the fully pooled limit.

**Documentation corrections in 0.13.0.** The docstring of
`MultiModeResult` said that `condition_number` is that of the weighted
design matrix when the error bars are scalars; the code has always
returned the condition number of the unweighted lever-arm matrix (the
value is unchanged; the docstring now says so). The docstring of
`fit_voigt` said its default start splits the width evenly between the
Gaussian and Lorentzian parts; it starts with a Gaussian FWHM equal to
the estimated width and a Lorentzian FWHM of half of it. The
`fitting` module docstring said the tests recover amplitude and
baseline to 1e-8; they check 1e-6. `SeparationResult.density` was
documented as cm^-2 by default; it is in the unit of the lever arms
(1e13 cm^-2 for the MoS2 sets).

**0.12.0 fixed three silent or misleading behaviours.**

- `fit_map` kept pixels whose fitted "peak" had a negative height (a
  dip fitted to noise). Such pixels are now masked (reason code 5).
  Good pixels give the same numbers as before; on the map of example
  12, the off-flake pixels kept went from 6 of 20 to 3 of 20 (and to 0
  with `min_snr=3`).
- `ThreeCauseModel.invert` refused a whole map, with a misleading
  "singular" message, when one pixel carried a NaN error bar, which is
  how `fit_map` marks masked pixels. Such pixels now come out as NaN,
  as in `MultiModeModel`.
- `SeparationModel.invert` and `compare_mode_sets` accepted negative
  shift errors and squared the sign away (for example `sigma1=-0.1`,
  `sigma2=0.1` with `mos2_a1_2la()` gave a strain error bar of
  0.0048 %). Both now refuse them. `compare_mode_sets` also returned
  NaN variances, ranked first, for a zero or NaN sigma, and silently
  left out every subset with an infinite lever arm; it now refuses
  those inputs too.

**0.11.1 fixed two silent failures.**

- `fit_map` treated a mistake in its own settings (for example
  `baseline="Linear"`, or overlapping windows) like a bad pixel: it
  masked every pixel and returned an all-NaN result without an error.
  It now refuses such settings before fitting, as `fit_two_modes`
  does. Valid settings give the same results as before. `fit_map` now
  also refuses a wavenumber axis that contains NaN or infinity; before,
  spectral points with a NaN wavenumber were silently left out of the
  fits.
- `bayesian_map_inversion` did not check that the lever-arm matrix
  can separate the two causes. With a rank-deficient matrix it
  returned NaN maps (at `lam = 0`) or arbitrary finite maps (at
  `lam > 0`). It now refuses, as `MultiModeModel` does.

**0.8.0:** `bayesian_map_inversion` now refuses NaN shifts; before, it
solved around them silently.

**Documentation corrections in 0.11.1.** The README of 0.11.0 said the
tests ran on Python 3.9-3.14, but the CI matrix skipped 3.10 (added
now, with an oldest-dependencies job). It called every test "pinned to
an exact result" at "machine precision"; many use stated tolerances,
now listed above. It said the tests reproduce published separations
for all three cited sets; they do so for `mos2_a1_2la` only. The
calibration check described as "an independent QR solve" uses
`numpy.linalg.lstsq`. It cited the E' + A'1 pair itself to Michail et
al. (2024); that paper supplies the pair's Grueneisen parameters, and
the pair's historical use is Michail et al. (2016), as the docstring
says.

The full history is in [CHANGELOG.md](CHANGELOG.md).

## Limits, and what is deliberately not included

- The model is linear: each shift is a fixed lever arm times each
  cause. Where a material responds nonlinearly (graphene doping, the
  W-based doping behaviour above), no linear density coefficient is
  shipped.
- No constants beyond the three cited sets. Response coefficients
  depend on material, mode pair, laser wavelength and substrate; a
  measurement tool that ships unverified constants propagates wrong
  results. For any other system you provide the numbers -- from the
  literature or from `calibrate_lever_arms` -- and the mandatory
  `reference` field makes their origin travel with the analysis.
- The smoothing weights of the Bayesian inversion are user-chosen,
  because estimating them would need noise assumptions this package
  does not make. Pixels it fills (`fill_masked=True`) are
  interpolations from their neighbours, not measurements.
- Peak-fit error bars are linearised estimates; they are reliable when
  the residuals are dominated by uncorrelated noise.
- By default `fit_two_modes` fits each peak in its own window, and the
  tail of the other peak can shift a fitted centre (by 0.08 cm^-1 in
  example 13; the round-trip test allows 5e-3 cm^-1). `joint=True`
  removes this, provided no third peak lies anywhere from the lower
  edge of the first window to the upper edge of the second; it fits
  Lorentzians only (no joint Voigt fit). The inversion then ignores the
  correlation between the two fitted centres (reported as
  `center_correlation`).
- `min_snr` and `fwhm_range` in `fit_map` are off by default and their
  values are your choice; the package does not decide what counts as
  "no peak" for your data.
- `load_map_table` groups stage positions by exact value (or after
  rounding with `decimals`) and needs one shared wavenumber axis; it
  does not resample spectra onto a common axis, and it does not read
  vendor binary formats or decimal commas. `design_references` plans
  for the through-origin calibration only; check a set chosen for
  `fit_offset=True` with `plan_calibration(..., fit_offset=True)`.
- `separation_with_calibration` and `multimode_with_calibration` are
  first-order in the calibration errors. The calibration part is
  common to all pixels, so it does not average down over a map. The
  chi-square check of `multimode_with_calibration` ignores the
  calibration uncertainty.

## Where it comes from

> T. M. Mahim and M. M. Rahman, "Two Raman phonons quantify the fixed
> edge charge left by patterning monolayer transition metal
> dichalcogenides" (under review). Code for the paper:
> https://github.com/Tanvir-Mahmud-Mahim/Width-scaling-in-monolayer-semiconductor-nanoribbon-transistors

This package is the general-purpose, material-agnostic tool; the
paper repository reproduces the specific published study.

## Citing, support and license

If `ramansep` helps your work, please cite it with the concept DOI
[10.5281/zenodo.22014913](https://doi.org/10.5281/zenodo.22014913),
which always resolves to the latest version; every release is
archived on Zenodo. [CITATION.cff](CITATION.cff) has the details.

Written and maintained by Tanvir Mahmud Mahim (Department of
Electrical and Electronic Engineering, BRAC University), who reviews
every change and takes the final decision on scope and releases.
Design questions are discussed in the open in issues and pull
requests, and the standing rule of
[CONTRIBUTING.md](CONTRIBUTING.md) binds the maintainer exactly as it
binds contributors: a change that touches physics arrives with a
test, and a constant arrives with its source.

Support runs through the
[issue tracker](https://github.com/TaN-MM-Org/ramansep/issues).
Usage questions are welcome alongside bug reports; a docstring that
left a unit or a sign convention unclear is treated as a
documentation bug, not user error. While the version is below 1.0
the API may still move between minor versions; such changes are
called out in the release notes.

Licensed under Apache-2.0.
