# -*-Python-*-
# Created by bgrierson at 23 Mar 2017 07:10
# Compare IDL-based EFIT mapping with rho_rz to python methods
# for same shot and time along Z=0 midplane
# Set the tree for holding the results
root['OUTPUTS'].setdefault('EFIT', OMFITtree())
root['OUTPUTS']['EFIT'].setdefault('MAPPING', OMFITtree())
# ------------------
# Run the IDL RHO_RZ mapping
# ------------------
idl_cmds = '''
g = READG(163303,3000,RUNID='EFIT01')
Rmin = 0.84
Rmax = 2.54
nR = 201
dR = (Rmax-Rmin)/(nR-1.0)
R = DINDGEN(nR)*dR + Rmin
Z = DBLARR(nR)
rho_linear = RHO_RZ(g,R,Z,/DO_LINEAR)
rho_bicubic = RHO_RZ(g,R,Z)
psi_linear = PSI_RZ(g,R,Z,0)
psi_bicubic_conv = PSI_RZ(g,R,Z,1)
psi_bicubic_spline = PSI_RZ(g,R,Z,2)
psin_linear = PSI_RZ(g,R,Z,0,/NORM)
psin_bicubic_conv = PSI_RZ(g,R,Z,1,/NORM)
psin_bicubic_spline = PSI_RZ(g,R,Z,2,/NORM)
SAVE,R,Z,rho_bicubic,rho_linear,psi_linear,psi_bicubic_conv,psi_bicubic_spline,psin_linear,psin_bicubic_conv,psin_bicubic_spline,FILENAME='rho_rz.sav'
WRITEG,163303,3000,RUNID='EFIT01'
exit
'''
runScript = OMFITascii('run_rho_rz.pro', fromString=idl_cmds)
root['SETTINGS']['REMOTE_SETUP']['serverPicker'] = 'iris'
OMFITx.executable(
root,
inputs=[runScript],
outputs=['rho_rz.sav', 'g163303.03000'],
executable=IDE[root['SETTINGS']['REMOTE_SETUP']['serverPicker']]['idl'] + ' run_rho_rz.pro',
)
root['OUTPUTS']['EFIT']['MAPPING']['rho_rz'] = OMFITidlSav('rho_rz.sav')
root['OUTPUTS']['EFIT']['MAPPING']['gEQDSK'] = OMFITgeqdsk('g163303.03000')
# ------------------
# Run the python mappings
# ------------------
# EFIT coordinates
gR = root['OUTPUTS']['EFIT']['MAPPING']['gEQDSK']['AuxQuantities']['R']
gZ = root['OUTPUTS']['EFIT']['MAPPING']['gEQDSK']['AuxQuantities']['Z']
gRHORZ = root['OUTPUTS']['EFIT']['MAPPING']['gEQDSK']['AuxQuantities']['RHORZ']
gPSIRZ = root['OUTPUTS']['EFIT']['MAPPING']['gEQDSK']['PSIRZ']
gPSINRZ = root['OUTPUTS']['EFIT']['MAPPING']['gEQDSK']['AuxQuantities']['PSIRZ_NORM']
# Measurement locations
R = root['OUTPUTS']['EFIT']['MAPPING']['rho_rz']['r']
Z = root['OUTPUTS']['EFIT']['MAPPING']['rho_rz']['z']
# Regular grid interpolation
rhoRGI = RegularGridInterpolator((gR, gZ), gRHORZ.T, method='linear', bounds_error=False, fill_value=nan)((R, Z))
psiRGI = RegularGridInterpolator((gR, gZ), gPSIRZ.T, method='linear', bounds_error=False, fill_value=nan)((R, Z))
psinRGI = RegularGridInterpolator((gR, gZ), gPSINRZ.T, method='linear', bounds_error=False, fill_value=nan)((R, Z))
spline = 1
Ri = interp1e(gR, arange(len(gR)), kind=spline)(R)
Zi = interp1e(gZ, arange(len(gZ)), kind=spline)(Z)
rhoSpline1 = ndimage.map_coordinates(gRHORZ.T, (Ri, Zi), order=3, mode='nearest')
psiSpline1 = ndimage.map_coordinates(gPSIRZ.T, (Ri, Zi), order=3, mode='nearest')
psinSpline1 = ndimage.map_coordinates(gPSINRZ.T, (Ri, Zi), order=3, mode='nearest')
spline = 2
Ri = interp1e(gR, arange(len(gR)), kind=spline)(R)
Zi = interp1e(gZ, arange(len(gZ)), kind=spline)(Z)
rhoSpline2 = ndimage.map_coordinates(gRHORZ.T, (Ri, Zi), order=3, mode='nearest')
psiSpline2 = ndimage.map_coordinates(gPSIRZ.T, (Ri, Zi), order=3, mode='nearest')
psinSpline2 = ndimage.map_coordinates(gPSINRZ.T, (Ri, Zi), order=3, mode='nearest')
spline = 3
Ri = interp1e(gR, arange(len(gR)), kind=spline)(R)
Zi = interp1e(gZ, arange(len(gZ)), kind=spline)(Z)
rhoSpline3 = ndimage.map_coordinates(gRHORZ.T, (Ri, Zi), order=3, mode='nearest')
psiSpline3 = ndimage.map_coordinates(gPSIRZ.T, (Ri, Zi), order=3, mode='nearest')
psinSpline3 = ndimage.map_coordinates(gPSINRZ.T, (Ri, Zi), order=3, mode='nearest')
fig, ax = plt.subplots()
ax.plot(R, root['OUTPUTS']['EFIT']['MAPPING']['rho_rz']['rho_linear'], label='IDL RHO_RZ linear', marker='o', alpha=0.5, ms=4.0)
ax.plot(R, root['OUTPUTS']['EFIT']['MAPPING']['rho_rz']['rho_bicubic'], label='IDL RHO_RZ bicubic', marker='o', alpha=0.5, ms=4.0)
ax.plot(R, rhoRGI, label='python RGI', marker='o', alpha=0.5, ms=4.0)
ax.plot(R, rhoSpline1, label='python spline 1', marker='o', alpha=0.5, ms=4.0)
ax.plot(R, rhoSpline2, label='python spline 2', marker='o', alpha=0.5, ms=4.0)
ax.plot(R, rhoSpline3, label='python spline 3', marker='o', alpha=0.5, ms=4.0)
ax.set_title('Normalized rho')
ax.set_xlabel('R (m)')
ax.set_ylabel('$\\rho$')
ax.legend().draggable()
fig, ax = plt.subplots()
ax.plot(R, root['OUTPUTS']['EFIT']['MAPPING']['rho_rz']['psi_linear'], label='IDL PSI_RZ linear', marker='o', alpha=0.5, ms=4.0)
ax.plot(
R, root['OUTPUTS']['EFIT']['MAPPING']['rho_rz']['psi_bicubic_conv'], label='IDL PSI_RZ bicubic conv.', marker='o', alpha=0.5, ms=4.0
)
ax.plot(
R, root['OUTPUTS']['EFIT']['MAPPING']['rho_rz']['psi_bicubic_spline'], label='IDL PSI_RZ bicubic spline.', marker='o', alpha=0.5, ms=4.0
)
ax.plot(R, psiRGI, label='python RGI', marker='o', alpha=0.5, ms=4.0)
ax.plot(R, psiSpline1, label='python spline 1', marker='o', alpha=0.5, ms=4.0)
ax.plot(R, psiSpline2, label='python spline 2', marker='o', alpha=0.5, ms=4.0)
ax.plot(R, psiSpline3, label='python spline 3', marker='o', alpha=0.5, ms=4.0)
ax.set_title('Polodial Flux')
ax.set_xlabel('R (m)')
ax.set_ylabel('$\\psi$')
ax.legend().draggable()
fig, ax = plt.subplots()
ax.plot(R, root['OUTPUTS']['EFIT']['MAPPING']['rho_rz']['psin_linear'], label='IDL PSI_RZ linear', marker='o', alpha=0.5, ms=4.0)
ax.plot(
R, root['OUTPUTS']['EFIT']['MAPPING']['rho_rz']['psin_bicubic_conv'], label='IDL PSI_RZ bicubic conv.', marker='o', alpha=0.5, ms=4.0
)
ax.plot(
R,
root['OUTPUTS']['EFIT']['MAPPING']['rho_rz']['psin_bicubic_spline'],
label='IDL PSI_RZ bicubic spline.',
marker='o',
alpha=0.5,
ms=4.0,
)
ax.plot(R, psinRGI, label='python RGI', marker='o', alpha=0.5, ms=4.0)
ax.plot(R, psinSpline1, label='python spline 1', marker='o', alpha=0.5, ms=4.0)
ax.plot(R, psinSpline2, label='python spline 2', marker='o', alpha=0.5, ms=4.0)
ax.plot(R, psinSpline3, label='python spline 3', marker='o', alpha=0.5, ms=4.0)
ax.set_title('Normalized Polodial Flux')
ax.set_xlabel('R (m)')
ax.set_ylabel('$\\psi_N$')
ax.legend().draggable()
