import numpy as np
#rom pygolem_lite import Shot
import pylab as pl
pl.ion()
import matplotlib.pyplot as plt
#plasma_start
#plasma_end
#For testing, the shot will be the following:
shot= 24867
#shot = 28305
# create data cache in the 'golem_cache' folder
ds = np.DataSource('golem_cache')
#Create a path to data and download and open the file
base_url = "http://golem.fjfi.cvut.cz/utils/data/"
##--------------------------------------------------
def Shot(number, identifier):
dat = ds.open(base_url+ str(number) + '/' + identifier)
sig = np.loadtxt(dat)
return sig
def get_mirnov(shotn):
'''Return signals from all Mirnov coils from given shot.
Call signature:
t, coils = get_mirnov(shot)'''
signal = Shot(shotn, "papouch_st")
return np.array(signal[:,0]), np.array(signal[:,1:7]).T
#return np.array(signal[0]), np.array(signal[1:]).T[0:6]
def get_vacuum_shot(current):
return int(21354)
for i in range(1,20):
vacuum = int (Shot(current-i, 'shotno'))
if (compare_shots(current, vacuum)):
return int(Shot(vacuum, 'shotno'))
return None
def compare_shots(current, vacuum):
if (Shot (vacuum, 'plasma')):
return 0
diagnostics = ('reversed_bfield', 'reversed_efield', 'tb', 'tcd', 'ub', 'ucd')
for diag in diagnostics:
if Shot(shot, diag)!=Shot(vacuum, diag):
return 0
return 1
##----------------------------------------------------
def integrate_coil(t, coil, coeff):
'''Return integrated signals from given coils from given shot.
Call signature:
t, coils = get_integrated_mirnov(shot, array_of_coil_signals)
Note: Unit of signal is T.'''
coil -= coil[0:4500].mean() #subtract drift
coil[4990:5060] = 0 #remove jump caused by thyristors
integrated_data = np.array([np.cumsum(coil)])
integrated_data *= (t[1] - t[0]) / coeff * (-1)
return integrated_data
def compute_plasma_position(mc):
#the following variables need to be added or substracted depending on the amount of MC working:
m1, m2, m3, m4, mSadd, m5 = mc
#vertical position
mirnDiam = 0.093 #radial position of Mirnov coils
denominator = m2 + m4
denominator[abs(denominator) < 1e-3] *= np.nan
x_vert = mirnDiam * (m2 - m4) / denominator
#horizontal position
Bz = m5 / -0.1052
a_L = 0.085
x_horiz = []
iterations = []
for i in np.arange(0, len(m1)):
a = 0.085
iteration = 0
if abs(m1[i] + m2[i] + m3[i] + m4[i]) < 1e-3:
x_horiz.append(np.nan)
iterations.append(0)
continue
while (a > 0) and (iteration<100):
Lambda = np.log(a/mirnDiam) - (2/(denominator[i]*0.2325))*(Bz[i] + (m1[i] - m3[i])/2) - 1
dr = (m3[i] - m1[i]) * mirnDiam / denominator[i] - (mirnDiam**2/(2*0.4)) * (np.log(a/mirnDiam) + (a**2/(mirnDiam**2) + 1)*(Lambda + 0.5) + 1)
epsilon = a - a_L + np.sqrt(dr**2 + x_vert[i]**2)
if abs(epsilon) < 1e-6:
break
a = a_L - np.sqrt(dr**2 + x_vert[i]**2)
iteration += 1
x_horiz.append(dr)
iterations.append(iteration)
return x_vert, x_horiz, iterations
def main():
coeff = (3.8e-3,)*4 + (1,)*2
t, mc_vac = get_mirnov(get_vacuum_shot(shot))
t, mc_pl = get_mirnov(shot)
mc = []
for i in range(0, mc_vac.shape[0]):
mc_vac[i] = integrate_coil(t, mc_vac[i], coeff[i])
mc_pl[i] = integrate_coil(t, mc_pl[i], coeff[i])
mc.append(mc_pl[i] - mc_vac[i])
np.savetxt('mc' + str(i) + '.txt', np.array((t, mc_vac[i], mc_pl[i], mc[i])).T, fmt='%s')
x_vert, x_horiz, iterations = compute_plasma_position(mc)
for i in range(len(x_vert)):
if (x_vert[i] > 0.1):
x_vert[i] = 0.1
if (x_vert[i] < -0.1):
x_vert[i] = -0.1
if (x_horiz[i] > 0.1):
x_horiz[i] = 0.1
if (x_horiz[i] < -0.1):
x_horiz[i] = -0.1
np.savetxt('plasma_position.txt', np.array((t, x_vert, x_horiz, iterations)).T, fmt='%s')
# #----------------------Plot----------------------------------------------
#
x_horiz= np.array(x_horiz)
plt.figure(1)
plt.subplot(211)
plt.plot(t, x_horiz, 'b-', label = 'R')
plt.subplot(212)
plt.plot(t, x_vert, 'r-', label = 'Z')
plt.legend()
plt.savefig('basicgraph.pdf')
plt.savefig('basicgraph.png')
plt.show()
'''t, x_vert, x_horiz = compute_plasma_position(0)'''
#header = 'time[s]\tvertical position [m]\thorizontal position [m]'
return
if __name__ == '__main__':
main()
