#!/usr/bin/env python3
"""Reanalyse published tabular facts. Python 3 standard library; no network."""
import csv
import hashlib
import json
from pathlib import Path

ROOT = Path(__file__).resolve().parent
H = 6.62607015e-34
C = 299792458
E = 1.602176634e-19

def photon_energy_ev(wavelength_nm):
    return H * C / (wavelength_nm * 1e-9) / E

def analyse(data):
    a, b = data['emitters']
    derived_devices = []
    for device in data['devices']:
        assert 0 < device['eqe1000Percent'] <= device['eqeMaxPercent'] <= 100
        retention = device['eqe1000Percent'] / device['eqeMaxPercent'] * 100
        derived_devices.append({
            'id': device['id'],
            'retainedEqe1000PercentOfPeak': retention,
            'rollOff1000PercentOfPeak': 100 - retention,
            'eqeLossPercentagePoints': device['eqeMaxPercent'] - device['eqe1000Percent'],
        })
    baseline = data['devices'][2]['lifetime']
    exponent = baseline['assumedAccelerationExponent']
    measured_luminance = baseline['initialLuminanceCdM2']
    measured_hours = baseline['measuredHours']
    prediction = lambda n, luminance: measured_hours * (measured_luminance / luminance) ** n
    result = {
        'schemaVersion': 1,
        'inputSha256': hashlib.sha256((ROOT/'data.json').read_bytes()).hexdigest(),
        'kind': 'Independent arithmetic on published measurements; no new experiment or fitted device model.',
        'formulas': {
            'retainedEqe1000': '100 * EQE(1000 cd m^-2) / EQEmax',
            'rollOff1000': '100 - retainedEqe1000',
            'photonEnergy': 'h*c/(lambda*e), exact SI h,c,e',
            'lifetimeExtrapolation': 'LT80(Ltarget) = LT80(Lmeasured) * (Lmeasured/Ltarget)^n',
        },
        'devices': derived_devices,
        'pairedMoleculeChanges': {
            'from': a['id'], 'to': b['id'],
            'plPeakShiftNm': b['plPeakNm'] - a['plPeakNm'],
            'plFwhmChangeNm': b['plFwhmNm'] - a['plFwhmNm'],
            'plqyChangePercentagePoints': 100 * (b['plqy'] - a['plqy']),
            'emissionPhotonEnergyEv': [photon_energy_ev(row['plPeakNm']) for row in (a,b)],
            'emissionPhotonEnergyChangeEv': photon_energy_ev(b['plPeakNm'])-photon_energy_ev(a['plPeakNm']),
            'delayedLifetimeChangePercent': 100 * (b['delayedLifetimeUs']/a['delayedLifetimeUs']-1),
            'riscRateRatio': b['ratesPerSecond']['risc']/a['ratesPerSecond']['risc'],
            'deltaESTChangeMeV': b['deltaESTMeV']-a['deltaESTMeV'],
            'photonEnergyCaveat': 'Energy of a photon at the PL maximum; not a measured HOMO-LUMO or transport gap.',
        },
        'architectureComparison': {
            'from': 'TADF-D2', 'to': 'PSF-D3',
            'sameTerminalEmitter': b['id'],
            'retentionChangePercentagePoints': derived_devices[2]['retainedEqe1000PercentOfPeak']-derived_devices[1]['retainedEqe1000PercentOfPeak'],
            'eqe1000ChangePercentagePoints': data['devices'][2]['eqe1000Percent']-data['devices'][1]['eqe1000Percent'],
            'caveat': 'Sensitizer, hosts, transport layers and emissive-layer thickness differ. This comparison does not isolate a causal sensitizer effect.',
        },
        'lifetimeExtrapolation': {
            'deviceId': 'PSF-D3', 'metric': 'LT80',
            'measuredHours': measured_hours, 'measuredInitialLuminanceCdM2': measured_luminance,
            'targetInitialLuminanceCdM2': 100, 'paperAssumedExponent': exponent,
            'recomputedHoursAt100': prediction(exponent,100), 'paperRoundedEstimateHoursAt100': baseline['estimatedHoursAt100'],
            'sensitivityLabel': 'Illustrative assumed exponents; not a confidence interval and not measured or refitted.',
            'sensitivity': [{'assumedExponent':n,'estimatedHoursAt100':prediction(n,100)} for n in [1.4,1.6,1.8,2.0,2.2]],
            'warning': 'Only 16 h at 1000 cd m^-2 is measured for this PSF device. The 1010 h value at 100 cd m^-2 is an estimate.',
        },
        'checks': {
            'numberOfMolecules': len(data['emitters']), 'numberOfDevices':len(data['devices']),
            'sourceLifetimeRoundingWithinOneHour': abs(prediction(exponent,100)-baseline['estimatedHoursAt100']) < 1,
            'allMeasuredDevicePointsPositive': all(p['eqePercent']>0 and p['luminanceCdM2']>0 for d in data['devices'] for p in d['eqeAtLuminance']),
        },
    }
    assert result['checks']['sourceLifetimeRoundingWithinOneHour']
    assert result['checks']['allMeasuredDevicePointsPositive']
    assert len(data['emitters'])==2 and len(data['devices'])==3
    return result

def write_csv(name,fields,rows):
    with (ROOT/name).open('w',newline='') as f:
        writer=csv.DictWriter(f,fieldnames=fields);writer.writeheader();writer.writerows(rows)

def main():
    data=json.loads((ROOT/'data.json').read_text())
    result=analyse(data)
    (ROOT/'analysis.json').write_text(json.dumps(result,ensure_ascii=False,indent=2)+'\n')
    fields=['id','name','formula','plPeakNm','plFwhmNm','plqy','promptYield','delayedYield','promptLifetimeNs','delayedLifetimeUs','deltaESTMeV','conditionId','source']
    write_csv('photophysics.csv',fields,[{k:e[k] for k in fields} for e in data['emitters']])
    fields=['id','emitterId','architecture','conditionId','emitterLoadingWtPercent','elPeakNm','elFwhmNm','eqeMaxPercent','eqe1000Percent','source']
    write_csv('devices.csv',fields,[{k:d[k] for k in fields} for d in data['devices']])
    fields=['deviceId','metric','measuredHours','initialLuminanceCdM2','estimatedHoursAt100','assumedAccelerationExponent']
    write_csv('lifetimes.csv',fields,[dict(deviceId=d['id'],**d['lifetime']) for d in data['devices']])
    print(json.dumps({'devices':result['devices'],'lifetime':result['lifetimeExtrapolation']['recomputedHoursAt100'],'checks':result['checks']},indent=2))

if __name__ == '__main__':main()