-
Matthias Reich authored
Change-Id: I3676986428936b7adb6746453f6652e723bbac1d Reviewed-on: http://gerrit.aug.ipp.mpg.de/c/5169 Reviewed-by:
Matthias Reich <matthias.reich@ipp.mpg.de> Tested-by:
Matthias Reich <matthias.reich@ipp.mpg.de>
Matthias Reich authoredChange-Id: I3676986428936b7adb6746453f6652e723bbac1d Reviewed-on: http://gerrit.aug.ipp.mpg.de/c/5169 Reviewed-by:
Matthias Reich <matthias.reich@ipp.mpg.de> Tested-by:
Matthias Reich <matthias.reich@ipp.mpg.de>
augdemo.py 3.45 KiB
#!/usr/bin/env python
import sys
from numpy import size
from ctbm import default_inbeams, readTopfile3, readPrdata, call_torbeam, plotEq, plotPr
from tbmdata_aug import getECRHanglesAUG
from tbmdata_aug import ecrh2tbm
import matplotlib.pyplot as plt
# set some defaults
ECsystem = 0
ECsystem = 5 # if running without dd, please comment out this line!
time = 5.43
shot = 30594
if size(sys.argv) > 1:
try:
shot = int(sys.argv[1])
except Exception:
print "Usage:"
print sys.argv[0], " [shot] [time] [ECsys]"
print " "
print "Without arguments, script loads data from dataset augdata"
sys.exit(0)
if size(sys.argv) > 2:
time = float(sys.argv[2])
if size(sys.argv) > 3:
ECsystem = int(sys.argv[3])
# Now build data from shotfiles:
import tbmdata_aug as tbmdata
dat = tbmdata.tbmd(shot)
dat.readAll()
m, n, eqdata = dat.eqdata(time=time)
k, l, prdata = dat.prdata(time=time)
intinbeam, dblinbeam = dat.default_inbeams()
dblinbeam = dat.beamdata(time=time, gyr=ECsystem, beamdbl=dblinbeam)
else:
# Build eqdata from input files
m, n, eqdata = readTopfile3('augdata/topfile')
k, l, prdata = readPrdata('augdata/ne.dat', 'augdata/Te.dat')
intinbeam, dblinbeam = default_inbeams()
# Set launcher position
dblinbeam[3] = 236.1
dblinbeam[5] = 32.
# AUG specifics
if ECsystem == 7 or ECsystem == 8:
dblinbeam[5] = -32.
# read angles from shotfile or skip
if ECsystem > 0:
pol, tor = getECRHanglesAUG(shot, ECsystem, time=time)
dblinbeam[2], dblinbeam[1] = ecrh2tbm(pol, tor, ECsystem, shot)
# Set beam power
dblinbeam[23] = 0.6806
# set working absorption routine (Westerhof=0,Farina=1)
intinbeam[11] = 1
# and profile calc (std=0, ala Maj=1)
intinbeam[13] = 0
# set temporary testing values:
if intinbeam[11] == 0:
abstr = 'Westerhof'
if intinbeam[11] == 1:
abstr = 'Farina'
if intinbeam[13] == 0:
prstr = 'standard'
if intinbeam[13] == 1:
prstr = 'a la Maj'
# activate current drive
intinbeam[4] = 1
# Call the library
rhoresult, t1data, t1tdata, t2data, t2n_data, icnt, ibgout, nv, volprof = \
call_torbeam(intinbeam, dblinbeam, m, n, eqdata, k, l, prdata)
# generate figure
plt.figure(figsize=(10, 8))
# Plot result
plt1 = plt.subplot(2, 2, 1)
plotEq(eqm=m, eqn=n, eqdata=eqdata, plotbt=0, plteq=plt1)
plt1.plot(t1data[0, :]/100., t1data[1, :]/100., 'b')
plt1.plot(t1data[2, :]/100., t1data[3, :]/100., 'b')
plt1.plot(t1data[4, :]/100., t1data[5, :]/100., 'b')
plt1.plot([rhoresult[1]/100.], [rhoresult[2]/100.], 'k+', mew=5, ms=14)
textstr = r"$\rho=%.2f$\n$\mathrm{R}=%.2f m$\n$\mathrm{z}=%.2f m$" % \
(round(rhoresult[0], 2), round(rhoresult[1]/100, 2), round(rhoresult[2]/100, 2))
props = dict(boxstyle='round', facecolor='white')
plt1.text(0.95, 0.05, textstr, transform=plt1.transAxes, fontsize=14, horizontalalignment='right', bbox=props)
plt2 = plt.subplot(2, 2, 3)
plt2.plot(t2n_data[0, :], t2n_data[1, :], 'b')
plt2.set_title('P_EC vs rho')
plt3 = plt.subplot(2, 2, 4)
plt3.plot(t2n_data[0, :], t2n_data[2, :], 'b')
plt3.set_title('I_ECCD vs rho')
plt4 = plt.subplot(2, 2, 2)
plt4.set_title('kinetic profiles vs rho')
plotPr(prdata=prdata, prk=k, prl=l, pltpr=plt4)
textstr = 'tor=%.2f pol=%.2f\nabsorb:%s\nprofiles:%s' % (round(dblinbeam[1], 2), round(dblinbeam[2], 2), abstr, prstr)
plt2.text(0.5, 0.77, textstr, transform=plt1.transAxes, fontsize=14, horizontalalignment='right', bbox=props)
plt.tight_layout()
plt.show()