Skip to content
Snippets Groups Projects
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()