Method: input.csv holds one energy axis (1.196-1.606 eV, 1024 points) and 15 spectra (5-100 uW), so the spectra are already merged. Each spectrum was fitted with a sum of three Lorentzians plus a constant (scipy curve_fit). The Lorentzians sit at about 1.33 eV (main peak), about 1.54 eV, and one at the high-energy edge near 1.6 eV. The edge peak is cut off by the data range and is poorly constrained, so I do not use it. Parameters per power are in peak_params.csv.

Main peak (~1.33 eV), Lorentzian fit:
- Intensity (peak height) goes from 123 counts at 5 uW to 1989 counts at 100 uW.
- Position xc goes from 1.3180 eV to 1.3350 eV.

The two plot datasets are intensity vs power and xc vs power. I could not render images because matplotlib is not installed, so only the fit values were produced. Allometric fit y = a * P^b on each:
- Intensity vs power: a = 21.3, b = 0.971 +/- 0.019 (unweighted least squares). A fit weighted by the Lorentzian-fit errors gives b = 0.927, but its reduced chi2 is 46, so the unweighted fit is the more honest one. The intensity is close to linear in power.
- Position vs power: a = 1.310 eV, b = 0.0042 +/- 0.0001. The peak blue-shifts only slightly, about 1.3% over the range.

Second peak (~1.54 eV), for reference: intensity b = 0.748 +/- 0.021, position b = 0.0018 +/- 0.0005.

The exponent reported below is b from the intensity-vs-power allometric fit of the main peak.

RESULT: b = 0.971 ± 0.019
