From ab69bd0619249503ec04371b60a6c65f10ee327d Mon Sep 17 00:00:00 2001 From: Amlakbekalu Achule Date: Thu, 18 Jun 2026 10:40:00 -0400 Subject: [PATCH 1/5] =?UTF-8?q?Added=20=CE=B4=5FFKR=20and=20=CE=B4=5Fvisco?= =?UTF-8?q?.?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- src/Tearing/InnerLayer/SLAYER/LayerThickness.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/Tearing/InnerLayer/SLAYER/LayerThickness.jl b/src/Tearing/InnerLayer/SLAYER/LayerThickness.jl index d2d649e0e..a169acdce 100644 --- a/src/Tearing/InnerLayer/SLAYER/LayerThickness.jl +++ b/src/Tearing/InnerLayer/SLAYER/LayerThickness.jl @@ -153,6 +153,9 @@ struct LayerWidths delta_s::ComplexF64 delta_s_m::Float64 d_beta::Float64 + δ_FKR::Float64 + δ_visco::Float64 + end """ @@ -168,6 +171,8 @@ Keyword arguments are forwarded to `riccati_del_s`. function slayer_layer_thickness(p::SLAYERParameters; kwargs...) dels_db = riccati_del_s(p; kwargs...) delta_s = dels_db * p.d_beta + δ_FKR = r_s · S^(−1/3) + δ_visco = r_s · S^(−1/3) · P_⊥^(1/6) return LayerWidths(p.ising, p.m, p.n, dels_db, delta_s, abs(delta_s), p.d_beta) From 6afe97e0e66266fadad6dbc03c0ca19fb24b8fc4 Mon Sep 17 00:00:00 2001 From: Amlakbekalu Achule Date: Fri, 24 Jul 2026 11:48:20 -0400 Subject: [PATCH 2/5] fixed delta_FKR and delta_visco in slayer_layer_thickness. --- src/Tearing/InnerLayer/SLAYER/LayerThickness.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/Tearing/InnerLayer/SLAYER/LayerThickness.jl b/src/Tearing/InnerLayer/SLAYER/LayerThickness.jl index a169acdce..3fcfb90c8 100644 --- a/src/Tearing/InnerLayer/SLAYER/LayerThickness.jl +++ b/src/Tearing/InnerLayer/SLAYER/LayerThickness.jl @@ -171,9 +171,9 @@ Keyword arguments are forwarded to `riccati_del_s`. function slayer_layer_thickness(p::SLAYERParameters; kwargs...) dels_db = riccati_del_s(p; kwargs...) delta_s = dels_db * p.d_beta - δ_FKR = r_s · S^(−1/3) - δ_visco = r_s · S^(−1/3) · P_⊥^(1/6) + δ_FKR = p.rs * p.lu^(-1.0/3.0) + δ_visco = δ_FKR * p.P_perp^(1.0/6.0) return LayerWidths(p.ising, p.m, p.n, dels_db, delta_s, abs(delta_s), - p.d_beta) + p.d_beta, δ_FKR, δ_visco) end From 064b627d4b7e2eae1497dc67b30386a3447dac40 Mon Sep 17 00:00:00 2001 From: Amlakbekalu Achule Date: Fri, 24 Jul 2026 11:55:21 -0400 Subject: [PATCH 3/5] added set_resistive_width_based_psihigh helper for resistive-layer-width-based psihigh truncation --- .../SLAYER/SetResistiveWidthPsihigh.jl | 137 ++++++++++++++++++ 1 file changed, 137 insertions(+) create mode 100644 src/Tearing/InnerLayer/SLAYER/SetResistiveWidthPsihigh.jl diff --git a/src/Tearing/InnerLayer/SLAYER/SetResistiveWidthPsihigh.jl b/src/Tearing/InnerLayer/SLAYER/SetResistiveWidthPsihigh.jl new file mode 100644 index 000000000..111dd8320 --- /dev/null +++ b/src/Tearing/InnerLayer/SLAYER/SetResistiveWidthPsihigh.jl @@ -0,0 +1,137 @@ +""" + set_resistive_width_based_psihigh(r; q_scan=2:1:40, n_tor=1, mu_i=2.0, zeff=1.0) + +Given a completed equilibrium/stability run `r`, computes the resistive-layer +based truncation boundary: finds rational surfaces via the iota spline, +computes each surface's real resistive layer width (SLAYER Riccati del_s +solve) and its visco-resistive comparison scale, and returns the more +conservative (smaller) psihigh at which neighboring surfaces first overlap +or a surface's own layer spans past psi=1. + +Returns `(psihigh, psihigh_dels, psihigh_visco)`. +""" +function set_resistive_width_based_psihigh(r; q_scan=2:1:40, n_tor=1, mu_i=2.0, zeff=1.0) + r_of_psi = r.equil.geometry.avg_r_spline + drdpsi = ψ -> r_of_psi(ψ; deriv=DerivOp(1)) + + x_knots = r.equil.profiles.xs + q_knots = r.equil.profiles.q_spline.y + iota_itp_real = cubic_interp(x_knots, 1.0 ./ q_knots; extrap=Extrap(:extend)) + q_of_psi_real = ψ -> 1.0 / iota_itp_real(ψ) + dqdpsi_real = ψ -> -iota_itp_real(ψ; deriv=DerivOp(1)) / iota_itp_real(ψ)^2 + + ψ1, ψ2 = x_knots[end-1], x_knots[end] + q1, q2 = q_of_psi_real(ψ1), q_of_psi_real(ψ2) + A = (q2 - q1) / (log(1-ψ1) - log(1-ψ2)) + ψ_last, q_last = ψ2, q2 + + q_of_psi_hybrid(ψ) = ψ <= ψ_last ? q_of_psi_real(ψ) : q_last - A*log(1-ψ) + A*log(1-ψ_last) + dqdpsi_hybrid(ψ) = ψ <= ψ_last ? dqdpsi_real(ψ) : A/(1-ψ) + + function find_psi_for_q(target_q) + if target_q <= q_last + ψgrid = range(x_knots[1], ψ_last, length=200_000) + vals = q_of_psi_real.(ψgrid) + for i in 1:length(ψgrid)-1 + y1, y2 = vals[i]-target_q, vals[i+1]-target_q + if y1*y2 <= 0 + lo, hi = ψgrid[i], ψgrid[i+1] + flo = y1 + for _ in 1:60 + mid = (lo+hi)/2 + fmid = q_of_psi_real(mid) - target_q + sign(fmid) == sign(flo) ? (lo,flo)=(mid,fmid) : (hi=mid) + end + return (lo+hi)/2, true + end + end + return nothing, false + else + x = (1-ψ_last) * exp((q_last - target_q)/A) + x < 1e-15 && return nothing, false + return 1.0 - x, true + end + end + + kfile_path = r.ctrl.kinetic_file + kf = h5open(kfile_path, "r") + psi_k = read(kf["psi"]); ne_k = read(kf["n_e"]); Ti_k = read(kf["T_i"]) + Te_k = read(kf["T_e"]); omE_k = read(kf["omega_E"]) + chie_k = read(kf["chi_e"]); chip_k = read(kf["chi_phi"]) + close(kf) + + mk_spline(y) = cubic_interp(psi_k, y; extrap=Extrap(:extend)) + ne_spl, Ti_spl, Te_spl = mk_spline(ne_k), mk_spline(Ti_k), mk_spline(Te_k) + omE_spl, chie_spl, chip_spl = mk_spline(omE_k), mk_spline(chie_k), mk_spline(chip_k) + omega_star(ne_val, dne_dr, T_eV, B, Z) = -(T_eV*E_CHG)/(Z*E_CHG*B) * (dne_dr/ne_val) + + bt = r.equil.params.bt0 + R0 = r.equil.params.rmean + + kept_psi = Float64[]; kept_w_dels = Float64[]; kept_w_visco = Float64[] + psihigh_dels, psihigh_visco = nothing, nothing + + for m in q_scan + ψ_m, found = find_psi_for_q(Float64(m)) + !found && continue + + rs_val = r_of_psi(ψ_m); q_val = q_of_psi_hybrid(ψ_m); q1_val = dqdpsi_hybrid(ψ_m) + dr_val_geo = drdpsi(ψ_m) + (!isfinite(dr_val_geo) || dr_val_geo == 0.0 || !isfinite(q1_val) || !isfinite(q_val)) && continue + + sval_r = try + r_based_shear(rs_val, q_val, q1_val, dr_val_geo) + catch + continue + end + + ne_val = ne_spl(ψ_m); Ti_val = Ti_spl(ψ_m); Te_val = Te_spl(ψ_m) + (!isfinite(ne_val) || !isfinite(Ti_val) || !isfinite(Te_val) || ne_val<=0 || Ti_val<=0 || Te_val<=0) && continue + + dne_dr = ne_spl(ψ_m; deriv=DerivOp(1)) / dr_val_geo + omega_e = omega_star(ne_val, dne_dr, Te_val, bt, -1.0) + omega_i = omega_star(ne_val, dne_dr, Ti_val, bt, 1.0) + omega = omE_spl(ψ_m) + chi_perp, chi_tor = chie_spl(ψ_m), chip_spl(ψ_m) + (!isfinite(chi_perp) || !isfinite(chi_tor) || chi_perp<=0 || chi_tor<=0) && continue + + p = try + slayer_parameters(; n_e=ne_val, t_e=Te_val, t_i=Ti_val, + omega=omega, omega_e=omega_e, omega_i=omega_i, + qval=q_val, sval_r=sval_r, bt=bt, rs=rs_val, R0=R0, + mu_i=mu_i, zeff=zeff, chi_perp=chi_perp, chi_tor=chi_tor, + m=Int(round(m)), n=n_tor, ising=length(kept_psi)+1) + catch + continue + end + + lw = slayer_layer_thickness(p) + isnan(lw.delta_s_m) && continue + + push!(kept_psi, ψ_m); push!(kept_w_dels, lw.delta_s_m); push!(kept_w_visco, lw.δ_visco) + + n_now = length(kept_psi) + if n_now >= 2 + right_prev_dels = kept_psi[n_now-1] + kept_w_dels[n_now-1]/2 + left_curr_dels = kept_psi[n_now] - kept_w_dels[n_now]/2 + right_prev_visco = kept_psi[n_now-1] + kept_w_visco[n_now-1]/2 + left_curr_visco = kept_psi[n_now] - kept_w_visco[n_now]/2 + + dels_stop = (left_curr_dels < right_prev_dels) || (right_prev_dels > 1.0) + visco_stop = (left_curr_visco < right_prev_visco) || (right_prev_visco > 1.0) + + if dels_stop && psihigh_dels === nothing + psihigh_dels = kept_psi[n_now-1] - kept_w_dels[n_now-1]/2 + end + if visco_stop && psihigh_visco === nothing + psihigh_visco = kept_psi[n_now-1] - kept_w_visco[n_now-1]/2 + end + dels_stop && visco_stop && break + end + end + + candidates = filter(!isnothing, [psihigh_dels, psihigh_visco]) + psihigh_final = isempty(candidates) ? nothing : minimum(candidates) + + return psihigh_final, psihigh_dels, psihigh_visco +end From 08d6150f71bc13804238d3da2efaefb21abf0bae Mon Sep 17 00:00:00 2001 From: Amlakbekalu Achule Date: Fri, 14 Aug 2026 11:13:12 -0400 Subject: [PATCH 4/5] Tearing - SLAYER - wire in set_resistive_width_based_psihigh, document delta_FKR/delta_visco fields --- src/InnerLayer/InnerLayer.jl | 2 ++ src/InnerLayer/SLAYER/LayerThickness.jl | 2 ++ src/InnerLayer/SLAYER/SLAYER.jl | 4 ++++ 3 files changed, 8 insertions(+) diff --git a/src/InnerLayer/InnerLayer.jl b/src/InnerLayer/InnerLayer.jl index bee2a8a5e..8c8c31012 100644 --- a/src/InnerLayer/InnerLayer.jl +++ b/src/InnerLayer/InnerLayer.jl @@ -24,6 +24,7 @@ import .GGJ: solve_inner_converged # experimental, not exported (reachable as a import .SLAYER: SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear import .SLAYER: riccati_del_s, slayer_layer_thickness, LayerWidths import .SLAYER: surface_minor_radius, surface_da_dpsi, build_slayer_inputs +import .SLAYER: set_resistive_width_based_psihigh, print_layer_width_table export InnerLayerModel, InnerLayerParameters, InnerLayerResponse, solve_inner export GGJ, GGJModel, GGJParameters @@ -34,5 +35,6 @@ export glasser_wang_2020_eq55 export SLAYER, SLAYERModel, SLAYERParameters, slayer_parameters, r_based_shear export riccati_del_s, slayer_layer_thickness, LayerWidths export surface_minor_radius, surface_da_dpsi, build_slayer_inputs +export set_resistive_width_based_psihigh, print_layer_width_table end # module InnerLayer diff --git a/src/InnerLayer/SLAYER/LayerThickness.jl b/src/InnerLayer/SLAYER/LayerThickness.jl index cfea1ada5..d59f5273e 100644 --- a/src/InnerLayer/SLAYER/LayerThickness.jl +++ b/src/InnerLayer/SLAYER/LayerThickness.jl @@ -137,6 +137,8 @@ is built from, retained as a drift-scale reference. - `delta_s` -- complex layer thickness `δ_s = dels_db · d_β` [m] - `delta_s_m` -- `|δ_s|`, the resistive layer thickness [m] (primary) - `d_beta` -- β-weighted ion scale `c_β·d_i` [m] (drift reference) + - `δ_FKR` -- FKR resistive-inertial length scale, `r_s · S^(-1/3)` [m] + - `δ_visco` -- visco-resistive length scale, `δ_FKR · P_perp^(1/6)` [m] `delta_s_m` should sit within a few orders of magnitude of `d_beta` for a well-posed surface (`dels_db` is O(1)); a large gap flags a normalisation diff --git a/src/InnerLayer/SLAYER/SLAYER.jl b/src/InnerLayer/SLAYER/SLAYER.jl index ece897c15..0e9a80dcd 100644 --- a/src/InnerLayer/SLAYER/SLAYER.jl +++ b/src/InnerLayer/SLAYER/SLAYER.jl @@ -21,6 +21,8 @@ module SLAYER using LinearAlgebra using StaticArrays +using FastInterpolations +using HDF5 import ..InnerLayerModel, ..InnerLayerResponse, ..solve_inner, ..InnerLayerParameters using ...Utilities.PhysicalConstants @@ -52,11 +54,13 @@ include("LayerParameters.jl") include("Riccati.jl") include("LayerThickness.jl") include("LayerInputs.jl") +include("SetResistiveWidthPsihigh.jl") export SLAYERModel, SLAYERParameters, slayer_parameters export r_based_shear export riccati_del_s, slayer_layer_thickness, LayerWidths export surface_minor_radius, surface_da_dpsi, build_slayer_inputs +export set_resistive_width_based_psihigh, print_layer_width_table export NeoResistivityModel, SpitzerModel, SpitzerHarmModel, SauterNeoModel, RedlNeoModel end # module SLAYER From 1ea7de32b9f51ec5f80bfe970dd73aaa2b400f73 Mon Sep 17 00:00:00 2001 From: Amlakbekalu Achule Date: Fri, 14 Aug 2026 11:43:40 -0400 Subject: [PATCH 5/5] Tearing - SLAYER - use real build_slayer_inputs helpers (surface_minor_radius, surface_da_dpsi) and KineticProfiles instead of manual reconstruction --- .../SLAYER/SetResistiveWidthPsihigh.jl | 115 +++++++++++------- 1 file changed, 68 insertions(+), 47 deletions(-) diff --git a/src/InnerLayer/SLAYER/SetResistiveWidthPsihigh.jl b/src/InnerLayer/SLAYER/SetResistiveWidthPsihigh.jl index 111dd8320..5ca01ae71 100644 --- a/src/InnerLayer/SLAYER/SetResistiveWidthPsihigh.jl +++ b/src/InnerLayer/SLAYER/SetResistiveWidthPsihigh.jl @@ -1,21 +1,8 @@ -""" - set_resistive_width_based_psihigh(r; q_scan=2:1:40, n_tor=1, mu_i=2.0, zeff=1.0) - -Given a completed equilibrium/stability run `r`, computes the resistive-layer -based truncation boundary: finds rational surfaces via the iota spline, -computes each surface's real resistive layer width (SLAYER Riccati del_s -solve) and its visco-resistive comparison scale, and returns the more -conservative (smaller) psihigh at which neighboring surfaces first overlap -or a surface's own layer spans past psi=1. - -Returns `(psihigh, psihigh_dels, psihigh_visco)`. -""" -function set_resistive_width_based_psihigh(r; q_scan=2:1:40, n_tor=1, mu_i=2.0, zeff=1.0) - r_of_psi = r.equil.geometry.avg_r_spline - drdpsi = ψ -> r_of_psi(ψ; deriv=DerivOp(1)) - - x_knots = r.equil.profiles.xs - q_knots = r.equil.profiles.q_spline.y +function set_resistive_width_based_psihigh(r; q_scan=2:1:40, n_tor=1, mu_i=2.0, zeff=1.0, z_i=1.0) + equil = r.equil + + x_knots = equil.profiles.xs + q_knots = equil.profiles.q_spline.y iota_itp_real = cubic_interp(x_knots, 1.0 ./ q_knots; extrap=Extrap(:extend)) q_of_psi_real = ψ -> 1.0 / iota_itp_real(ψ) dqdpsi_real = ψ -> -iota_itp_real(ψ; deriv=DerivOp(1)) / iota_itp_real(ψ)^2 @@ -53,54 +40,68 @@ function set_resistive_width_based_psihigh(r; q_scan=2:1:40, n_tor=1, mu_i=2.0, end end - kfile_path = r.ctrl.kinetic_file - kf = h5open(kfile_path, "r") - psi_k = read(kf["psi"]); ne_k = read(kf["n_e"]); Ti_k = read(kf["T_i"]) - Te_k = read(kf["T_e"]); omE_k = read(kf["omega_E"]) - chie_k = read(kf["chi_e"]); chip_k = read(kf["chi_phi"]) - close(kf) + kfile = joinpath("examples", "DIIID-like_ideal_example", "TkMkr_D3Dlike_Hmode_kinetic.h5") + psi_k, ne_k, Ti_k, Te_k, omE_k, chie_k, chip_k = h5open(kfile, "r") do kf + (read(kf["psi"]), read(kf["n_e"]), read(kf["T_i"]), read(kf["T_e"]), + read(kf["omega_E"]), read(kf["chi_e"]), read(kf["chi_phi"])) + end - mk_spline(y) = cubic_interp(psi_k, y; extrap=Extrap(:extend)) - ne_spl, Ti_spl, Te_spl = mk_spline(ne_k), mk_spline(Ti_k), mk_spline(Te_k) - omE_spl, chie_spl, chip_spl = mk_spline(omE_k), mk_spline(chie_k), mk_spline(chip_k) - omega_star(ne_val, dne_dr, T_eV, B, Z) = -(T_eV*E_CHG)/(Z*E_CHG*B) * (dne_dr/ne_val) + kp = KineticProfiles( + psi = psi_k, + n_e = ne_k, + T_e = Te_k, + T_i = Ti_k, + omega = omE_k, + omega_e = zeros(length(psi_k)), + omega_i = zeros(length(psi_k)), + ) + chie_spl = cubic_interp(psi_k, chie_k; extrap=Extrap(:extend)) + chip_spl = cubic_interp(psi_k, chip_k; extrap=Extrap(:extend)) + + chi1 = 2π * equil.psio + function omega_star_at(ψ) + n_e_val = kp.n_e(ψ); dn_e = kp.n_e(ψ; deriv=DerivOp(1)) + T_e_val = kp.T_e(ψ); dT_e = kp.T_e(ψ; deriv=DerivOp(1)) + T_i_val = kp.T_i(ψ); dT_i = kp.T_i(ψ; deriv=DerivOp(1)) + ω_star_e = (2π / chi1) * (T_e_val * dn_e / n_e_val + dT_e) + ω_star_i = -(2π / (z_i * chi1)) * (T_i_val * dn_e / n_e_val + dT_i) + return (ω_star_e, ω_star_i) + end - bt = r.equil.params.bt0 - R0 = r.equil.params.rmean + R0 = equil.ro + bt_at(ψ) = Float64(equil.profiles.F_spline(ψ)) / (2π * R0) - kept_psi = Float64[]; kept_w_dels = Float64[]; kept_w_visco = Float64[] + kept_r, kept_q, kept_psi = Float64[], Float64[], Float64[] + kept_w_dels, kept_w_visco = Float64[], Float64[] psihigh_dels, psihigh_visco = nothing, nothing for m in q_scan ψ_m, found = find_psi_for_q(Float64(m)) !found && continue - rs_val = r_of_psi(ψ_m); q_val = q_of_psi_hybrid(ψ_m); q1_val = dqdpsi_hybrid(ψ_m) - dr_val_geo = drdpsi(ψ_m) - (!isfinite(dr_val_geo) || dr_val_geo == 0.0 || !isfinite(q1_val) || !isfinite(q_val)) && continue + rs_val = surface_minor_radius(equil, ψ_m) + q_val = q_of_psi_hybrid(ψ_m) + q1_val = dqdpsi_hybrid(ψ_m) + da_dpsi = surface_da_dpsi(equil, ψ_m) + (!isfinite(da_dpsi) || da_dpsi == 0.0 || !isfinite(q1_val) || !isfinite(q_val)) && continue sval_r = try - r_based_shear(rs_val, q_val, q1_val, dr_val_geo) + r_based_shear(rs_val, q_val, q1_val, da_dpsi) catch continue end - ne_val = ne_spl(ψ_m); Ti_val = Ti_spl(ψ_m); Te_val = Te_spl(ψ_m) - (!isfinite(ne_val) || !isfinite(Ti_val) || !isfinite(Te_val) || ne_val<=0 || Ti_val<=0 || Te_val<=0) && continue - - dne_dr = ne_spl(ψ_m; deriv=DerivOp(1)) / dr_val_geo - omega_e = omega_star(ne_val, dne_dr, Te_val, bt, -1.0) - omega_i = omega_star(ne_val, dne_dr, Ti_val, bt, 1.0) - omega = omE_spl(ψ_m) + prof = kp(ψ_m) + ω_e_use, ω_i_use = omega_star_at(ψ_m) chi_perp, chi_tor = chie_spl(ψ_m), chip_spl(ψ_m) (!isfinite(chi_perp) || !isfinite(chi_tor) || chi_perp<=0 || chi_tor<=0) && continue p = try - slayer_parameters(; n_e=ne_val, t_e=Te_val, t_i=Ti_val, - omega=omega, omega_e=omega_e, omega_i=omega_i, - qval=q_val, sval_r=sval_r, bt=bt, rs=rs_val, R0=R0, + slayer_parameters(; n_e=prof.n_e, t_e=prof.T_e, t_i=prof.T_i, + omega=prof.omega, omega_e=ω_e_use, omega_i=ω_i_use, + qval=q_val, sval_r=sval_r, bt=bt_at(ψ_m), rs=rs_val, R0=R0, mu_i=mu_i, zeff=zeff, chi_perp=chi_perp, chi_tor=chi_tor, - m=Int(round(m)), n=n_tor, ising=length(kept_psi)+1) + m=Int(round(m)), n=n_tor, ising=length(kept_q)+1) catch continue end @@ -108,7 +109,8 @@ function set_resistive_width_based_psihigh(r; q_scan=2:1:40, n_tor=1, mu_i=2.0, lw = slayer_layer_thickness(p) isnan(lw.delta_s_m) && continue - push!(kept_psi, ψ_m); push!(kept_w_dels, lw.delta_s_m); push!(kept_w_visco, lw.δ_visco) + push!(kept_r, rs_val); push!(kept_q, m); push!(kept_psi, ψ_m) + push!(kept_w_dels, lw.delta_s_m); push!(kept_w_visco, lw.δ_visco) n_now = length(kept_psi) if n_now >= 2 @@ -133,5 +135,24 @@ function set_resistive_width_based_psihigh(r; q_scan=2:1:40, n_tor=1, mu_i=2.0, candidates = filter(!isnothing, [psihigh_dels, psihigh_visco]) psihigh_final = isempty(candidates) ? nothing : minimum(candidates) + print_layer_width_table(kept_q, kept_psi, kept_w_dels, "Resistive layer width (Riccati ODE)") + print_layer_width_table(kept_q, kept_psi, kept_w_visco, "Visco-resistive comparison scale") + return psihigh_final, psihigh_dels, psihigh_visco end + +function print_layer_width_table(kept_q, kept_psi, kept_w, label) + println("\n$label") + println(rpad("q", 6), rpad("psi", 14), rpad("left edge", 14), rpad("right edge", 14), "width") + for i in eachindex(kept_q) + left_edge = kept_psi[i] - kept_w[i]/2 + right_edge = kept_psi[i] + kept_w[i]/2 + println( + rpad(string(kept_q[i]), 6), + rpad(string(round(kept_psi[i], digits=8)), 14), + rpad(string(round(left_edge, digits=8)), 14), + rpad(string(round(right_edge, digits=8)), 14), + string(round(kept_w[i], digits=10)) + ) + end +end