1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166 | 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()
|
