Examples Overview

This page collects short, self-contained recipes for the most common GASFIR use cases. Longer worked examples live on the Basic Usage, Advanced Usage, Parameter Fitting, and Performance Optimization pages.

Minimal end-to-end calculation

Create a pulse, load parameters, compute the ionization probability:

import numpy as np
from gasfir import (
    create_pulse,
    get_parameters,
    get_diabatic_ionization_probability,
    get_quasi_static_rate_for_field,
)

laser  = create_pulse(800, 1e14, 0, 30)     # 800 nm, 1e14 W/cm², CEP 0, 30 cycles
params = get_parameters("He_Hacc5_QS")

# Non-adiabatic (GASFIR) probability
prob = get_diabatic_ionization_probability(pulse=laser, param_dict=params)
print(f"P_NA = {prob:.4e}")

# Quasi-static rate over a range of field amplitudes
fields = np.linspace(0.01, 0.15, 100)       # atomic units
rates  = get_quasi_static_rate_for_field(fields, params)

cos8 pulse (the default envelope)

gasfir.create_pulse() builds a linearly polarized pulse with a cosN intensity envelope; N defaults to 8. The duration is given in optical cycles:

from gasfir import create_pulse, get_parameters, get_diabatic_ionization_probability

# 800 nm, 1e14 W/cm², CEP = 0, FWHM = 6 cycles, cos^8 envelope (N=8 default)
laser  = create_pulse(800, 1e14, 0, 6, N=8)
params = get_parameters("Hydrogen_SFA")

prob = get_diabatic_ionization_probability(pulse=laser, param_dict=params)

The time grid and fields are available through the pulse object:

t = laser.get_tgrid(dt=0.25)        # time grid (a.u.)
E = laser.get_electric_field(t)     # electric field E(t)
A = laser.get_vector_potential(t)   # vector potential A(t)

Two-colour field with SumOfPulses

gasfir.pulse.SumOfPulses superposes any number of pulses (e.g. a fundamental and its second harmonic). The result behaves like any other Pulse:

from gasfir import create_pulse, get_parameters, get_diabatic_ionization_probability
from gasfir.pulse import SumOfPulses

fundamental = create_pulse(800, 1e14, 0, 6)     # 800 nm
harmonic    = create_pulse(400, 2e13, 0, 6)     # 400 nm (SHG)
two_colour  = SumOfPulses([fundamental, harmonic])

prob = get_diabatic_ionization_probability(
    pulse=two_colour, param_dict=get_parameters("H_SFA")
)

Pump–probe scan with TIPTOE (H_diabatic)

gasfir.TIPTOE holds a strong pump fixed and scans a weak probe in delay, producing a SumOfPulses at each delay. return_delay_array() returns a Nyquist-sampled delay grid:

import numpy as np
from gasfir import (
    TIPTOE,
    create_pulse,
    get_parameters,
    get_diabatic_ionization_probability,
)

pump  = create_pulse(800, 3e14, 0, 3)   # strong pump
probe = create_pulse(800, 1e13, 0, 1)   # weak probe
scanner = TIPTOE(pump, probe)

delays = scanner.return_delay_array()   # Nyquist-sampled (oversample=4 default)
params = get_parameters("H_diabatic")

signal = np.array([
    get_diabatic_ionization_probability(
        pulse=scanner.at_delay(tau), param_dict=params
    )
    for tau in delays
])
# `signal` vs `delays` is the simulated TIPTOE trace

Increase oversample for finer delay resolution:

delays_fine = scanner.return_delay_array(oversample=8)

Numerical field with DataPulse

gasfir.DataPulse wraps a numerically sampled field so it can be used anywhere a Pulse is accepted. Supply either the electric field E(t) or the vector potential A(t) on a time grid in atomic units. A(t) must vanish at both ends of the window:

import numpy as np
from gasfir import DataPulse, get_parameters, get_diabatic_ionization_probability

t      = np.linspace(-400, 400, 8001)            # atomic units
omega  = 0.057                                    # ~800 nm carrier (a.u.)
env    = np.exp(-(t / 150.0) ** 2)                # Gaussian, vanishes at ends
A      = env * np.sin(omega * t) / omega          # vector potential A(t)

pulse = DataPulse(t=t, A=A)                       # E(t) derived as -dA/dt

prob = get_diabatic_ionization_probability(
    pulse=pulse, param_dict=get_parameters("H_SFA")
)

If you only have the electric field, pass E= instead and GASFIR integrates it to obtain A (make sure the field switches on/off smoothly so the resulting A returns to zero):

pulse = DataPulse(t=t, E=my_field_array)

Retrieval strategies and pipelines (gasfir[retrieval])

The three-phase pipeline is global search → least-squares polish → emcee MCMC, driven by gasfir.retrieval.retrieve() and configured with gasfir.retrieval.RetrievalConfig. The training data is a DataFrame with "pulses" and "Y" (target probabilities) columns.

Full pipeline, known medium — initial guess and bounds from the stored reference parameters:

from gasfir.retrieval import RetrievalConfig, retrieve

cfg = RetrievalConfig(medium_name="Hydrogen_SFA", emcee_nsteps=3000)
result = retrieve(data_NA=data, config=cfg)
# outputs (JSON, HDF5 chain, corner/trace plots, .tex) saved to ./Hydrogen_SFA/

Switch the global search strategy — CMA-ES (default) or differential evolution:

cfg = RetrievalConfig(
    medium_name="Hydrogen_SFA",
    global_method="differential_evolution",   # default is "cma"
    de_maxiter=2000,
)

Run only selected phases via run_phases=(global, polish, mcmc). For example, a quick CMA-ES + least-squares fit with no MCMC:

result = retrieve(data_NA=data, config=cfg, run_phases=(True, True, False))

…or extend an existing MCMC chain (skip global search and polish, resume emcee from the saved HDF5 backend):

cfg = RetrievalConfig(medium_name="Hydrogen_SFA", emcee_nsteps=6000)
result = retrieve(data_NA=data, config=cfg, run_phases=(False, False, True))

Quick sanity check with minimal iterations:

cfg = RetrievalConfig(medium_name="Hydrogen_SFA", trial_run=True)

Simultaneous fit to non-adiabatic and quasi-static data — pass a data_QS DataFrame (columns "field" and "Y" = log QS rate):

cfg = RetrievalConfig(medium_name="Hydrogen_SFA", simultaneous=True)
result = retrieve(data_NA=data, data_QS=qs_data, config=cfg)

Cold start for an atomE_g is looked up from the NIST database, no other guess required:

cfg = RetrievalConfig(medium_name="Kr", cma_maxiter=5000)
result = retrieve(data_NA=data, config=cfg)

Cold start for a crystalE_g must be supplied; fit the electron density instead of a probability:

cfg = RetrievalConfig(
    medium_name="MyMaterial",
    initial_params={"E_g": 0.30},
    is_crystal=True,
    ret_electron_density=True,
    cma_maxiter=5000,
)
result = retrieve(data_NA=data, config=cfg)

See Parameter Fitting for fully worked retrieval examples, post-processing of saved chains, and the low-level residual API.