SCRIPTS EFIT volint_volavg

# -*-Python-*-
# Created by grierson at 14 Aug 2018  07:32

"""
This script demonstrates a volume integral and volume average of a quantity
based on an EFIT equilibrium reconstruction

defaultVars parameters
----------------------
:param var: Variable to process. Can be ne, Te, Ti, omega, nimp, ni, Zeff, Qrad
"""

defaultVars(var='ne')

# For DIII-D, get the volume integrated, volume aveage and peaking factor for the electron density
device = 'DIII-D'
shot = 163303
time = 3500
SNAPfile = 'EFIT01'

if 'EFIT' not in root:
    root['EFIT'] = OMFIT.loadModule('EFIT')

if 'ZIPFITprofiles' not in root:
    root['ZIPFITprofiles'] = OMFIT.loadModule('ZIPFITprofiles')

# Get the EFIT
root['EFIT'].experiment(device=device, shot=shot, time=time)
root['EFIT']['SETTINGS']['PHYSICS']['SNAPfile'] = SNAPfile
root['EFIT']['SCRIPTS']['fromMDSplus'].run()

# Get the ZIPFITs
root['ZIPFITprofiles'].experiment(device=device, shot=shot, time=time)
root['ZIPFITprofiles']['SETTINGS']['SETUP']['MDS_ZipFit'] = True
root['ZIPFITprofiles']['SCRIPTS']['extractProfiles'].run()

# ####
# Now use the EFIT and profile to form our quantities
# ####
# Get gEQDSK from EFIT
gEQDSK = root['EFIT']['FILES']['gEQDSK']

# Get the quantity and it's x-coordinate rho
quant = root['ZIPFITprofiles']['PROFILES'][var]
rho = root['ZIPFITprofiles']['PROFILES']['rho']

# Now map our quantity in rho to the EFIT rho
quant_on_RHO = interpolate.interp1d(rho, quant)(gEQDSK['AuxQuantities']['RHO'])

# Perform the volume integral of our quantity
quant_vint = gEQDSK.volume_integral(quant_on_RHO)

# Perform the volume integral of unity to get the accumulated plasma volume
ones_vint = gEQDSK.volume_integral(np.ones(len(gEQDSK['AuxQuantities']['RHO'])))

# Form the volume average as the volume integral to that surface divided by the accumulated volume to that surface
quant_vavg = quant_vint / ones_vint

# Form the peaking factor var[0] / total VOLAVG(var)
quant_pf = quant[0] / quant_vavg[-1]

print('VOLINT({}): Total number of electrons = {}'.format(var, quant_vint[-1]))
print('VOLAVG({}): {}'.format(var, quant_vint[-1] / ones_vint[-1]))
print('Peaking Factor({}): {}'.format(var, quant_pf))
fig, ax = plt.subplots(nrows=3)

ax[0].plot(rho, quant, marker='o')
ax[0].set_xlabel('rho')
ax[0].set_ylabel(var)

ax[1].plot(gEQDSK['AuxQuantities']['RHO'], quant_vint, marker='o')
ax[1].set_xlabel('EFIT RHO')
ax[1].set_ylabel('VOLINT({})'.format(var))

ax[2].plot(gEQDSK['AuxQuantities']['RHO'], quant_vint / ones_vint, marker='o')
ax[2].set_xlabel('EFIT RHO')
ax[2].set_ylabel('VOLAVG({})'.format(var))

fig.tight_layout()