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
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283 | ####################################################################################################################
####################################################################################################################
#### ####
#### GOLEM interferometer algorithm ####
#### ####
#### Lukas Matena 2015 ####
#### ####
####################################################################################################################
####################################################################################################################
# #
# The algorithm evaluates phase shift between two signals - a sawtooth used to modulate microwave generator and #
# a signal detected after the diagnostic and reference waves interfere. Plasma density can then be calculated. #
# The algorithm has to deal with situations when the signal is momentarily lost. The modulating sawtooth can be #
# either ascending or descending (so that the generator can be easily changed). Any 2pi fringes are eliminated. #
# #
# Limitations: the longer #
# sawtooth edge is not shorter than 1400 ns. The input channels have to be sampled simultaneously #
# so that its time values are exactly the same. #
# #
# Calibration: Assuming the microwave generator frequency is 75 GHz and the path through the plasma is 0.17 m #
# (as should be the case for GOLEM), phase shift change of 2pi means that average density changed by #
# about 3.28E18 m^(-3). At least as long as the interferometer hardware is properly configured. #
# #
####################################################################################################################
####################################################################################################################
from numpy import *
import matplotlib
matplotlib.rcParams['backend'] = 'Agg'
from matplotlib.pylab import *
from scipy import *
from time import time
from pygolem_lite import save_adv,load_adv,saveconst, Shot
from pygolem_lite.modules import multiplot, get_data, paralel_multiplot
import os
##################################################################
### INPUT DATA LOADING ###
##################################################################
def LoadData():
Data = Shot()
Bt_trigger = Data['Tb']
gd = Shot().get_data
if Shot()['shotno'] > 18674:
tvec, density2 = gd('any', 'interframpsign')
tvec, density1 = gd('any', 'interfdiodeoutput')
else:
tvec, density1 = gd('any', 'density1')
tvec, density2 = gd('any', 'density2')
start = Data['plasma_start']
end = Data['plasma_end']
density1 = asarray(density1)
density2 = asarray(density2)
return tvec, start, end, density2,density1,Bt_trigger
def main():
for path in ['graphs', 'results' ]:
if not os.path.exists(path):
os.mkdir(path)
#NOTE never store these files as ASCI!! it is huge and slow
#NOTE these files should be loaded from GOLEM database
if sys.argv[1] == "analysis":
t, start, end, pila,sig,Bt_trigger = LoadData()
##################################################################
### BASIC DATA ANALYSIS ###
##################################################################
dt = (t[-1]-t[0])/len(t) # calculating sampling period
t = linspace(t[0],t[-1],len(t))
# at higher sampling rates the time values
# are truncated - let's replace it with more precise numbers (equidistant sampling assumed)
# now determine whether the sawtooth is ascending or descending
OldGen = -1 if sum(np.diff(pila[:2000])) < 0 else 1
pila*= OldGen
pila-= mean(pila)
sig-= mean(sig)
##################################################################
### FINDING MARKS OF SAWTOOTH AND SIGNAL PERIODS ###
##################################################################
pila_znacky = [] # to store time marks of the sawtooth
sig_znacky = [] # to store time marks of the signal
## let's go through all the data and find marks of the sawtooth and signal periods
posledni_znacka = 0
#BUG the slowest part!!!
for i in xrange(int(2e-6/dt)+1,len(t)):
# first the sawtooth: if it crosses its mean in given direction and it is not just noise...
if pila[i]<0. and pila[i-1]>0. and pila[i-int(.7e-6/dt)]>0. and posledni_znacka+1e-6 < t[i]:
# ...fit a 1400 ns window with a line...
ind = slice(i-.7e-6/dt+1,i+.7e-6/dt-1)
par = polyfit(t[ind], pila[ind],1)
#NOTE it is much better to use a linear least square method than a nonlinear method
posledni_znacka = -par[1]/par[0]
# and find its intersection with the mean - that is the mark we are looking for
pila_znacky.append(posledni_znacka)
# now the diode signal - zero crossing up defines the mark, the lookback condition eliminates noise influence
# the precise mark position is calculated by linear interpolation
if sig[i] >= 0. and sig[i-1] < 0. and sig[i-int(.5e-6/dt)] < 0. :
sig_znacky.append(t[i-1]-sig[i-1]*(t[i]-t[i-1])/(sig[i]-sig[i-1]))
pila_znacky = asarray(pila_znacky)
sig_znacky = asarray(sig_znacky)
# sawtooth frequency calculation
f = 1/((pila_znacky[-1]-pila_znacky[0])/(len(pila_znacky)-1))
##################################################################
### MAKING SAWTOOTH MARKS EQUIDISTANT ###
##################################################################
# compensates the noise introduced by root finding
# improves the result a little bit
from scipy.signal import butter,filtfilt
def butter_lowpass(cutOff, fs, order=5):
nyq = 0.5 * fs
normalCutoff = cutOff / nyq
b, a = butter(order, normalCutoff, btype='low', analog = False)
return b, a
def butter_lowpass_filter(data, cutOff, fs, order=4):
b, a = butter_lowpass(cutOff, fs, order=order)
#print abs(roots( a))
y = filtfilt(b, a, data)
return y
#estimate the slow variation in the carrier frequency
cutOff = 1e5#Hz
dpila = np.diff(pila_znacky)-1/f
dpila[0] = 0
y = butter_lowpass_filter(dpila,cutOff, 1/dt)
pila_znacky = cumsum(y+1/f)+pila_znacky[0]
##################################################################
### CALCULATING PHASE SHIFT ###
##################################################################
# the marks are compared and the phase shift is evaluated
# the diode signal can disappear for a while, the algorithm has to deal with it
i = j = 0
cas = [] # to store time...
shift = [] # ...and phase shift data
while i<len(pila_znacky) and j<len(sig_znacky):
while pila_znacky[i]<sig_znacky[j] and i<len(pila_znacky)-1: i+=1
phase_shift = 2*pi*f*(pila_znacky[i]-sig_znacky[j])
if phase_shift <= 2*pi:
cas.append(pila_znacky[i])
shift.append(phase_shift)
j+=1
cas = asarray(cas)
shift = asarray(shift)
##################################################################
### REMOVING FRINGES AND CALIBRATION ###
##################################################################
# The fringes are removed in two phases forward and backward
lim = 0.5
i=0
if start > cas[0]: # if plasma appeares after the beginning
i = 1
fringes = 0
lastshift=shift[0]
offset=shift[0]
shift[0]=0
while i < len(cas) and cas[i]-cas[i-1] < 6E-6 : # walk through from the beginning
diff=shift[i]-lastshift
if diff > 2*pi-lim and diff < 2*pi+lim : fringes+=1 # increment fringes counter
elif -diff > 2*pi-lim and -diff < 2*pi+lim : fringes-=1 # decrement fringes counter
elif fabs(diff) > lim: # failure
break
lastshift = shift[i] # remember the uncorrected value
shift[i]-=offset+fringes*2*pi # remove the fringe and introduce correct offset
i+=1
j = len(shift)-2
if cas[-1]+0.001 > end and i!=len(cas): # now the same from the end (in case there is a gap)
fringes = 0
lastshift = shift[-1]
offset = shift[-1]
shift[-1] = 0
while cas[j]-cas[j+1] < 6E-6 and j >= i:
diff=shift[j]-lastshift
if diff > 2*pi-lim and diff < 2*pi+lim : fringes+=1
elif -diff > 2*pi-lim and -diff < 2*pi+lim : fringes-=1
elif fabs(diff) > lim: break
lastshift = shift[j]
shift[j]-=offset+fringes*2*pi
j-=1
if j > i-1 : # remove the parts where the offset couldn't be determined
cas = r_[cas[:i],cas[j:]]
shift = r_[shift[:i],shift[j:]]
shift*= 3.29E18/(2*pi)
##################################################################
### SMOOTHING THE OUTPUT ###
##################################################################
cut_off = 5e6 #Hz
shift = butter_lowpass_filter( shift, cut_off, 1/dt)
shift-= mean(shift[(cas<start)|(cas>end)]) #remove offset
##################################################################
### CREATING OUTPUT FILE ###
##################################################################
save_adv('results/electron_density', cas, shift)
saveconst('results/carrier_freq', abs(f))
#savetxt('out_shift.txt', c_[cas,shift],fmt=('%.7f','%.3e'))
if sys.argv[1] == "plots":
tvec, n_e = load_adv('results/electron_density')
data = get_data('electron_density', 'Average electron density', '$<n_e>$ [$10^{19}\,m^{-3}$]', data_rescale=1e-19)
multiplot(data, '' , 'graphs/electron_density', (9,3) )
paralel_multiplot(data, '', 'icon', (4,3), 40)
saveconst('status', 0)
#print 'done',time()-T
#plot(cas,shift)
#savefig('./graphs/electron_density.png')
#show()
if __name__ == "__main__":
main()
|
