#################################################################################################################### #################################################################################################################### #### #### #### 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.28E18 m^(-3). At least as long as the interferometer hardware is properly configured. # # # #################################################################################################################### #################################################################################################################### from scipy.optimize import curve_fit def linf(x,a,b): # linear function definition for curve fitting return a*x+b ################################################################## ### INPUT DATA LOADING ### ################################################################## zdroj = open('sig.txt','r') t,sig = loadtxt(zdroj,usecols=(0,1),unpack=True) zdroj.close() zdroj = open('pila.txt','r') t,pila = loadtxt(zdroj,usecols=(0,1),unpack=True) zdroj.close() start = 0.0178 # plasma appearence time end = 0.0261 # plasma disappearence time tcd = 0.015 # current drive trigger time ################################################################## ### 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 ################################################################## ### 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 zmereno = 0 # counter for identified sawtooth periods i = int(2E-6/dt)+1 # starting at the beginning would make a lookback more complicated # 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-700E-9/dt)+1]>OldGen*mean_pila and (len(pila_znacky)==0 or pila_znacky[len(pila_znacky)-1]+1E-6 < t[i])): # ...fit a 1400 ns window with a line... par,cov = curve_fit(linf, t[int(i-(700E-9)/dt+1):int(i+(700E-9)/dt-1)], pila[int(i-(700E-9)/dt+1):int(i+(700E-9)/dt-1)],[524000,0]) # 524 kHz - frequency estimate pila_znacky.append((mean_pila-par[1])/par[0]) # 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 ): sig_znacky.append(t[i-1]-sig[i-1]*(t[i]-t[i-1])/(sig[i]-sig[i-1])) sig_dole=0 if (sig[i] < -0.05): sig_dole=1 i+=1 # sawtooth frequency calculation f = 1/((pila_znacky[len(pila_znacky)-1]-pila_znacky[0])/(zmereno-1)) ################################################################## ### MAKING SAWTOOTH MARKS EQUIDISTANT ### ################################################################## # compensates for sawtooth generator frequency instability # not necessary, improves the result a little bit i = j = 0 while True: j+=50 if (j>=len(pila_znacky)): break par,cov = curve_fit(linf, range(i,j), pila_znacky[i:j],[1/f,0]) # 1/f - initial estimate for a in range(i,j): pila_znacky[a] = linf(a,par[0],par[1]) i = j+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 = [] # to store time... shift = [] # ...and phase shift data i = j = 0 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 ( 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 = len(shift)-2 if ( cas[len(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[len(shift)-1] offset = shift[len(shift)-1] shift[len(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 ( abs(diff) > lim): break lastshift = shift[j] shift[j]=shift[j]-offset-fringes*2*pi j-=1 if ( j > i-1): # remove the parts where the offset couldn't be determined for a in range(0,j-i+1): cas.pop(j-a) shift.pop(j-a) for i in range(0,len(shift)): shift[i] *= (3.29E18)/(2*pi) # calibration (assuming 75 GHz wave and 0.17m limiter diameter) ################################################################## ### SMOOTHING THE OUTPUT ### ################################################################## p = 3 for i in range(int(p/2+1),len(shift)-int(p/2+1)-1): souc=0 for a in range(i-int(p/2),i-int(p/2)+p): souc+=shift[a] shift[i]=souc/p ################################################################## ### CREATING OUTPUT FILE ### ################################################################## outshift_f = open("out_shift.txt","w") for i in range(0,len(shift)-1): outshift_f.write(str(cas[i])) outshift_f.write(" ") outshift_f.write(str(shift[i])) outshift_f.write("\n") outshift_f.close()