# -*-Python-*-
# Created by grierson at 26 Jul 2017 08:55
"""
This script demos the reading of an ADAS ADF12 file
"""
try:
from adaslib import read_adf12
except ImportError as _excp:
printe(repr(_excp))
OMFITx.End()
printi('ADASHOME: ' + os.environ['ADASHOME'])
printi('ADASCENT: ' + os.environ['ADASCENT'])
# The ADF12 file for H beam interaction with C+6
file = os.environ['ADASCENT'] + '/adf12/qef93#h/qef93#h_c6.dat'
# Choose data block for C-VI (8-7) at 5290.5 A
block = 5
#%%%%%%%%%%%%%%%%%%%%%%%%
# Beam neutral properties
#%%%%%%%%%%%%%%%%%%%%%%%%
# Set beam mass to deuterium
ab = 2.0
# Set beam injection energy to 81 kV
vaccel = 81e3
# Get beam energy in eV/amu
ein = vaccel / ab
#%%%%%%%%%%%%%%%%%%%%%%%%
# Plasma properties
#%%%%%%%%%%%%%%%%%%%%%%%%
# Background ion is deuterium
ai = 2.0
Zi = 1.0
# Impurity ion is carbon
az = 12.0
Zz = 6.0
# Electron density (cm**-3) and temperature (eV)
ne = 2e13
Te = 1e3
# Impurity density (cm**-3)
nz = 0.1e13
# Ion density and temperature
ni = ne - Zz * nz
Ti = 1e3
# Zeff
Zeff = (ni * Zi ** 2 + nz * Zz ** 2) / ne
# Magnetic field strengh (T)
Bmag = 2.0
#%%%%%%%%%%%%%%%%%%%%%%%%
# Read the data for single set of parameters
#%%%%%%%%%%%%%%%%%%%%%%%%
qeff = read_adf12(file=file, block=block, ein=ein, tion=Ti, dion=ni, zeff=Zeff, bmag=Bmag)
print('ADAS qeff = {} ph-cm**3/s'.format(qeff[0]))
#%%%%%%%%%%%%%%%%%%%%%%%%
# Read the data for a range of energy
#%%%%%%%%%%%%%%%%%%%%%%%%
vaccel_min = 20e3
vaccel_max = 100e3
nvaccel = 101
vaccel = linspace(vaccel_min, vaccel_max, nvaccel)
ein = vaccel / ab
qeff = read_adf12(file=file, block=block, ein=ein, tion=Ti, dion=ni, zeff=Zeff, bmag=Bmag)
#%%%%%%%%%%%%%%%%%%%%%%%%
# Repeat exercise from IDL
#%%%%%%%%%%%%%%%%%%%%%%%%
idl_cmd = '''
file = '{}'
block = {}
ab = {}
vaccel_min = {}
vaccel_max = {}
nvaccel = {}
dvaccel = (vaccel_max-vaccel_min)/(nvaccel-1)
vaccel = FINDGEN(nvaccel)*dvaccel + vaccel_min
ein = vaccel/ab
dion = {}
tion = {}
zeff = {}
bmag = {}
read_adf12,file=file,block=block,ein=ein,dion=dion,tion=tion,zeff=zeff,bmag=bmag,/ENERGY,DATA=qeff
OPENW,lun,'qeff.dat',/GET_LUN
PRINTF,lun,'data'
PRINTF,lun,transpose(qeff)
PRINTF,lun
FREE_LUN,lun
exit
'''.format(
file, block, ab, vaccel_min, vaccel_max, nvaccel, ni, Ti, Zeff, Bmag
)
# IDL batch file
runName = 'run_read_adf12.pro'
runScript = OMFITascii(runName, fromString=idl_cmd)
# Set to run on iris
serverPicker = 'iris'
serverOptions = SERVER[serverPicker]
idl = serverOptions['idl']
# Run the script
OMFITx.executable(root, inputs=[runScript], outputs=['qeff.dat'], executable='module load adas; ' + idl + ' ' + runName + '; ls')
# Gather the output
qeff_idl = OMFITasciitable('qeff.dat')
# Plot the comparison
fig, ax = plt.subplots()
ax.plot(ein * 1e-3, qeff, label='Python', lw=2.0, color='black')
ax.plot(ein * 1e-3, qeff_idl['data'], label='IDL', ls='dashed', lw=2.0, color='red')
ax.set_xlabel('Energy (keV/amu)')
ax.set_ylabel('Effective Emission Coefficient (ph-cm**3/s)')
ax.text(0.05, 0.9, file, transform=ax.transAxes)
ax.legend()
