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 | '''Module for computing plasma position.'''
import numpy as np
from pygolem_lite import Shot
import pylab as pl
pl.ion()
#plasma_start
#plasma_end
def get_mirnov(shot):
'''Return signals from all Mirnov coils from given shot.
Call signature:
t, coils = get_mirnov(shot)'''
signal = Shot(shot)['papouch_st']
return np.array(signal[0]), np.array(signal[1]).T[0:6]
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 get_vacuum_shot(current):
return 21354
for i in range(1,20):
vacuum = Shot(int(current['shotno'] - i))
if (compare_shots(current, vacuum)):
return current['shotno'] - i
return None
def compare_shots(current, vacuum):
if (vacuum['plasma']):
return 0
diagnostics = ('reversed_bfield', 'reversed_efield', 'tb', 'tcd', 'ub', 'ucd')
for diag in diagnostics:
if (current[diag]!=vacuum[diag]):
return 0
return 1
def get_magnetic_field(shot):
'''Return magnetic field from given shot. Units in T, positive.
Call signature:
t, coils = get_magnetic_field(shot)'''
t, coils_plasma = integrate_coils(*get_mirnov(shot))
t, coils_vacuum = integrate_coils(*get_mirnov(get_vacuum_shot(shot)))
return t, (coils_plasma - coils_vacuum)
def compute_plasma_position(mc):
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(0)))
t, mc_pl = get_mirnov(0)
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='%f')
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='%f')
'''t, x_vert, x_horiz = compute_plasma_position(0)'''
#header = 'time[s]\tvertical position [m]\thorizontal position [m]'
return
if __name__ == '__main__':
main()
|
