Method: input.csv holds one merged energy axis (1.196-1.606 eV, 1024 points) and 15 spectra (5-100 uW). The dominant PL peak sits near 1.33 eV. Each spectrum was fitted with a Lorentzian plus linear background, y = y0 + m(E-1.33) + (2A/pi) w / (4(E-xc)^2 + w^2), over 1.25-1.42 eV. (A single Lorentzian over the full range fails, because of the other features near 1.55 eV and above 1.6 eV.) Peak intensity is the fitted height 2A/(pi w).

Per-power results (power uW: xc eV, intensity counts):
5: 1.31684, 125.1 | 10: 1.32191, 230.0 | 15: 1.32411, 325.5 | 20: 1.32592, 413.9 | 25: 1.32734, 505.7
30: 1.32831, 605.1 | 35: 1.32891, 680.4 | 40: 1.33027, 805.9 | 45: 1.33049, 867.3 | 50: 1.33158, 965.6
60: 1.33209, 1158.1 | 70: 1.33341, 1354.8 | 80: 1.33360, 1542.8 | 90: 1.33430, 1772.0 | 100: 1.33494, 2046.8
(FWHM about 0.065-0.074 eV throughout.)

Allometric fit y = a * P^b, least squares, uncertainty from the covariance matrix:
- Intensity vs power: a = 19.6, b = 1.002 +/- 0.019 (unweighted). The fit weighted by the Lorentzian-fit errors gives b = 0.967 +/- 0.008. Changing the fit window (1.22-1.45 or 1.27-1.40 eV) gives b between 0.96 and 1.00. So the intensity is linear in power within about 3%.
- Peak position vs power: a = 1.308, b = 0.0045 +/- 0.0001 (weighted: 0.0042 +/- 0.00005). This is a small blue shift of about 18 meV over 5-100 uW.

Note: the plots (intensity vs power and xc vs power, each with the power-law curve) were not rendered as image files; the fit tables above hold the data for them.

RESULT: b = 1.002 ± 0.019
