Convert DESC equilibrium into EFIT inputs (at least profiles and Psi)

[dpanici]

Would be interesting to also do the reverse transformation, where we take a DESC equilibrium and convert the profile to $p'(\chi)$ and $FF'(\chi)$

[dpanici]

This seems reasonable

import sys
sys.path.insert(1, "../../efit2desc")
from efit2desc import (compute_Bp_probe_signals_from_DESC, get_coilset_for_shot,
                        compute_flux_loop_signals_from_DESC, convert_EFIT_to_DESC)

import numpy as np
from desc.io import load 
from desc.coils import MixedCoilSet
from desc.grid import LinearGrid, QuadratureGrid
from desc.magnetic_fields import ToroidalMagneticField, SumMagneticField
# from desc.objectives import (FluxLoop, PointBMeasurement, ObjectiveFunction, 
#                              BoundaryError, ForceBalance, FixBoundaryR,
#                                FixBoundaryZ, FixPressure, FixCurrent, FixPsi)

gfile = "g166439.04400"
eq, efit = convert_EFIT_to_DESC(gfile, L=16, M=16)

fluxsurf = efit["fluxSurfaces"]
efit_chi = fluxsurf["geo"]["psi"]
fluxsurf.keys()


# compute FF', compare to efit

import matplotlib.pyplot as plt
from desc.io import load
# eq = load("rev_C_tok_DESC_output.h5")
from desc.integrals import surface_averages
from desc.grid import LinearGrid
g = LinearGrid(L=eq.L_grid,M=eq.M_grid,N=0,axis=False)
data = eq.compute(["R","B_phi","sqrt(g)","chi"],grid=g)
F = surface_averages(g,q=data["R"]*data["B_phi"],
                     sqrt_g = data["sqrt(g)"],expand_out=False)

chi = g.compress(data["chi"])
from desc.profiles import PowerSeriesProfile
F_fxn_chi = PowerSeriesProfile.from_values(x=chi,y=F,sym=False, order=24)
FFp_desc = F_fxn_chi(chi)*F_fxn_chi(chi,dr=1)

plt.figure()
plt.plot(chi, F_fxn_chi(chi))
plt.scatter(chi, F)

# F should be even I think, so FF' should be odd
plt.figure()
# plt.scatter(efit_chi, fluxsurf["avg"]["FFPRIM"])
efit_chi_shift = efit_chi - (np.min(efit_chi))
plt.scatter(efit_chi_shift, fluxsurf["avg"]["FFPRIM"],label="EFIT")
plt.plot(chi, FFp_desc,"k",label="DESC")
plt.legend()

# pressure
from desc.profiles import PowerSeriesProfile
p_DESC = eq.pressure(g.nodes[g.unique_rho_idx,0])
p_fxn_chi = PowerSeriesProfile.from_values(x=chi,y=p_DESC,sym=False, order=24)

plt.figure()
plt.plot(chi, p_fxn_chi(chi))
plt.scatter(chi, p_DESC)

plt.figure()
efit_chi_shift = efit_chi - (np.min(efit_chi))
plt.scatter(efit_chi_shift, fluxsurf["avg"]["PPRIME"],label="EFIT")
plt.plot(chi, p_fxn_chi(chi,dr=1),"k",label="DESC")
plt.legend()

FF'

p'