Data: input.csv has one energy axis (1.196-1.606 eV, 1024 points) and 15 spectra (5-100 uW). Each column is one spectrum, so no merging was needed. Each spectrum has a dominant broad peak near 1.33 eV, a weaker broad band near 1.54 eV, and a rising edge at the high-energy limit (>1.6 eV, truncated).

Lorentzian fit: for every power I fitted y = y0 + A / (1 + ((E - xc)/(w/2))^2) to the main peak (window 1.22-1.45 eV) with scipy curve_fit. Intensity A grows from about 114 counts at 5 uW to about 1980 counts at 100 uW. xc shifts from 1.3193 eV to 1.3357 eV, a blue shift. The FWHM stays at about 0.062-0.066 eV.

Allometric fit (y = a * P^b, curve_fit, unweighted):
- Intensity vs power: a = 18.4, b = 1.009 +/- 0.019. This is linear in power. Weighting by the Lorentzian fit errors gives b = 0.966 +/- 0.002; other fit windows give b = 1.01-1.05.
- Peak position vs power: a = 1.312 eV, b = 0.0039 +/- 0.0001. This is a very weak power-law blue shift.

A 3-Lorentzian-plus-baseline fit of the full spectrum (second peak near 1.54 eV) gave main-peak b = 0.97 and second-peak b = 0.75. The third feature is truncated at the edge of the range and is poorly constrained.

The plots (intensity vs power and xc vs power, each with the allometric fit) could not be rendered because matplotlib is not installed. The fitted values behind them are listed above.

RESULT: b = 1.009 ± 0.019
