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