AthermalElectrons <: SimulationTypes
    Enabled::Bool

    AthermalElectron_ElectronCoupling::Bool # Enables coupling to an electronic bath
    AthermalElectron_PhononCoupling::Bool # Enables coupling to a phononic bath
    Conductivity::Bool # Provides conductivity of a ballistic nature using velocity given by v_g
    EmbeddedAthEM::Bool #WIP!! Do not set to true

    ElectronicRelaxation::Symbol # Implementations are Fermi Liquid Theory (:FLT) or constant (:constant)
    PhononicRelaxation::Symbol # Implementations are constant (:constant) or quasiparticle scattering (:quasi)
    ExcitationMatrixElements::Symbol # Implementation is only match internal energy (:unity)
    Conductive_Velocity::Symbol # Implementation of how gorup velocity is calculated, :constant, :fermigas or :effectiveoneband
    
    FE::Union{Float64,Vector{Float64}} # Shifted Fermi energy to the bottom of the valence band for FLT relaxation and group velocity
    τ::Union{Float64,Vector{Float64}} # Material dependent scale-factor for :FLT relaxation time or the constant value for :constant
    τep::Union{Float64,Vector{Float64}} # Constant relaxation time for phonons
    v_g::Union{Vector{Float64},Matrix{Float64}} # Group velocity of electrons calculated assuming a Fermi liquid with μ = FE
end

Struct that defines and holds all values for the propagation of athermal electrons
Enabling this struct assumes an AthEM like system (https://arxiv.org/abs/2503.09479) so can be coupled to electronic
and phononic thermal baths. Coupling implicitly assumes the other system is enabled.

WIP!!!
DensityMatrix <: SimulationTypes
    Enabled::Bool = false
end 

Struct that defines and holds all values for the density matrix propagation.
This Simulation object doesn't function or couple with the others due to the difference in propagation from coupled 
ODE to a von-Neumann equation.

Dimension <: SimulationTypes
    length::Int # The length of the grid, not the depth of the slab
    grid::AbstractArray{Float64} # The grid the simulation is solved over
    spacing::Union{Float64, Vector{Float64}} #The spacing between grid points
    InterfaceHeight::Union{Float64, Vector{Float64}} # Height sorted list of the interfaces between materials
end

Struct that contains all information regarding the spatial grid that the simulation is performed on.

struct ElectronicDistribution <: SimulationTypes
    Enabled::Bool = false

    Electron_PhononCoupling::Bool = false
end
Struct that defines and holds all values for the propagation of an electronic distribution

ElectronicTemperature <: SimulationTypes
    Enabled::Bool = false

    AthermalElectron_ElectronCoupling::Bool = false # Enables coupling to athermal electrons
    Electron_PhononCoupling::Bool = false # Enables coupling to a phonon thermostat
    Conductivity::Bool = false # Provides diffusive thermal conductivity

    ElectronicHeatCapacity::Symbol = :linear # Whether to use linear (:linear), non-linear (:nonlinear) or constant :(constant)
                                             # electronic heat capacity
    ElectronPhononCouplingValue::Symbol = :constant # Whether to use constant (:constant) or variable (:variable)
                                                    # electron phonon coupling

    γ::Union{Float64,Vector{Float64}} = 1.0 # Specific heat capacity of electrons at room temperature for linear heat capacity
    κ::Union{Float64,Vector{Float64}} = 1.0 # Thermal conductivity of electrons at room temperature
    λ::Union{Float64,Vector{Float64}} = 1.0 # Electron-phonon mass enhancement factor for non-linear electron-phonon coupling
    ω::Union{Float64,Vector{Float64}} = 1.0 # Second moment of phonon spectral function for non-linear electron-phonon coupling
    g::Union{Float64,Vector{Float64}} = 1.0 # Constant electron-phonon coupling value 
end

Struct that defines and holds all values for the propagation of an electronic temperature
This can be coupled solely to a thermal phonon bath for a Two-Temperature Model simulation or to athermal electrons
for AthEM propagation with relaxation. Coupling implicitly assumes the other system is enabled.

Fields <: SimulationTypes
    Electric::Expr # Expression for the magnitude of the electric field in the simulation
    Nagnetic::Expr # Expression for the magnitude of the magnetic field in the simulation
end

Struct that contains all information regarding electromagnetic fields in the simulation.

Laser <: SimulationTypes
    envelope::Symbol = :Gaussian # Currently implemented are :Gaussian, :HyperbolicSecant, :Lorentzian and :Rectangular
    Transport::Symbol # :optical, :ballistic and :combined are the options for how the laser decays into a slab

    FWHM::Float64 # The Full-Width Half-Maximum of the laser, for rectnagular half the length
    ϕ::Float64 # The unabsorbed fluence of the laser
    hv::Float64 # The photon frequency of the laser
    ϵ::Union{Float64,Vector{Float64}} # The inverse of the absorption coefficient
    R::Float64 # The reflectivity of the sample
    δb::Union{Float64,Vector{Float64}} # The ballistic length of electrons
end

Struct that contains all laser parameters and any material parameters that affect laser absorption

W.I.P!!!
struct PhononicDistribution <: SimulationTypes
    Enabled::Bool = false

    Electron_PhononCoupling::Bool = false
end
Struct that defines and holds all values for the propagation of a phononic distribution

struct PhononicTemperature <: SimulationTypes
    Enabled::Bool = false

    AthermalElectron_PhononCoupling::Bool = false # Enables coupling to athermal electrons
    Electron_PhononCoupling::Bool = false # Enables coupling to an electron thermostat
    Conductivity::Bool = false # Provides diffusive thermal conductivity

    PhononicHeatCapacity::Symbol = :constant # Whether to use constant (:constant) or non-linear/Simpson's Rule (:nonlinear) 
                                             # Phononic Heat Capacity
    
    θ::Union{Float64,Vector{Float64}} = 1.0 # Debye temperature for non-linear phonon heat capacity
    n::Union{Float64,Vector{Float64}} = 1.0 # Atomic density for non-linear phonon heat capacity
    Cph::Union{Float64,Vector{Float64}} = 1.0 # Constant phonon heat capacity
    κ::Union{Float64,Vector{Float64}} = 1.0 # Constant phonon thermal conductivity
end
Struct that defines and holds all values for the propagation of a phononic temperature
This can be coupled solely to a thermal electronic bath for a Two-Temperature Model simulation or to athermal electrons
for AthEM propagation with phonon-relaxation. Coupling implicitly assumes the other system is enabled.

struct Simulation <: SimulationTypes
    densitymatrix::DensityMatrix
    electronictemperature::ElectronicTemperature
    phononictemperature::PhononicTemperature
    athermalelectrons::AthermalElectrons
    electronicdistribution::ElectronicDistribution
    phononicdistribution::PhononicDistribution
    structure::Structure
    laser::Laser
end
This struct contains all the others and is the main simulation object both in assembling a simulation and during it

SimulationTypes

Parent type of all subtypes in LightMatter.jl

Structure <: SimulationTypes
    Spatial_DOS::Bool # Whether to vary the DOS with height - if so the DOS becomes a vector

    Elemental_System::Int # The number of elemental systems, if > 1 then each constant and vector
                        # of material parameters needs to become a vector of length=Elemental_System

    DOS::Union{spl,Vector{spl},Vector{Vector{spl}}} # The density of states of the simulation
    egrid::Vector{Float64} # An energy grid for electronic or phononic distributions to be solved on

    dimension::Union{Dimension} # A struct holding all spatial grid structure (0D or 1D)
    bandstructure::Union{Vector{<:DataInterpolations.AkimaInterpolation}, Vector{<:Vector{DataInterpolations.AkimaInterpolation}},Nothing} 
                    # The band structure of the simulation both in terms of k->E and E-> k
end

Struct that contains any spatial information including the DOS, the spatial grid to solve the simulation on and
the elemental composition of the simulation (e.g. an antenna-reactor system would contain two elemental systems)

TotalFields <: SimulationTypes
    laser::Fields # The laser fields in the simulation
    external::Fields # The external fields in the simulation, e.g. a magnetic field
end

Struct that contains both the field generated by the laser and any external fields.

spl=Interpolations.AbstractExtrapolation

A convenience type definition to make type specificity easier throughout the code

build_AthermalElectrons(;structure::Structure, Enabled = false, AthermalElectron_ElectronCoupling = false, 
                        AthermalElectron_PhononCoupling = false, Conductivity = false, ElectronicRelaxation = :FLT, 
                        PhononicRelaxation = :constant, ExcitationMatrixElements = :unity, FE=0.0, τ=1.0, τep = 1000.0, 
                        v_g = nothing, Conductive_Velocity = :constant, EmbeddedAthEM = false)

Outer constructor function to assemble the AthermalElectrons struct. Unit conversion is detected on all parameters.
The function will build the group veolcity if one isn't provided by the user.
Defaults allow any unneccessary parameters for users simulation to be ignored.

# Arguments
- 'Enabled': Bool for enabling an athermal electron subssystem
- 'structure': Structure struct, provide if you want the group velocity calculated for you
- 'AthermalElectron_ElectronCoupling': Enables athermal electron - thermal electron coupling
- 'AthermalElectron_PhononCoupling': Enables athermal electron - thermal phonon coupling
- 'Conductivity': Bool for enabling athermal electron ballistic transport
- 'ElectronicRelaxation': Method for calculating athermal electron lifetime due to e-e collisions, :FLT or :constant
- 'PhononicRelaxation': Method for calculating athermal electron lifetime due to e-p collisions, :quasi or :constant
- 'ExcitationMatrixElements': Decides method to calculate excitation matrix elements, only :unity is currently implemented
- 'FE': unit = eV: The Fermi energy defined as the difference between the bottom of the valence band in the DOS and 0.0
- 'τ': unit = fs: A material dependent scalar for the :FLT lifetime or the constant value for :constant e-e lifetime
- 'τep': unit = fs: The constant lifetime for the athermal electrons due to electron-phonon coupling
- 'v_g': unit = nm/fs: The group velocity of the ballistic electrons, for the user to define their own group velocity and will overwrite the 
         one calculated by build_group_velocity. Also the value used for a constant velocity.
- 'Conductive_Velocity': Define the group velocity that build_group_velocity should use, :constant :fermigas, :effectiveoneband
- 'EmbeddedAthEM': Bool for setting only the surface layer to AthEM and the rest a TTM. Can't be used alongside athermal electron
                   transport

# Returns
- The AthermalElectrons struct with the users settings and parameters with any neccessary unit conversion.

WIP!!!
build_DensityMatrix(; Enabled = false)

Once implemented will build a density matrix and store Hamiltonian for propagation via the vonNeumann equation.

build_Dimension(grid=[0.0]::AbstractArray{Float64}, cutoff=0.0::Union{Vector{Float64},Float64})

Outer constructor function to assemble the Dimension struct. The user provides an evenly spaced grid 
and sorted list of interface heights for antenna-reactor complexes. The user must ensure the length of 
cutoff = Elemental_System - 1. No unit conversion is performed when assembling this struct.
Defaults allow any unneccessary parameters for users simulation to be ignored.

# Arguments
- 'grid': unit = nm: Vector representing spatial grid. If [0.0] then homogenous (0D) calculation is performed
- 'cutoff': unit = nm:  Sorted list of all interface heights. Only used when Elemental_System > 1. 

# Returns
- The Dimension struct with the users grid and interface heights

build_ElectronicDistribution(;Enabled = false, Electron_PhononCoupling = false, Ω=1.0, me = Constants.me)

Boltzmann clacluations are currently WIP.
Outer constructor function to assemble the ElectronicDistribution struct. Unit conversion is detected on all parameters.
Defaults allow any unneccessary parameters for users simulation to be ignored. Recommended to use Unitful Quantity
for me due to the unusual mass units in LightMatter.jl

# Arguments
- 'Enabled': Bool for enabling the full electronic distriubtion propagation
- 'Electron_PhononCoupling': Enables electronic distribution - phononic distribution bath coupling
- 'me': unit = eVm : Mass of conduction electron in LightMatter.jl mass units

# Returns
- The ElectronicDistribution struct with the users settings and parameters with any neccessary unit conversion.

build_ElectronicTemperature(; Enabled = false, AthermalElectron_ElectronCoupling = false, Electron_PhononCoupling = false, Conductivity = false,
                           ElectronicHeatCapacity = :linear, ElectronPhononCouplingValue = :constant, γ = 1.0, κ = 1.0, λ = 1.0, ω = 1.0, g = 1.0)

Outer constructor function to assemble the ElectronicTemperature struct. Unit conversion is detected on all parameters.
Defaults allow any unneccessary parameters for users simulation to be ignored.

# Arguments
- 'Enabled': Bool for enabling an thermal electronic bath
- 'AthermalElectron_ElectronCoupling': Enables athermal electron - thermal electron coupling
- 'Electron_PhononCoupling': Enables electron bath - phonon bath coupling
- 'Conductivity': Bool for enabling thermnal electron diffusive transport
- 'ElectronicHeatCapacity': Method for calculating the electronic heat capacity, :linear or :nonlinear
- 'ElectronPhononCouplingValue': Method for calculating the electron phonon coupling value, either :constant or :variable
- 'γ': unit = eV/nm³/K²: Specific heat capacity of electronic bath for :linear ElectronicHeatCapacity
- 'κ': unit = eV/fs/nm/K: Thermal conductivity of electrons at room temperature
- 'λ': unit = unitless: Electron-phonon mass enhancement parameter
- 'ω': unit = eV^2: The second moment of the phonon spectrum
- 'g': unit = eV/fs/nm³/K: Constant value for the electron-phonon coupling if using :constant

# Returns
- The ElectronicTemperature struct with the users settings and parameters with any neccessary unit conversion.

build_laser(;envelope=:Gaussian, FWHM=10.0, ϕ=10.0, hv=5.0, Transport=:optical, ϵ=1.0, R=0.0, δb=1.0)

Outer constructor function to assemble the Laser struct. Can handle unit conversions if the user provides a 
unitful quantity to the laser. 
Defaults allow any unneccessary parameters for users simulation to be ignored.

# Arguments
- 'envelope': Symbol representing the shape of the envelope, :Gaussian, :HyperbolicSecant, :Lorentzian, :Rectangular
- 'FWHM': unit = fs:  Full-Width Half-Maximum of the laser pulse or half the duration of the Rectangular laser ± 0.0
- 'ϕ': unit = eV/nm²: The fluence of the laser 
- 'hv': unit = eV: The photon energy of the laser
- 'Transport': The method of spatial transport of the laser, :optical, :ballistic, :combined
- 'ϵ': unit = nm: The penetration depth of the material (1/α), for the :Optical & :Combined transport
- 'R': unit = unitless:The reflectivity of the sample surface, leave at 0.0 if your provided fluence is the absorbed fluence
- 'δb': unit = nm:The ballistic length of electrons, for :ballistic & :combined transport

# Returns
- The Laser struct with the user settings and neccessary values converted to the correct units

build_PhononicDistribution(;Enabled = false, Electron_PhononCoupling = false, Ω=1.0, me = Constants.me)

Boltzmann clacluations are currently WIP.
Outer constructor function to assemble the PhononicDistribution struct. Unit conversion is detected on all parameters.
Defaults allow any unneccessary parameters for users simulation to be ignored. Recommended to use Unitful Quantity
for me due to the unusual mass units in LightMatter.jl

# Arguments
- 'Enabled': Bool for enabling the full phononic distriubtion propagation
- 'Electron_PhononCoupling': Enables electronic distribution - phononic distribution bath coupling
- 'cs': unit = nm/fs : Speed of sound of longitudinal mode
= 'ED': unit = eV : Deybe energy of the material

# Returns
- The PhononicDistribution struct with the users settings and parameters with any neccessary unit conversion.

build_PhononicTemperature(;Enabled = false, AthermalElectron_PhononCoupling = false, Electron_PhononCoupling = false, 
                           Conductivity = false, PhononicHeatCapacity = :linear, θ = 1.0, n = 1.0, Cph = 1.0, κ = 1.0)

Outer constructor function to assemble the PhononicTemperature struct. Unit conversion is detected on all parameters.
Defaults allow any unneccessary parameters for users simulation to be ignored.

# Arguments
- 'Enabled': Bool for enabling an thermal electronic bath
- 'AthermalElectron_PhononCoupling': Enables athermal electron - thermal electron coupling
- 'Electron_PhononCoupling': Enables electron bath - phonon bath coupling
- 'Conductivity': Bool for enabling thermal phonon diffusive transport
- 'PhononicHeatCapacity': Method for calculating the phononic heat capacity, :constant or :nonlinear
- 'θ': unit = K: Debye temperature of the material
- 'n': unit = atoms/nm³: Float64 of atoms per nm³
- 'Cph': unit = eV/nm³/K: Constant heat capacity for :constant
- 'κ': unit = eV/nm: Constant thermal conductivity of phonons

# Returns
- The PhononicTemperature struct with the users settings and parameters with any neccessary unit conversion.

build_Simulation(;densitymatrix::Union{DensityMatrix,NamedTuple,Nothing}=nothing, electronictemperature::Union{ElectronicTemperature,NamedTuple,Nothing}=nothing,
                       phononictemperature::Union{PhononicTemperature,NamedTuple,Nothing}=nothing, athermalelectrons::Union{AthermalElectrons,NamedTuple,Nothing}=nothing,
                       electronicdistribution::Union{ElectronicDistribution,NamedTuple,Nothing}=nothing, phononicdistribution::Union{PhononicDistribution,NamedTuple,Nothing}=nothing,
                       structure::Union{Structure,NamedTuple,Nothing}=nothing, laser::Union{Laser,NamedTuple,Nothing}=nothing)

Assembles the full Simulation struct from the requested components. Any systems not provided to the function are disabled by default.
The user can send either a completed Struct of the correct type or a NamedTuple with the correct key-word arguments to assemble the struct directly within the function
using the representative build_x function where x is the subsystem.

# Arguments
- 'densitymatrix': The DensityMatrix subsystem
- 'electronictemperature': The ElectronicTemperature subsystem
- 'phononictemperature': The PhononicTemperature subsystem
- 'athermalelectrons': The AthermalElectrons subsystem
- 'electronicdistribution': The ElectronicDistribution subsystem
- 'phononicdistribution': The PhononicDistribution subsystem
- 'structure': The Structure subsystem
- 'laser': The Laser subsystem

# Returns
- The PhononicTemperature struct with the users settings and parameters with any neccessary unit conversion.

build_Structure(; las::Laser=build_Laser(), Spatial_DOS::Bool = false, Elemental_System::Int = 1, dimension::Dimension = build_Dimension(),
                bulk_DOS::Union{String,Vector{String},Nothing} = nothing, DOS_folder::Union{String,Vector{String},Nothing} = nothing, 
                bulk_geometry::Union{String,Vector{String},Nothing} = nothing, slab_geometry::Union{String,Vector{String},Nothing} = nothing, 
                atomic_layer_tolerance::Union{Float64,Vector{Float64}} = 0.1, DOS::Union{spl,Vector{spl},Nothing} = nothing, 
                egrid::Union{Vector{Float64},Nothing} = nothing)

Outer constructor function to assemble the Structure struct. No unit conversion is performed.
All DOS files must be in the format |energy (eV), states (eV⁻¹Vᵤ⁻¹)|. Comment lines (#) are ignored
Defaults allow any unneccessary parameters for users simulation to be ignored.

# Arguments
- 'las': Laser struct, provide if not providing a pre-made energy grid
- 'Spatial_DOS': Bool for determening whether the DOS is spatially resolved or bulk
- 'Elemental_System': Float64 of different crystal systems in the structure
- 'dimension': Dimension struct, provide if not wanting a 0D calculation
- 'bulk_DOS': File location of the bulk DOS file
- 'DOS_folder': Location of a folder containing atom projected DOS. These must be in units of (eV⁻¹atom⁻¹) and be .dat files
- 'bulk_geometry': File location of the bulk DOS' geometry.in file
- 'slab_geometry': File location of the geometry.in to create the atom projected DOS' found in DOS_folder
- 'atomic_layer_tolerance': Minimum height spacing in Å for two atoms to be apart to be considered different layers
- 'DOS': A user made DOS or vector of DOS' for the simulation if not wanting to use the DOS_initialization function. Will overwrite this function.
         Must be of the type DataInterpolations.LinearInterpolation or a vector of them. It is advised to have extrapolation enabled just in case. 
- 'egrid': A user made energy grid if not wanting to use build_egrid function. Will overwrite this function. Must be an evenly spaced grid that
           has length(egrid) % 4 == 1 for the numerical integration algorithm to work. 
- 'BandStructure': Contains splines of [k->E, E->k] for the evaluation of magnetotransport and group velocity

# Returns
- The Structure struct with the DOS and egrid assembled or provided by the user

build_egrid(hv)

Builds an energy grid with spacing of 0.01 eV, limits of -2*hv to 2*hv and ensures that length(egrid) % 4 == 1

# Arguments
- 'hv': Photon energy of the laser

# Returns
- Evenly spaced energy grid that is suitable for the numerical integration algorithm and of sufficient accuracy/discretisation for accurate dynamics

FE_initalization(bulk_DOS::Union{String, Vector{String}})

Initializes Fermi energy based on bulk density of states. This is defined as the difference between the bottom of the
valence band and 0.0. This is required for calculations such as the ballistic velocity or athermal electron lifetime due
to electronic interactions.

# Arguments
- `bulk_DOS`: (Array of) file path to bulk density-of-states

# Returns
- Fermi energy/energies for each DOS input.

athem_electronelectroninteraction(sim::Simulation)

Returns an expression representing the electron-electron scattering term in the AthEM model.

# Arguments
- 'sim': settings of the desired simulation

# Returns
- Symbolic expression for relaxation or 0.0 if disabled.

athem_electronelectronscattering(fdis::VectorTel::Float64,μ::Float64,sim::Simulation,fneq::Vector{Float64},DOS::spl,n::Float64,τee::Union{Float64,Vector{Float64}})

Calculates the electron-electron scattering contribution using a modified relaxation time approximation.

# Arguments
- 'Tel': Electronic temperature of the system
- 'μ': Chemical potential at the current temperature
- 'sim': Simulation settings and parameters
- 'fneq': Current non-equilibrium electron distribution
- 'DOS': Spline of the density-of-states
- 'n': The number of electrons in the thermal system
- 'τee': The athermal electron lifetime

# Returns
- Change in the non-equilibrium distribution due to scattering with a thermal electronic system

athem_electrongeneration!(tmp::Vector, egrid::Vector{Float64},DOS::spl,ftotspl::spl,hv::Float64,M::Union{Float64,Vector{Float64}})

Computes the shape of the distribution change due to electron generation.

# Arguments
- `tmp`: Temporary vector for in-place update
- `egrid`: Energy grid.
- `DOS`: Density of states spline.
- `ftotspl::spl`: Spline of the total distribution.
- `hv`: Photon energy.
- `M`: Matrix elements.

# Returns
- Change in distribution due to electron generation

athem_electronelectroninteraction(sim::Simulation)

Returns an expression representing the athermal electron-phonon scattering term in the AthEM model.

# Arguments
- 'sim': settings of the desired simulation

# Returns
- Symbolic expression for relaxation or 0.0 if disabled.

athem_holegeneration!(tmp::Vector, egrid::Vector{Float64},DOS::spl,ftotspl::spl,hv::Float64,M::Union{Float64,Vector{Float64}})

Computes the shape of the distribution change due to hole generation.

# Arguments
- `tmp`: Temporary vector for in-place update
- `egrid`: Energy grid.
- `DOS`: Density of states spline.
- `ftotspl::spl`: Spline of the total distribution.
- `hv`: Photon energy.
- `M`: Matrix elements.

# Returns
- Change in distribution due to hole generation

athem_thermalelectronparticlechange(sim::Simulation)

Returns an expression for the total thermal electronum number due to relaxation and optional conductivity.

# Arguments
- 'sim': settings of the desired simulation

# Returns
- Expr for the time dependence of the thermal electron number.

athemdistribution_factory(sim, laser_expr)

Constructs the time evolution of the non-equilibrium distribution in the AthEM model
Find out more at https://arxiv.org/abs/2503.09479

# Arguments
- `sim`::Simulation
       Simulation struct containing physical and model parameters.
- `laser_expr`::Expr
              Expression representing the laser excitation term.

# Returns
- The total expression for the evolution of the athermal electron distribution, combining excitation and scattering terms.

athemexcitation!(Δfneqe::DiffCache, Δfneqh::DiffCache, ftot::Vector{Float64}, egrid::Vector{Float64}, DOS::spl, hv::Float64, M::Union{Float64,Vector{Float64}})

Computes the net non-equilibrium excitation using Fermi's Golden Rule.

# Arguments
- 'Δfneqe': A pre-initialized vector to reduce allocations. Can be DiffCache or Vector. 
            This vector contains the final result
- 'Δfneqh': A pre-initialized vector to reduce allocations. Can be DiffCache or Vector.
- `ftot`: Total distribution (f_eq + f_neq).
- `egrid`: Energy grid.
- `DOS`: Density of states spline.
- `hv`: Photon energy.
- `M`: Matrix elements.

# Returns
- In-place normalized change in particle distribution shape. Multiply by inputted laser energy at time t to get 
  full excitation shape change

athemexcitation_matrixelements(sim::Simulation)

Returns the symbolic matrix element expression for photo-excitation.
Currently implemented:
- :unity : No specified matrix elements - Returns 1.0

# Arguments
- 'sim': settings of the desired simulation

# Returns
- Matrix element expression.

build_athemdistribution(sim, athemexcite, Elecelec, Elecphon)

Combines excitation and scattering contributions into a single broadcasted sum expression in the AthEM model

# Arguments
- `sim`::Simulation
       The simulation object.
- `athemexcite`::Expr
               Expression for athermal excitation.
- `Elecelec`::Union{Expr, Float64}
            Electron-electron scattering Expr or 0.0 if disabled
- `Elecphon`::Union{Expr, Float64}
            Electron-phonon scattering Expr or 0.0 if disabled

# Returns
- A broadcasted summation of all interaction terms.

electron_distribution_transport!(v_g::Vector{Float64},f::AbstractArray{Float64},Δf::AbstractArray{Float64},dz::Float64)

Computes ballistic transport of an electronic distribution using second-order finite differences via a kinetic like model.
Uses forward(reverse) difference for the boundaries. 

# Arguments
- `v_g`: Group velocity vector or vector of vectors(spatially-resolved DOS)
- `f`: Distribution function.
- `Δf`: Output array to store result.
- `dz`: Spatial resolution.

# Returns
- In-place change to Δf

electron_relaxationtime(sim::Simulation)
Returns the requested expression for the electronic relaxation time. 
Currently implemented: 
- :constant : A constant time defined by sim.athermalelectrons.τ
- :FLT : Fermi Liquid lifetime defined in https://arxiv.org/abs/2503.09479

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Expression for the athermal electron lifetime due to electron-electron interactions.

find_relaxeddistribution(egrid::Vector{Float64},goal::Float64,n::Float64,DOS::spl,kB::Float64)

Solves for a Fermi distribution with the same internal energy as a given target.

# Arguments
- 'egrid': Energy grid distributions are solved on
- 'goal': The internal energy of the total electronic system
- 'n': Float64 of electrons in the thermal system
- 'DOS': Spline of the density-of-states

# Returns
- Fermi-Dirac distribution with same internal energy as the goal.

find_temperatureandμ(Tel::Float64,n::Float64,DOS::spl,egrid::Vector{Float64})

Given a temperature guess, computes chemical potential and internal energy.

# Arguments
- 'Tel': Guessed electronic temperature to match the goal
- 'n': Float64 of electrons in the thermal system
- 'egrid': Energy grid distributions are solved on
- 'DOS': Spline of the density-of-states

# Returns
- Internal energy of the current temperature guess.

phonon_relaxationtime(sim::Simulation)

Determines the requested expression for the electron-phonon relaxation time. 
Currently implemented:
- :constant : a constant lifetime given by simulation.athermalelectrons.τep
- :quasi : quasiparticle relaxation time defined in E.Carpene, Phys. Rev. B., 2006, 74, 024301.

# Arguments
- 'sim': settings of the desired simulation

# Returns
- Expression for the athermal electron lifetime due to phonon interactions.

build_dm(N::Int, occ::Vector{Float64})

Builds a density matrix of NxN states with diagonal occupations occ

# Arguments
- 'N': Number of states
- 'occ': Initial occupations of the density matrix 

# Returns
- An assembled density matrix with diagonal given by occ

discretize_DOS(sim::Simulation)

Constructs an expression for the propagation of a Hamiltonian's reaction to light
in the dipole approximation. It takes the dipole matrix from sim.densitymatrix.DipoleMatrix

# Arguments
- 'sim': Settings of the simulation containing the dipole matrix and electric fields 

# Returns
- Expression for the propagation of a density matrix under illumination in the dipole approximation

discretize_DOS(dos_file::String, no_states::Int, egrid::Vector{Float64})

Discretizes a DOS and creates a vector of repeating energy states of length no_states.
The number of these repeats should mimic the shape of the original DOS. The discretization is 
only performed along the vector given by egrid. 

# Arguments
- 'dos_file': File path to DOS
- 'no_states': Approximate length of the final vector (there may be slight differences in the final result)
- 'egrid': Vector of energy values to discretize the DOS onto 

# Returns
- An approximation of a DOS to use as the diagonal of a Hamiltonian

thermal_bath_densitymatrix(H::Matrix{Float64}, β::Float64, ne::Int)

Builds a thermal density matrix assuming Fermi-Dirac statistics 

# Arguments
- 'H': Hamiltonian of the system, no external forces should be included here
- 'Β': kB*T in the same units as the Hamiltonian
- 'ne': The total number of electrons the density matrix represents

# Returns
- An assembled density matrix with a thermal occupation

discretize_DOS(DOS::spl, egrid::Vector{Float64})

Calculates the total number of states within the range of the egrid of the DOS
provided by dos_file in discretize_DOS.

# Arguments
- 'DOS': Spline of the DOS generated from dos_file
- 'egrid': Vector of energy values to discretize the DOS onto 

# Returns
- Total number of states in the energy window (egrid) of the DOS

vonNeumann(dρ::Matrix{Float64},ρ::Matrix{Float64},H::Matrix{Float64})

Equation of motion for a density matrix

# Arguments
- 'dρ': In-place derivative of the density matrix
- 'ρ': Density matrix of the system
- 'H': Hamiltonian matrix of the system

# Returns
- In-place function calculating dρ

FermiDirac(Tel::Float64, μ::Float64, E::Union{Vector{Float64},Float64}) 

Returns the thermal occupation of electrons at a given temperature, chemical potential and energy.
If a vector of energies is given then it will return the distribution across that range.

# Arguments
- 'Tel': Thermal electronic temperature
- 'μ': Chemical Potential
- 'E': Energy value or range

# Returns
- Value or vector of the Fermi-Dirac distribution

dFDdE(Tel::Float64, μ::Float64, E::Union{Vector{Float64},Float64})

Returns the change in the Fermi distribution with respect to energy at the given energy value or range.

# Arguments
- 'Tel': Thermal electronic temperature
- 'μ': Chemical Potential
- 'E': Energy value or range

# Returns
- Value or vector of the partial derivative of the Fermi distribution with respect to energy

dFDdT(Tel::Float64, μ::Float64, E::Union{Vector{Float64},Float64})

Returns the change in the Fermi distribution with respect to temperature at the given energy value or range.

# Arguments
- 'Tel': Thermal electronic temperature
- 'μ': Chemical Potential
- 'E': Energy value or range

# Returns
- Value or vector of the partial derivative of the Fermi distribution with respect to temperature

dFDdμ(Tel::Float64, μ::Float64, E::Union{Vector{Float64},Float64})

Returns the change in the Fermi distribution with respect to chemical potential at the given energy value or range.

# Arguments
- 'Tel': Thermal electronic temperature
- 'μ': Chemical Potential
- 'E': Energy value or range

# Returns
- Value or vector of the partial derivative of the Fermi distribution with respect to chemical potential

athem_electempenergychange(sim::Simulation)

Builds an expression for the change in internal energy of an AthEM electronic bath

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Expression for the change of internal energy of an AthEM electronic bath

build_athemelectron(Δu::Expr)

Builds the differential equation (expression) for an AthEM thermal bath. The difference to the TTM
is due to accounting for the changing particle number

# Arguments
- 'Δu': Expression for the change in internal energy of the electronic bath

# Returns
- Expression for the time evolution of a, AthEM electronic temperature

build_electronTTM(sim::Simulation, Source::Expr, ElecPhon::Expr, HeatCapacity::Expr)

Builds the differential equation (expression) for a Two-Temperature like system.

# Arguments
- 'sim': Simulation settings and parameters
- 'Source': Expression for the incoming energy source (typically a laser)
- 'ElecPhon': Expression for the electron-phonon coupling
- 'HeatCapacity': Expression for calculating the heat capacity

# Returns
- Expression for the time evolution of a two-temperature model electronic temperature

depthderivative!(vec::Vector{Float64}, dz::Float64, Diff::Vector{Float64})

Calculates a 2nd order finite difference method of a vector along a grid with spacing dz.
Uses central difference in the middle and forward(reverse) difference at the top(end) of the vector

# Arguments
- 'vec': Vector the finite difference is being performed over
- 'dz': z-grid spacing
- 'Diff': Vector to store the finite difference result

# Returns
- Updates the diff vector with the finite difference result

elec_energychange(egrid::Vector{Float64}, relax_dis::Vector{Float64}, DOS::spl)

Calculates the energy change of the thermal bath due to non-equilibrium-equilibrium
electron scattering.

# Arguments
- 'egrid: Energy grid the simulation is solved over
- 'relax_dis': The change in distribution due to the electron-electron scattering
- 'DOS': Density-of-states of the system

# Returns
- The change in the internal energy of the thermal system due to e-e scattering

electronphonon_coupling(sim::Simulation)

Determines the expression for the coupling between an electronic and phononic thermal bath.
Currently implemented:
- :variable : Calculates the electron-phonon coupling parameter from the density-of-states of the system
- :constant : Uses a constant value for the electron-phonon coupling parameter (g)

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Expression for the electron-phonon thermal bath energy transfer

electrontemperature_conductivity!(Tel::Vector{Float64}, κ::Union{Float64,Vector{Float64}}, dz::Float64, Tph::Vector{Float64}, cond::Vector{Float64})

Calculates the change in temperature due to diffusive transport of energy through the system. Uses the boundary conditions of
setting dTel / dz to 0.0.

# Arguments
- 'Tel': Temperature of the electronic bath
- 'κ': Electronic thermal conductivity at room temperature
- 'dz': z-grid spacing
- 'Tph': Temperature of the phonon bath
- 'cond': Vector to store the change in electronic temperature

# Returns
- Updates the cond vector with the change in electronic temperature at each grid point

electrontemperature_factory(sim::Simulation, laser::Expr)

Assembles expresssion for how the electronic temperature should be propgated through time.

# Arguments
- 'sim': Simulation settings and parameters
- 'laser': Expression for laser energy as a funciton of time

# Returns
- Expression for the time-propagation of the electronic temperature subsystem

electrontemperature_heatcapacity(sim::Simulation)

Determines the expression for the electronic temperature's heat capacity.
Currently implemented:
- :linear : Specific heat of electrons(γ) multiplied by the current temperature
- :nonlinear : Calculated from the density-of-states of the system

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Expression for the electronic temperature's heat capacity

nonlinear_electronheatcapacity(Tel::Float64, μ::Float64, DOS::spl, egrid::Vector{Float64})

Calculates non-linear electronic bath heat capacity. A more accurate method than the
linear form.

# Arguments
- 'Tel': Temperature of the electronic bath
- 'μ': Chemical potential of the electronic bath
- 'DOS': Density-of-states of the system
- 'egrid': Energy grid the simulation is solved over

# Returns
- The current heat capacity of the electronic thermal bath

variable_electronphononcoupling(λ::Float64, ω::Float64, DOS::spl, Tel::Float64, μ::Float64, Tph::Float64, egrid::Vector{Float64})

Calculates the non-linear electron phonon coupling parameter and subsequent energy flow from the density-of-states
of the system. More accurate than a constant value. 
The expression can be found in Z. Lin, L. V. Zhigilei and V. Celli, Phys. Rev. B, 2008, 77, 075133.

# Arguments
- 'λ': Electron-phonon mass enhancement parameter
- 'ω': Second moment of the phonon spectrum
- 'DOS': Density-of-states of the system
- 'Tel': Temperature of the electronic bath
- 'μ': Chemical potential of the electronic bath
- 'Tph': Temperature of the phonon bath
- 'egrid': Energy grid the simulation is solved over

# Returns
- Energy flow between an electronic and phononic bath with a calculate g parameter

E_magnitude(las_field::Vector{Expr}, ext_field::Vector{Expr})

# Arguments
- 'las_field': Electric field vector from the laser
- 'ext_field': Electric field vector from external sources

# Returns
- Expression for the total magntitude of electric fields

get_laser_fields(las::Laser)

Returns an expression for the electric field induced by the laser
Currently implemented:
- :Gaussian
- :Rectangular

# Arguments
- 'las': Laser settings 

# Returns
- Expression for the electric field of the laser

laser_factory(sim::Simulation)

Assembles the user desired laser expression

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Expression for the power of the laser as a function of time and depth

photon_energytofrequency(hv::Real)

# Arguments
- 'hv': Photon energy in eV

# Returns
- Photon frequency in 1 / fs

spatial_laser(sim::Simulation)

Assembles the expression for the laser penetration depth. Currently only works in the z-axis

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Expression for the spatial shape of the laser

WIP: spatial_xy_laser(sim::Simulation)

Defines how the laser energy changes with distance from the spot centre. Do not use as simulations
aren't constructed for more than 1D.

Currently implemented:
- :Cylindrical : Circular shaped sample
- :Cubic : Cubic shaped sample

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Expression for the spatial shape of the laser in the xy plane

spatial_z_laser(sim::Simulation)

Defines how the laser energy changes with depth into the sample
Currently implemented:
- :ballistic : Ballistic depth of electrons defines spatial change, controlled by δb
- :optical : Absorbtion depth of the laser controls the spatial change, controlled by ϵ (1/α)
- :combined : Both ballistic and optical parameters enabled, sums the two effects together

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Expression for the spatial shape of the laser

temporal_laser(las::Laser)

Assembles the expression for the user desired temporal shape of the laser
Currently implemented:
- :Gaussian : Gaussian laser shape
- :Lorentzian : Lorentzian laser shape
- :HyperbolicSecant : Hyperbolic secant laser shape
- :Rectangular : Constant illumination style method. The laser is on for 2*FWHM at either side of 0.0 fs

# Arguments
- 'las': Settings and parameters of the laser

# Returns
- Expression for the temporal shape of the laser

build_phonontemperature(sim::Simulation, Source::Expr, ElecPhon::Expr, HeatCapacity::Expr)

Builds the differential equation (expression) for the phonon bath

# Arguments
- 'sim': Simulation settings and parameters
- 'Source': Expression for energy input from athermal electrons
- 'ElecPhon': Expression for the thermal electron-phonon coupling
- 'HeatCapacity': Expression for calculating the heat capacity

# Returns
- Expression for the time evolution of a phnonic thermal bath

neqelectron_phonontransfer(fneq::Vector{Float64}, egrid::Vector{Float64}, τep::Float64, DOS::spl)

Calculates energy input into the phonon bath due to non-equilibrium electron-phonon scattering

# Arguments
- 'fneq': Non-equilibrium electron distribution
- 'egrid': Energy grid the distributions are solved on
- 'τep': Lifetime of athermal electrons due to electron-phonon scattering
- 'DOS': Density-of-states of the system

# Returns
- Value of the change in the phonon internal energy

phonontemperature_conductivity!(Tph::Vector{Float64}, κ::Union{Float64,Vector{Float64}}, dz::Float64, cond::Vector{Float64})

Calculates thermal phonon conductivity due to diffusive transport

# Arguments
- 'Tph': Phonon bath temperature
- 'κ': Phonon thermal conductivity (assumed constant)
- 'dz': Spacing of the spatial grid
- 'cond': Vector to store the temperature change

# Returns
- Updates cond with the change in temperature at each z-grid point

phonontemperature_factory(sim::Simulation)

Assembles expresssion for how the phononic temperature should be propgated through time.

# Arguments
- 'sim': Simulation settings and parameters
# Returns
- Expression for the time-propagation of the phononic temperature subsystem

phonontemperature_heatcapacity(sim::Simulation)

Determines the expression for the phononic temperature's heat capacity.
Currently implemented:
- :constant : Constant phonon heat capacity (Cph)
- :variable : Calculated from Simpson's rule

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Expression for the phononic temperature's heat capacity

phonontemperature_source(sim::Simulation)

Determines the expresision of any additional source terms into the phonon bath. Currently this
is just athermal electron-phonon coupling
Currently implemented:

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Expression for any source term into the phonon bath

variable_phononheatcapacity(Tph::Float64, n::Float64, θ::Float64)

Calculates non-linear phononic bath heat capacity. A more accurate method than the
constant form.

# Arguments
- 'Tph': Temperature of the phononic bath
- 'n': Float64 of atoms per nm³
- 'θ': Debye temperature of the system

# Returns
- The current heat capacity of the phononic thermal bath

Bode_rule(y::Vector{Float64}, x::Vector{Float64})

Performs numerical integration on a grid using the higher order Boole's method.
Will integrate from end to end of the x vector

# Arguments
- 'y': The spectrum on a grid to be integrated
- 'x': The grid the spectrum is on w

# Returns
- The integration value across the range

c_T(μ::Float64, Tel::Float64, DOS::spl, egrid::Vector{Float64})

Calculates the change in internal energy of a Fermi Dirac distribution with respect to temperature

# Arguments
- 'μ': Electronic distribution
- 'Tel': The electronic bath temperature
- 'DOS': Density-of-states of the system
- 'egrid': Energy grid the simulation is solved over

# Returns
- The value of dU/dT

c_μ(μ::Float64, Tel::Float64, DOS::spl, egrid::Vector{Float64})

Calculates the change in internal energy of a Fermi Dirac distribution with respect to chemical potential

# Arguments
- 'μ': Electronic distribution
- 'Tel': The electronic bath temperature
- 'DOS': Density-of-states of the system
- 'egrid': Energy grid the simulation is solved over

# Returns
- The value of dU/dμ

find_chemicalpotential(no_part::Float64, Tel::Float64, DOS::spl, egrid::Vector{Float64})

Determines the chemical potential at the current temperature

# Arguments
- 'no_part': Float64 of particles in the thermal electronic system
- 'Tel': Temperature of the electronic bath
- 'DOS': Density-of-states of the system
- 'egrid': Energy grid the simulation is solved over

# Returns
- The value of the chemical potential

get_internalenergy(Dis::Vector{Float64}, DOS::spl, egrid::Vector{Float64})

Calculates the internal energy of the given electronic distribution

# Arguments
- 'Dis': Electronic distribution
- 'DOS': Density-of-states of the system
- 'egrid': Energy grid the simulation is solved over

# Returns
- The internal energy of the given distribution

get_noparticles(Dis::Vector{Float64}, DOS::spl, egrid::Vector{Float64})

Calculates number of electrons in the given distribution

# Arguments
- 'Dis': Electronic distribution
- 'DOS': Density-of-states of the system
- 'egrid': Energy grid the simulation is solved over

# Returns
- The number of electrons

get_thermalparticles(μ::Float64, Tel::Float64, DOS::spl, egrid::Vector{Float64})

Calculates number of electrons assuming a Fermi Dirac distribution

# Arguments
- 'μ': Chemical potential
- 'Tel': Temperature of the electronic bath
- 'DOS': Density-of-states of the system
- 'egrid': Energy grid the simulation is solved over

# Returns
- The number of electrons

p_T(μ::Float64, Tel::Float64, DOS::spl, egrid::Vector{Float64})

Calculates the change in the number of particles of a Fermi Dirac distribution with respect to temperature

# Arguments
- 'μ': Electronic distribution
- 'Tel': The electronic bath temperature
- 'DOS': Density-of-states of the system
- 'egrid': Energy grid the simulation is solved over

# Returns
- The value of dn/dT

p_μ(μ::Float64, Tel::Float64, DOS::spl, egrid::Vector{Float64})

Calculates the change in the number of particles of a Fermi Dirac distribution with respect to chemical potential

# Arguments
- 'μ': Electronic distribution
- 'Tel': The electronic bath temperature
- 'DOS': Density-of-states of the system
- 'egrid': Energy grid the simulation is solved over

# Returns
- The value of dn/dμ

run_simulation(sys::Dict{String,Union{Expr,Vector{Expr}}}, initialtemps::Dict{String, Float64},
tspan::Tuple{Float64,Float64}, sim::Simulation; 
save, tolerance, max_step, min_step, callbacks)

Generates the problem the dynamics will solve and then solves the coupled system of ODE's.
Currently always uses Tsit5 for the integration routine but in the future that may be user-defined

# Arguments
- 'sys': Dictionary of ODE equations to be propagated
- 'initialtemps': Dictionary of initial temperatures of the bath
- 'tspan': Tuple of values for the dynamics to run between (the laser is centred on 0.0)
- 'sim': Simulation settings and parameters

# KWARGS
- Any key-word arguemnts from DiffEq that work with ODEProblems can be including in a namedtuple here

# Returns
- The solution of the dynamics calculation

convert_symbols_to_strings(dict::Dict{Any,Any})

Converts any Symbol key values to their respecitve string for saving

# Arguments
- 'dict': The dictionary which needs keys converting

# Returns
- Dictionary with all Symbol keys changed to strings

create_datafile_and_structure(file_name::String)

Creates a file with the standard data structure:
- Athermal Electrons
- Density Matrix
- Electronic Temperature
- Phononic Temperature
- Electronic Distribution
- Phononic Distribution
- Laser
- Structure

# Arguments
- 'file_name': The name the user wishes the file to be saved as 

# Returns
- A created HDF5 file

dict_to_hdf5(f,d::Dict{String,Any})

Unpacks a dictionary and saves it to the designated file location (f)

# Arguments
- 'f': Location in the HDF5 file the dictionary information will be saved to
- 'd': The dictionary to be unpacked and saved

# Returns
- Nothing returned but dictionary saved to location

extract_structure(f, structure::Structure)

Specific way of writing structure settings to the file as both DOS and Dimension need special handling

# Arguments
- 'f': File to be written to
- 'structure': Structure settings and parameters

# Returns
- Nothing returned but structure parameters and settings wrote to file

generate_valuedict(sol, sim::Simulation, fields::Vector{Symbol})

Generates a dictionary of all propagated systems with placeholder arrays

# Arguments
- 'sol': The solution of the simulation
- 'sim': Simulation settings and parameters
- 'fields': Symbolic name of the propgated subsystems

# Returns
- Dictionary of values of each of the propgated subsystems

particlenumber_values(sim::Simulation)

Calculates the number of thermal electrons for a system when not propagated

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Array of the number of thermal electrons at each spatial point

populate_unpropagatedvalues!(sol, initial_temps::Dict{String,Float64}, fields::Vector{Symbol}, sim::Simulation, vals::Dict{String,AbstractArray{Float64}})

Adds placeholder information to any unpropagated fields

# Arguments
- 'sol': The solution of the simulation
- 'initial_temps': Initial temperature the baths have been set to
- 'sim': Simulation settings and parameters
- 'fields': Symbolic name of the propgated subsystems
- 'vals': Dictionary of the propgated subsytems with placeholder arrays

# Returns
- vals dictionary with the added unpropagated subsystems

populate_value_dict!(sol ,fields::Vector{Symbol}, vals::Dict{String,AbstractArray{Float64}})

Populates the subsystem dictionary with the reuslting vlaues from the simulation

# Arguments
- 'sol': The solution of the simulation
- 'fields': Symbolic name of the propgated subsystems
- 'vals': Dictionary of the propgated subsytems with placeholder arrays

# Returns
- The vals dictionary with the actual results inside

post_production(sol,file_name::String,initial_temps::Dict{String,Float64},output::Symbol,sim::Simulation)

Handles processing, saving the simulation after it has completed. Uses HDF5 file format
Currently the only output setting supported is :minimum which only saves the parameters, chemical potential and
progated systems

# Arguments
- 'sol': The solution of the simulation
- 'file_name': The name the user wishes the file to be saved as 
- 'initial_temps': Initial temperature the baths have been set to
- 'output': The output method required - may change in future
- 'sim': Simulation settings and parameters

# Returns
- Nothing is returned but a file is created

seperate_results(sol, initial_temps::Dict{String,Float64}, sim::Simulation)

Seperates the results held inside of the solution object
Also fills all unpropagated subsystems with parameter/temp information where neccessary

# Arguments
- 'sol': The solution of the simulation
- 'initial_temps': Initial temperature the baths have been set to
- 'sim': Simulation settings and parameters

# Returns
- Dictionary of values of each of the seperated systems

write_DOS(structure::Structure)

Handles converting the DOS splines into a writable format

# Arguments
- 'structure': Structure settings and parameters

# Returns
- Dictionary of the DOS

write_simulation(f,sim::Simulation)

Writes all settings and parameters to their respective location in the file (f)

# Arguments
- 'f': File to be written to
- 'sim': Simulation settings and parameters

# Returns
- Nothing returned but all parameters and settings wrote to file

DOSScale!(Temp::Matrix{Float64},bulk::Vector{Float64},Energies::Vector{Float64})

Ensures that all DOS are scaled to the same number of particles as the bulk

# Arguments
- 'Temp': Matrix of number of states in height X Energy
- 'bulk': The bulk DOS on the same energy grid as Temp
- 'Energies': The energy grid used for Temp and bulk

# Returns
- Temp with the rescaled DOS'

DOS_initialization(bulk_DOS::Union{String,Vector{String}}, bulk_geometry::String, DOS_folder::String, slab_geometry::String,
                   atomic_layer_tolerance::Float64, dimension::Dimension, zDOS::Bool, DOS::Union{Nothing, spl})

Determines the desired DOS configuration and assembles it accordingly

# Arguments
- 'bulk_DOS': The bulk DOS file location
- 'bulk_geometry': The bulk DOS geometry location
- 'DOS_folder': The folder where the atom-projected DOS' are present
- 'slab_geometry': The geometry.in file used in the creation of the atom-projected DOS'
- 'atomic_layer_tolerance': The minimum height 2 atoms need to be apart to be considered seperate layers (default = 0.1 nm)
- 'dim': The Dimension struct which holds the z-grid
- 'zDOS': Bool for enabling a spatially dependent DOS
' 'DOS': Allows the user to use their own splined DOS, for other regions of the code it must be the same type as normal DOS

# Returns
- A spline or vector of splines for the desired DOS structure

build_DOS(dos_file::String, geometry_file::String)

A convenient constructor for building a DOS from the DOS and geometry file.

# Arguments
- 'dos_file': Path to the total DOS file.
- 'geometry_file': Path to the geometry.in file

# Returns
- An interpolation object representing the DOS.

build_group_velocity(v_g::Union{Vector{Float64},Nothing}, FE::Union{Float64,Vector{Float64}}, Conductivity::Bool, conductive_velocity::Symbol, structure::Structure)

Creates a vector or array of vectors (spatial DOS) for the group veolcity for ballistic electron transport. Users can also provide a constant value in the form
of v_g, they must also set conductive_veolcity to constant.
Currently Implemented:
- :fermigas : Assumes a free electron gas solution therefore is an analytical form of the group velocity
- :effectiveoneband : Uses the effective one band model to convert a DOS into a group velocity, for more details see Mueller & Rethfeld, Phys. Rev. B 87, 035139.
- :constant : Uses the v_g argument to set a constant group velocity for all energy ranges

# Arguments
- 'v_g': A constant group velocity value if :constant is requested
- 'FE': The Fermi energy, calculated from get_FermiEnergy
- 'Conductivity': Sets whether ballistic transport should be enabled
- 'conductive_velocity': The form the user wants the group velocity to take
- 'structure': Contains all structural information including DOS and number of elemental systems

# Returns
- The group velocity vector or array of vectors as requested by the user for ballistic electron transport

build_zDOSArray(egrid::Vector{Float64},folder::String,files::Vector{String},heights::Vector{Float64})

Builds a matrix of the DOS as a function of height and energy for the individual layers. 

# Arguments
- 'egrid': Energy grid the DOS is calculated on
- 'folder': The folder where the atom-projected DOS' are present
- 'files': Vector of file names 
- 'heights': Vector of each file names height

# Returns
- A matrix of states as a function of height and energy

effective_one_band_velocity(DOS::spl, egrid::Vector{Float64}, FE::Float64)

Calculates the group velocity from the effective one band model.
For more details see Mueller & Rethfeld, Phys. Rev. B 87, 035139.

# Arguments
- 'DOS': The density-of-states of the system
- 'egrid': Energy grid all distributions are solved on
- 'FE': The Fermi energy, calculated from get_FermiEnergy

# Returns
- The effective one band model group velocity as a vector or vector of vectors depending on the structure
of the system

effective_onebandmodel(DOS::spl, egrid::Vector{Float64}, FE::Float64)

Calculates the dispersion relation within the effective one band model.
For more details see Mueller & Rethfeld, Phys. Rev. B 87, 035139.

# Arguments
- 'DOS': The density-of-states of the system
- 'egrid': Energy grid all distributions are solved on
- 'FE': The Fermi energy, calculated from get_FermiEnergy

# Returns
- The effective one band model dispersion relation

generate_DOS(File::String, unit_scalar::Float64)

Generates a spline of a DOS from a file. Assumes the structure of the DOS is column 1 is Energy in eV
and column 2 is States in eV⁻¹V⁻¹ (volume of unit cell) 

# Arguments
- 'File': Path to the total DOS file.
- 'unit_scalar': Scalar to convert the units (1/V in nm⁻³).

# Returns
- An interpolation object representing the DOS.

get_FermiEnergy(File::String)

Calculates the Fermi energy defined as the difference between 0.0 and the
bottom of the valence band. Assumes the DOS provided has the Fermi energy
set to 0.0.

# Arguments
- 'File': Path to a file containing the density of states data.

# Returns
- The Fermi energy calculated at the bottom of the valence band.

get_atomiclayers(atoms::Matrix{Float64},tolerance::Float64)

Seperats the atoms into their layers and selects a single atom from each layer. To remove degeneracy for 
larger supercell structures.

# Arguments
- 'atoms': Matrix of the atom number and it's repseictve coordinates
- 'tolerance': The minimum height 2 atoms need to be apart to be considered seperate layers (default = 0.1Å)

# Returns
- A trimmed matrix of atoms now containing one atom per layer

get_fermigas_dos(egrid, FE)

Function for calculating a free electron gas.

# Arguments
- 'egrid': Energy grid all distributions are solved on
- 'FE': The Fermi energy, calculated from get_FermiEnergy

# Returns
- The free electron gas denisty-of-states

get_fermigas_velocity(egrid::Vector{Float64}, EF::Float64)

The analytical free electron gas group velocity, requested by conductive_velocity = :fermigas

# Arguments
- 'egrid': Energy grid all distributions are solved on
- 'FE': The Fermi energy, calculated from get_FermiEnergy

# Returns
- The free electron gas group velocity

get_files_heights_forDOS(folder::String,geometry::String,tolerance::Float64)

Extracts atom from geometry, removes all bar one from each layer defined by tolerance, then connects the atom
to the corresponding file in the folder of DOS'. Readjusts the heights to set the top layer to 0.0 and the
lower layers to increase from there. 

# Arguments
- 'folder': The folder where the atom-projected DOS' are present
- 'geometry': The geometry.in file used in the creation of the atom-projected DOS'
- 'tolerance': The minimum height 2 atoms need to be apart to be considered seperate layers (default = 0.1Å)

# Returns
- A vector of DOS files and their respective heights

get_interpolant(xvals::Vector{Float64},yvals::Vector{Float64})

Generates a linear spline of any two vectors with a constant extrapolation applied to the boundaries.

# Arguments
- 'xvals': x-axis of the desired spline
- 'yvals': y-axis of the desired spline

# Returns
- Spline of yvals vs xvals

get_slabgeometry(file_path::String)

Extracts the atoms and their coordinates from a FHI-aims slab geometry.in file. It ignores any 
atoms that have their relaxation constrained. 

# Arguments
- 'file_path': The geometry.in file used in the creation of the atom-projected DOS'

# Returns
- A matrix of the atom number and their coordinates

get_unitcellvolume(geometry_file::String)

Calculates the volume of the unit cell for DOS unit conversion

# Arguments
- 'geometry_file': Path to the geometry.in file

# Returns
- Volume of the unit cell in nm⁻³

spatial_DOS(folder::String,geometry::String,bulk::String,Vbulk::Float64,dim::Dimension,tolerance::Float64)

Creates a spline of a DOS at each z-grid point in the simulation. Reads a folder of atom projected DOS's and the 
respective geomwtry.in file to determine the height of each DOS and interpolates between them to create the final
DOS'

# Arguments
- 'folder': The folder where the atom-projected DOS' are present
- 'geometry': The geometry.in file used in the creation of the atom-projected DOS'
- 'bulk': The bulk DOS file location
- 'Vbulk': The volume of the bulk unit cell
- 'dim': The Dimension struct which holds the z-grid
- 'tolerance': The minimum height 2 atoms need to be apart to be considered seperate layers (default = 0.1 nm)

# Returns
- A vector of splines corresponding to the DOS at each z-height

convert_units(value::Union{Quantity, AbstractArray{<:Quantity}, Float64, AbstractArray{Float64}})

Converts any user-given parameters that they have attached units to, to the correct units for LightMatter.jl
All values without Unitful.jl units attached are assumed to be in the correct units given by LightMatter_units

# Arguments
- 'value': The value(array of values) of the parameter, either as a Unitful.jl Quantity or a Float64

# Returns
- The Quantity converted into LightMatter.jl's preferred units, or the Float64 left as is

BaseUnits = 

Global Tuple for conversion factors from base SI units to base LightMatter units

Constants = (ħ = 0.6582 eVfs, kB = 8.617e-5 eV/K, me = 5.686 eVm)

Global named tuple for accessing constant physical values during a Simulation

LightMatter_units

A list of the units used in LightMatter.jl: Please convert all units to this

build_backgroundTTM(sim::Simulation)

Builds the structs for the TTM to be used in the sub-surface layers of the embedding AthEM nmethod

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Simulation struct setup to perform the background TTM for the embedded AthEM method

build_loopbody(sys, sim::Simulation)

Builds the multi-threaded section of the ODE problem

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Quote block for the multithreaded section of the ODE problem

conductivity_expressions(sim::Simulation)

Creates vector of expressions for the different conductivities occuring during the simulation

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Vector of expression for the conductivity of each subsytem if they are enabled

function_builder(sim::Simulation)

Assembles the correct dictionary of equations for the subsystems that are propagated. Mainly used for
seperating embedded and un-embedded methods

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Dictionary of subsystem names and their respective equation / expression

generate_expressions(sim::Simulation, laser::Expr)

Calculates and groups each subsystems expression in turn into a dictionary

# Arguments
- 'sim': Simulation settings and parameters
- 'laser': Expression for the temporal evolution and spatial decay of the laser

# Returns
- Dictionary of subssytems and their respective expressions

generate_initialconditions(sim::Simulation, initialtemps::Dict{String, Float64})

Generates the initial conditions (u0) NamedArrayPartition for the ODE 

# Arguments
- 'sim': Simulation settings and parameters
- 'initialtemps': Dictionary containing initial temepratures for electronic and phononic baths

# Returns
- NamedArrayPartition containing the initial conditions of the simulation

generate_initialconditions(sim::Simulation, initialtemps::Dict{String, Float64})

Generates the parameters as a NamedTuple for the ODE 

# Arguments
- 'sim': Simulation settings and parameters
- 'initialtemps': Dictionary containing initial temepratures for electronic and phononic baths

# Returns
- NamedTuple containing the parameters of the simulation

monometallic_system(sys, sim::Simulation)

Creates the expression block for the entire ODE function including multithreading.
This function is sepcfically for a monometallic system, Elemental_System == 1

# Arguments
- 'sys': Dictionary of expressions for each subssytem propagated
- 'sim': Simulation settings and parameters

# Returns
- Quote block for the monometallic problem

simulation_construction(sys, sim::Simulation)

Creates the expression block for the entire ODE function including multithreading.
This function calls monometallic_system or antenna_reactor_system respectively.

# Arguments
- 'sys': Dictionary of expressions for each subssytem propagated
- 'sim': Simulation settings and parameters

# Returns
- Quote block for the entire ODE problem

variable_renaming(sim::Simulation)

Renames variables such as u.x or p.y to just x or y so that during expression evaluation they 
are named correctly. 

# Arguments
- 'sim': Simulation settings and parameters

# Returns
- Expression for the variabble renaming to enter the top of the multithreaded loop

antenna_reactor_system(sys, sim)

Constructs the main expression for an antenna reactor simulation system.

# Arguments
- `sys`::Dict{String, Union{Expr, Vector{Expr}}}
       A dictionary or configuration object containing the expressions for each system propagated
- `sim`::Simulation
       A simulation object containing all necessary input data for the simulation.

# Returns
- An expression block containing the time output, conductivity expressions, and a threaded simulation loop.

ar_variable_renaming(sim)

Generates variable renaming expressions for translation between the variibale names in DiffEq.jl and LightMatter.jl

# Arguments
- `sim`::Simulation
       The Simulation struct containing all information about the simulation

# Returns
- An expression block assigning simulation-specific variable names.

mat_picker(height, cutoffs)

Selects an index based on material interface height and given cutoffs.

# Arguments
- `height`::Float64
          The height value to compare.
- `cutoffs`::Union{Float64,Vector{Float64}}
           A vector of cutoff heights defining material regions.

# Returns
- The index of the region in which `height` lies.

sim_seperation(sim)

Splits a composite `Simulation` object into separate simulations for each elemental subsystem.

# Arguments
- `sim`::Simulation
       The composite simulation containing multiple subsystems.

# Returns
- A vector of `Simulation` objects, each corresponding to a single elemental subsystem.

split_grid(grid, cutoffs)

Splits a numerical grid into regions that connect to each material

# Arguments
- `grid`::Vector{Float64}
        A vector of Float64 numbers representing the full z-grid
- `cutoffs`::Union{Float64, Vector{Float64}}
           A single value or vector of values defining the interfaces between each material 

# Returns
- A vector of sub-vectors representing segments of the original z-grid.

split_struct(data, number)

Splits fields of a composite object with vector fields into a vector of scalar instances.

# Arguments
- `data`::SimulationTypes 
        A subsystem of the simulation e.g. AthermalElectrons or ElectronicTemperature
- `number`::Int 
          Number of different elemental systems.

# Returns
- A vector of the subsystem with scalar data extracted from the original vector fields.

split_structure(structure)

Splits a `Structure` object into multiple structures if `DOS` is a vector. Different to split_struct due to the possibility 
of spatially resolved DOS' though currently that isn't implemented

# Arguments
- `structure`::Structure
             The struct containing the information of the simulatkion structure

# Returns
- A vector of `Structure` objects, one per element if applicable.