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 | ####################################################################################################################
####################################################################################################################
#### ####
#### 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 algorithm assumes that the sawtooth and signal frequency is around 524 kHz and that 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.29E18 m^(-3). At least as long as the interferometer hardware is properly configured. #
# #
####################################################################################################################
####################################################################################################################
from numpy import *
from scipy.signal import butter,filtfilt
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
import sys
from matplotlib.pylab import *
OldGen = 1 # 1..descending sawtooth, 2..ascending sawtooth
cut_Off = 13e5 # smoothing frequency (2e5-20e5)
def LoadData():
Data = Shot()
Bt_trigger = Data['Tb']
gd = Shot().get_data
if Shot()['shotno'] > 18674:
tvec, density1 = gd('any', 'interframpsign')
tvec, density2 = gd('any', 'interfdiodeoutput')
else:
tvec, density1 = gd('any', 'density1')
tvec, density2 = gd('any', 'density2')
start = Data['plasma_start']
end = Data['plasma_end']
return tvec, start, end, density1,density2,Bt_trigger
def main():
for path in ['graphs', 'results' ]:
if not os.path.exists(path):
os.mkdir(path)
if sys.argv[1] == "analysis":
print 'analysis'
t, start, end, pila ,sig, Bt_trigger = LoadData()
pila = pila.signal
sig = sig.signal
print("Nacteno")
zac = time()
##################################################################
### BASIC DATA ANALYSIS ###
##################################################################
t = asarray(t)
pila = asarray(pila)
sig = asarray(sig)
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)
mean_pila = mean(pila) # mean sawtooth voltage
mean_sig = mean(sig) # mean sawtooth voltage
pila-=mean_pila
sig-=mean_sig
pila*=OldGen # make the sawtooth descending
##################################################################
### FINDING MARKS OF SAWTOOTH AND SIGNAL PERIODS ###
##################################################################
pila_marks = np.empty((t[-1]-t[0])/1.5E-6) # to store time marks of the sawtooth
sig_marks = np.empty((t[-1]-t[0])/1.5E-6) # to store time marks of the signal
zmereno = 0 # counter for identified sawtooth periods
i = int(2E-6/dt)+1 # starting at the beginning would make a lookback more complicated
sig_index = 0
ns700=(700E-9)/dt
# let's go through all the data and find marks of the sawtooth and signal periods
while i < len(t) - int(1E-6/dt+1):
# 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[int(i-ns700+1)] > 0. and (zmereno==0 or pila_marks[zmereno-1]+1.5E-6 < t[i])):
# ...fit a 1400 ns window with a line...
ind1=int(i-ns700+1)
ind2=int(i+ns700)
sx=sum(t[ind1:ind2])
sy=sum(pila[ind1:ind2])
sxy=sum(t[ind1:ind2]*pila[ind1:ind2])
sxx=sum(t[ind1:ind2]*t[ind1:ind2])
slope = (-sx*sy+(ind2-ind1)*sxy)/((ind2-ind1)*sxx-sx*sx)
offset = (sy-slope*sx)/(ind2-ind1)
pila_marks[zmereno]=((-offset)/slope) # and find its intersection with the mean - that is the mark we are looking for
zmereno += 1 # identified periods counter
# 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[int(i-(500E-9)/dt)] < 0. and (sig_index==0. or sig_marks[sig_index-1]+1.5E-6 < t[i] )):
sig_marks[sig_index] = (t[i-1]-sig[i-1]*(t[i]-t[i-1])/(sig[i]-sig[i-1]))
sig_index+=1
i+=1
# sawtooth frequency calculation
f = 1/((pila_marks[zmereno-1]-pila_marks[0])/(zmereno-1))
##################################################################
### 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
cas = np.zeros(zmereno) # to store time...
shift = np.zeros(zmereno) # ...and phase shift data
i = j = 0
index=0
while i<min(zmereno,sig_index) and j<min(zmereno,sig_index):
while (pila_marks[i]<sig_marks[j] and i<zmereno-1): i+=1
phase_shift = 2*pi*f*(pila_marks[i]-sig_marks[j])
if ( phase_shift <= 2*pi):
cas[index] = pila_marks[i]
shift[index] = phase_shift
index+=1
i+=1
j+=1
##################################################################
### REMOVING FRINGES AND CALIBRATION ###
##################################################################
# The fringes are removed in two phases.
lim = pi/2
time_lim = 10E-6
lasttime = cas[0]
i=0
if ( start > cas[0]): # if plasma appeares after the beginning
i = 1
fringes = 0
while ( i < index and cas[i]-lasttime < time_lim ): # walk through from the beginning
diff=shift[i]+2*pi*fringes-shift[i-1]
if ( diff < -2*pi+lim ): fringes+=1 # increment fringes counter
elif ( diff > 2*pi-lim ): fringes-=1 # decrement fringes counter
if (abs(shift[i]+2*pi*fringes-shift[i-1]) < lim or (i-1+int(time_lim*f)>0 and min(abs(shift[i]+2*pi*fringes-shift[i-k]) for k in range(1,int(time_lim*f))) < lim)):
lasttime = cas[i]
shift[i]=shift[i]+fringes*2*pi # remove the fringe and introduce correct offset
else: shift[i]=shift[i-1]
i+=1
j = index-2
lasttime=cas[j]
if ( cas[index-1] > end+0.001 and i!=index): # now the same from the end (in case there is a gap)
fringes = 0
while ( lasttime-cas[j] < time_lim and j >= i ):
diff=shift[j]+2*pi*fringes-shift[j+1]
if ( diff > 2*pi-lim and diff < 2*pi+lim ): fringes+=1
elif ( -diff > 2*pi-lim and -diff < 2*pi+lim ): fringes-=1
if (abs(shift[j]+2*pi*fringes-shift[j+1]) < lim or (j+1+int(time_lim*f)<index and min(abs(shift[j]+2*pi*fringes-shift[j+k]) for k in range(1,int(time_lim*f))) < lim)):
lasttime = cas[j]
shift[j]=shift[j]+fringes*2*pi # remove the fringe and introduce correct offset
else: shift[j]=shift[j+1]
j-=1
a=index
if ( j > i-1): # remove the parts where the offset couldn't be determined
for a in range(i,index-j+i+1):
cas[a] = cas[a+j-i+1]
shift[a] = shift[a+j-i+1]
cas=resize(cas,a-1)
shift=resize(shift,a-1)
shift *= (3.29E18)/(2*pi) # calibration (assuming 75 GHz wave and 0.17m limiter diameter)
shift-= mean(shift[(cas<start)|(cas>end+0.002)])
##################################################################
### SMOOTHING THE OUTPUT ###
##################################################################
b, a = butter(4, 2*dt*cut_Off, btype='low', analog = False)
shift[0:i] = filtfilt(b, a, shift[0:i])
if (j>i-1): shift[i:-1] = filtfilt(b, a, shift[i:-1])
shift-= mean(shift[(cas<start)|(cas>end+0.002)]) #remove offset
print("Dokonceno v case ",time()-zac)
save_adv('results/electron_density', cas, shift)
if sys.argv[1] == "plots":
tvec, n_e = load_adv('results/electron_density')
data = get_data('electron_density_LM', 'Average electron density', '$<n_e>$ [$10^{18}\,m^{-3}$]', data_rescale=1e-18)
multiplot(data, '' , 'graphs/electron_density', (9,3) )
paralel_multiplot(data, '', 'icon', (4,3), 40)
saveconst('status', 0)
if __name__ == "__main__":
main()
|
