####################################################################################################################
####################################################################################################################
#### ####
#### 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 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 *
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 graphs():
#tvec, phase_pila = load_adv('results/phase_saw')
#tvec, phase_sinus = load_adv('results/phase_sinus')
#tvec, phase = load_adv('results/phase_substracted')
#tvec, phase_corr = load_adv('results/phase_corrected')
#tvec, amplitude = load_adv('results/amplitude_sinus')
tvec, n_e = load_adv('results/electron_density_LM')
#print mean(n_e)
#data = [[get_data([tvec,-phase_pila+mean(phase_pila)], 'phase 1', 'phase [rad]', xlim=[0,40], fmt="--"),
# get_data([tvec,-phase_sinus+mean(phase_sinus)], 'phase 2', 'phase [rad]', xlim=[0,40], fmt="--" ),
# get_data([tvec,phase], 'substracted phase', 'phase [rad]', xlim=[0,40], fmt="k" ),
# get_data([tvec,phase_corr], 'corrected phase', 'phase [rad]', xlim=[0,40], fmt="k:" )],
# get_data([tvec,amplitude], 'amplitude', 'amplitude [a.u.]', xlim=[0,40],ylim=[0,None] )]
# multiplot(data, '' , 'graphs/demodulation', (10,6) )
data = get_data('electron_density_LM', 'Average electron density', '$<n_e>$ [$10^{18}\,m^{-3}$]', data_rescale=1e-18)
multiplot(data, '' , 'graphs/electron_density_LM', (9,3) )
paralel_multiplot(data, '', 'icon', (4,3), 40)
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
#start = start.signal
#end = end.signal
print("Nacteno")
##################################################################
##################################################################
##################################################################
##################################################################
### BASIC DATA ANALYSIS ###
##################################################################
dt = (t[-1]-t[0])/len(t) # calculating sampling period
mean_pila = mean(pila) # mean sawtooth voltage
for i in range(1,len(t)-1): t[i]=t[0]+i*dt # 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
asc = 0
desc = 0
for i in range(1,min(len(pila),2000)):
if (pila[i]-pila[i-1]>0): asc+=1
if (pila[i]-pila[i-1]<0): desc+=1
OldGen = 1
if (asc>desc): OldGen=-1
# descending sawtooth => old generator is used
print("Uvod hotov")
##################################################################
### FINDING MARKS OF SAWTOOTH AND SIGNAL PERIODS ###
##################################################################
pila_znacky = np.empty((t[-1]-t[0])/1.5E-6) # to store time marks of the sawtooth
sig_znacky = 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 (OldGen*pila[i]<OldGen*mean_pila and OldGen*pila[i-1]>OldGen*mean_pila and OldGen*pila[int(i-ns700+1)]>OldGen*mean_pila and (zmereno==0 or pila_znacky[zmereno-1]+1E-6 < t[i])):
# ...fit a 1400 ns window with a line...
ind1=int(i-ns700+1)
ind2=int(i+ns700-1)+1
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_znacky[zmereno]=((mean_pila-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_znacky[sig_index-1]+1E-6 < t[i] )):
sig_znacky[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_znacky[zmereno-1]-pila_znacky[0])/(zmereno-1))
print("Znacky nalezeny")
##################################################################
### 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
print(zmereno,sig_index)
while i<min(zmereno,sig_index) and j<min(zmereno,sig_index):
while (pila_znacky[i]<sig_znacky[j] and i<zmereno-1): i+=1
phase_shift = 2*pi*f*(pila_znacky[i]-sig_znacky[j])
if ( phase_shift <= 2*pi):
cas[index] = pila_znacky[i]
shift[index] = phase_shift
index+=1
j+=1
##################################################################
### REMOVING FRINGES AND CALIBRATION ###
##################################################################
# The fringes are removed in two phases.
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 < index 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 ( abs(diff) > lim): # failure
break
lastshift = shift[i] # remember the uncorrected value
shift[i]=shift[i]-offset-fringes*2*pi # remove the fringe and introduce correct offset
i+=1
j = index-2
if ( cas[index-1]+0.001 > end and i!=index): # now the same from the end (in case there is a gap)
fringes = 0
lastshift = shift[index-1]
offset = shift[index-1]
shift[index-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 ( abs(diff) > lim): break
lastshift = shift[j]
shift[j]=shift[j]-offset-fringes*2*pi
j-=1
if ( j > i-1): # zero the parts where the offset couldn't be determined
for a in range(0,j-i+1):
shift[j-a] = 0
for i in range(0,index):
shift[i] *= (3.29E18)/(2*pi) # calibration (assuming 75 GHz wave and 0.17m limiter diameter)
cas.resize(index)
shift.resize(index)
##################################################################
### SMOOTHING THE OUTPUT ###
##################################################################
p = 5
for i in range(int(p/2+1),index-int(p/2+1)-1):
souc=0
for a in range(i-int(p/2),i-int(p/2)+p):
if (shift[a]==0):
souc=p*shift[i]
break
souc+=shift[a]
shift[i]=souc/p
##################################################################
##################################################################
##################################################################
#save_adv('results/phase_saw', tvec, phase_pila)
#save_adv('results/phase_sinus', tvec, phase_sinus)
#save_adv('results/phase_substracted', tvec, phase)
#save_adv('results/amplitude_sinus', tvec, amplitude)
#save_adv('results/phase_corrected', tvec, phase_new)
#save_adv('results/electron_density_line', tvec, n_e)
#saveconst('results/electron_density_mean', mean(n_e[ind_plasma]))
save_adv('results/electron_density_LM', cas, shift)
#phase_skip = abs(phase[0] - phase[-1])/2*pi
#negativity = 1-sum(phase[ind_plasma])/sum(abs(phase[ind_plasma]))
#cum_var = mean(abs(diff(phase[ind_plasma]))) / mean(abs(phase[ind_plasma]))
#saveconst('results/carrier_freq', abs(f_carrier))
#saveconst('results/harmonics_distortion', k)
#saveconst('results/norm_ampl', norm_ampl)
#saveconst('results/probability', p)
#saveconst('results/reliability', 0 )
#print 'carrier_freq ', f_carrier
#print 'harmonics_distortion ',k
#print 'norm_ampl ',norm_ampl
#print 'probability ',p
#print "cum_var", cum_var
#print "negativity", negativity
#print "phase_skip", phase_skip
if sys.argv[1] == "plots":
graphs()
saveconst('status', 0)
if __name__ == "__main__":
main()
