Add MXH to IMAS.jl

src/physics/fluxsurfaces.jl

@@ -1650,7 +1650,7 @@ function flux_surfaces(eqt::equilibrium__time_slice{T1}, wall_r::AbstractVector{
 
     # trace flux surfaces
     Br, Bz = Br_Bz(eqt2d)
-    surfaces = trace_surfaces(eqt1d.psi, eqt1d.f, r, z, eqt2d.psi, Br, Bz, PSI_interpolant, RA, ZA, wall_r, wall_z)
+    surfaces = trace_surfaces(eqt1d.psi, eqt1d.f, r, z, eqt2d.psi, Br, Bz, PSI_interpolant, RA, ZA, wall_r, wall_z;refine_extrema=true)
 
     # calculate flux surface averaged and geometric quantities
     N = length(eqt1d.psi)
@@ -2342,3 +2342,62 @@ end
 
 @compat public areal_elongation
 push!(document[Symbol("Physics flux-surfaces")], :areal_elongation)
+
+
+
+"""
+    rec_to_mxh_coeffs(dd::IMAS.dd)
+
+fits rectangular flux surface mesh in dd to MXH
+"""
+function rec_to_mxh_coeffs(dd; MXH_modes=4, λ=10.0)
+    return rec_to_mxh_coeffs(dd.equilibrium.time_slice[], dd.wall; MXH_modes, λ)
+end
+
+function rec_to_mxh_coeffs(eqt::IMAS.equilibrium__time_slice, wall::IMAS.wall; MXH_modes=4, λ=10.0)
+
+    function extrap(x::Vector, u::Vector)
+        m = (u[3] - u[2]) / (x[3] - x[2])
+        return u[3] - m * x[3]
+    end
+    fw = first_wall(wall)
+
+    pnorm = eqt.profiles_1d.psi_norm
+    surfaces = trace_surfaces(eqt, fw.r, fw.z)
+    ns = length(surfaces)
+
+    nimas = 5 + 2*MXH_modes
+    coeffs = zeros(nimas, ns)
+    mxh_fits = Vector{Any}(nothing, ns)
+
+    for isurface in 1:ns
+        r,z  = surfaces[isurface].r, surfaces[isurface].z
+        # Inside-out with extrapolated prior
+        mxh_prev = if isurface <3
+            nothing
+        else
+            # Extrapolate from two already-fitted inner surfaces
+            x0, x1, x2 = pnorm[isurface], pnorm[isurface-1], pnorm[isurface-2]
+            mxh1, mxh2 = mxh_fits[isurface-1], mxh_fits[isurface-2]
+            c0_ext = mxh1.c0 + (mxh1.c0 - mxh2.c0) / (x1 - x2) * (x0 - x1)
+            c_ext  = mxh1.c  .+ (mxh1.c  .- mxh2.c)  ./ (x1 - x2) .* (x0 - x1)
+            s_ext  = mxh1.s  .+ (mxh1.s  .- mxh2.s)  ./ (x1 - x2) .* (x0 - x1)
+            MXH(mxh1.R0, mxh1.Z0, mxh1.ϵ, mxh1.κ, c0_ext, c_ext, s_ext)
+        end
+        mxh_fits[isurface] = MXH(r, z, MXH_modes; spline=false, optimize_fit=false, mxh_prev, λ)
+        coeffs[:, isurface]   = flat_coeffs(mxh_fits[isurface])
+    end
+
+    rnorm = sqrt.(pnorm)
+    coeffs[1, 1] = extrap(pnorm, coeffs[1, :])  # R0: linear in pnorm
+    coeffs[2, 1] = extrap(rnorm, coeffs[2, :])  # Z0: linear in rnorm
+    coeffs[3, 1] = 0.0                               # ϵ = 0 at axis
+    coeffs[4, 1] = extrap(rnorm, coeffs[4, :])  # κ: linear in rnorm
+    coeffs[5, 1] = extrap(pnorm, coeffs[5, :])  # c0: linear in pnorm
+    coeffs[6:end, 1] .= 0.0                          # Fourier terms = 0
+
+    return coeffs
+end
+
+@compat public rec_to_mxh_coeffs
+push!(document[Symbol("Physics flux-surfaces")], :rec_to_mxh_coeffs)

src/physics/interferometer.jl

@@ -61,7 +61,7 @@ function line_average(
     @assert ϕ1 == ϕ2 "line_average cannot handle different ϕ yet"
 
     # Create 1D interpolant for the quantity q
-    q_interp = cubic_interp1d(rho_tor_norm, q)
+    q_interp = linear_interp1d(rho_tor_norm, q)
 
     # Find intersections between line and LCFS
     intersections = intersection([r1, r2], [z1, z2], lcfs_r, lcfs_z)
@@ -83,7 +83,7 @@ function line_average(
         rho_seg = rho_interp.RHO_interpolant.(r_seg, z_seg)
 
         # Interpolate quantity q along the segment
-        q_seg = q_interp.(rho_seg)
+        q_seg = max.(q_interp.(rho_seg), zero(T))
 
         # Calculate path length inside LCFS
         path_length_inside_lcfs = sqrt((r_exit - r_entry)^2 + (z_exit - z_entry)^2)