diff --git a/Project.toml b/Project.toml
index ca7c4fbe..17551b22 100644
--- a/Project.toml
+++ b/Project.toml
@@ -20,6 +20,7 @@ Interpolations = "a98d9a8b-a2ab-59e6-89dd-64a1c18fca59"
 Jedis = "b89ccfe0-2c5f-46f6-b89b-da3e1c2e286f"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 Measurements = "eff96d63-e80a-5855-80a2-b1b0885c5ab7"
 Measures = "442fdcdd-2543-5da2-b0f3-8c86c306513e"
 Memoize = "c03570c3-d221-55d1-a50c-7939bbd78826"
@@ -52,6 +53,7 @@ IMASutils = "1.3"
 Interpolations = "0.14.7, 0.15"
 Jedis = "0.3"
 LaTeXStrings = "1"
+MacroTools = "0.5.13"
 Measurements = "2"
 Measures = "0.3"
 Memoize = "0.4"
diff --git a/src/IMAS.jl b/src/IMAS.jl
index 53ab1ff7..49f228ed 100644
--- a/src/IMAS.jl
+++ b/src/IMAS.jl
@@ -3,6 +3,7 @@ module IMAS
 using Printf
 using Compat:@compat
 import OrderedCollections
+import MacroTools
 const document = OrderedCollections.OrderedDict()
 
 #= ====== =#
@@ -46,6 +47,7 @@ include("physics.jl")
 #= =========== =#
 #= EXPRESSIONS =#
 #= =========== =#
+include(joinpath("expressions", "expr_info.jl"))
 include(joinpath("expressions", "onetime.jl"))
 include(joinpath("expressions", "dynamic.jl"))
 
@@ -57,6 +59,6 @@ include("plot.jl")
 #= ====== =#
 #= EXPORT =#
 #= ====== =#
-export @ddtime, @findall
+export @ddtime, @findall, @explain
 
 end # module
diff --git a/src/expressions/dynamic.jl b/src/expressions/dynamic.jl
index d9548aa2..be5b0b17 100644
--- a/src/expressions/dynamic.jl
+++ b/src/expressions/dynamic.jl
@@ -24,26 +24,26 @@ dyexp = dynamic_expressions
 
 #  core_profiles.profiles_1d  #
 
-dyexp["core_profiles.profiles_1d[:].electrons.density"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].electrons.density"] =
     (rho_tor_norm; electrons, _...) -> electrons.density_thermal .+ electrons.density_fast
 
-dyexp["core_profiles.profiles_1d[:].electrons.density_fast"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].electrons.density_fast"] =
     (rho_tor_norm; _...) -> zero(rho_tor_norm)
 
-dyexp["core_profiles.profiles_1d[:].electrons.pressure_thermal"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].electrons.pressure_thermal"] =
     (rho_tor_norm; electrons, _...) -> pressure_thermal(electrons)
 
-dyexp["core_profiles.profiles_1d[:].electrons.pressure_fast_parallel"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].electrons.pressure_fast_parallel"] =
     (rho_tor_norm; _...) -> zero(rho_tor_norm)
 
-dyexp["core_profiles.profiles_1d[:].electrons.pressure_fast_perpendicular"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].electrons.pressure_fast_perpendicular"] =
     (rho_tor_norm; _...) -> zero(rho_tor_norm)
 
-dyexp["core_profiles.profiles_1d[:].electrons.pressure"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].electrons.pressure"] =
     (rho_tor_norm; electrons, _...) -> electrons.pressure_thermal .+ electrons.pressure_fast_parallel .+ electrons.pressure_fast_perpendicular .* 2.0
 
 
-dyexp["core_profiles.profiles_1d[:].ion[:].z_ion"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].ion[:].z_ion"] =
     (; ion, _...) -> begin
         if length(ion.element) == 1
             return ion.element[1].z_n
@@ -52,10 +52,10 @@ dyexp["core_profiles.profiles_1d[:].ion[:].z_ion"] =
         end
     end
 
-dyexp["core_profiles.profiles_1d[:].t_i_average"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].t_i_average"] =
     (rho_tor_norm; profiles_1d, _...) -> t_i_average(profiles_1d)
 
-dyexp["core_profiles.profiles_1d[:].ion[:].density"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].ion[:].density"] =
     (rho_tor_norm; ion, _...) -> begin
         if ismissing(ion, :density_thermal) && ismissing(ion, :density_fast)
             return zero(rho_tor_norm)
@@ -68,70 +68,70 @@ dyexp["core_profiles.profiles_1d[:].ion[:].density"] =
         end
     end
 
-dyexp["core_profiles.profiles_1d[:].ion[:].density_fast"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].ion[:].density_fast"] =
     (rho_tor_norm; ion, _...) -> ion.density .- ion.density_thermal
 
-dyexp["core_profiles.profiles_1d[:].ion[:].density_thermal"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].ion[:].density_thermal"] =
     (rho_tor_norm; ion, _...) -> ion.density .- ion.density_fast
 
-dyexp["core_profiles.profiles_1d[:].ion[:].pressure_thermal"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].ion[:].pressure_thermal"] =
     (rho_tor_norm; ion, _...) -> pressure_thermal(ion)
 
-dyexp["core_profiles.profiles_1d[:].ion[:].pressure_fast_parallel"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].ion[:].pressure_fast_parallel"] =
     (rho_tor_norm; _...) -> zero(rho_tor_norm)
 
-dyexp["core_profiles.profiles_1d[:].ion[:].pressure_fast_perpendicular"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].ion[:].pressure_fast_perpendicular"] =
     (rho_tor_norm; _...) -> zero(rho_tor_norm)
 
-dyexp["core_profiles.profiles_1d[:].ion[:].pressure"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].ion[:].pressure"] =
     (rho_tor_norm; ion, _...) -> ion.pressure_thermal .+ ion.pressure_fast_parallel .+ ion.pressure_fast_perpendicular .* 2.0
 
 
-dyexp["core_profiles.profiles_1d[:].pressure_ion_total"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].pressure_ion_total"] =
     (rho_tor_norm; profiles_1d, _...) -> pressure_thermal(profiles_1d.ion)
 
-dyexp["core_profiles.profiles_1d[:].pressure_thermal"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].pressure_thermal"] =
     (rho_tor_norm; profiles_1d, _...) -> pressure_thermal(profiles_1d)
 
-dyexp["core_profiles.profiles_1d[:].pressure_parallel"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].pressure_parallel"] =
     (rho_tor_norm; profiles_1d, _...) ->
         profiles_1d.pressure_thermal ./ 3.0 .+ profiles_1d.electrons.pressure_fast_parallel .+ sum(ion.pressure_fast_parallel for ion in profiles_1d.ion)
 
-dyexp["core_profiles.profiles_1d[:].pressure_perpendicular"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].pressure_perpendicular"] =
     (rho_tor_norm; profiles_1d, _...) ->
         profiles_1d.pressure_thermal ./ 3.0 .+ profiles_1d.electrons.pressure_fast_perpendicular .+ sum(ion.pressure_fast_perpendicular for ion in profiles_1d.ion)
 
-dyexp["core_profiles.profiles_1d[:].pressure"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].pressure"] =
     (rho_tor_norm; profiles_1d, _...) -> profiles_1d.pressure_perpendicular .* 2.0 .+ profiles_1d.pressure_parallel
 
 
-dyexp["core_profiles.profiles_1d[:].conductivity_parallel"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].conductivity_parallel"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> neo_conductivity(dd.equilibrium.time_slice[Float64(profiles_1d.time)], profiles_1d)
 
-dyexp["core_profiles.profiles_1d[:].j_bootstrap"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].j_bootstrap"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> Sauter_neo2021_bootstrap(dd.equilibrium.time_slice[Float64(profiles_1d.time)], profiles_1d)
 
-dyexp["core_profiles.profiles_1d[:].j_ohmic"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].j_ohmic"] =
     (rho_tor_norm; profiles_1d, _...) -> profiles_1d.j_total .- profiles_1d.j_non_inductive
 
-dyexp["core_profiles.profiles_1d[:].j_non_inductive"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].j_non_inductive"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> total_sources(dd.core_sources, profiles_1d; exclude_indexes=[7, 13], fields=[:j_parallel]).j_parallel .+ profiles_1d.j_bootstrap
 
-dyexp["core_profiles.profiles_1d[:].j_total"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].j_total"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> profiles_1d.j_ohmic .+ profiles_1d.j_non_inductive
 
-dyexp["core_profiles.profiles_1d[:].j_tor"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].j_tor"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> begin
         eqt = dd.equilibrium.time_slice[Float64(profiles_1d.time)]
         Jpar_2_Jtor(rho_tor_norm, profiles_1d.j_total, true, eqt)
     end
 
-dyexp["core_profiles.profiles_1d[:].zeff"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].zeff"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> zeff(profiles_1d)
 
 #  core_profiles.global_quantities  #
 
-dyexp["core_profiles.global_quantities.current_non_inductive"] =
+@store_expr dyexp["core_profiles.global_quantities.current_non_inductive"] =
     (time; dd, core_profiles, _...) ->
         return [
             begin
@@ -142,7 +142,7 @@ dyexp["core_profiles.global_quantities.current_non_inductive"] =
             for time0 in time
         ]
 
-dyexp["core_profiles.global_quantities.current_bootstrap"] =
+@store_expr dyexp["core_profiles.global_quantities.current_bootstrap"] =
     (time; dd, core_profiles, _...) ->
         return [
             begin
@@ -153,7 +153,7 @@ dyexp["core_profiles.global_quantities.current_bootstrap"] =
             for time0 in time
         ]
 
-dyexp["core_profiles.global_quantities.current_ohmic"] =
+@store_expr dyexp["core_profiles.global_quantities.current_ohmic"] =
     (time; dd, core_profiles, _...) ->
         return [
             begin
@@ -164,16 +164,16 @@ dyexp["core_profiles.global_quantities.current_ohmic"] =
             for time0 in time
         ]
 
-dyexp["core_profiles.global_quantities.ip"] =
+@store_expr dyexp["core_profiles.global_quantities.ip"] =
     (time; core_profiles, _...) -> [Ip(core_profiles.profiles_1d[Float64(time0)]) for time0 in time]
 
-dyexp["core_profiles.global_quantities.beta_tor_norm"] =
+@store_expr dyexp["core_profiles.global_quantities.beta_tor_norm"] =
     (time; dd, _...) -> [beta_tor_norm(dd.equilibrium, dd.core_profiles.profiles_1d[Float64(time0)]) for time0 in time]
 
-dyexp["core_profiles.global_quantities.v_loop"] =
+@store_expr dyexp["core_profiles.global_quantities.v_loop"] =
     (time; dd, _...) -> [v_loop(core_profiles.profiles_1d[Float64(time0)]) for time0 in time]
 
-dyexp["core_profiles.profiles_1d[:].time"] =
+@store_expr dyexp["core_profiles.profiles_1d[:].time"] =
     (; core_profiles, profiles_1d_index, _...) -> begin
         return core_profiles.time[profiles_1d_index]
     end
@@ -181,7 +181,7 @@ dyexp["core_profiles.profiles_1d[:].time"] =
 #= ============ =#
 # core_transport #
 #= ============ =#
-dyexp["core_transport.model[:].profiles_1d[:].time"] =
+@store_expr dyexp["core_transport.model[:].profiles_1d[:].time"] =
     (; core_transport, profiles_1d_index, _...) -> begin
         return core_transport.time[profiles_1d_index]
     end
@@ -189,47 +189,47 @@ dyexp["core_transport.model[:].profiles_1d[:].time"] =
 #= ========= =#
 # equilibrium #
 #= ========= =#
-dyexp["equilibrium.time_slice[:].global_quantities.energy_mhd"] =
+@store_expr dyexp["equilibrium.time_slice[:].global_quantities.energy_mhd"] =
     (; time_slice, _...) -> 3 / 2 * trapz(time_slice.profiles_1d.volume, time_slice.profiles_1d.pressure)
 
-dyexp["equilibrium.time_slice[:].global_quantities.q_95"] =
+@store_expr dyexp["equilibrium.time_slice[:].global_quantities.q_95"] =
     (; time_slice, _...) -> Interpolations.linear_interpolation(norm01(time_slice.profiles_1d.psi), time_slice.profiles_1d.q)(0.95)
 
-dyexp["equilibrium.time_slice[:].global_quantities.q_axis"] =
+@store_expr dyexp["equilibrium.time_slice[:].global_quantities.q_axis"] =
     (; time_slice, _...) -> time_slice.profiles_1d.q[1]
 
-dyexp["equilibrium.time_slice[:].global_quantities.q_min"] =
+@store_expr dyexp["equilibrium.time_slice[:].global_quantities.q_min"] =
     (; time_slice, _...) -> minimum(time_slice.profiles_1d.q)
 
-dyexp["equilibrium.time_slice[:].global_quantities.psi_axis"] =
+@store_expr dyexp["equilibrium.time_slice[:].global_quantities.psi_axis"] =
     (; time_slice, _...) -> time_slice.profiles_1d.psi[1]
 
-dyexp["equilibrium.time_slice[:].global_quantities.psi_boundary"] =
+@store_expr dyexp["equilibrium.time_slice[:].global_quantities.psi_boundary"] =
     (; time_slice, _...) -> time_slice.profiles_1d.psi[end]
 
-dyexp["equilibrium.time_slice[:].global_quantities.magnetic_axis.r"] =
+@store_expr dyexp["equilibrium.time_slice[:].global_quantities.magnetic_axis.r"] =
     (; time_slice, _...) -> time_slice.profiles_1d.geometric_axis.r[1]
 
-dyexp["equilibrium.time_slice[:].global_quantities.magnetic_axis.z"] =
+@store_expr dyexp["equilibrium.time_slice[:].global_quantities.magnetic_axis.z"] =
     (; time_slice, _...) -> time_slice.profiles_1d.geometric_axis.z[1]
 
 
-dyexp["equilibrium.time_slice[:].global_quantities.vacuum_toroidal_field.b0"] =
+@store_expr dyexp["equilibrium.time_slice[:].global_quantities.vacuum_toroidal_field.b0"] =
     (; dd, time_slice, _...) -> get_time_array(dd.pulse_schedule.tf.b_field_tor_vacuum, :reference, time_slice.time, :linear)
 
-dyexp["equilibrium.time_slice[:].global_quantities.vacuum_toroidal_field.r0"] =
+@store_expr dyexp["equilibrium.time_slice[:].global_quantities.vacuum_toroidal_field.r0"] =
     (; dd, _...) -> dd.pulse_schedule.tf.r0
 
-dyexp["equilibrium.time_slice[:].global_quantities.magnetic_axis.b_field_tor"] =
+@store_expr dyexp["equilibrium.time_slice[:].global_quantities.magnetic_axis.b_field_tor"] =
     (; time_slice, _...) -> time_slice.global_quantities.vacuum_toroidal_field.b0 * time_slice.global_quantities.vacuum_toroidal_field.r0 / time_slice.boundary.geometric_axis.r
 
-dyexp["equilibrium.time_slice[:].profiles_1d.geometric_axis.r"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_1d.geometric_axis.r"] =
     (psi; time_slice, _...) -> (time_slice.profiles_1d.r_outboard .+ time_slice.profiles_1d.r_inboard) .* 0.5
 
-dyexp["equilibrium.time_slice[:].profiles_1d.geometric_axis.z"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_1d.geometric_axis.z"] =
     (psi; time_slice, _...) -> psi .* 0.0 .+ time_slice.global_quantities.magnetic_axis.z
 
-dyexp["equilibrium.time_slice[:].boundary.geometric_axis.r"] =
+@store_expr dyexp["equilibrium.time_slice[:].boundary.geometric_axis.r"] =
     (; time_slice, _...) -> begin
         if !ismissing(time_slice.profiles_1d.geometric_axis, :r)
             return time_slice.profiles_1d.geometric_axis.r[end]
@@ -239,7 +239,7 @@ dyexp["equilibrium.time_slice[:].boundary.geometric_axis.r"] =
         end
     end
 
-dyexp["equilibrium.time_slice[:].boundary.geometric_axis.z"] =
+@store_expr dyexp["equilibrium.time_slice[:].boundary.geometric_axis.z"] =
     (; time_slice, _...) -> begin
         if !ismissing(time_slice.profiles_1d.geometric_axis, :z)
             return time_slice.profiles_1d.geometric_axis.z[end]
@@ -250,87 +250,87 @@ dyexp["equilibrium.time_slice[:].boundary.geometric_axis.z"] =
     end
 
 
-dyexp["equilibrium.time_slice[:].boundary.minor_radius"] =
+@store_expr dyexp["equilibrium.time_slice[:].boundary.minor_radius"] =
     (; time_slice, _...) -> (time_slice.profiles_1d.r_outboard[end] - time_slice.profiles_1d.r_inboard[end]) * 0.5
 
-dyexp["equilibrium.time_slice[:].boundary.elongation"] =
+@store_expr dyexp["equilibrium.time_slice[:].boundary.elongation"] =
     (; time_slice, _...) -> (time_slice.boundary.elongation_lower + time_slice.boundary.elongation_upper) * 0.5
 
-dyexp["equilibrium.time_slice[:].boundary.elongation_lower"] =
+@store_expr dyexp["equilibrium.time_slice[:].boundary.elongation_lower"] =
     (; time_slice, _...) -> time_slice.profiles_1d.elongation[end] # <======= IMAS limitation
 
-dyexp["equilibrium.time_slice[:].boundary.elongation_upper"] =
+@store_expr dyexp["equilibrium.time_slice[:].boundary.elongation_upper"] =
     (; time_slice, _...) -> time_slice.profiles_1d.elongation[end] # <======= IMAS limitation
 
-dyexp["equilibrium.time_slice[:].boundary.triangularity"] =
+@store_expr dyexp["equilibrium.time_slice[:].boundary.triangularity"] =
     (; time_slice, _...) -> (time_slice.boundary.triangularity_lower + time_slice.boundary.triangularity_upper) * 0.5
 
-dyexp["equilibrium.time_slice[:].boundary.triangularity_lower"] =
+@store_expr dyexp["equilibrium.time_slice[:].boundary.triangularity_lower"] =
     (; time_slice, _...) -> time_slice.profiles_1d.triangularity_lower[end]
 
-dyexp["equilibrium.time_slice[:].boundary.triangularity_upper"] =
+@store_expr dyexp["equilibrium.time_slice[:].boundary.triangularity_upper"] =
     (; time_slice, _...) -> time_slice.profiles_1d.triangularity_upper[end]
 
-dyexp["equilibrium.time_slice[:].boundary.squareness_lower_inner"] =
+@store_expr dyexp["equilibrium.time_slice[:].boundary.squareness_lower_inner"] =
     (; time_slice, _...) -> time_slice.profiles_1d.squareness_lower_inner[end]
 
-dyexp["equilibrium.time_slice[:].boundary.squareness_upper_inner"] =
+@store_expr dyexp["equilibrium.time_slice[:].boundary.squareness_upper_inner"] =
     (; time_slice, _...) -> time_slice.profiles_1d.squareness_upper_inner[end]
 
-dyexp["equilibrium.time_slice[:].boundary.squareness_lower_outer"] =
+@store_expr dyexp["equilibrium.time_slice[:].boundary.squareness_lower_outer"] =
     (; time_slice, _...) -> time_slice.profiles_1d.squareness_lower_outer[end]
 
-dyexp["equilibrium.time_slice[:].boundary.squareness_upper_outer"] =
+@store_expr dyexp["equilibrium.time_slice[:].boundary.squareness_upper_outer"] =
     (; time_slice, _...) -> time_slice.profiles_1d.squareness_upper_outer[end]
 
-dyexp["equilibrium.time_slice[:].boundary.squareness"] =
+@store_expr dyexp["equilibrium.time_slice[:].boundary.squareness"] =
     (; time_slice, _...) -> begin
         mxh = MXH(time_slice.boundary.outline.r, time_slice.boundary.outline.z, 2)
         return -mxh.s[2]
     end
 
-dyexp["equilibrium.time_slice[:].profiles_1d.j_tor"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_1d.j_tor"] =
     (psi; profiles_1d, _...) -> begin
         j_parallel = profiles_1d.j_parallel
         Jpar_2_Jtor(profiles_1d.rho_tor_norm, j_parallel, true, time_slice)
     end
 
-dyexp["equilibrium.time_slice[:].profiles_1d.j_parallel"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_1d.j_parallel"] =
     (psi; time_slice, profiles_1d, _...) -> begin
         j_tor = profiles_1d.j_tor
         Jtor_2_Jpar(profiles_1d.rho_tor_norm, j_tor, true, time_slice)
     end
 
 
-dyexp["equilibrium.time_slice[:].profiles_1d.dpressure_dpsi"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_1d.dpressure_dpsi"] =
     (psi; time_slice, profiles_1d, _...) -> gradient(psi, profiles_1d.pressure)
 
-dyexp["equilibrium.time_slice[:].profiles_1d.dvolume_drho_tor"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_1d.dvolume_drho_tor"] =
     (psi; profiles_1d, _...) -> gradient(profiles_1d.rho_tor, profiles_1d.volume)
 
-dyexp["equilibrium.time_slice[:].profiles_1d.dvolume_dpsi"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_1d.dvolume_dpsi"] =
     (psi; profiles_1d, _...) -> gradient(psi, profiles_1d.volume)
 
-dyexp["equilibrium.time_slice[:].profiles_1d.darea_drho_tor"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_1d.darea_drho_tor"] =
     (psi; profiles_1d, _...) -> gradient(profiles_1d.rho_tor, profiles_1d.area)
 
-dyexp["equilibrium.time_slice[:].profiles_1d.darea_dpsi"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_1d.darea_dpsi"] =
     (psi; profiles_1d, _...) -> gradient(psi, profiles_1d.area)
 
-dyexp["equilibrium.time_slice[:].profiles_1d.dpsi_drho_tor"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_1d.dpsi_drho_tor"] =
     (psi; profiles_1d, _...) -> gradient(profiles_1d.rho_tor, profiles_1d.psi)
 
-dyexp["equilibrium.time_slice[:].profiles_1d.psi_norm"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_1d.psi_norm"] =
     (psi; _...) -> norm01(psi)
 
 # 2D
-dyexp["equilibrium.time_slice[:].profiles_2d[:].r"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_2d[:].r"] =
     (dim1, dim2; _...) -> ones(length(dim2))' .* dim1
 
-dyexp["equilibrium.time_slice[:].profiles_2d[:].z"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_2d[:].z"] =
     (dim1, dim2; _...) -> dim2' .* ones(length(dim1))
 
-dyexp["equilibrium.time_slice[:].profiles_2d[:].b_field_tor"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_2d[:].b_field_tor"] =
     (dim1, dim2; time_slice, profiles_2d, _...) -> begin
         psi = time_slice.profiles_1d.psi
         fpol = time_slice.profiles_1d.f
@@ -338,17 +338,17 @@ dyexp["equilibrium.time_slice[:].profiles_2d[:].b_field_tor"] =
         return FPOL ./ profiles_2d.r
     end
 
-dyexp["equilibrium.time_slice[:].profiles_2d[:].b_field_r"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_2d[:].b_field_r"] =
     (dim1, dim2; profiles_2d, _...) -> begin
         return Br_Bz(profiles_2d).Br
     end
 
-dyexp["equilibrium.time_slice[:].profiles_2d[:].b_field_z"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_2d[:].b_field_z"] =
     (dim1, dim2; profiles_2d, _...) -> begin
         return Br_Bz(profiles_2d).Bz
     end
 
-dyexp["equilibrium.time_slice[:].profiles_2d[:].j_tor"] =
+@store_expr dyexp["equilibrium.time_slice[:].profiles_2d[:].j_tor"] =
     (dim1, dim2; profiles_2d, _...) -> begin
         dBzdR = gradient(dim1, dim2, profiles_2d.b_field_z, 1)
         dBrdZ = gradient(dim1, dim2, profiles_2d.b_field_r, 2)
@@ -356,7 +356,7 @@ dyexp["equilibrium.time_slice[:].profiles_2d[:].j_tor"] =
     end
 
 
-dyexp["equilibrium.time_slice[:].time"] =
+@store_expr dyexp["equilibrium.time_slice[:].time"] =
     (; equilibrium, time_slice_index, _...) -> begin
         return equilibrium.time[time_slice_index]
     end
@@ -364,85 +364,85 @@ dyexp["equilibrium.time_slice[:].time"] =
 #= ============ =#
 #  core_sources  #
 #= ============ =#
-dyexp["core_sources.source[:].profiles_1d[:].electrons.power_inside"] =
+@store_expr dyexp["core_sources.source[:].profiles_1d[:].electrons.power_inside"] =
     (rho_tor_norm; profiles_1d, _...) -> begin
         energy = profiles_1d.electrons.energy
         cumtrapz(profiles_1d.grid.volume, energy)
     end
 
-dyexp["core_sources.source[:].profiles_1d[:].electrons.energy"] =
+@store_expr dyexp["core_sources.source[:].profiles_1d[:].electrons.energy"] =
     (rho_tor_norm; profiles_1d, _...) -> begin
         power_inside = profiles_1d.electrons.power_inside
         gradient(profiles_1d.grid.volume, power_inside)
     end
 
 
-dyexp["core_sources.source[:].profiles_1d[:].total_ion_power_inside"] =
+@store_expr dyexp["core_sources.source[:].profiles_1d[:].total_ion_power_inside"] =
     (rho_tor_norm; profiles_1d, _...) -> begin
         total_ion_energy = profiles_1d.total_ion_energy
         cumtrapz(profiles_1d.grid.volume, total_ion_energy)
     end
 
-dyexp["core_sources.source[:].profiles_1d[:].total_ion_energy"] =
+@store_expr dyexp["core_sources.source[:].profiles_1d[:].total_ion_energy"] =
     (rho_tor_norm; profiles_1d, _...) -> begin
         total_ion_power_inside = profiles_1d.total_ion_power_inside
         gradient(profiles_1d.grid.volume, total_ion_power_inside)
     end
 
 
-dyexp["core_sources.source[:].profiles_1d[:].electrons.particles_inside"] =
+@store_expr dyexp["core_sources.source[:].profiles_1d[:].electrons.particles_inside"] =
     (rho_tor_norm; profiles_1d, _...) -> begin
         particles = profiles_1d.electrons.particles
         cumtrapz(profiles_1d.grid.volume, particles)
     end
 
-dyexp["core_sources.source[:].profiles_1d[:].electrons.particles"] =
+@store_expr dyexp["core_sources.source[:].profiles_1d[:].electrons.particles"] =
     (rho_tor_norm; profiles_1d, _...) -> begin
         particles_inside = profiles_1d.electrons.particles_inside
         gradient(profiles_1d.grid.volume, particles_inside)
     end
 
 
-dyexp["core_sources.source[:].profiles_1d[:].current_parallel_inside"] =
+@store_expr dyexp["core_sources.source[:].profiles_1d[:].current_parallel_inside"] =
     (rho_tor_norm; profiles_1d, _...) -> begin
         j_parallel = profiles_1d.j_parallel
         cumtrapz(profiles_1d.grid.area, j_parallel)
     end
 
-dyexp["core_sources.source[:].profiles_1d[:].j_parallel"] =
+@store_expr dyexp["core_sources.source[:].profiles_1d[:].j_parallel"] =
     (rho_tor_norm; profiles_1d, _...) -> begin
         current_parallel_inside = profiles_1d.current_parallel_inside
         gradient(profiles_1d.grid.area, current_parallel_inside)
     end
 
 
-dyexp["core_sources.source[:].profiles_1d[:].torque_tor_inside"] =
+@store_expr dyexp["core_sources.source[:].profiles_1d[:].torque_tor_inside"] =
     (rho_tor_norm; profiles_1d, _...) -> begin
         momentum_tor = profiles_1d.momentum_tor
         cumtrapz(profiles_1d.grid.volume, momentum_tor)
     end
 
-dyexp["core_sources.source[:].profiles_1d[:].momentum_tor"] =
+@store_expr dyexp["core_sources.source[:].profiles_1d[:].momentum_tor"] =
     (rho_tor_norm; profiles_1d, _...) -> begin
         torque_tor_inside = profiles_1d.torque_tor_inside
         gradient(profiles_1d.grid.volume, torque_tor_inside)
     end
 
 
-dyexp["core_sources.source[:].profiles_1d[:].ion[:].particles_inside"] =
+@store_expr dyexp["core_sources.source[:].profiles_1d[:].ion[:].particles_inside"] =
     (rho_tor_norm; profiles_1d, ion, _...) -> begin
         particles = ion.particles
         cumtrapz(profiles_1d.grid.volume, particles)
     end
 
-dyexp["core_sources.source[:].profiles_1d[:].ion[:].particles"] =
+@store_expr dyexp["core_sources.source[:].profiles_1d[:].ion[:].particles"] =
     (rho_tor_norm; profiles_1d, ion, _...) -> begin
         particles_inside = ion.particles_inside
         gradient(profiles_1d.grid.volume, ion.particles_inside)
     end
 
 
-dyexp["core_sources.source[:].profiles_1d[:].time"] =
+@store_expr dyexp["core_sources.source[:].profiles_1d[:].time"] =
     (; core_sources, profiles_1d_index, _...) -> begin
         return core_sources.time[profiles_1d_index]
     end
@@ -450,145 +450,145 @@ dyexp["core_sources.source[:].profiles_1d[:].time"] =
 #= ===== =#
 #  build  #
 #= ===== =#
-dyexp["build.layer[:].identifier"] =
+@store_expr dyexp["build.layer[:].identifier"] =
     (; build, layer, _...) -> begin
         plasma_index = index(get_build_layer(build.layer; type=_plasma_))
         in_len = length(get_build_layers(build.layer; fs=_in_))
         return -(index(layer) - plasma_index) * layer.side
     end
 
-dyexp["build.layer[:].outline.r"] =
+@store_expr dyexp["build.layer[:].outline.r"] =
     (x; layer, _...) -> opposite_side_layer(layer).outline.r
 
-dyexp["build.layer[:].outline.z"] =
+@store_expr dyexp["build.layer[:].outline.z"] =
     (x; layer, _...) -> opposite_side_layer(layer).outline.z
 
-dyexp["build.layer[:].shape"] =
+@store_expr dyexp["build.layer[:].shape"] =
     (; layer, _...) -> opposite_side_layer(layer).shape
 
-dyexp["build.layer[:].shape_parameters"] =
+@store_expr dyexp["build.layer[:].shape_parameters"] =
     (; layer, _...) -> opposite_side_layer(layer).shape_parameters
 
-dyexp["build.layer[:].start_radius"] =
+@store_expr dyexp["build.layer[:].start_radius"] =
     (; build, layer_index, _...) -> build_radii(build.layer)[1:end-1][layer_index]
 
-dyexp["build.layer[:].end_radius"] =
+@store_expr dyexp["build.layer[:].end_radius"] =
     (; build, layer_index, _...) -> build_radii(build.layer)[2:end][layer_index]
 
-dyexp["build.layer[:].area"] =
+@store_expr dyexp["build.layer[:].area"] =
     (; layer, _...) -> area(layer)
 
-dyexp["build.layer[:].volume"] =
+@store_expr dyexp["build.layer[:].volume"] =
     (; layer, _...) -> volume(layer)
 
-dyexp["build.tf.ripple"] =
+@store_expr dyexp["build.tf.ripple"] =
     (; build, _...) -> tf_ripple(get_build_layer(build.layer; type=_plasma_).end_radius, get_build_layer(build.layer; type=_tf_, fs=_lfs_).start_radius, build.tf.coils_n)
 
-dyexp["build.tf.wedge_thickness"] =
+@store_expr dyexp["build.tf.wedge_thickness"] =
     (; build, _...) -> 2π * get_build_layer(build.layer; type=_tf_, fs=_hfs_).end_radius / build.tf.coils_n
 
 #= ======= =#
 #  costing  #
 #= ======= =#
-dyexp["costing.cost_direct_capital.system[:].cost"] =
+@store_expr dyexp["costing.cost_direct_capital.system[:].cost"] =
     (; system, _...) -> isempty(system.subsystem) ? error("no subsystem") : sum(sub.cost for sub in system.subsystem if !ismissing(sub, :cost))
 
-dyexp["costing.cost_direct_capital.cost"] =
+@store_expr dyexp["costing.cost_direct_capital.cost"] =
     (; cost_direct_capital, _...) -> isempty(cost_direct_capital.system) ? error("no system") : sum(sys.cost for sys in cost_direct_capital.system if !ismissing(sys, :cost))
 
-dyexp["costing.cost_operations.system[:].yearly_cost"] =
+@store_expr dyexp["costing.cost_operations.system[:].yearly_cost"] =
     (; system, _...) -> isempty(system.subsystem) ? error("no subsystem") : sum(sub.yearly_cost for sub in system.subsystem if !ismissing(sub, :yearly_cost))
 
-dyexp["costing.cost_operations.yearly_cost"] =
+@store_expr dyexp["costing.cost_operations.yearly_cost"] =
     (; cost_operations, _...) -> isempty(cost_operations.system) ? error("no system") : sum(sys.yearly_cost for sys in cost_operations.system if !ismissing(sys, :yearly_cost))
 
-dyexp["costing.cost_decommissioning.system[:].cost"] =
+@store_expr dyexp["costing.cost_decommissioning.system[:].cost"] =
     (; system, _...) -> isempty(system.subsystem) ? error("no subsystem") : sum(sub.cost for sub in system.subsystem if !ismissing(sub, :cost))
 
-dyexp["costing.cost_decommissioning.cost"] =
+@store_expr dyexp["costing.cost_decommissioning.cost"] =
     (; cost_decommissioning, _...) -> isempty(cost_decommissioning.system) ? error("no system") : sum(sys.cost for sys in cost_decommissioning.system if !ismissing(sys, :cost))
 
 #= ============== =#
 #  BalanceOfPlant  #
 #= ============== =#
-dyexp["balance_of_plant.Q_plant"] =
+@store_expr dyexp["balance_of_plant.Q_plant"] =
     (time; balance_of_plant, _...) -> balance_of_plant.power_plant.power_electric_generated ./ balance_of_plant.power_electric_plant_operation.total_power
 
-dyexp["balance_of_plant.power_electric_net"] =
+@store_expr dyexp["balance_of_plant.power_electric_net"] =
     (time; balance_of_plant, _...) -> balance_of_plant.power_plant.power_electric_generated .- balance_of_plant.power_electric_plant_operation.total_power
 
-dyexp["balance_of_plant.power_electric_plant_operation.total_power"] =
+@store_expr dyexp["balance_of_plant.power_electric_plant_operation.total_power"] =
     (time; power_electric_plant_operation, _...) -> sum(sys.power for sys in power_electric_plant_operation.system)
 
 #= ========= =#
 #  divertors  #
 #= ========= =#
-dyexp["divertors.divertor[:].power_black_body.time"] =
+@store_expr dyexp["divertors.divertor[:].power_black_body.time"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_black_body)[1]
 
-dyexp["divertors.divertor[:].power_black_body.data"] =
+@store_expr dyexp["divertors.divertor[:].power_black_body.data"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_black_body)[2]
 
 
-dyexp["divertors.divertor[:].power_conducted.time"] =
+@store_expr dyexp["divertors.divertor[:].power_conducted.time"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_conducted)[1]
 
-dyexp["divertors.divertor[:].power_conducted.data"] =
+@store_expr dyexp["divertors.divertor[:].power_conducted.data"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_conducted)[2]
 
 
-dyexp["divertors.divertor[:].power_convected.time"] =
+@store_expr dyexp["divertors.divertor[:].power_convected.time"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_convected)[1]
 
-dyexp["divertors.divertor[:].power_convected.data"] =
+@store_expr dyexp["divertors.divertor[:].power_convected.data"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_convected)[2]
 
 
-dyexp["divertors.divertor[:].power_currents.time"] =
+@store_expr dyexp["divertors.divertor[:].power_currents.time"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_currents)[1]
 
-dyexp["divertors.divertor[:].power_currents.data"] =
+@store_expr dyexp["divertors.divertor[:].power_currents.data"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_currents)[2]
 
 
-dyexp["divertors.divertor[:].power_incident.time"] =
+@store_expr dyexp["divertors.divertor[:].power_incident.time"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_incident)[1]
 
-dyexp["divertors.divertor[:].power_incident.data"] =
+@store_expr dyexp["divertors.divertor[:].power_incident.data"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_incident)[2]
 
 
-dyexp["divertors.divertor[:].power_neutrals.data"] =
+@store_expr dyexp["divertors.divertor[:].power_neutrals.data"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_neutrals)[1]
 
-dyexp["divertors.divertor[:].power_neutrals.data"] =
+@store_expr dyexp["divertors.divertor[:].power_neutrals.data"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_neutrals)[2]
 
 
-dyexp["divertors.divertor[:].power_radiated.time"] =
+@store_expr dyexp["divertors.divertor[:].power_radiated.time"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_radiated)[1]
 
-dyexp["divertors.divertor[:].power_radiated.data"] =
+@store_expr dyexp["divertors.divertor[:].power_radiated.data"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_radiated)[2]
 
 
-dyexp["divertors.divertor[:].power_recombination_neutrals.time"] =
+@store_expr dyexp["divertors.divertor[:].power_recombination_neutrals.time"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_recombination_neutrals)[1]
 
-dyexp["divertors.divertor[:].power_recombination_neutrals.data"] =
+@store_expr dyexp["divertors.divertor[:].power_recombination_neutrals.data"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_recombination_neutrals)[2]
 
 
-dyexp["divertors.divertor[:].power_recombination_plasma.time"] =
+@store_expr dyexp["divertors.divertor[:].power_recombination_plasma.time"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_recombination_plasma)[1]
 
-dyexp["divertors.divertor[:].power_recombination_plasma.data"] =
+@store_expr dyexp["divertors.divertor[:].power_recombination_plasma.data"] =
     (time; divertor, _...) -> divertor_totals_from_targets(divertor, :power_recombination_plasma)[2]
 
 #= ============== =#
 #  pulse_schedule  #
 #= ============== =#
-dyexp["pulse_schedule.time"] =
+@store_expr dyexp["pulse_schedule.time"] =
     (time; pulse_schedule, _...) -> begin
         all_times = Float64[]
         for item in keys(pulse_schedule)
@@ -602,25 +602,25 @@ dyexp["pulse_schedule.time"] =
         return sort!(unique(all_times))
     end
 
-dyexp["pulse_schedule.tf.b_field_tor_vacuum_r.reference"] =
+@store_expr dyexp["pulse_schedule.tf.b_field_tor_vacuum_r.reference"] =
     (time; tf, _...) -> tf.r0 .* tf.b_field_tor_vacuum.reference
 
-dyexp["pulse_schedule.tf.r0"] =
+@store_expr dyexp["pulse_schedule.tf.r0"] =
     (; dd, _...) -> dd.equilibrium.vacuum_toroidal_field.r0
 
-dyexp["pulse_schedule.tf.b_field_tor_vacuum.reference"] =
+@store_expr dyexp["pulse_schedule.tf.b_field_tor_vacuum.reference"] =
     (time; dd, _...) -> dd.equilibrium.vacuum_toroidal_field.b0
 
-dyexp["pulse_schedule.tf.time"] =
+@store_expr dyexp["pulse_schedule.tf.time"] =
     (time; dd, _...) -> dd.equilibrium.time
 
 #= ========= =#
 #  stability  #
 #= ========= =#
-dyexp["stability.model[:].cleared"] =
+@store_expr dyexp["stability.model[:].cleared"] =
     (time; model, _...) -> Int.(model.fraction .<= 1.0)
 
-dyexp["stability.all_cleared"] =
+@store_expr dyexp["stability.all_cleared"] =
     (time; stability, _...) -> begin
         all_cleared = ones(Int, length(time))
         for model in stability.model
@@ -632,7 +632,7 @@ dyexp["stability.all_cleared"] =
 #= ======= =#
 #  summary  #
 #= ======= =#
-dyexp["summary.fusion.power.value"] = # NOTE: This is the fusion power that is coupled to the plasma
+@store_expr dyexp["summary.fusion.power.value"] = # NOTE: This is the fusion power that is coupled to the plasma
     (time; dd, summary, _...) -> begin
         tmp = eltype(summary)[]
         for time in summary.time
@@ -641,10 +641,10 @@ dyexp["summary.fusion.power.value"] = # NOTE: This is the fusion power that is c
         return tmp
     end
 
-dyexp["summary.global_quantities.ip.value"] =
+@store_expr dyexp["summary.global_quantities.ip.value"] =
     (time; dd, summary, _...) -> [dd.equilibrium.time_slice[Float64(time0)].global_quantities.ip for time0 in time]
 
-dyexp["summary.global_quantities.current_bootstrap.value"] =
+@store_expr dyexp["summary.global_quantities.current_bootstrap.value"] =
     (time; dd, summary, _...) -> begin
         tmp = eltype(summary)[]
         for time in summary.time
@@ -655,7 +655,7 @@ dyexp["summary.global_quantities.current_bootstrap.value"] =
         return tmp
     end
 
-dyexp["summary.global_quantities.current_non_inductive.value"] =
+@store_expr dyexp["summary.global_quantities.current_non_inductive.value"] =
     (time; dd, summary, _...) -> begin
         tmp = eltype(summary)[]
         for time in summary.time
@@ -666,7 +666,7 @@ dyexp["summary.global_quantities.current_non_inductive.value"] =
         return tmp
     end
 
-dyexp["summary.global_quantities.current_ohm.value"] =
+@store_expr dyexp["summary.global_quantities.current_ohm.value"] =
     (time; dd, summary, _...) -> begin
         tmp = eltype(summary)[]
         for time in summary.time
@@ -678,50 +678,50 @@ dyexp["summary.global_quantities.current_ohm.value"] =
     end
 
 
-dyexp["summary.global_quantities.beta_pol_mhd.value"] =
+@store_expr dyexp["summary.global_quantities.beta_pol_mhd.value"] =
     (time; dd, summary, _...) -> [dd.equilibrium.time_slice[Float64(time0)].global_quantities.beta_pol for time0 in time]
 
-dyexp["summary.global_quantities.beta_tor.value"] =
+@store_expr dyexp["summary.global_quantities.beta_tor.value"] =
     (time; dd, summary, _...) -> [beta_tor(dd.equilibrium, dd.core_profiles.profiles_1d[Float64(time0)]) for time0 in time]
 
-dyexp["summary.global_quantities.beta_tor_mhd.value"] =
+@store_expr dyexp["summary.global_quantities.beta_tor_mhd.value"] =
     (time; dd, summary, _...) -> [dd.equilibrium.time_slice[Float64(time0)].global_quantities.beta_tor for time0 in time]
 
-dyexp["summary.global_quantities.beta_tor_norm_mhd.value"] =
+@store_expr dyexp["summary.global_quantities.beta_tor_norm_mhd.value"] =
     (time; dd, summary, _...) -> [dd.equilibrium.time_slice[Float64(time0)].global_quantities.beta_normal for time0 in time]
 
-dyexp["summary.global_quantities.beta_tor_norm.value"] =
+@store_expr dyexp["summary.global_quantities.beta_tor_norm.value"] =
     (time; dd, summary, _...) -> [beta_tor_norm(dd.equilibrium, dd.core_profiles.profiles_1d[Float64(time0)]) for time0 in time]
 
-dyexp["summary.global_quantities.beta_tor_thermal_norm.value"] =
+@store_expr dyexp["summary.global_quantities.beta_tor_thermal_norm.value"] =
     (time; dd, summary, _...) -> [beta_tor_thermal_norm(dd.equilibrium, dd.core_profiles.profiles_1d[Float64(time0)]) for time0 in time]
 
-dyexp["summary.global_quantities.energy_thermal.value"] =
+@store_expr dyexp["summary.global_quantities.energy_thermal.value"] =
     (time; dd, summary, _...) -> [energy_thermal(dd.core_profiles.profiles_1d[Float64(time0)]) for time0 in time]
 
-dyexp["summary.global_quantities.tau_energy.value"] =
+@store_expr dyexp["summary.global_quantities.tau_energy.value"] =
     (time; dd, summary, _...) -> [tau_e_thermal(dd.core_profiles.profiles_1d[Float64(time0)], dd.core_sources; subtract_radiation_losses=false) for time0 in time]
 
-dyexp["summary.global_quantities.tau_energy_98.value"] =
+@store_expr dyexp["summary.global_quantities.tau_energy_98.value"] =
     (time; dd, summary, _...) -> [tau_e_h98(dd; time0, subtract_radiation_losses=false) for time0 in time]
 
-dyexp["summary.global_quantities.h_98.value"] =
+@store_expr dyexp["summary.global_quantities.h_98.value"] =
     (time; dd, summary, _...) -> summary.global_quantities.tau_energy.value ./ summary.global_quantities.tau_energy_98.value
 
 
-dyexp["summary.heating_current_drive.power_launched_ec.value"] =
+@store_expr dyexp["summary.heating_current_drive.power_launched_ec.value"] =
     (time; dd, summary, _...) -> sum(interp1d(beam.power_launched.time, beam.power_launched.data, :constant).(summary.time) for beam in dd.ec_launchers.beam)
 
-dyexp["summary.heating_current_drive.power_launched_ic.value"] =
+@store_expr dyexp["summary.heating_current_drive.power_launched_ic.value"] =
     (time; dd, summary, _...) -> sum(interp1d(antenna.power_launched.time, antenna.power_launched.data, :constant).(summary.time) for antenna in dd.ic_antennas.antenna)
 
-dyexp["summary.heating_current_drive.power_launched_lh.value"] =
+@store_expr dyexp["summary.heating_current_drive.power_launched_lh.value"] =
     (time; dd, summary, _...) -> sum(interp1d(antenna.power_launched.time, antenna.power_launched.data, :constant).(summary.time) for antenna in dd.lh_antennas.antenna)
 
-dyexp["summary.heating_current_drive.power_launched_nbi.value"] =
+@store_expr dyexp["summary.heating_current_drive.power_launched_nbi.value"] =
     (time; dd, summary, _...) -> sum(interp1d(unit.power_launched.time, unit.power_launched.data, :constant).(summary.time) for unit in dd.nbi.unit)
 
-dyexp["summary.heating_current_drive.power_launched_total.value"] =
+@store_expr dyexp["summary.heating_current_drive.power_launched_total.value"] =
     (time; dd, summary, _...) ->
         getproperty(dd.summary.heating_current_drive.power_launched_nbi, :value, zeros(length(summary.time))) .+
         getproperty(dd.summary.heating_current_drive.power_launched_ec, :value, zeros(length(summary.time))) .+
@@ -729,33 +729,33 @@ dyexp["summary.heating_current_drive.power_launched_total.value"] =
         getproperty(dd.summary.heating_current_drive.power_launched_lh, :value, zeros(length(summary.time)))
 
 
-dyexp["summary.local.magnetic_axis.t_e.value"] =
+@store_expr dyexp["summary.local.magnetic_axis.t_e.value"] =
     (time; dd, summary, _...) -> [dd.core_profiles.profiles_1d[Float64(time0)].electrons.temperature[1] for time0 in time]
 
-dyexp["summary.local.magnetic_axis.n_e.value"] =
+@store_expr dyexp["summary.local.magnetic_axis.n_e.value"] =
     (time; dd, summary, _...) -> [dd.core_profiles.profiles_1d[Float64(time0)].electrons.density[1] for time0 in time]
 
-dyexp["summary.local.magnetic_axis.t_i_average.value"] =
+@store_expr dyexp["summary.local.magnetic_axis.t_i_average.value"] =
     (time; dd, summary, _...) -> [dd.core_profiles.profiles_1d[Float64(time0)].t_i_average[1] for time0 in time]
 
-dyexp["summary.local.magnetic_axis.zeff.value"] =
+@store_expr dyexp["summary.local.magnetic_axis.zeff.value"] =
     (time; dd, summary, _...) -> [dd.core_profiles.profiles_1d[Float64(time0)].zeff[1] for time0 in time]
 
 
-dyexp["summary.local.separatrix.t_e.value"] =
+@store_expr dyexp["summary.local.separatrix.t_e.value"] =
     (time; dd, summary, _...) -> [dd.core_profiles.profiles_1d[Float64(time0)].electrons.temperature[end] for time0 in time]
 
-dyexp["summary.local.separatrix.n_e.value"] =
+@store_expr dyexp["summary.local.separatrix.n_e.value"] =
     (time; dd, summary, _...) -> [dd.core_profiles.profiles_1d[Float64(time0)].electrons.density[end] for time0 in time]
 
-dyexp["summary.local.separatrix.t_i_average.value"] =
+@store_expr dyexp["summary.local.separatrix.t_i_average.value"] =
     (time; dd, summary, _...) -> [dd.core_profiles.profiles_1d[Float64(time0)].t_i_average[end] for time0 in time]
 
-dyexp["summary.local.separatrix.zeff.value"] =
+@store_expr dyexp["summary.local.separatrix.zeff.value"] =
     (time; dd, summary, _...) -> [dd.core_profiles.profiles_1d[Float64(time0)].zeff[end] for time0 in time]
 
 
-dyexp["summary.volume_average.t_e.value"] =
+@store_expr dyexp["summary.volume_average.t_e.value"] =
     (time; dd, summary, _...) -> begin
         tmp = eltype(summary)[]
         for time in summary.time
@@ -765,7 +765,7 @@ dyexp["summary.volume_average.t_e.value"] =
         return tmp
     end
 
-dyexp["summary.volume_average.n_e.value"] =
+@store_expr dyexp["summary.volume_average.n_e.value"] =
     (time; dd, summary, _...) -> begin
         tmp = eltype(summary)[]
         for time in summary.time
@@ -775,7 +775,7 @@ dyexp["summary.volume_average.n_e.value"] =
         return tmp
     end
 
-dyexp["summary.volume_average.t_i_average.value"] =
+@store_expr dyexp["summary.volume_average.t_i_average.value"] =
     (time; dd, summary, _...) -> begin
         tmp = eltype(summary)[]
         for time in summary.time
@@ -785,7 +785,7 @@ dyexp["summary.volume_average.t_i_average.value"] =
         return tmp
     end
 
-dyexp["summary.volume_average.zeff.value"] =
+@store_expr dyexp["summary.volume_average.zeff.value"] =
     (time; dd, summary, _...) -> begin
         tmp = eltype(summary)[]
         for time in summary.time
diff --git a/src/expressions/expr_info.jl b/src/expressions/expr_info.jl
new file mode 100644
index 00000000..361b3a99
--- /dev/null
+++ b/src/expressions/expr_info.jl
@@ -0,0 +1,69 @@
+function IMASdd.get_expr_info_dict(::Type{Val{:onetime_and_dynamic}})
+    return expr_info_dict
+end
+
+struct ExprInfo
+    name::String
+    args::String
+    body::String
+end
+
+const expr_info_dict = Dict{String,ExprInfo}()
+
+function Base.show(io::IO, ::MIME"text/plain", expr_info::ExprInfo)
+    printstyled(io, "-"^length(expr_info.name) * "\n"; color=:red)
+    printstyled(io, "$(expr_info.name)\n"; color=:red, bold=true)
+    printstyled(io, "-"^length(expr_info.name) * "\n"; color=:red)
+    printstyled(io, "[Args]\n"; color=:blue)
+    print(io, "$(expr_info.args)\n")
+    printstyled(io, "[Body]\n"; color=:blue)
+    return print(io, "$(expr_info.body)\n")
+end
+
+"""
+    @store_expr expr
+
+Takes an assignment expression of the form `dict["key"] = (args...) -> body` and:
+
+ 1. Executes the original assignment.
+ 2. Extracts the argument list and body from the given anonymous function.
+ 3. Stores this information as an `ExprInfo` object in `expr_info_dict` under the specified key.
+"""
+macro store_expr(expr)
+    @assert expr.head == :(=)
+    lhs = expr.args[1]
+    rhs = expr.args[2]
+
+    @assert lhs.head == :ref && length(lhs.args) == 2
+    key_expr = lhs.args[2]
+
+    expr_args = sprint(show, MacroTools.prettify(rhs.args[1]))
+    expr_body = sprint(show, MacroTools.prettify(rhs.args[2]))
+
+    expr_info_dict = :expr_info_dict
+
+    return quote
+        # execute the given original expression
+        $(esc(expr))
+        # store the ExprInfo into expr_info_dict
+        $expr_info_dict[$(key_expr)] = ExprInfo($key_expr, $expr_args, $expr_body)
+    end
+end
+
+"""
+    @explain expr
+
+Temporarily sets a flag to explain the execution of `expr`. When `expr` runs under this macro,
+`IMASdd.FLAG_EXPLAIN_EXPRESSION` is enabled, allowing additional diagnostic or explanatory
+information to be displayed. After `expr` completes (successfully or not), the flag is reset.
+"""
+macro explain(expr)
+    return quote
+        IMASdd.FLAG_EXPLAIN_EXPRESSION[] = true
+        try
+            $(esc(expr))
+        finally
+            IMASdd.FLAG_EXPLAIN_EXPRESSION[] = false
+        end
+    end
+end
\ No newline at end of file
diff --git a/src/expressions/onetime.jl b/src/expressions/onetime.jl
index 5ad52004..48819ff7 100644
--- a/src/expressions/onetime.jl
+++ b/src/expressions/onetime.jl
@@ -27,24 +27,24 @@ otexp = onetime_expressions
 #= =========== =#
 # core_profiles #
 #= =========== =#
-otexp["core_profiles.profiles_1d[:].grid.psi_norm"] =
+@store_expr otexp["core_profiles.profiles_1d[:].grid.psi_norm"] =
     (rho_tor_norm; grid, _...) -> norm01(grid.psi)
 
-otexp["core_profiles.profiles_1d[:].grid.volume"] =
+@store_expr otexp["core_profiles.profiles_1d[:].grid.volume"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> begin
         eqt = dd.equilibrium.time_slice[Float64(profiles_1d.time)]
         volume = eqt.profiles_1d.volume
         return interp1d(eqt.profiles_1d.rho_tor_norm, sqrt.(volume), :cubic).(rho_tor_norm) .^ 2
     end
 
-otexp["core_profiles.profiles_1d[:].grid.area"] =
+@store_expr otexp["core_profiles.profiles_1d[:].grid.area"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> begin
         eqt = dd.equilibrium.time_slice[Float64(profiles_1d.time)]
         area = eqt.profiles_1d.area
         return interp1d(eqt.profiles_1d.rho_tor_norm, sqrt.(area), :cubic).(rho_tor_norm) .^ 2
     end
 
-otexp["core_profiles.profiles_1d[:].grid.surface"] =
+@store_expr otexp["core_profiles.profiles_1d[:].grid.surface"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> begin
         eqt = dd.equilibrium.time_slice[Float64(profiles_1d.time)]
         surface = eqt.profiles_1d.surface
@@ -52,7 +52,7 @@ otexp["core_profiles.profiles_1d[:].grid.surface"] =
     end
 
 
-otexp["core_profiles.profiles_1d[:].grid.psi"] =
+@store_expr otexp["core_profiles.profiles_1d[:].grid.psi"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> begin
         eqt = dd.equilibrium.time_slice[Float64(profiles_1d.time)]
         psi = eqt.profiles_1d.psi
@@ -63,37 +63,37 @@ otexp["core_profiles.profiles_1d[:].grid.psi"] =
 #= ============ =#
 # core_transport #
 #= ============ =#
-otexp["core_transport.model[:].profiles_1d[:].grid_flux.rho_tor_norm"] =
+@store_expr otexp["core_transport.model[:].profiles_1d[:].grid_flux.rho_tor_norm"] =
     (rho_tor_norm; profiles_1d, _...) -> profiles_1d.grid_d.rho_tor_norm
 
-otexp["core_transport.model[:].profiles_1d[:].grid_d.rho_tor_norm"] =
+@store_expr otexp["core_transport.model[:].profiles_1d[:].grid_d.rho_tor_norm"] =
     (rho_tor_norm; profiles_1d, _...) -> profiles_1d.grid_flux.rho_tor_norm
 
-otexp["core_transport.model[:].profiles_1d[:].grid_flux.psi_norm"] =
+@store_expr otexp["core_transport.model[:].profiles_1d[:].grid_flux.psi_norm"] =
     (rho_tor_norm; grid_flux, _...) -> norm01(grid_flux.psi)
 
-otexp["core_transport.model[:].profiles_1d[:].grid_flux.volume"] =
+@store_expr otexp["core_transport.model[:].profiles_1d[:].grid_flux.volume"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> begin
         eqt = dd.equilibrium.time_slice[Float64(profiles_1d.time)]
         volume = eqt.profiles_1d.volume
         return interp1d(eqt.profiles_1d.rho_tor_norm, volume, :cubic).(rho_tor_norm)
     end
 
-otexp["core_transport.model[:].profiles_1d[:].grid_flux.area"] =
+@store_expr otexp["core_transport.model[:].profiles_1d[:].grid_flux.area"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> begin
         eqt = dd.equilibrium.time_slice[Float64(profiles_1d.time)]
         area = eqt.profiles_1d.area
         return interp1d(eqt.profiles_1d.rho_tor_norm, area, :cubic).(rho_tor_norm)
     end
 
-otexp["core_transport.model[:].profiles_1d[:].grid_flux.surface"] =
+@store_expr otexp["core_transport.model[:].profiles_1d[:].grid_flux.surface"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> begin
         eqt = dd.equilibrium.time_slice[Float64(profiles_1d.time)]
         surface = eqt.profiles_1d.surface
         return interp1d(eqt.profiles_1d.rho_tor_norm, surface, :cubic).(rho_tor_norm)
     end
 
-otexp["core_transport.model[:].profiles_1d[:].grid_flux.psi"] =
+@store_expr otexp["core_transport.model[:].profiles_1d[:].grid_flux.psi"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> begin
         eqt = dd.equilibrium.time_slice[Float64(profiles_1d.time)]
         psi = eqt.profiles_1d.psi
@@ -103,31 +103,31 @@ otexp["core_transport.model[:].profiles_1d[:].grid_flux.psi"] =
 #= ============ =#
 #  core_sources  #
 #= ============ =#
-otexp["core_sources.source[:].profiles_1d[:].grid.psi_norm"] =
+@store_expr otexp["core_sources.source[:].profiles_1d[:].grid.psi_norm"] =
     (rho_tor_norm; grid, _...) -> norm01(grid.psi)
 
-otexp["core_sources.source[:].profiles_1d[:].grid.volume"] =
+@store_expr otexp["core_sources.source[:].profiles_1d[:].grid.volume"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> begin
         eqt = dd.equilibrium.time_slice[Float64(profiles_1d.time)]
         volume = eqt.profiles_1d.volume
         return interp1d(eqt.profiles_1d.rho_tor_norm, volume, :cubic).(rho_tor_norm)
     end
 
-otexp["core_sources.source[:].profiles_1d[:].grid.area"] =
+@store_expr otexp["core_sources.source[:].profiles_1d[:].grid.area"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> begin
         eqt = dd.equilibrium.time_slice[Float64(profiles_1d.time)]
         area = eqt.profiles_1d.area
         return interp1d(eqt.profiles_1d.rho_tor_norm, area, :cubic).(rho_tor_norm)
     end
 
-otexp["core_sources.source[:].profiles_1d[:].grid.surface"] =
+@store_expr otexp["core_sources.source[:].profiles_1d[:].grid.surface"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> begin
         eqt = dd.equilibrium.time_slice[Float64(profiles_1d.time)]
         surface = eqt.profiles_1d.surface
         return interp1d(eqt.profiles_1d.rho_tor_norm, surface, :cubic).(rho_tor_norm)
     end
 
-otexp["core_sources.source[:].profiles_1d[:].grid.psi"] =
+@store_expr otexp["core_sources.source[:].profiles_1d[:].grid.psi"] =
     (rho_tor_norm; dd, profiles_1d, _...) -> begin
         eqt = dd.equilibrium.time_slice[Float64(profiles_1d.time)]
         psi = eqt.profiles_1d.psi