Other results ...
Procedure (This notebook to download)
In [1]:
import os
import numpy as np
import matplotlib.pyplot as plt
from scipy import integrate, signal, interpolate
import pandas as pd
import holoviews as hv
hv.extension('bokeh')
import hvplot.pandas
import requests
from IPython.display import Markdown
#execute time
from timeit import default_timer as timer
start_execute = timer()
In [2]:
# data_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/DASs/LimiterMirnovCoils/{identifier}.csv" #Mirnov coils and quadrupole
BDdata_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/DASs/StandardDAS/{identifier}.csv" #BD = basic diagnostic
# data_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/DASs/022Daniela0InnerStabilization/{identifier}.csv"
data_URL = "http://golem.fjfi.cvut.cz/shots/{shot_no}/DASs/PlasmaPosition/Nidatap_6133.lvm" #NI Standart (DAS)
shot_no = 33986 # to be replaced by the actual discharge number
# shot_no = 33947
vacuum_shot = 33980 # to be replaced by the discharge command line paramater "vacuum_shot"
# vacuum_shot = 33944 #number of the vacuum shot or 'False'
ds = np.DataSource(destpath='')
In [3]:
def open_remote(shot_no, identifier, url_template=data_URL):
return ds.open(url_template.format(shot_no=shot_no, identifier=identifier))
def read_signal(shot_no, identifier, url = data_URL, data_type='csv'):
file = open_remote(shot_no, identifier, url)
if data_type == 'lvm':
channels = ['Time', 'mc1', 'mc5', 'mc9', 'mc13', '5', '6', '7', '8', '9', '10', '11', '12']
return pd.read_table(file, sep = '\t', names=channels)
else:
return pd.read_csv(file, names=['Time', identifier],
index_col = 'Time', squeeze=True)
Data integration and $B_t$ elimination¶
In [4]:
def elimination (shot_no, identifier, vacuum_shot = False):
#load data
mc = (read_signal(shot_no, identifier, data_URL, 'lvm'))
mc = mc.set_index('Time')
mc = mc.replace([np.inf, -np.inf, np.nan], value = 0)
mc = mc[identifier]
konst = 1/(3.8e-03)
if vacuum_shot == False:
signal_start = mc.index[0]
length = len(mc)
Bt = read_signal(shot_no, 'BtCoil_integrated', BDdata_URL).loc[signal_start:signal_start+length*1e-06]
if len(Bt)>len(mc):
Bt = Bt.iloc[:length]
if len(Bt)<len(mc):
mc = mc.iloc[:len(Bt)]
if identifier == 'mc1':
k=300
elif identifier == 'mc5':
k= 14
elif identifier == 'mc9':
k = 31
elif identifier == 'mc13':
k = -100
mc_vacuum = Bt/k
else:
mc_vacuum = (read_signal(vacuum_shot, identifier, data_URL, 'lvm'))
mc_vacuum = mc_vacuum.set_index('Time')
mc_vacuum = mc_vacuum[identifier]
mc_vacuum -= mc_vacuum.loc[:0.9e-3].mean() #remove offset
mc_vacuum = mc_vacuum.replace([np.inf, -np.inf], value = 0)
mc_vacuum = pd.Series(integrate.cumtrapz(mc_vacuum, x=mc_vacuum.index, initial=0) * konst,
index=mc_vacuum.index*1000, name= identifier) #integration
mc -= mc.loc[:0.9e-3].mean() #remove offset
mc = pd.Series(integrate.cumtrapz(mc, x=mc.index, initial=0) * konst,
index=mc.index*1000, name= identifier) #integration
#Bt elimination
mc_vacuum = np.array(mc_vacuum)
mc_elim = mc - mc_vacuum
return mc_elim
Plasma life time¶
In [5]:
loop_voltage = read_signal(shot_no, 'LoopVoltageCoil_raw', BDdata_URL)
dIpch = read_signal(shot_no, 'RogowskiCoil_raw', BDdata_URL)
dIpch -= dIpch.loc[:0.9e-3].mean()
Ipch = pd.Series(integrate.cumtrapz(dIpch, x=dIpch.index, initial=0) * (-5.3*1e06),
index=dIpch.index, name='Ipch')
U_l_func = interpolate.interp1d(loop_voltage.index, loop_voltage)
def dIch_dt(t, Ich):
return (U_l_func(t) - 0.0097 * Ich) / (1.2e-6/2)
t_span = loop_voltage.index[[0, -1]]
solution = integrate.solve_ivp(dIch_dt, t_span, [0], t_eval=loop_voltage.index, )
Ich = pd.Series(solution.y[0], index=loop_voltage.index, name='Ich')
Ip = Ipch - Ich
Ip.name = 'Ip'
Ip_detect = Ip.loc[0.0025:]
dt = (Ip_detect.index[-1] - Ip_detect.index[0]) / (Ip_detect.index.size)
window_length = int(0.5e-3/dt)
if window_length % 2 == 0:
window_length += 1
dIp = pd.Series(signal.savgol_filter(Ip_detect, window_length, polyorder=3, deriv=1, delta=dt),
name='dIp', index=Ip_detect.index) / 1e6
threshold = 0.05
CD = requests.get("http://golem.fjfi.cvut.cz/shots/%i/Production/Parameters/CD_orientation" % shot_no)
CD_orientation = CD.text
if "ACW" in CD_orientation:
plasma_start = dIp[dIp < dIp.min()*threshold].index[0]*1e3
plasma_end = dIp.idxmax()*1e3
else:
plasma_start = dIp[dIp > dIp.max()*threshold].index[0]*1e3
plasma_end = dIp.idxmin()*1e3 + 0.5
print ('Plasma start =', round(plasma_start, 3), 'ms')
print ('Plasma end =', round(plasma_end, 3), 'ms')
# print (CD_orientation)
Horizontal plasma position $\Delta r$ calculation¶
In [6]:
def horpol(shot_no, vacuum_shot=False):
mc1 = elimination(shot_no, 'mc1', vacuum_shot)
mc9 = elimination (shot_no, 'mc9', vacuum_shot)
b = 93
r = ((mc1-mc9)/(mc1+mc9))*b
r = r.replace([np.nan], value = 0)
# return r.loc[plasma_start:]
return r.loc[plasma_start:plasma_end]
# return r
In [7]:
r = horpol(shot_no, vacuum_shot)
ax = r.plot()
ax.set(ylim=(-85,85), xlim=(plasma_start,plasma_end), xlabel= 'Time [ms]', ylabel = '$\Delta$r [mm]', title = 'Horizontal plasma position #{}'.format(shot_no))
ax.axhline(y=0, color='k', ls='--', lw=1, alpha=0.4)
Out[7]:
Vertical plasma position $\Delta z$ calculation¶
In [8]:
def vertpol(shot_no, vacuum_shot = False):
mc5 = elimination(shot_no, 'mc5', vacuum_shot)
mc13 = elimination (shot_no, 'mc13', vacuum_shot)
b = 93
z = ((mc5-mc13)/(mc5+mc13))*b
z = z.replace([np.nan], value = 0)
# return z.loc[plasma_start:]
return z.loc[plasma_start:plasma_end]
# return z
In [9]:
z = vertpol (shot_no, vacuum_shot)
ax = z.plot()
ax.set(ylim=(-85, 85), xlim=(plasma_start,plasma_end), xlabel= 'Time [ms]', ylabel = '$\Delta$z [mm]', title = 'Vertical plasma position #{}'.format(shot_no))
ax.axhline(y=0, color='k', ls='--', lw=1, alpha=0.4)
Out[9]:
Plasma column radius $a$ calculation¶
In [10]:
def plasma_radius(shot_no, vacuum_shot=False):
r = horpol(shot_no, vacuum_shot)
z = vertpol(shot_no, vacuum_shot)
a0 = 85
a = a0 - np.sqrt((r**2)+(z**2))
a = a.replace([np.nan], value = 0)
# return a.loc[plasma_start:]
return a.loc[plasma_start:plasma_end]
# return a
In [11]:
a = plasma_radius(shot_no,vacuum_shot)
ax = a.plot()
ax.set(ylim=(0,85), xlim=(plasma_start,plasma_end), xlabel= 'Time [ms]', ylabel = '$a$ [mm]', title = 'Plasma column radius #{}'.format(shot_no))
Out[11]:
In [12]:
plasma_time = []
t = 0
for i in a:
if i>85 or i <0:
a = a.replace(i, value = 0)
else:
plasma_time.append(a.index[t])
t+=1
start = plasma_time[0]-1e-03
end = plasma_time[-1]-1e-03
print('start =', round(start, 3), 'ms')
print('end =', round(end, 3), 'ms')
Graphs¶
In [13]:
r_cut = r.loc[start:end]
a_cut = a.loc[start:end]
z_cut = z.loc[start:end]
df_processed = pd.concat(
[r_cut.rename('r'), z_cut.rename('z'), a_cut.rename('a')], axis= 'columns')
df_processed
Out[13]:
In [14]:
savedata = 'plasma_position.csv'
df_processed.to_csv(savedata)
Data to download¶
In [15]:
Markdown("[Plasma position data - r, z, a ](./{})".format(savedata))
Out[15]:
In [16]:
hline = hv.HLine(0)
hline.opts(
color='k',
line_dash='dashed',
alpha = 0.4,
line_width=1.0)
layout = hv.Layout([df_processed[v].hvplot.line(
xlabel='', ylabel=l,ylim=(-85,85), xlim=(start,end),legend=False, title='', grid=True, group_label=v)
for (v, l) in [('r', ' r [mm]'), ('z', 'z [mm]'), ('a', 'a [mm]')] ])*hline
plot=layout.cols(1).opts(hv.opts.Curve(width=600, height=200),
hv.opts.Curve('a', xlabel='time [ms]'))
plot
Out[16]:
In [17]:
fig, axs = plt.subplots(3, 1, sharex=True, dpi=200)
r.plot(grid=True, ax=axs[0])
z.plot(grid=True, ax=axs[1])
a.plot(grid=True, ax=axs[2])
axs[2].set(ylim=(0,85), xlim=(start,end), xlabel= 'Time [ms]', ylabel = '$a$ [mm]')
axs[1].set(ylim=(-85,85), xlim=(start,end), xlabel= 'Time [ms]', ylabel = '$\Delta$z [mm]')
axs[0].set(ylim=(-85,85), xlim=(start,end), xlabel= 'Time [ms]', ylabel = '$\Delta$r [mm]', title = 'Horizontal, vertical plasma position and radius #{}'.format(shot_no))
plt.savefig('icon-fig')
Video¶
In [18]:
from matplotlib.patches import FancyArrowPatch
# from moviepy.editor import VideoClip
# from moviepy.video.io.bindings import mplfig_to_npimage
import matplotlib.animation as animation
from matplotlib.animation import PillowWriter
In [19]:
a = plasma_radius(shot_no,vacuum_shot).loc[plasma_start:plasma_end]
atime=np.array(a.index)
a = np.array(a)
r = np.array(r)
z = np.array(z)
n = 100 #the number of frames that will be skipped
FPS = 40
DURATION = len(a)/(FPS*n)
index = 0
x=0
y=0
Uloop = read_signal(shot_no, 'LoopVoltageCoil_raw', BDdata_URL).loc[plasma_start*1e-03:plasma_end*1e-03]
BDtime = np.array(Uloop.index)*1e03
Uloop = np.array(Uloop)
dBt = read_signal(shot_no, 'BtCoil_raw', BDdata_URL)
dBt -= dBt.loc[:0.9e-3].mean()
Btor = pd.Series(integrate.cumtrapz(dBt, x=dBt.index, initial=0) * 70.42,
index=dBt.index, name='Bt').loc[plasma_start*1e-03:plasma_end*1e-03]
# Btor = read_signal(shot_no, 'BtCoil_integrated', BDdata_URL).loc[plasma_start*1e-03:plasma_end*1e-03]
Bt = np.array(Btor)
responseBt = requests.get("http://golem.fjfi.cvut.cz/shots/%i/Production/Parameters/Bt_orientation" % shot_no)
Bt_orientation = responseBt.text
if 'ACW' in Bt_orientation:
Bt *= -1
dIpch = read_signal(shot_no, 'RogowskiCoil_raw', BDdata_URL)
dIpch -= dIpch.loc[:0.9e-3].mean()
Ipch = pd.Series(integrate.cumtrapz(dIpch, x=dIpch.index, initial=0) * (-5.3*1e06),
index=dIpch.index, name='Ipch')
Ip = Ipch - loop_voltage/0.0097
Ip = Ip.loc[plasma_start*1e-03:plasma_end*1e-03]
Ip = np.array(Ip)/1e03
if 'ACW' in CD_orientation:
Uloop *= -1
Ip *= -1
In [20]:
%%time
fig = plt.figure(figsize = (12,6), constrained_layout=True)
responsestab = requests.get("http://golem.fjfi.cvut.cz/shots/%i/Production/Parameters/PowSup4InnerStab" % shot_no)
stab = responsestab.text
if '-1' in stab or 'Not Found' in stab:
gs = plt.GridSpec(3, 7, figure=fig)
ax1 = fig.add_subplot(gs[0:,:4])
ax2 = fig.add_subplot(gs[0,4:])
ax3 = fig.add_subplot(gs[1,4:])
ax4 = fig.add_subplot(gs[2,4:])
else:
gs = plt.GridSpec(4, 8, figure=fig)
ax1 = fig.add_subplot(gs[0:,:5])
ax2 = fig.add_subplot(gs[0,5:])
ax3 = fig.add_subplot(gs[1,5:])
ax4 = fig.add_subplot(gs[2,5:])
ax5 = fig.add_subplot(gs[3,5:])
Irog = read_signal(shot_no, 'RogQuadr', data_URL).loc[start*1e-03:end*1e-03]
if abs(Irog.max()) > abs(Irog.min()): #orientation of the vertical magnetic field
IrogOr = 'up'
else:
IrogOr = 'down'
Irog = abs(Irog)
Irogmax = Irog.idxmax()*1e03
Irogtime = np.array(Irog.index)*1e03
Irog = np.array(Irog)*1/5e-03
alpha = 0
def make_frame(t):
global index
global x
global y
global alpha
time_a = atime[index]
ax1.clear()
ax2.clear()
ax3.clear()
ax4.clear()
ax1.set(xlim=(-100,100), ylim=(-100,100), xlabel = 'r [mm]', ylabel = 'z [mm]',
title = '#%i $t$ = %.3f ms'%(shot_no, time_a))
ax1.grid(linestyle='--')
ax1.set_aspect(1)
circle1 = plt.Circle((0, 0), 85, color='k', fill = False) #plot limiter
ax1.add_artist(circle1)
if a[index] < 0 or a[index] > 85:
radius = 0
x = 2000
y = 2000
else:
radius = a[index]
x = r[index]
y = z[index]
circle = plt.Circle(((0 + x), (0 + y)), radius, color='r', fill = False) #plot plasma column
ax1.add_artist(circle)
ax1.plot(0+x, 0+y, 'x', c = 'k')
ax1.legend([circle1,circle], ['limiter, $a_{0}$ = 85 mm',
'plasma column, $a$ = %i mm' %(radius)],
bbox_to_anchor=(1,0.8), loc="lower right")
time_i = np.array(BDtime[:index])
Uloop_i = np.array(Uloop[:index])
ax2.set(xlim=(plasma_start, plasma_end), ylim=(0,20), xlabel = 'time [ms]', ylabel = 'U$_{loop}$ [V]')
ax2.plot(time_i, Uloop_i, c='b')
ax2.grid(True)
Ip_i = np.array(Ip[:index])
Ip_max = Ip.max()
ax3.set(xlim=(plasma_start, plasma_end), ylim=(0,Ip_max), xlabel = 'time [ms]', ylabel = 'I$_{p}$ [kA]')
ax3.plot(time_i, Ip_i, c='r')
ax3.grid(True)
Bt_i = np.array(Bt[:index])
ax4.set(xlim=(plasma_start, plasma_end), ylim=(0,0.55), xlabel = 'time [ms]', ylabel = 'B$_{t}$ [T]')
ax4.plot(time_i, Bt_i, c='g')
ax4.grid(True)
#inner quadrupole
if not '-1' in stab and not 'Not Found' in stab:
ax5.clear()
Irog_i = np.array(Irog[:index])
Irogtime_i = np.array(Irogtime[:index])
ax5.set(xlim=(plasma_start, plasma_end), ylim=(0,500), xlabel = 'time [ms]', ylabel = 'I$_{quadr}$ [A]')
ax5.plot(Irogtime_i, Irog_i, c='k')
ax5.grid(True)
color = 'k'
if Irog[index] > 10 and Irogtime[index] < Irogmax:
alpha += ((4*n*1e-06)/(2*np.pi*np.sqrt(13)*1e-03))*2
#part of the quarter the period T for quadrupole; if index += 5 -> 5e-06
if alpha > 1:
alpha = 1
elif Irog[index]> 10 and Irogtime[index]>Irogmax:
alpha -= (4*n*1e-06)/(2*np.pi*np.sqrt(13)*1e-03)
if alpha < 0:
alpha = 0
else:
alpha = 0
if IrogOr == 'up':
A = (83, 83, 'x', 71, 70, 91)
B = (83, -83, 'x', -95, 70, -75)
C = (-83, -83, 'o', -95, -70, -75)
D = (-83, 83, 'o', 71, -70, 91)
if IrogOr == 'down':
A = (83, 83, 'o', 95, 70, 75)
B = (83, -83, 'o', -71, 70, -91)
C = (-83, -83, 'x', -71, -70, -91)
D = (-83, 83, 'x', 95, -70, 75)
for i in (A, B,C,D):
ax1.plot(i[0], i[1], i[2], c = color , alpha = alpha )
c = plt.Circle((i), 5, color= color, fill = False, alpha = alpha)
ax1.add_artist(c)
e = FancyArrowPatch((i[0], i[3]), (i[4], i[5]),
arrowstyle='->',
connectionstyle='angle3,angleA=0, angleB=90',
mutation_scale=10.0,
linewidth=2,
color=color,
alpha = alpha)
ax1.add_patch(e)
index += n
# return mplfig_to_npimage(fig) #have to be commented when matplotlib is used to make the animation
# animation in moviepy
# try:
# animation = VideoClip(make_frame, duration=DURATION)
# index = 0
# animation.write_videofile("video.mp4", fps=FPS)
# except IndexError:
# pass
#animation in matplotlib
try:
anim = animation.FuncAnimation(fig, make_frame, frames = 500)
anim.save('video.mp4', fps = FPS) #writer = ffmpeg
# anim.save('video.gif', writer = PillowWriter(fps = FPS))
except IndexError:
pass
Video to watch¶
In [21]:
from IPython.display import Video, Image
Video("video.mp4")
# Image("video.gif")
Out[21]:
In [22]:
end_execute = timer()
time = round(end_execute-start_execute,2)
print(time)
value = time
file = open('ExecutionTime.txt', 'w')
file.write(str(value)+' seconds'+'\n')
file.close()
Speed Camera¶
In [23]:
from PIL import Image
from io import BytesIO
from IPython.display import Image as DispImage
SpeedCamera = True
try:
url = "http://golem.fjfi.cvut.cz/shots/%i/Diagnostics/FastExCameras/plasma_film.png"%(shot_no)
image = requests.get(url)
img = Image.open(BytesIO(image.content)).convert('P') #'L','1'
data = np.asanyarray(img)
r = []
for i in range(data.shape[1]):
a=0
b=0
for j in range(336):
a += data[j,i]*j
b += data[j,i]
r.append((a/b))
camera_time = np.linspace(start, end, len(r))
r_camera = pd.Series(r, index = camera_time)-150
r_mirnov = horpol(shot_no, vacuum_shot).loc[start:end]
print(start, end, len(r))
fig = plt.figure()
ax = r_camera.plot(label = 'Speed Camera')
ax = r_mirnov.plot(label = 'Mirnov coils')
plt.legend()
ax.set(ylim=(-85, 85), title='#%i'%shot_no,ylabel='$\Delta$r visible radiation')
ax.axhline(y=0, color='k', ls='--', lw=1, alpha=0.4)
fig.savefig('SpeedCamera_deltaR')
except OSError:
SpeedCamera = False
In [24]:
if SpeedCamera == True:
plt.imshow(img)
In [ ]: