Base de données des spectres optiques

#!/usr/bin/env python3
"""
make_laser_spectrum.py - generate an OghmaNano spectrum file for a laser source.

Output is an OghmaNano spectra directory:
    <outdir>/data.json    spectra item metadata
    <outdir>/spectra.csv  #oghma_csv {json header}*
                          <wavelength in m>\t<spectral irradiance in W m^-2 m^-1>

If <outdir> exists it is reused; only data.json and spectra.csv are
overwritten, nothing is deleted.

The spectrum is normalised so that the integral over wavelength (in metres)
equals --power (W/m^2). Default is 10 W/m^2.

Models
------
gaussian : single line, Gaussian (or Lorentzian) of FWHM --fwhm.
fp       : Fabry-Perot multimode laser. Longitudinal modes spaced by
           d_lambda = lambda^2 / (2 n_g L) under a Gaussian gain envelope of
           FWHM --fwhm. Each mode has width --mode-width (effectively the
           OSA resolution; true FP mode linewidths are ~1e-4 nm).
dfb      : DFB laser. Main line of FWHM --fwhm plus two side modes at the
           grating stopband edges (lambda0 +/- stopband/2), suppressed by
           --smsr dB, with optional extra suppression of one side
           (--side-asym) since real DFBs usually lase with one dominant
           side mode.

Any model can add a broad amplified-spontaneous-emission (ASE) pedestal
--ase-db dB below the peak with FWHM --ase-fwhm.

Examples
--------
  # no arguments: 1550 nm, 2 nm FWHM Gaussian, 10 W/m^2 -> ./laser_1550nm/
  python3 make_laser_spectrum.py

  # Fabry-Perot, 300 um InGaAsP cavity, 2 nm envelope, ASE 40 dB down
  python3 make_laser_spectrum.py fp_1550 --model fp --cavity-length 300 --ase-db 40

  # DFB, 0.1 nm line, 40 dB SMSR, 3 nm stopband, preview plot
  python3 make_laser_spectrum.py dfb_1550 --model dfb --fwhm 0.1 --smsr 40 --plot
"""

import argparse
import json
import math
import os
import sys

import numpy as np

LN2 = math.log(2.0)
_trapz = getattr(np, "trapezoid", None) or np.trapz

# data.json for an OghmaNano spectra item, written verbatim.
DATA_JSON = """{
\t"item_type":"spectra",
\t"status":"public",
\t"changelog":"",
\t"spectra_import": {
\t\t"import_file_path":"none",
\t\t"import_x_combo_pos":9,
\t\t"import_data_combo_pos":5,
\t\t"import_x_spin":0,
\t\t"import_data_spin":2,
\t\t"import_title":"Voltage - J",
\t\t"import_xlabel":"Voltage",
\t\t"import_data_label":"J",
\t\t"import_area":0.104,
\t\t"import_data_invert":"False",
\t\t"import_x_invert":"False",
\t\t"data_file":"none"
\t\t}
\t}
"""


def gaussian(x, x0, fwhm):
    return np.exp(-4.0 * LN2 * ((x - x0) / fwhm) ** 2)


def lorentzian(x, x0, fwhm):
    hw = 0.5 * fwhm
    return hw * hw / ((x - x0) ** 2 + hw * hw)


def line(x, x0, fwhm, shape):
    return gaussian(x, x0, fwhm) if shape == "gaussian" else lorentzian(x, x0, fwhm)


def mode_spacing_nm(a):
    if a.mode_spacing is not None:
        return a.mode_spacing
    lam = a.wavelength * 1e-9
    L = a.cavity_length * 1e-6
    return lam * lam / (2.0 * a.group_index * L) * 1e9


def build(a):
    """Return (wavelength_nm, relative_intensity) before power normalisation."""
    lam0, fwhm = a.wavelength, a.fwhm

    # Window half-width: enough to make the line tails negligible.
    tail = 5.0 if a.shape == "gaussian" else 50.0
    half = tail * fwhm

    if a.model == "fp":
        half = max(half, 3.0 * fwhm + tail * a.mode_width)

    if a.model == "dfb":
        half = max(half, 0.5 * a.stopband + tail * fwhm)

    if a.ase_db is not None:
        half = max(half, 2.0 * a.ase_fwhm)

    if a.window is not None:
        half = 0.5 * a.window

    # Step: resolve the narrowest feature with ~20 points per FWHM.
    narrowest = a.mode_width if a.model == "fp" else fwhm
    step = a.step if a.step is not None else narrowest / 20.0

    n = int(round(2.0 * half / step)) + 1

    if n > 500000:
        sys.exit(f"error: {n} points needed; increase --step or reduce --window")

    lam = np.linspace(lam0 - half, lam0 + half, n)

    if a.model == "gaussian":
        I = line(lam, lam0, fwhm, a.shape)

    elif a.model == "fp":
        dl = mode_spacing_nm(a)
        kmax = int(half / dl) + 1
        I = np.zeros_like(lam)

        for k in range(-kmax, kmax + 1):
            lk = lam0 + k * dl
            I += gaussian(lk, lam0, fwhm) * line(
                lam, lk, a.mode_width, a.shape
            )

    elif a.model == "dfb":
        s = 10.0 ** (-a.smsr / 10.0)
        s_short = s * 10.0 ** (-a.side_asym / 10.0)

        I = (
            line(lam, lam0, fwhm, a.shape)
            + s * line(
                lam, lam0 + 0.5 * a.stopband, fwhm, a.shape
            )
            + s_short * line(
                lam, lam0 - 0.5 * a.stopband, fwhm, a.shape
            )
        )

    if a.ase_db is not None:
        I = I + I.max() * 10.0 ** (-a.ase_db / 10.0) * gaussian(
            lam, lam0, a.ase_fwhm
        )

    # Zero anchors just outside the window so interpolation beyond it gives 0.
    lam = np.concatenate(([lam[0] - step], lam, [lam[-1] + step]))
    I = np.concatenate(([0.0], I, [0.0]))

    return lam, I


def measured_fwhm(lam, I):
    """FWHM of the outer half-maximum crossings (envelope width for fp)."""
    above = np.where(I >= 0.5 * I.max())[0]
    i0, i1 = above[0], above[-1]

    l = np.interp(
        0.5 * I.max(),
        [I[i0 - 1], I[i0]],
        [lam[i0 - 1], lam[i0]]
    )

    r = np.interp(
        0.5 * I.max(),
        [I[i1 + 1], I[i1]],
        [lam[i1 + 1], lam[i1]]
    )

    return r - l


def write_oghma(path, lam_m, I, title):
    header = {
        "title": title,
        "type": "xy",
        "y_label": "Wavelength",
        "data_label": "Intensity",
        "y_units": "nm",
        "y_mul": 1000000000.0,
        "data_units": "m^{-1}.W.m(^-2)",
        "icon": "spectra_file",
        "time ": 0.0,
        "Vexternal": 0.0,
        "x_len": 1,
        "y_len": int(len(lam_m)),
        "z_len": 1,
        "cols": "yd",
    }

    with open(path, "w") as f:
        f.write(
            "#oghma_csv "
            + json.dumps(header, separators=(",", ":"))
            + "*\n"
        )

        for l, v in zip(lam_m, I):
            f.write(f"{l:e}\t{v:e}\n")


def plot(path, lam_nm, I):
    import matplotlib
    matplotlib.use("Agg")
    import matplotlib.pyplot as plt

    fig, (ax1, ax2) = plt.subplots(
        2, 1, figsize=(7, 6), sharex=True
    )

    ax1.plot(lam_nm, I)
    ax1.set_ylabel("W m$^{-2}$ m$^{-1}$")

    db = 10.0 * np.log10(
        np.maximum(I, I.max() * 1e-8) / I.max()
    )

    ax2.plot(lam_nm, db)
    ax2.set_ylabel("dB rel. peak")
    ax2.set_xlabel("Wavelength (nm)")
    ax2.set_ylim(-80, 3)

    fig.tight_layout()
    fig.savefig(path, dpi=120)


def main():
    p = argparse.ArgumentParser(
        description=__doc__,
        formatter_class=argparse.RawDescriptionHelpFormatter
    )

    p.add_argument(
        "outdir", nargs="?", default="laser_1550nm",
        help="output directory (default: laser_1550nm)"
    )

    p.add_argument(
        "--model",
        choices=["gaussian", "fp", "dfb"],
        default="gaussian"
    )

    p.add_argument(
        "--wavelength",
        type=float,
        default=1550.0,
        help="centre, nm"
    )

    p.add_argument(
        "--fwhm",
        type=float,
        default=2.0,
        help="line FWHM (gaussian/dfb) or gain-envelope FWHM (fp), nm"
    )

    p.add_argument(
        "--power",
        type=float,
        default=10.0,
        help="integrated irradiance, W/m^2"
    )

    p.add_argument(
        "--shape",
        choices=["gaussian", "lorentzian"],
        default="gaussian",
        help="lineshape of each line/mode"
    )

    # Fabry-Perot
    p.add_argument(
        "--cavity-length",
        type=float,
        default=300.0,
        help="fp: cavity length, um"
    )

    p.add_argument(
        "--group-index",
        type=float,
        default=3.6,
        help="fp: group index"
    )

    p.add_argument(
        "--mode-spacing",
        type=float,
        default=None,
        help="fp: override mode spacing, nm"
    )

    p.add_argument(
        "--mode-width",
        type=float,
        default=0.05,
        help="fp: width of each mode, nm"
    )

    # DFB
    p.add_argument(
        "--smsr",
        type=float,
        default=40.0,
        help="dfb: side-mode suppression, dB"
    )

    p.add_argument(
        "--stopband",
        type=float,
        default=3.0,
        help="dfb: stopband width, nm"
    )

    p.add_argument(
        "--side-asym",
        type=float,
        default=6.0,
        help="dfb: extra suppression of the short-wavelength side mode, dB"
    )

    # ASE pedestal
    p.add_argument(
        "--ase-db",
        type=float,
        default=None,
        help="add ASE pedestal this many dB below peak (off by default)"
    )

    p.add_argument(
        "--ase-fwhm",
        type=float,
        default=40.0,
        help="ASE pedestal FWHM, nm"
    )

    # Grid
    p.add_argument(
        "--step",
        type=float,
        default=None,
        help="grid step, nm (auto)"
    )

    p.add_argument(
        "--window",
        type=float,
        default=None,
        help="total grid width, nm (auto)"
    )

    p.add_argument(
        "--plot",
        action="store_true",
        help="also save a PNG preview"
    )

    a = p.parse_args()

    lam_nm, I = build(a)
    lam_m = lam_nm * 1e-9

    I = I * a.power / _trapz(I, lam_m)

    title = (
        f"{a.model} laser {a.wavelength:g} nm, "
        f"FWHM {a.fwhm:g} nm, {a.power:g} W/m2"
    )

    os.makedirs(a.outdir, exist_ok=True)

    csv_path = os.path.join(a.outdir, "spectra.csv")

    write_oghma(csv_path, lam_m, I, title)

    with open(os.path.join(a.outdir, "data.json"), "w") as f:
        f.write(DATA_JSON)

    print(f"wrote {a.outdir}/data.json")
    print(
        f"wrote {csv_path}: {len(lam_m)} points, "
        f"{lam_nm[0]:.3f}-{lam_nm[-1]:.3f} nm"
    )

    print(
        f"  integral    = {_trapz(I, lam_m):.6g} W/m^2"
    )

    print(
        f"  peak        = {lam_nm[np.argmax(I)]:.4f} nm, "
        f"{I.max():.4e} W/m^2/m"
    )

    if a.model == "fp":
        print(
            f"  mode spacing= {mode_spacing_nm(a):.4f} nm"
        )
    else:
        print(
            f"  meas. FWHM  = {measured_fwhm(lam_nm, I):.4f} nm"
        )

    if a.plot:
        png = os.path.join(a.outdir, "spectra.png")
        plot(png, lam_nm, I)
        print(f"  plot        = {png}")


if __name__ == "__main__":
    main()

#!/usr/bin/env python3
"""Generate an independent library of phenomenological LED emission spectra.

The output is intended for OghmaNano illumination examples. Spectra are not
manufacturer measurements or fits to specific commercial devices.

A specified, illustrative external optical power (mW) is converted into the
on-axis irradiance of a point-like Lambertian LED, 1 cm from its emitting face:

    E [W/m^2] = (P_mW * 1e-3) / (pi * distance_m**2)

Each spectra.csv stores wavelength (metres) and spectral irradiance (W/m^3).
Trapezoidal integration over wavelength in metres returns E in W/m^2.

Dependencies: numpy (no SciPy needed).
Usage: python3 oghma_led_material_spectra_generator.py [output_directory]
"""

import csv
import json
import math
import sys
from pathlib import Path

import numpy as np

DISTANCE_M = 0.01
TEMPERATURE_K = 300.0
K_B_EV_K = 8.617333262145e-5
HC_EV_NM = 1239.841984

# The wavelengths, FWHMs, and powers are selected independently as useful
# representative examples of material families. They are ILLUSTRATIVE, not
# measurements or specifications of particular commercial LED devices.
# 'band_edge' is an asymmetric, broadened bulk-like phenomenological curve.
# 'localized' is a mildly asymmetric broadened emission-centre curve, used for
# GaP:N and GaAs:Si, which are not well represented by direct-gap band-edge EL.
LEDS = [
    # peak_nm, material,             emission label,             FWHM_nm, power_mW, shape
    (275,  'AlGaN',                  'UV-C',                       12,      0.2, 'band_edge'),
    (310,  'AlGaN',                  'UV-B',                       14,      0.4, 'band_edge'),
    (365,  'GaN/AlGaN',              'Near ultraviolet',           13,      2.0, 'band_edge'),
    (405,  'InGaN/GaN',              'Violet',                     17,      4.0, 'band_edge'),
    (450,  'InGaN/GaN',              'Blue',                       22,     12.0, 'band_edge'),
    (480,  'InGaN/GaN',              'Cyan',                       27,      7.0, 'band_edge'),
    (525,  'InGaN/GaN',              'Green',                      34,      3.0, 'band_edge'),
    (565,  'GaP:N',                  'Green (isoelectronic)',      48,      0.8, 'localized'),
    (590,  'AlGaInP',                'Amber',                      22,      3.0, 'band_edge'),
    (610,  'AlGaInP',                'Orange',                     22,      5.0, 'band_edge'),
    (630,  'AlGaInP',                'Red',                        22,      8.0, 'band_edge'),
    (660,  'AlGaInP',                'Deep red',                   25,      7.0, 'band_edge'),
    (740,  'AlGaAs',                 'Far red',                    30,      5.0, 'band_edge'),
    (780,  'AlGaAs',                 'Near infrared',             34,      8.0, 'band_edge'),
    (850,  'AlGaAs/GaAs',            'Near infrared',             38,     12.0, 'band_edge'),
    (870,  'GaAs/Al0.3Ga0.7As',      'GaAs double heterostructure',40,      8.0, 'band_edge'),
    (940,  'GaAs:Si',                'Near infrared',             55,      8.0, 'localized'),
    (1050, 'InGaAs/GaAs',            'Near infrared',             65,      4.0, 'band_edge'),
    (1300, 'InGaAsP/InP',            'Short-wave infrared',      100,      1.0, 'band_edge'),
    (1550, 'InGaAsP/InP',            'Short-wave infrared',      120,      0.5, 'band_edge'),
]

# This example is the exact user-specified layer stack, top to bottom.
# It documents an LED device, not a fully parameterised OghmaNano sim.json.
GAAS_870_DEVICE = {
    'example_name': 'GaAs/AlGaAs double-heterostructure infrared LED',
    'nominal_peak_wavelength_nm': 870,
    'orientation': 'top_to_bottom',
    'device_total_thickness_nm': 800,
    'layers': [
        {
            'material': 'Al0.3Ga0.7As',
            'type': 'p',
            'thickness_nm': 200,
            'acceptor_doping_m^-3': 5e23,
            'donor_doping_m^-3': 0,
            'aluminium_fraction': 0.30,
        },
        {
            'material': 'GaAs',
            'type': 'intrinsic',
            'role': 'active_emitting_layer',
            'thickness_nm': 300,
            'acceptor_doping_m^-3': 0,
            'donor_doping_m^-3': 0,
        },
        {
            'material': 'Al0.3Ga0.7As',
            'type': 'n',
            'thickness_nm': 200,
            'acceptor_doping_m^-3': 0,
            'donor_doping_m^-3': 5e23,
            'aluminium_fraction': 0.30,
        },
        {
            'material': 'GaAs',
            'type': 'n+',
            'role': 'contact_layer',
            'thickness_nm': 100,
            'acceptor_doping_m^-3': 0,
            'donor_doping_m^-3': 1e24,
        },
    ],
    'note': 'Structure-only example; optical spectrum is phenomenological. '
            'An OghmaNano electrical simulation still needs its full model input.',
}


def wavelength_grid(peak_nm, fwhm_nm):
    """Cover each emission band with plenty of margin but no huge zero tails."""
    margin = max(130.0, 8.0 * fwhm_nm)
    lo = max(210.0, peak_nm - margin)
    hi = peak_nm + margin
    return np.arange(lo, hi + 0.125, 0.25)


def energy_profile(shape_type, gaussian_width_ev):
    """Return a spectral-power shape tabulated against relative photon energy."""
    step_ev = 0.001
    relative_ev = np.arange(-1.5, 1.5 + step_ev / 2, step_ev)

    if shape_type == 'band_edge':
        above_edge = np.maximum(relative_ev, 0.0)
        e_thermal = K_B_EV_K * TEMPERATURE_K
        ideal = np.sqrt(above_edge) * np.exp(-above_edge / e_thermal)

        # Convolution with a Gaussian in energy. Fourier convolution is fast
        # and the grid is wide enough for negligible boundary wraparound.
        frequencies = np.fft.rfftfreq(len(relative_ev), step_ev)
        gaussian_fft = np.exp(-0.5 * (2.0 * np.pi * frequencies * gaussian_width_ev) ** 2)
        profile = np.fft.irfft(np.fft.rfft(ideal) * gaussian_fft, n=len(relative_ev))
        profile = np.maximum(profile, 0.0)
    elif shape_type == 'localized':
        # A broadened impurity/isoelectronic-centre line, with a somewhat wider
        # low-photon-energy tail; this is an empirical shape, not a band solver.
        local_sigma = np.where(relative_ev < 0.0,
                               gaussian_width_ev * 1.25,
                               gaussian_width_ev * 0.90)
        profile = np.exp(-0.5 * (relative_ev / local_sigma) ** 2)
    else:
        raise ValueError(f'Unknown shape type: {shape_type}')
    return relative_ev, profile


def spectrum_at_center(wavelength_nm, center_ev, relative_ev, profile):
    """Convert energy spectral density to wavelength spectral density."""
    photon_energy_ev = HC_EV_NM / wavelength_nm
    spectral_ev = np.interp(photon_energy_ev - center_ev, relative_ev, profile,
                            left=0.0, right=0.0)
    return spectral_ev * HC_EV_NM / (wavelength_nm * wavelength_nm)


def fwhm_nm_of(y, wavelength_nm):
    indices = np.flatnonzero(y >= y.max() / 2.0)
    if not len(indices):
        raise ValueError('No valid emission band found')
    return wavelength_nm[indices[-1]] - wavelength_nm[indices[0]]


def tune_peak(peak_nm, wavelength_nm, relative_ev, profile):
    """Shift the chosen energy-domain profile until its wavelength peak matches."""
    target_ev = HC_EV_NM / peak_nm
    lo, hi = target_ev - 0.5, target_ev + 0.3
    for _ in range(20):
        center = (lo + hi) / 2.0
        y = spectrum_at_center(wavelength_nm, center, relative_ev, profile)
        measured_peak_nm = wavelength_nm[np.argmax(y)]
        if measured_peak_nm > peak_nm:
            lo = center  # larger energy shifts the peak toward shorter wavelengths
        else:
            hi = center
    return spectrum_at_center(wavelength_nm, (lo + hi) / 2.0,
                              relative_ev, profile)


def spectrum_for_peak(peak_nm, desired_fwhm_nm, shape_type):
    wavelength_nm = wavelength_grid(peak_nm, desired_fwhm_nm)
    lo, hi = 0.0005, 0.30  # Gaussian energy broadening (eV)
    for _ in range(18):
        mid = (lo + hi) / 2.0
        relative_ev, profile = energy_profile(shape_type, mid)
        y = tune_peak(peak_nm, wavelength_nm, relative_ev, profile)
        if fwhm_nm_of(y, wavelength_nm) < desired_fwhm_nm:
            lo = mid
        else:
            hi = mid
    relative_ev, profile = energy_profile(shape_type, (lo + hi) / 2.0)
    y = tune_peak(peak_nm, wavelength_nm, relative_ev, profile)
    return wavelength_nm, y


def on_axis_irradiance_w_m2(external_power_mw, distance_m=DISTANCE_M):
    if distance_m <= 0.0:
        raise ValueError('distance must be positive')
    return external_power_mw * 1e-3 / (math.pi * distance_m * distance_m)


def write_oghma_csv(path, title, wavelength_m, spectral_w_m3):
    header = {
        'title': title, 'type': 'xy', 'y_label': 'Wavelength',
        'data_label': 'Intensity', 'y_units': 'nm', 'y_mul': 1e9,
        'data_units': 'm^{-1}.W.m(^-2)', 'icon': 'spectra_file',
        'time ': 0.0, 'Vexternal': 0.0,
        'x_len': 1, 'y_len': len(wavelength_m), 'z_len': 1,
        'cols': 'yd',
    }
    with path.open('w', encoding='utf-8') as out:
        out.write('#oghma_csv ' + json.dumps(header, separators=(',', ':')) + '*\n')
        for wl, value in zip(wavelength_m, spectral_w_m3):
            out.write(f'{wl:.12e}\t{value:.12e}\n')


def write_data_json(path):
    # Preserve the supplied OghmaNano spectra-item metadata format.
    metadata = {
        'item_type': 'spectra', 'status': 'public', 'changelog': '',
        'spectra_import': {
            'import_file_path': 'none', 'import_x_combo_pos': 9,
            'import_data_combo_pos': 5, 'import_x_spin': 0,
            'import_data_spin': 2, 'import_title': 'Voltage - J',
            'import_xlabel': 'Voltage', 'import_data_label': 'J',
            'import_area': 0.104, 'import_data_invert': 'False',
            'import_x_invert': 'False', 'data_file': 'none',
        },
    }
    path.write_text(json.dumps(metadata, indent='\t') + '\n', encoding='utf-8')


def generate(root):
    root.mkdir(parents=True, exist_ok=True)
    manifest = []
    wavelengths_seen = set()

    for peak_nm, material, label, width_nm, power_mw, shape_type in LEDS:
        if peak_nm in wavelengths_seen:
            raise ValueError(f'Duplicate peak wavelength/folder: {peak_nm}nm')
        wavelengths_seen.add(peak_nm)

        folder = root / f'{peak_nm}nm'
        folder.mkdir(parents=True, exist_ok=True)
        wave_nm, shape = spectrum_for_peak(peak_nm, width_nm, shape_type)
        wave_m = wave_nm * 1e-9
        area = np.trapezoid(shape, wave_m) if hasattr(np, 'trapezoid') else np.trapz(shape, wave_m)
        if area <= 0.0:
            raise ValueError(f'Unable to normalise {peak_nm} nm spectrum')
        target_w_m2 = on_axis_irradiance_w_m2(power_mw)
        spectral_w_m3 = shape * (target_w_m2 / area)
        actual_area = (np.trapezoid(spectral_w_m3, wave_m)
                       if hasattr(np, 'trapezoid') else np.trapz(spectral_w_m3, wave_m))

        title = (f'{material} {label} LED {peak_nm} nm '
                 f'at {DISTANCE_M * 100:g} cm {target_w_m2:.4f} Wm^{{-2}}')
        write_oghma_csv(folder / 'spectra.csv', title, wave_m, spectral_w_m3)
        write_data_json(folder / 'data.json')

        if peak_nm == 870:
            (folder / 'device_structure.json').write_text(
                json.dumps(GAAS_870_DEVICE, indent=2) + '\n', encoding='utf-8')
            (folder / 'DEVICE_STRUCTURE.txt').write_text(
                'GaAs/AlGaAs double-heterostructure infrared LED (~870 nm)\n'
                'Top to bottom:\n'
                '1. p-type Al0.3Ga0.7As, 200 nm, NA = 5e23 m^-3, Al = 0.30\n'
                '2. Intrinsic GaAs (emitting active layer), 300 nm, undoped\n'
                '3. n-type Al0.3Ga0.7As, 200 nm, ND = 5e23 m^-3, Al = 0.30\n'
                '4. n+ GaAs (contact layer), 100 nm, ND = 1e24 m^-3\n'
                'Total thickness: 800 nm.\n'
                'This is a layer specification, not a complete OghmaNano sim.json.\n',
                encoding='utf-8',
            )

        manifest.append({
            'folder': folder.name, 'peak_wavelength_nm': peak_nm,
            'material_system': material, 'emission': label,
            'nominal_fwhm_nm': width_nm,
            'lineshape': shape_type,
            'assumed_external_optical_power_mw': power_mw,
            'reference_distance_m': DISTANCE_M,
            'angular_emission_assumption': 'point Lambertian, forward hemisphere',
            'integrated_irradiance_w_m2': float(actual_area),
            'device_structure_file': 'device_structure.json' if peak_nm == 870 else None,
        })
        print(f'{peak_nm:4d} nm  {material:20s} {power_mw:5.1f} mW'
              f' -> {actual_area:8.3f} W/m^2'
              f'  FWHM {fwhm_nm_of(shape, wave_nm):6.1f} nm')

    (root / 'manifest.json').write_text(json.dumps(manifest, indent=2) + '\n', encoding='utf-8')

    with (root / 'catalogue.csv').open('w', newline='', encoding='utf-8') as out:
        writer = csv.writer(out)
        writer.writerow(['folder', 'material_system', 'emission', 'peak_nm',
                         'fwhm_nm', 'external_power_mw', 'irradiance_at_1cm_w_m2', 'shape'])
        for m in manifest:
            writer.writerow([m['folder'], m['material_system'], m['emission'],
                             m['peak_wavelength_nm'], m['nominal_fwhm_nm'],
                             m['assumed_external_optical_power_mw'],
                             f"{m['integrated_irradiance_w_m2']:.8f}", m['lineshape']])

    (root / 'README.md').write_text(
        '# Independent phenomenological LED emission examples\n\n'
        '20 representative material systems / emitting bands from 275 to 1550 nm.\n'
        'These spectra are analytically generated, not copied or fitted from any\n'
        'manufacturer spectrum. Chosen wavelengths, FWHMs and powers are\n'
        'illustrative example operating points, not device specifications.\n\n'
        'Each `xxxnm/` contains `spectra.csv` and `data.json`.\n'
        '`870nm/` additionally contains the specified 800 nm-thick GaAs/AlGaAs\n'
        'double-heterostructure device layer stack. That JSON describes its\n'
        'layers and doping, not a complete electrical simulation input.\n\n'
        'The spectral file uses wavelength in **metres** and spectral irradiance\n'
        'in **W/m^3**. Numerically integrating with respect to wavelength in\n'
        'metres gives irradiance in **W/m^2**.\n\n'
        '## Reference geometry\n\n'
        'A point Lambertian emitter radiates its assumed external optical power\n'
        'into the forward hemisphere. At 1 cm directly on-axis:\n\n'
        '`E [W/m^2] = P_external [W] / (pi * (0.01 m)^2)`.\n\n'
        'This ignores packaging optics, near-field geometry, finite chip size,\n'
        'source orientation, drive current and temperature dependence. In\n'
        'particular, edge-emitting infrared sources may not be Lambertian;\n'
        'these spectra are example incident-source conventions only.\n\n'
        '## Phenomenological lineshapes\n\n'
        'Most bands: `(E-Eg)^0.5 exp[-(E-Eg)/kBT]` above the effective edge,\n'
        'Gaussian-broadened in photon energy and converted to wavelength.\n'
        'The GaP:N and GaAs:Si bands use a phenomenological broadened, mildly\n'
        'asymmetric localized-emission-centre shape instead. Width is tuned\n'
        'numerically to the per-example target FWHM and all curves are scaled\n'
        'to the listed irradiance. These are not predictions from alloy\n'
        'composition or band-structure calculations.\n',
        encoding='utf-8',
    )


if __name__ == '__main__':
    generate(Path(sys.argv[1]) if len(sys.argv) > 1 else Path('led_material_spectra_20'))