####################################################################################################################
####################################################################################################################
#### ####
#### 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[i-1]>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<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
# the data acquisition for the interferometer is triggered with current drive
# it is neccessary to shift the result in time
for i in range(0,len(shift)):
cas[i]+=tcd
##################################################################
### 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 < 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()
