#################################################################################################################### #################################################################################################################### #### #### #### 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', '$$ [$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 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[(casend)]) #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', '$$ [$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()