####################################################################################################################
####################################################################################################################
#### ####
#### 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 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)
#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))
save_adv('results/electron_density_LM', cas, shift)
#savetxt('out_shift.txt', c_[cas,shift],fmt=('%.7f','%.3e'))
if sys.argv[1] == "plots":
graphs()
saveconst('status', 0)
#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()
