#!/usr/bin/env python3 """ csdb.py - load the T-matrix Cross Section Database (binary file) into Python dictionaries. REQUIREMENTS Python 3.8+ and NumPy. Nothing else. The file is a plain NumPy ".npz" archive: it holds arrays and one JSON text, no pickled objects, so opening it never executes code. QUICK START import csdb db = csdb.load("csdb_v0.1_full.npz") # the file you downloaded csdb.show(db) # print an overview (a good first check) csdb.show(db, "forsterite", "compact") # print every curve of one species and shape c = csdb.get(db, "forsterite", "compact", 1.0) # one curve: species, shape, radius in micron lam, ext = c["wavelength_um"], c["ext_m2"] # plain NumPy arrays, ready to plot From the command line: python csdb.py csdb_v0.1_full.npz (prints the overview and runs the checks) UNITS (also stored in db["units"]) wavelength_um micron radius_um micron, volume-equivalent radius of the whole grain monomer_radius_um micron, radius of one monomer (aggregates only; NaN for a single sphere) ext_m2, abs_m2, sca_m2 m^2, extinction, absorption and scattering cross sections (ext = abs + sca) (the JSON files of the website use cm^2: 1 m^2 = 1e4 cm^2) positions in geometry micron, centres of the monomers; the grain is centred on its own centre of mass STRUCTURE OF THE DICTIONARY (db = csdb.load(...)) db["release"] {"version": "v0.1", "doi": ..., "released": "2026-11-03"} db["variant"] "full" (everything) or "public" (what is open to everybody) db["units"] the units above db["license"] data licence db["cite"] list of dictionaries {"paper", "doi", "arxiv", "when"}: what to cite, and when db["species"] dictionary keyed by the species slug, for example "forsterite", "spinel", "silicon_monoxide": db["species"]["forsterite"] = { "name": "forsterite", "formula": "Mg2SiO4", "category": "silicate", "nk_source": "Jaeger et al. (2003), ...", # where the refractive index comes from: cite it too "access": "public", "shapes": { # keyed by the shape slug "sphere" | "compact" | "elongated" | "porous": { "label": "Compact aggregate", "kind": "aggregate", "description": "...", "n_curves": 28, # number of sizes "radius_um": array, shape (n,) # one value per size, increasing "n_monomers": array, shape (n,) # 1 for a single sphere "monomer_radius_um": array, shape (n,) "wavelength_um": array, shape (L,) # shared by all the sizes of the group "ext_m2": array, shape (n, L) # row i = size i "abs_m2": array, shape (n, L) "sca_m2": array, shape (n, L) "geometry": list of n items # None, or array (n_monomers, 3) in micron "code_version": list of n strings # code that computed each curve } } } So the extinction of the 1 micron compact aggregate of forsterite is g = db["species"]["forsterite"]["shapes"]["compact"] i = list(g["radius_um"]).index(1.0) # or use csdb.get(...), which finds the nearest size for you g["ext_m2"][i] # array of L values, one per wavelength FUNCTIONS load(path) -> the dictionary described above show(db) -> print the overview and run the consistency checks (what is inside, ranges, problems) show(db, species, shape) -> print one table with every size of that species and shape get(db, species, shape, radius_um, tol=0.01) -> one curve as a flat dictionary (radius_um, wavelength_um, ext_m2, abs_m2, sca_m2, n_monomers, monomer_radius_um, geometry, code_version); the size nearest to radius_um, raising an error if none is within a relative tolerance tol curves(db) -> iterate over every curve as (species, shape, flat dictionary) CITATION If you use these data, please cite the papers in db["cite"] (print them with csdb.show_cite(db)), the refractive-index source of each species (db["species"][...]["nk_source"]) and the release you used (db["release"]). """ import json import sys import numpy as np FORMAT = "csdb-npz" __all__ = ["load", "show", "show_cite", "get", "curves"] def load(path): """Read the .npz file and return the dictionary described at the top of this file.""" try: z = np.load(path, allow_pickle=False) meta = json.loads(bytes(z["meta_json"]).decode("utf-8")) except (OSError, ValueError, EOFError, KeyError) as e: raise ValueError("%s is not a database file of this site (%s: %s); download the .npz file from the Downloads page" % (path, type(e).__name__, e)) from None with z: if meta.get("format") != FORMAT: raise ValueError("%s is not a %s file" % (path, FORMAT)) if meta.get("format_version", 0) > 1: raise ValueError("this file was written by a newer version of the database (format %s): download the latest csdb.py" % meta["format_version"]) db = {k: meta[k] for k in ("release", "variant", "units", "license", "cite") if k in meta} db["species"] = {slug: dict(info, shapes={}) for slug, info in meta["species"].items()} for g in meta["groups"]: p = g["prefix"] n = int(g["n"]) geo_n = z[p + "_geo_n"] geo_pos = z[p + "_geo_pos"] ends = np.cumsum(geo_n) geometry = [geo_pos[e - m:e].copy() if m else None for m, e in zip(geo_n, ends)] shape = dict(meta["shapes"][g["shape"]], n_curves=n) shape.update( radius_um=z[p + "_radius_um"], n_monomers=z[p + "_n_monomers"], monomer_radius_um=z[p + "_monomer_radius_um"], wavelength_um=z[p + "_wavelength_um"], ext_m2=z[p + "_ext_m2"], abs_m2=z[p + "_abs_m2"], sca_m2=z[p + "_sca_m2"], geometry=geometry, code_version=g["code_version"]) db["species"][g["species"]]["shapes"][g["shape"]] = shape return db def curves(db): """Iterate over every curve: (species slug, shape slug, flat dictionary).""" for sp, s in db["species"].items(): for sh, g in s["shapes"].items(): for i in range(g["n_curves"]): yield sp, sh, _curve(g, i) def _curve(g, i): return {"radius_um": float(g["radius_um"][i]), "n_monomers": int(g["n_monomers"][i]), "monomer_radius_um": float(g["monomer_radius_um"][i]), "wavelength_um": g["wavelength_um"], "ext_m2": g["ext_m2"][i], "abs_m2": g["abs_m2"][i], "sca_m2": g["sca_m2"][i], "geometry": g["geometry"][i], "code_version": g["code_version"][i]} def get(db, species, shape, radius_um, tol=0.01): """One curve (flat dictionary) for a species, a shape and a radius in micron: the nearest size, within a relative tolerance.""" try: g = db["species"][species]["shapes"][shape] except KeyError: have = {sp: sorted(s["shapes"]) for sp, s in db["species"].items()} raise KeyError("no %r / %r in this file; available: %s" % (species, shape, have)) from None r = g["radius_um"] i = int(np.argmin(np.abs(r - radius_um))) if abs(r[i] - radius_um) > tol * radius_um: raise ValueError("no size within %.0f%% of %g um for %s / %s; available radii (um): %s" % (tol * 100, radius_um, species, shape, ", ".join("%g" % x for x in r))) return _curve(g, i) def _problems(db): """Consistency checks of the whole file; returns a list of messages (empty = all good).""" bad = [] for sp, sh, c in curves(db): tag = "%s/%s r=%g" % (sp, sh, c["radius_um"]) arrs = (c["ext_m2"], c["abs_m2"], c["sca_m2"]) if not all(np.all(np.isfinite(a)) for a in arrs): bad.append("%s: non-finite values" % tag) elif min(a.min() for a in arrs) < 0: bad.append("%s: negative cross section" % tag) else: dev = np.max(np.abs(c["ext_m2"] - c["abs_m2"] - c["sca_m2"]) / c["ext_m2"]) if dev > 1e-3: bad.append("%s: extinction differs from absorption + scattering by up to %.1e" % (tag, dev)) if c["geometry"] is not None and len(c["geometry"]) != c["n_monomers"]: bad.append("%s: %d monomer positions for %d monomers" % (tag, len(c["geometry"]), c["n_monomers"])) return bad def show(db, species=None, shape=None): """Print an overview of the file (and run the checks); with species and shape, print one table of its sizes.""" rel = db["release"] if species is not None: g = db["species"][species]["shapes"][shape] if shape else None if g is None: raise ValueError("give a shape as well, for example show(db, %r, 'compact')" % species) s = db["species"][species] print("%s (%s) - %s, n,k from: %s" % (s["name"], s["formula"], g["label"], s["nk_source"])) print("%d sizes; wavelength %g-%g um, %d points" % (g["n_curves"], g["wavelength_um"][0], g["wavelength_um"][-1], len(g["wavelength_um"]))) print("%10s %9s %13s %14s %14s %s" % ("radius_um", "monomers", "monomer_um", "max ext (m2)", "max sca/ext", "geometry")) for i in range(g["n_curves"]): e = g["ext_m2"][i] print("%10g %9d %13s %14.3e %14.3f %s" % (g["radius_um"][i], g["n_monomers"][i], "-" if np.isnan(g["monomer_radius_um"][i]) else "%g" % g["monomer_radius_um"][i], e.max(), np.max(g["sca_m2"][i] / e), "yes" if g["geometry"][i] is not None else "-")) return print("T-matrix Cross Section Database, release %s (%s), %s data, released %s" % (rel["version"], rel.get("doi") or "DOI pending", db["variant"], rel.get("released"))) print("licence: %s" % db["license"]) total = 0 print("%-18s %-10s %-20s %6s %-19s %-15s" % ("species", "formula", "shape", "curves", "radius_um", "wavelength_um")) for sp, s in db["species"].items(): for sh, g in s["shapes"].items(): total += g["n_curves"] print("%-18s %-10s %-20s %6d %-19s %-15s" % (sp, s["formula"], sh, g["n_curves"], "%g - %g" % (g["radius_um"][0], g["radius_um"][-1]), "%g - %g" % (g["wavelength_um"][0], g["wavelength_um"][-1]))) print("%d species, %d curves in total" % (len(db["species"]), total)) bad = _problems(db) print("checks: %s" % ("all good (finite, positive, ext = abs + sca within 0.1%, geometry consistent)" if not bad else "%d problem(s)" % len(bad))) for b in bad[:20]: print(" -", b) def show_cite(db): """Print what to cite.""" print("Please cite:") for c in db["cite"]: print(" - %s" % c["paper"]) print(" %s" % c["when"]) print(" %s" % ("https://doi.org/" + c["doi"] if c.get("doi") else "arXiv:" + c["arxiv"])) print(" - the refractive-index source of each species you use (db['species'][slug]['nk_source'])") print(" - the release: %s%s" % (db["release"]["version"], " (DOI %s)" % db["release"]["doi"] if db["release"].get("doi") else "")) if __name__ == "__main__": if len(sys.argv) < 2: sys.exit(__doc__) d = load(sys.argv[1]) show(d) print() show_cite(d) # ---------------------------------------------------------------------------------------------------------------------- # QUICK EXAMPLE - copy these lines into your own script (next to this file), remove the leading "# " and run it: # the absorption cross section of forsterite, for a 2.5 micron particle, at a wavelength of 5 micron, for every shape. # # import numpy as np # import csdb # # db = csdb.load("csdb_v0.1_full.npz") # the file you downloaded # for shape in db["species"]["forsterite"]["shapes"]: # sphere, compact, elongated, porous # try: # c = csdb.get(db, "forsterite", shape, 2.5) # the curve of the 2.5 micron particle # except ValueError: # continue # this shape has no 2.5 micron size # i = int(np.argmin(np.abs(c["wavelength_um"] - 5.0))) # index of the wavelength nearest to 5 micron # print(f"{shape:10s} lambda = {c['wavelength_um'][i]:.2f} um abs = {c['abs_m2'][i]:.4e} m^2") # # For absorption use c["abs_m2"], for scattering c["sca_m2"], for extinction c["ext_m2"]. # ----------------------------------------------------------------------------------------------------------------------