#!/usr/bin/env python3
"""Reproduce Mol2Mat's oligophenyl analysis, offline, using only Python's stdlib.

Download this script, manifest.json and the three raw/ files into one folder.
Run: python3 analyse.py
Source: PhotochemCAD, Lindsey laboratory, NC State University.
No fitting, smoothing, baseline correction or spectral simulation is applied.
"""
import csv
import hashlib
import json
from pathlib import Path

ROOT = Path(__file__).resolve().parent
# Exact SI constants. Photon energy in eV for wavelength in nm.
HC_EV_NM = 6.62607015e-34 * 299792458 / 1.602176634e-19 * 1e9


def analyse():
    manifest = json.loads((ROOT / 'manifest.json').read_text())
    result = dict(manifest, hc_ev_nm=HC_EV_NM, records=[])
    for record in manifest['records']:
        raw = (ROOT / 'raw' / f"{record['id']}.absorption.txt").read_bytes()
        if hashlib.sha256(raw).hexdigest() != record['sha256']:
            raise ValueError(f"Source checksum mismatch: {record['id']}")
        lines = raw.decode().splitlines()
        # B01 includes a normalised middle column. Always use the final epsilon column.
        points = [(float(row[0]), float(row[-1])) for row in csv.reader(lines[1:], delimiter='\t') if row]
        if any(b[0] <= a[0] for a, b in zip(points, points[1:])):
            raise ValueError('Wavelengths must increase strictly')
        low, high = manifest['peak_window_nm']
        peak = max((p for p in points if low <= p[0] <= high), key=lambda p: p[1])
        start, end = manifest['chart_window_nm']
        # Keep every supplied sample within the plot window (0.25 nm spacing).
        chart = [p for p in points if start <= p[0] <= end]
        result['records'].append(dict(record, sample_count=len(points),
                                      peak_nm=peak[0], peak_epsilon=peak[1],
                                      photon_energy_ev=HC_EV_NM / peak[0], spectrum=chart))
    first, middle, last = result['records']
    result['comparison'] = {
        'shift_nm': last['peak_nm'] - first['peak_nm'],
        'energy_drop_ev': first['photon_energy_ev'] - last['photon_energy_ev'],
        'step_shifts_nm': [middle['peak_nm'] - first['peak_nm'], last['peak_nm'] - middle['peak_nm']],
    }
    return result


if __name__ == '__main__':
    result = analyse()
    (ROOT / 'analysis.json').write_text(json.dumps(result, separators=(',', ':')) + '\n')
    with (ROOT / 'peaks.csv').open('w', newline='') as stream:
        columns = ['id', 'name', 'rings', 'solvent', 'peak_nm', 'peak_epsilon', 'photon_energy_ev',
                   'temperature', 'concentration', 'source_url', 'raw_url', 'sha256']
        writer = csv.DictWriter(stream, fieldnames=columns, extrasaction='ignore')
        writer.writeheader()
        writer.writerows(result['records'])
    for record in result['records']:
        print(f"{record['name']}: {record['peak_nm']:.2f} nm; {record['photon_energy_ev']:.3f} eV")