%pylab notebook
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
# %matplotlib inline
from importlib import reload
import pandas as pd
from scipy.optimize import curve_fit
import scipy.ndimage as flt
from scipy.stats import linregress
from matplotlib import mlab as mlab
import xarray as xr
import json
import scipy
from astropy.table import Table, Column, MaskedColumn
from astropy.io import ascii
from scipy.signal import butter, filtfilt, lfilter
from urllib.request import urlopen
def linear_fit(x, a, b):
return a*x+b
def my_butter_filter(cutoff, fs, order = 32):
nyq = 0.5 * fs
normal_cutoff = cutoff / nyq
b, a = butter(order, normal_cutoff, btype='low', analog=False)
return b, a
def my_lowpass_filter(data, cutoff, fs, order = 32):
#scipy filtfilt
b, a = my_butter_filter(cutoff, fs, order = order)
return filtfilt(b, a, data)
#FS = 1 / dsp.get_sampling_step(current),... tj. 1/T
def lowpasss(t, s, cutoff, vb=False, order = 3):
f_sample = (1./(t[1]-t[0]))
sigma = f_sample/(2.*np.pi*cutoff)
return flt.gaussian_filter1d(s, sigma)
def my_bandstop_filter(data, lowcut, highcut, fs, order = 32):
nyq = 0.5 * fs
low = lowcut / nyq
high = highcut /nyq
b, a = butter(order, [low, high], btype = 'bandstop')
return filtfilt(b, a, data)
def bin_average(current, voltage, binsize, vb=False):
min_voltage = np.nanmin(voltage)
max_voltage = np.nanmax(voltage)
if binsize == 0.: return {'I': current, 'V': voltage} #binsize zero will returned the input values
voltage_bins = np.arange(min_voltage + 0.5*binsize, max_voltage - 0.5*binsize, binsize)
current_bins_mean = np.empty_like(voltage_bins)
voltage_bins_mean = np.empty_like(voltage_bins)
current_bins_std = np.empty_like(voltage_bins)
voltage_bins_std = np.empty_like(voltage_bins)
for i, voltage_bin in enumerate(voltage_bins):
bin_condition = (voltage_bin + 0.5 * binsize > voltage) & (voltage > voltage_bin - 0.5 * binsize)
current_bins_mean[i] = np.mean(current[bin_condition])
voltage_bins_mean[i] = np.mean(voltage[bin_condition])
current_bins_std[i] = np.std(current[bin_condition])
voltage_bins_std[i] = np.std(voltage[bin_condition])
current_bins_std[current_bins_std == 0.] = np.nanmean(current_bins_std) #zero error can happen but is unphysical.
voltage_bins_std[voltage_bins_std == 0.] = np.nanmean(voltage_bins_std)
return {'I': current_bins_mean, 'V': voltage_bins_mean, 'I_err': current_bins_std, 'V_err': voltage_bins_std}
class Interpolator(object):
def __init__(self, x, y, kind = 'linear'):
"""
Construct 1 dimensional interpolator for function of kind f(x)=y
:param x: x coordinates, array
:param y: y coordinates, array
:param kind: method of interpolation, default linear
"""
self.interpolator = scipy.interpolate.interp1d(x, y, kind = kind)
def interpolate(self, x):
"""
Interpolates new values of y for values of x for f(x)=y
:param x: x coordinates, array
:retun: intepolated values of y, array
"""
return self.interpolator(x)
def mean_of_interpolate(self, x):
"""
Interpolates new values of y for values of x for f(x)=y and calculates their mean.
:param x: x coordinates, array
:return: mean value of interpolates values of y
"""
interval = self.interpolate(x)
return np.mean( interval )
class Intervals(object):
def __init__(self, interval, x, y, kind = 'linear'):
"""
Costructs interval with 1 dimensional interpolator for function of kind f(x)=y, where interval is x.
:param interval: interval of all possible values, array
:param x: x coordinates, array
:param y: y coordinates, array
:param kind: method of interpolation, default linear
"""
self._interpolator = Interpolator(x, y, kind = kind)
self._interval = interval
def interval(self, begin, end):
"""
Returns interval of values between begin and end.
:param begin: Smallest possible value of interval.
:param end: Largest possible value of interval.
"""
return self._interval[
np.where(
np.logical_and(
begin <= self._interval, self._interval <= end
)
)
]
def interpolate_for_interval(self, begin, end):
"""
Interpolates new values of y for values of x for f(x)=y
for interval between begin and end.
:param begin: Smallest possible value of interval.
:param end: Largest possible value of interval.
:retun: intepolated values of y within interval, array
"""
return self._interpolator.interpolate(self.interval(begin, end))
def interpolate_mean_for_interval(self, begin, end):
"""
Interpolates new values of y for values of x for f(x)=y and calculates their mean
for interval between begin and end.
:param begin: Smallest possible value of interval.
:param end: Largest possible value of interval.
:return: mean value of interpolates values of y within interval
"""
return self._interpolator.mean_of_interpolate(self.interval(begin, end))
def nearest(array, value):
array = np.asarray(array)
idx = (np.abs(array - value)).argmin()
return idx
##############FOURIER TRANSFORMATION SHOT CHECK################
#1. normal method
def freq_gen(t, Fs): #FS default value change to freq of data aquisition
n = len(t)
half = int(n/2)
k = np.arange(n)
T = n/Fs
f = k/T
f_half = f[0:half]
return f_half
#2. averaged method
from scipy.signal import welch
def my_welch(signals, nperseg, fs):
#FS sampling frq, nperseg is window lenght
#(The bigger the box the lower frequency/wavelenght we can see)
frequencies, spectras = welch(signals, fs = fs, nperseg = nperseg, scaling = 'spectrum')
return frequencies, spectras
def running_mean(l, N):
#he running mean is a case of the mathematical operation of convolution.
#For the running mean, you slide a window along the input and compute the mean of the window's contents.
# Also works for the(strictly invalid) cases when N is even.
if (N//2)*2 == N:
N = N - 1
front = np.zeros(N//2)
back = np.zeros(N//2)
for i in range(1, (N//2)*2, 2):
front[i//2] = np.convolve(l[:i], np.ones((i,))/i, mode = 'valid')
for i in range(1, (N//2)*2, 2):
back[i//2] = np.convolve(l[-i:], np.ones((i,))/i, mode = 'valid')
return np.concatenate([front, np.convolve(l, np.ones((N,))/N, mode = 'valid'), back[::-1]])
def ivchar_fit(v, Ti, R, Vf, Isat):
return np.exp(alpha)* Isat * (1 + R*(v-Vf)) - Isat * np.exp((Vf - v) / Ti)
#makes a 3-param fit with fixed Vp (eats potential gives 3p fit)###########
def Afterburner(Vf):
def ivchar_fit(v, Ti, R, Isat):
return np.exp(alpha)* Isat * (1 + R*(v-Vf)) - Isat * np.exp((Vf - v) / Ti)
return ivchar_fit
# for applications (eats Voltage and params gives I)
def ivchar_fit_Weighted(v, Ti, R, Isat):
return np.exp(alpha)* Isat * (1 + R*( v-rrpopt2[2])) - Isat * np.exp((rrpopt2[2]-v ) / Ti)
def ivchar_fit_Unweighted(v, Ti, R, Isat):
return np.exp(alpha)* Isat * (1 + R*(v-rrpopt[2] )) - Isat * np.exp(( rrpopt[2]-v) / Ti)
####condition = (150+rrpopt2[2]) < 3*crrpopt2[0] ####
def Afterburner_1(Vf):
def ivchar_fit_2p(v, Ti, Isat):
return np.exp(alpha)* Isat * (1) - Isat * np.exp(( Vf-v) / Ti)
return ivchar_fit_2p
def Afterburner_2(Vf):
def ivchar_fit_3p(v, Ti, R, Isat):
return np.exp(alpha)* Isat * (1 + R*(v-Vf )) - Isat * np.exp((Vf-v) / Ti)
return ivchar_fit_3p
def Rzero(condition):
if condition:
# print(alpha)
return Afterburner_1
else:
# print('4-param')
# print(alpha)
return Afterburner_2
def q_curve_fit(boundss, *args, **kwargs):
_len = len(boundss[0])
if _len == 3:
res_11, res_22 = curve_fit(*args, **kwargs)
return np.asarray([res_11[0], res_11[1], res_11[2]]) , np.asarray([np.sqrt(res_22[0][0]), np.sqrt(res_22[1][1]) ,np.sqrt(res_22[2][2])])
elif _len == 2:
res_1, res_2 = curve_fit(*args, **kwargs)
#return np.asarray([res_1[0], 0, res_1[1]]), res_2
return np.asarray([res_1[0], 0, res_1[1]]) , np.asarray([np.sqrt(res_2[0][0]), 0 ,np.sqrt(res_2[1][1] )])
def signal_cleaner(problem):
for i in range(len(problem)):
if np.isnan(problem[i]) == True:
problem[i]=0.0
print(i)
elif problem[i] == -inf:
problem[i] = 0.0
print('inf = ' +str(i))
elif problem[i] == inf:
problem[i] = 0.0
print('inf = ' +str(i))
return problem
Populating the interactive namespace from numpy and matplotlib
shot_number = 0
radial_probe_position = pd.read_csv(urlopen('http://golem.fjfi.cvut.cz/shots/'+str(shot_number)+'/Diagnostics/PetiProbe/Parameters/r_lp_tip'),names = ['R'])['R'][0]
sweep_frequency = pd.read_csv(urlopen('http://golem.fjfi.cvut.cz/shots/'+str(shot_number)+'/Diagnostics/PetiProbe/Parameters/f_fg')
,names = ['f_fg'])['f_fg'][0]/1e3 # [kHz]
print('sweep_frequency =' +str(sweep_frequency))
resistor = pd.read_csv(urlopen('http://golem.fjfi.cvut.cz/shots/'+str(shot_number)+'/Diagnostics/PetiProbe/Parameters/r_i'),names = ['r_i'])['r_i'][0] # [Ohm]
print('resistor =' +str(resistor))
data_file = pd.read_csv(urlopen('http://golem.fjfi.cvut.cz/shots/'+str(shot_number)+'/Diagnostics/PetiProbe/DAS_raw_data_dir/TektrMSO64_ALL.csv'), skiprows=10)
#BPP
current_BPP = signal_cleaner(data_file['CH2'])
bpp_time = 1e3*(data_file['TIME'])
current = current_BPP/resistor*1e3 ##odpor bol 47 *1000 to [mA]
voltage_BPP = (data_file['CH1'])
voltage_time = 1e3*data_file['TIME']
voltage = -1*signal_cleaner(voltage_BPP)
time_ax = 1e3*data_file['TIME']
##BPP float
# voltage_BPP_float = pd.read_csv(urlopen('http://golem.fjfi.cvut.cz/shots/' + shot_number_float +'/Diagnostics/PetiProbe/DAS_raw_data_dir/ch5.csv'),names = ['t','V'])
# BPP_time_float = 1e3*voltage_BPP_float.t
# Vp_BPP_float = 1*voltage_BPP_float.V
# OffsetBPP_float = np.mean(Vp_BPP_float[10:nearest(time_ax, tt)])
# Vp_BPP = Vp_BPP_float - OffsetBPP_float
# print('ofset_BPP_float: '+str(OffsetBPP_float))
f_sample = (1./(bpp_time[1]-bpp_time[0]))
print('DAS freqency = '+ str(f_sample) + ' MHz') ### = 1MHz
FS = f_sample
# zakladni parametry plazmatu
Bt = pd.read_csv(urlopen('http://golem.fjfi.cvut.cz/shots/' + str(shot_number) +'/Diagnostics/BasicDiagnostics/U_IntBtCoil.csv'),names = ['t','B'])
Bt.t = 1000*Bt.t
Ip = pd.read_csv(urlopen('http://golem.fjfi.cvut.cz/shots/' + str(shot_number) +'/Diagnostics/BasicDiagnostics/U_IntRogCoil.csv'),names = ['t','I'])
Ip.t = 1000*Ip.t
Ip.I = (Ip.I)
Uloop = pd.read_csv(urlopen('http://golem.fjfi.cvut.cz/shots/' + str(shot_number) +'/Diagnostics/BasicDiagnostics/U_Loop.csv'),names = ['t','V'])
Uloop.t = 1000*Uloop.t
sweep_frequency =50.0 resistor =47 DAS freqency = 12500.000000004846 MHz
fig,ax = plt.subplots(3)
fromm = Bt.t[nearest(Bt.B, 0.30)] ### usefull signal start est.
untill = Bt.t[nearest(Bt.B, 0.40)]
ax[0].plot(Uloop.t,Uloop.V)
ax[0].set_xticks([])
ax[0].set_ylabel('Uloop [V]')
ax[0].axvline(x = fromm)
ax[0].axvline(x = untill)
ax[1].plot(Ip.t,Ip.I,color = 'red')
ax[1].set_xticks([])
ax[1].set_ylabel('Ip [kA]')
ax[1].axvline(x = fromm)
ax[1].axvline(x = untill)
st = nearest(Bt.t , fromm)
ed = nearest(Bt.t , untill)
ax[2].plot(Bt.t,Bt.B, color = 'green', label = '$B_t$ interval aver. ~ '+str('{0:4.2f}'.format((Bt.B[st] + Bt.B[ed])/2))+' T')
ax[2].legend()
ax[2].set_xlabel('t [ms]')
ax[2].set_ylabel('B[T]')
ax[2].axvline(x = fromm)
ax[2].axvline(x = untill)
ax[2].grid()
B_tor_avg = (Bt.B[st] + Bt.B[ed])/2
#alpha_calc_avg = -2.735 * B_tor_avg + 2.041 ### calib_koef from lin fit of results alpha_bpp on B [P.Macha]
alpha = 0.25
print('alpha used = ' + str(alpha))
if Bt.B[st]< 0.22:
print('WARNING !!! ' +'B min calculated = ' +str(Bt.B[st]) +' WARNING !!!' )
else:
print('B min calculated = ' + str(Bt.B[st]))
st = nearest(time_ax , fromm)
ed = nearest(time_ax , untill)
B_interpol = scipy.interpolate.interp1d(Bt.t, Bt.B)
B_interpolated = B_interpol(time_ax[st:ed])
ax[2].plot(time_ax[st:ed] , B_interpolated)
plt.savefig('Results/plasma_params'+str(shot_number)+'.png')
alpha used = 0.25 B min calculated = 0.299995035
fig,ax = plt.subplots(2 , figsize = (8,8) ,sharex= False)
ax[0].plot(bpp_time, current , color= 'steelblue', label = 'BPP Current')
# ax[0].plot(bpp_time[:-1], reconstructed_shifted, 'b', label='reconstructed')
ax[0].set_ylabel('I [mA]')
#ax[0].set_ylim(-0.02,0.02)
ax[0].legend(loc= 'lower right' ,fontsize = 12)
ax[1].plot(bpp_time, voltage, color= 'lightcoral', label = 'BPP sweept voltage')
ax[1].set_xlabel('t [ms]')
ax[1].set_ylabel('U [V]')
# ax[0].set_xlim(0.75,0.77)
# ax[1].set_xlim(0.75,0.77)
# ax[1].set_ylim(-100,100)
ax[1].legend(loc= 'lower right' ,fontsize = 12)
# ax[0].axhline(x = 0)
# ax[1].axhline(x = 0)
#ax[1].plot(bpp_time, lp_I , color= 'Darkkhaki', label = 'lp Current original')
# ax[2].plot(time_ax, Vfl_LP , color= 'green', label = 'Lp_floating')
# ax[2].set_ylabel('U [V]')
# #ax[2].set_ylim(-6,6)
# ax[2].legend(loc= 'lower right' ,fontsize = 12)
plt.savefig('Results/raw_probes'+str(shot_number)+'.png')
# ------------------------------------------------------------------------------------------------------------
# ---------------------------------- Change this -------------------------------------------------------------
shift0 = 0 #shift the array to find the relative phase
shift1 = 0
tcap = 1 # end of the unperturbed voltage signal
lowpass= sweep_frequency*10 # up to 10x frequency estimate
fig, ax =plt.subplots(3,1,figsize=(8, 8))
plt.grid()
artificial = np.diff(my_lowpass_filter(voltage[nearest(bpp_time, 0.1):nearest(bpp_time, tcap)], lowpass, FS, order = 2))
real = my_lowpass_filter(current[nearest(bpp_time, 0.1):nearest(bpp_time, tcap)][:-1], lowpass, FS, order = 2)
ax[0].plot(np.roll(a=artificial, shift=shift0), real, linestyle='None', Marker='*', color='g', label='reference')
ax[0].plot(np.roll(a=artificial, shift=shift1), real, linestyle='None', Marker='*', color='b', label='shifted')
popt, pcov = curve_fit(linear_fit, artificial, real)
ax[0].plot(artificial, linear_fit(artificial, popt[0], popt[1]), linewidth=5, color='orange', label = 'linear fit')
#plt.ylim(0., 0.008)
ax[0].set_xlabel(r'Smooth $\frac{{dV}}{{dt}}$', fontsize=12)
ax[0].set_ylabel('Real current', fontsize=12)
ax[0].legend(fontsize=12)
x = np.diff(my_lowpass_filter(voltage, lowpass, FS, order = 5))
###############SHIFTING ##############################################
reconstr_shift_right = list(linear_fit(x, popt[0], popt[1]))
reconstr_shift_right.insert(0, 0.0)
reconstr_shift_right.insert(0, 0.0)
reconstructed_shifted = reconstr_shift_right[:-2]
# reconstr_shift_left = list(linear_fit(x, popt[0], popt[1]))
# reconstr_shift_left.append(0.0)
# reconstr_shift_left.append(0.0)
# reconstr_shift_left.append(0.0)
# reconstructed_shifted = reconstr_shift_left[3:]
reconstructed_noShift = linear_fit(x, popt[0], popt[1])
###############SHIFTING ##############################################
current_corrected_1 = current[:-1]-reconstructed_shifted
ax[1].grid()
ax[1].plot(bpp_time, current, 'k', label='Current raw')
ax[1].plot(bpp_time[:-1], current_corrected_1, 'red', label='Current corrected 1')
ax[1].plot(bpp_time[:-1], reconstructed_shifted, 'b', label='reconstructed 1')
###############SHIFTING RESULT COMPARATION##############################################
# ax[1].plot(bpp_time[:-1], current[:-1]-reconstructed_noShift, 'red', label='Current corrected')
# ax[1].plot(bpp_time[:-1], current[:-1]-reconstructed_shifted, 'green', label='Current corrected SHIFTED')
########################################################################################
ax[1].legend()
ax[1].set_xlabel('time [ms]', fontsize=12)
ax[1].set_ylabel('Current [mA]', fontsize=12)
ax[2].grid()
ax[2].plot(bpp_time, current, 'k',linewidth = 2, label='Current raw')
ax[2].plot(bpp_time[:-1], current_corrected_1, 'red', label='Current corrected 1')
ax[2].plot(bpp_time[:-1], reconstructed_shifted, 'b', label='reconstructed 1')
ax[2].set_xlim(0.1 , 0.3)
ax[2].set_ylim(-0.0015*1e3, 0.0015*1e3)
ax[2].legend()
ax[2].set_xlabel('time [ms]', fontsize=12)
ax[2].set_ylabel('Current [mA]', fontsize=12)
plt.tight_layout()
plt.savefig('Results/cleaning1A_'+str(shot_number)+'.png')
# ------------------------------------------------------------------------------------------------------------
# SECONDARY REMOVAL OF LEFTOVER STRAY CURRENT
# ------------------------------------------------------------------------------------------------------------
N_points = ((1/sweep_frequency)/4) // (1/(FS)) #### input in kHz!!!. This stray current is shifted by pi/2 , calculating the num. of points which represent pi/2 shift
print('second shift is '+ str(N_points))
shift0 = 0 #shift the array to find the relative phase
shift1 = 0
fig, ax =plt.subplots(3,1,figsize=(8, 8))
plt.grid()
artificial = (my_lowpass_filter(voltage[nearest(bpp_time, 0.1):nearest(bpp_time, tcap)], lowpass, FS, order = 5))
real = my_lowpass_filter(current_corrected_1[nearest(bpp_time, 0.1):nearest(bpp_time, tcap)], lowpass, FS, order = 5)
artificial1 = artificial[(np.where(artificial>0 ))]
artificial2= artificial1[(np.where(artificial1<100 ))]
real1 = real[(np.where(artificial>0 ))]
real2= real1[(np.where(artificial1<100 ))]
ax[0].plot(np.roll(a=artificial, shift=shift0), real, linestyle='None', Marker='*', color='g', label='reference')
ax[0].plot(np.roll(a=artificial, shift=shift1), real, linestyle='None', Marker='*', color='b', label='shifted')
popt, pcov = curve_fit(linear_fit, np.roll(a=artificial2, shift=shift1), real2)
ax[0].plot(artificial2, linear_fit(artificial2, popt[0], popt[1]), linewidth=5, color='orange', label = 'linear fit')
ax[0].set_xlabel(r'Smooth Voltage [V]', fontsize=12)
ax[0].set_ylabel('Real current', fontsize=12)
ax[0].legend(fontsize=12)
x = np.diff(my_lowpass_filter(voltage, lowpass, FS, order = 5))
volt_derived = my_lowpass_filter(voltage, lowpass, FS, order = 5)
reconstructed_shifted_2 = np.roll(a=volt_derived, shift=shift1)* popt[0] + popt[1]
###############SHIFTING ##############################################
current_corrected_2 = current_corrected_1 -reconstructed_shifted_2[:-1]
ax[1].grid()
ax[1].plot(bpp_time[:-1], current_corrected_1, 'k', label='Current corrected 1')
ax[1].plot(bpp_time[:-1], current_corrected_2, 'red', label='Current corrected 2')
ax[1].plot(bpp_time[:-1], reconstructed_shifted_2[:-1], 'b', label='reconstructed 2 ')
###############SHIFTING RESULT COMPARATION##############################################
# ax[1].plot(bpp_time[:-1], current[:-1]-reconstructed_noShift, 'red', label='Current corrected')
# ax[1].plot(bpp_time[:-1], current[:-1]-reconstructed_shifted, 'green', label='Current corrected SHIFTED')
########################################################################################
ax[1].legend()
ax[1].set_xlabel('time [ms]', fontsize=12)
ax[1].set_ylabel('Current [mA]', fontsize=12)
ax[1].set_ylim(-20, 0.01*1e3)
ax[2].plot(bpp_time[:-1], current_corrected_1, 'k', label='Current corrected 1')
ax[2].plot(bpp_time[:-1], current_corrected_2, 'red', label='Current corrected 2')
ax[2].plot(bpp_time[:-1], reconstructed_shifted_2[:-1], 'b', label='reconstructed 2 ')
ax[2].grid()
ax[2].set_xlim(0.1 , 0.3)
ax[2].set_ylim(-1.5, 1.5)
ax[2].legend()
ax[2].set_xlabel('time [ms]', fontsize=12)
ax[2].set_ylabel('Current [mA]', fontsize=12)
plt.tight_layout()
plt.savefig('Results/cleaning1B_'+str(shot_number)+'.png')
<ipython-input-1-debd3c57a73c>:18: MatplotlibDeprecationWarning: Case-insensitive properties were deprecated in 3.3 and support will be removed two minor releases later ax[0].plot(np.roll(a=artificial, shift=shift0), real, linestyle='None', Marker='*', color='g', label='reference') <ipython-input-1-debd3c57a73c>:19: MatplotlibDeprecationWarning: Case-insensitive properties were deprecated in 3.3 and support will be removed two minor releases later ax[0].plot(np.roll(a=artificial, shift=shift1), real, linestyle='None', Marker='*', color='b', label='shifted')
second shift is 62.0
<ipython-input-1-debd3c57a73c>:93: MatplotlibDeprecationWarning: Case-insensitive properties were deprecated in 3.3 and support will be removed two minor releases later ax[0].plot(np.roll(a=artificial, shift=shift0), real, linestyle='None', Marker='*', color='g', label='reference') <ipython-input-1-debd3c57a73c>:94: MatplotlibDeprecationWarning: Case-insensitive properties were deprecated in 3.3 and support will be removed two minor releases later ax[0].plot(np.roll(a=artificial, shift=shift1), real, linestyle='None', Marker='*', color='b', label='shifted')
lowpass= sweep_frequency*7.4
smooth_voltage = my_lowpass_filter(voltage, lowpass, FS, order = 10)
V = smooth_voltage[:-1]
Id = current_corrected_2
trt = bpp_time[:-1]
t = np.array(trt)
fourier_from = nearest(trt,fromm)
fourier_until = nearest(trt,untill )
I = my_lowpass_filter(Id, lowpass, FS, order = 10 ) # order32 like default for compass(only current), For 1kHz -> 200 50kHz->500
#PLOTING:
fig, ax = plt.subplots(3,1, sharex= False, figsize =(10,10))
ax[0].plot(bpp_time, current, label=r'raw', color='orange', alpha=.5)
ax[0].plot(t, current_corrected_2, label='I after stray current removal', color='gray', alpha=.7)
# ax[0].plot(t, Id, label='I before lowpass', color='black', alpha=.7)
ax[0].plot(t, I, label='I - offset, after lowpass', color='blue')
ax[0].plot(t[fourier_from:fourier_until], Id[fourier_from:fourier_until], label=r'Fourier analysed', color='red', alpha=1)
ax[0].legend(fontsize = 12)
ax[0].set_xlabel('t [ms]', fontsize=14)
ax[0].set_ylabel('I [mA]', fontsize=14)
ax[0].axhline(y=0)
ax[1].plot(bpp_time, current, label=r'raw', color='orange', alpha=.5)
ax[1].plot(t, current_corrected_2, label='I after stray current removal', color='gray', alpha=.7)
ax[1].plot(t, I, label='I - offset, after lowpass', color='blue')
ax[1].plot(t[fourier_from:fourier_until], Id[fourier_from:fourier_until], label=r'Fourier analysed', color='red', alpha=1)
ax[1].legend(fontsize = 12)
ax[1].set_xlabel('t [ms]', fontsize=14)
ax[1].set_ylabel('I [mA]', fontsize=14)
ax[1].set_xlim(0.1 , 4)
ax[1].set_ylim(-0.002*1e3, 0.002*1e3)
# ax[0].set_ylim(-0.04*1e3, 0.01*1e3)
fffff, sssss = my_welch(Id[fourier_from:fourier_until], nperseg=3500, fs = FS)
fff, sss = my_welch(I[fourier_from:fourier_until], nperseg=3500, fs = FS)
ax[2].loglog(fffff, sssss, color= 'gray', alpha = 0.5, label= 'before lowpass')
ax[2].loglog(fff, sss, label= 'after lowpass' )
# ax[1].set_ylim(10**(-13), 10**(-6))
ax[2].set_ylim(10**(-12), 10**(1))
# ax[1].set_xlim(0, 5*10**(6))
ax[2].axvline(lowpass, lineStyle= 'dashed' , color = "gray", label = 'cut-off freq = ' +str(lowpass)+' kHz')
ax[2].set_xlabel('frequency [kHz]', fontsize=14)
ax[2].set_ylabel('amplitude', fontsize=14)
plt.legend(fontsize=12)
plt.tight_layout()
plt.savefig('Results/cleaning2_'+str(shot_number)+'.png')
amplitude_of_cap_curr_upbound = np.std(I[nearest(bpp_time, 1.5):nearest(bpp_time, 2.5)])*3
amplitude_of_cap_curr_lobound = np.std(I[nearest(bpp_time, 1.5):nearest(bpp_time, 2.5)])*3
amplitude_of_cap_curr=np.std(I[nearest(bpp_time, 1.5):nearest(bpp_time, 2.5)])
print('Capacitive current max amplitude after all cleaning is '+str(amplitude_of_cap_curr_upbound))
<ipython-input-1-f897446c4c3a>:42: MatplotlibDeprecationWarning: Case-insensitive properties were deprecated in 3.3 and support will be removed two minor releases later ax[2].axvline(lowpass, lineStyle= 'dashed' , color = "gray", label = 'cut-off freq = ' +str(lowpass)+' kHz')
Capacitive current max amplitude after all cleaning is 0.2807192608600938
VAchar_voltage_smoothed=running_mean(V,5)
VAchar_current_smoothed=running_mean(I,5)
data=VAchar_voltage_smoothed
extrems_tmp = (np.diff(np.sign(np.diff(data))).nonzero()[0] + 1) # local min+max
minima = ((np.diff(np.sign(np.diff(data))) > 0).nonzero()[0] + 1) # local min
maxima = ((np.diff(np.sign(np.diff(data))) < 0).nonzero()[0] + 1) # local max
# Sometimes double extrems appear, removal:
#print extrems
extrems=[]
for i in range(len(extrems_tmp)-1):
if abs(extrems_tmp[i]-extrems_tmp[i+1])>10:
extrems.append(extrems_tmp[i])
# graphical output...
from pylab import *
plt.figure(figsize=(10,6), dpi= 80, facecolor='w', edgecolor='k')
plt.plot(t,data , label = 'smoothed voltage')
plt.plot(t,V, '-+', label = 'voltage')
plt.plot(t[maxima], data[maxima], "o", label="max",markersize=12)
plt.plot(t[minima], data[minima], "o", label="min",markersize=12)
plt.ylabel('Voltage [V]', fontsize = 12)
plt.xlabel('Time [ms]',fontsize = 12)
#plt.title('Maxima and minima identification')
plt.legend(fontsize = 10)
plt.axvline(x=fromm)
plt.axvline(x=untill)
# plt.ylim( -10.,110)
interval_pivots =np.concatenate([minima, maxima])
interval_pivots.sort()
interval_pivots = t[interval_pivots]
#plt.xlim(interval_pivots[192],interval_pivots[196])
sweep_amplitude = np.mean(data[maxima])-5
print('sweep_amplitude = '+str(sweep_amplitude))
sweep_amplitude = 131.51762322269064
iv_start = 1158
# iv_start = Index_filter[2]
index_to_time = interval_pivots[iv_start]
# print('For current B_tor = '+str('{0:4.2f}'.format(B_interpolated[nearest(time_ax[st:ed] , index_to_time)]))+ ' T ' + ' ln(alpha_BPP) is '+str('{0:4.2f}'.format(alpha)))
time_bin = 15
Ti_lobound= 0
Ti_upboud = 50
Vp_lobound = -20 ### The Vp can also be neg on golem
Vp_upbound = 40
R_lobound = 0
R_upbound = 1
Isat_lobound = -10## mA
Isat_upbound = 0
MAX_voltage = sweep_amplitude
potential_shift = 0.1 ### shifting the obtained Vp before cutting to the left by a factor: potetial_shift*rrpopt[2]
lobound = [Ti_lobound, R_lobound, Vp_lobound, Isat_lobound] # fit lower bound for Ti, R, Vf, Isat
upbound = [Ti_upboud, R_upbound, Vp_upbound, Isat_upbound] # fit upper bound for Ti, R, Vf, Isat
binsize = 5
# ------------------------------------------------------------------------------------------------------------
iv_stop = iv_start + 1 # +1 perioda # LOL
cut_err = (interval_pivots[iv_start-time_bin//2] < t) & (t < interval_pivots[iv_start+time_bin//2])
cut_dat = (interval_pivots[iv_start] < t) & (t < interval_pivots[iv_stop])
I_for_float_potential_index = nearest(I[cut_dat],0.0)
V_for_float_potential = V[cut_dat][I_for_float_potential_index]
cut_V = (-20 < V) & (V < MAX_voltage)
cut_I = True
cut_I = (I < 0.0)
I_short = I[cut_dat & cut_I & cut_V]
V_short = V[cut_dat & cut_I & cut_V]
I_long = I[cut_err & cut_I & cut_V]
V_long = V[cut_err & cut_I & cut_V]
if -np.nanmin(I_short) < abs(amplitude_of_cap_curr_lobound)*5:
print('WARNING: EFFECTIVE CURRENT LESS THAN 5 TIMES THE CAPACITIVE')
binned_data = bin_average(I_long, V_long, binsize) # Compute the standard deviations for each bin
verrs = binned_data['V_err'] # voltage stds
ierrs = binned_data['I_err'] # current stds
vbins = binned_data['V'] # un-center the bins
ibins = binned_data['I']
current_err = np.ones_like(I_short) # create an array of ones with the same size as I
voltage_err = np.ones_like(V_short) # which is also the size of V, so both are the same size
for i, vbin in enumerate(vbins): # iterate on the bins
current_err[V_short >= vbin] = ierrs[i]
voltage_err[V_short >= vbin] = verrs[i]
# This is the unweighted fit
popt, pcov = curve_fit(ivchar_fit, V_short, I_short, bounds=(lobound, upbound))
# This is the weighted fit
popt2, pcov2 = curve_fit(ivchar_fit, V_short, I_short, bounds=(lobound, upbound), sigma=current_err, absolute_sigma=False)
################# repeat the fit for better Vfl evaluation ##########################################################
# This is the unweighted fit
rpopt, rpcov = curve_fit(ivchar_fit, V_short, I_short, bounds=([Ti_lobound, R_lobound, popt[2]-10, Isat_lobound], [Ti_upboud, R_upbound, popt[2]+10, Isat_upbound]))
# This is the weighted fit
rpopt2, rpcov2 = curve_fit(ivchar_fit, V_short, I_short, bounds=([Ti_lobound, R_lobound, popt2[2]-10, Isat_lobound], [Ti_upboud, R_upbound, popt2[2]+10, Isat_upbound]), sigma=current_err, absolute_sigma=False)
################# 2x repeat the fit for better Vfl evaluation ######################################################
# This is the unweighted fit
rrpopt, rrpcov = curve_fit(ivchar_fit, V_short, I_short, bounds=([Ti_lobound, R_lobound, rpopt[2]-5, Isat_lobound], [Ti_upboud, R_upbound, rpopt[2]+5, Isat_upbound]))
# This is the weighted fit
rrpopt2, rrpcov2 = curve_fit(ivchar_fit, V_short, I_short, bounds=([Ti_lobound, R_lobound, rpopt2[2]-5, Isat_lobound], [Ti_upboud, R_upbound, rpopt2[2]+5, Isat_upbound]), sigma=current_err, absolute_sigma=False)
rcut_V = ((rrpopt[2] - potential_shift*rrpopt[2]) < V) & (V < MAX_voltage)
rI_short = I[cut_dat & cut_I & rcut_V]
rV_short = V[cut_dat & cut_I & rcut_V]
rI_long = I[cut_err & cut_I & rcut_V]
rV_long = V[cut_err & cut_I & rcut_V]
binned_data = bin_average(rI_long, rV_long, binsize) # Compute the standard deviations for each bin
rverrs = binned_data['V_err'] # voltage stds
rierrs = binned_data['I_err'] # current stds
rvbins = binned_data['V']
ribins = binned_data['I']
r_current_err = np.ones_like(rI_short) # create an array of ones with the same size as I
r_voltage_err = np.ones_like(rV_short) # which is also the size of V, so both are the same size
for i, vbin in enumerate(rvbins): # iterate on the bins
r_current_err[rV_short >= vbin] = rierrs[i]
r_voltage_err[rV_short >= vbin] = rverrs[i]
# This is the unweighted fit
crrpopt, crrpcov = curve_fit(Afterburner(rrpopt[2]), rV_short, rI_short, bounds=([Ti_lobound, R_lobound, Isat_lobound], [Ti_upboud, R_upbound, Isat_upbound]))
#This is the weighted fit
crrpopt2, crrpcov2 = curve_fit(Afterburner(rrpopt2[2]), rV_short, rI_short, bounds=([Ti_lobound, R_lobound, Isat_lobound], [Ti_upboud, R_upbound, Isat_upbound]), sigma=r_current_err, absolute_sigma=False)
conditions = [(MAX_voltage-rrpopt[2]) < 3*crrpopt[0] , np.sqrt(crrpcov[0][0])/crrpopt[0] < 0.6]
print(conditions)
condition = all(conditions)
if condition :
boundss = ([Ti_lobound, Isat_lobound], [Ti_upboud, Isat_upbound])
else:
boundss = ([Ti_lobound, R_lobound, Isat_lobound], [Ti_upboud, R_upbound, Isat_upbound])
#()() means call chaining... FUN()(xyz) means calling the second function returned by the first one
Ccrrpopt, Ccrrpcov = q_curve_fit(boundss, Rzero(condition)(rrpopt[2]), rV_short, rI_short, bounds=boundss)
Ccrrpopt2, Ccrrpcov2 = q_curve_fit(boundss, Rzero(condition)(rrpopt2[2]), rV_short, rI_short, bounds=boundss, sigma= r_current_err, absolute_sigma=False)
fig, ax = plt.subplots(figsize = (8,6))
ax.plot([],[],ls ='None', label ='#' +str(shot_number)+', time: '+str('{0:4.3f}'.format(index_to_time))+' s; '+ ' B_tor = '+'{0:4.2f}'.format(B_interpolated[nearest(time_ax[st:ed] , index_to_time)])+ ' T; ' +' alpha_BPP = '+str('{0:4.2f}'.format(alpha)))
ax.plot(V_short, -I_short, color='k', Marker = 'o',alpha = 0.3, linestyle='None')
ax.plot(rV_short, -rI_short, color='k', label='data', Marker = 'o', linestyle='None')
# ax.plot(vbins , -ibins , '*',color ='gray', label = 'bins through ' +str(time_bin)+' IVs')
#ax.errorbar(vbins , -ibins ,ierrs)
xx = np.linspace(rrpopt[2]- 0.1*rrpopt[2], np.nanmax(V_short), num=100) ## cutted IV
#ax.errorbar(V_short, -I_short, yerr=current_err, color='grey', label='sigma', Marker = '', linestyle='None')
##### FIRST Fits ###########
x = np.linspace(np.nanmin(V_short) , np.nanmax(rV_short), num=100) ## noncutted IV
ax.plot(x, -ivchar_fit(x, *rrpopt), color='gray', linewidth=4, alpha=.53 ,label='UNweighted_FF Ti = {0:4.2f} ± {1:4.2f} eV, Vfl = {2:4.2f} V, Isat = {3:4.2f} A'.format(rrpopt[0], np.sqrt(rrpcov[0][0]),rrpopt[2], -rrpopt[3]))
# ax.plot(x, -ivchar_fit(x, *rrpopt2), color='green', linewidth=4, alpha=.3, label='Weighted_FF Ti = {0:4.2f} ± {1:4.2f} eV, Vfl = {2:4.2f} V, Isat = {3:4.2f} A'.format(rrpopt2[0], np.sqrt(rrpcov2[0][0]),rrpopt2[2], -rrpopt2[3]))
###################### Results of normal fits ########
# ax.plot(xx, -ivchar_fit_Weighted(xx, *Ccrrpopt2), color='red', linewidth=4, alpha=.5,
# label='Weighted Ti Rzero= {0:4.2f} ± {1:4.2f} eV, Vfl = {2:4.2f} ± {3:4.2f} V, Isat = {4:4.2f} A'.format(Ccrrpopt2[0], Ccrrpcov2[0], rrpopt2[2], np.sqrt(rrpcov2[2][2]) ,-Ccrrpopt2[2]))
ax.plot(xx, -ivchar_fit_Unweighted(xx, *Ccrrpopt), color='blue', linewidth=4, alpha=.5,
label='fit: Ti= {0:4.1f} ± {1:4.1f} eV, V_pl = {2:4.1f} ± {3:4.1f} V, Isat = {4:4.3f} A'.format(Ccrrpopt[0], Ccrrpcov[0], rrpopt[2],np.sqrt(rrpcov[2][2]) , -Ccrrpopt[2]))
# ax.set_xlim(0,)
# ax.set_ylim(0,)
ax.axvline(x=0, color = 'black')
ax.axvline(x=V_for_float_potential, color = 'gray', ls = '--', label = 'BPP floating potential')
ax.axhline(y=0, color = 'black')
ax.axhline(y=-amplitude_of_cap_curr_lobound, label='stray current max. amplitude')
ax.axhline(y=-amplitude_of_cap_curr_upbound)
ax.set_xlabel('Voltage [V]', fontsize=12)
ax.set_ylabel('Current [mA]', fontsize=12)
ax.legend(loc='lower right', fontsize=10)
ax.legend(loc='upper left', fontsize=10)
plt.savefig('Results/Example_'+str(index_to_time)+'_shot'+str(shot_number)+'.png')
WARNING: EFFECTIVE CURRENT LESS THAN 5 TIMES THE CAPACITIVE [False, True]
<ipython-input-1-d607b0275828>:123: MatplotlibDeprecationWarning: Case-insensitive properties were deprecated in 3.3 and support will be removed two minor releases later ax.plot(V_short, -I_short, color='k', Marker = 'o',alpha = 0.3, linestyle='None') <ipython-input-1-d607b0275828>:124: MatplotlibDeprecationWarning: Case-insensitive properties were deprecated in 3.3 and support will be removed two minor releases later ax.plot(rV_short, -rI_short, color='k', label='data', Marker = 'o', linestyle='None')
time_bin = 15
Ti_lobound= 0
Ti_upboud = 50
Vp_lobound = -20 ### The Vp can also be neg on golem
Vp_upbound = 40
R_lobound = 0
R_upbound = 1
Isat_lobound = -10## mA
Isat_upbound = 0
MAX_voltage = sweep_amplitude
potential_shift = 0.1 ### shifting the obtained Vp before cutting to the left by a factor: potetial_shift*rrpopt[2]
min_float_potential = -50
max_float_potential = 40
time_start = fromm
time_end = untill
iv_start = nearest(interval_pivots,time_start)
iv_terminate = nearest(interval_pivots,time_end)
# iv_start = 666
# iv_terminate = 667
j = iv_start
cutoff_resolution = 5
resolution = 1 # radial temperature resolution. It is a multiplication constant! #############################1
###############################################################################################################
Ti = []
Vfl = []
R = []
I_sat = []
Ti_err = []
R_err = []
Vfl_err = []
I_sat_err= []
index = []
recognition = []
iv_time = []
V_BPP_floating_pot = []
B_IV = []
alpha_for_Te = []
for j in range(iv_start,iv_terminate, 1):
iv_start = j
lobound = [Ti_lobound, R_lobound, Vp_lobound, Isat_lobound] # fit lower bound for Ti, R, Vf, Isat
upbound = [Ti_upboud, R_upbound, Vp_upbound, Isat_upbound] # fit upper bound for Ti, R, Vf, Isat
# ------------------------------------------------------------------------------------------------------------
iv_stop = iv_start + resolution #+ 1 perioda # LOL
cut_dat = (interval_pivots[iv_start] < t) & (t < interval_pivots[iv_stop])
####### estimation of BPP floating potential##############
cut_V_for_float = (min_float_potential < V) & (V < max_float_potential)
I_for_float_potential_index = nearest(I[cut_dat & cut_V_for_float],0.0)
V_for_float_potential = V[cut_dat & cut_V_for_float][I_for_float_potential_index]
V_BPP_floating_pot.append(V_for_float_potential)
#############################################################
index_to_time = interval_pivots[iv_start]
alpha_for_Te.append(1.89 * B_interpolated[nearest(time_ax[st:ed] , index_to_time)] + 1.85)
B_IV.append(B_interpolated[nearest(time_ax[st:ed] , index_to_time)])
# cut_V = (V_for_float_potential < V) & (V < MAX_voltage)
# cut_V = (-30 < V) & (V < MAX_voltage)
cut_V = (-20 < V) & (V < MAX_voltage)
cut_I = (I < 0.0)
I_short = I[cut_dat & cut_I & cut_V]
V_short = V[cut_dat & cut_I & cut_V]
V_check = V[cut_dat & cut_V ]
if ( len(V_check) < 1 ):
Ti.append(-1000)
Ti_err.append(-1000)
Vfl.append(-1000)
R.append(-1000)
I_sat.append(-1000)
R_err.append(-1000)
Vfl_err.append(-1000)
I_sat_err.append(-1000)
index.append(iv_start)
recognition.append('unfitable V_check')
iv_time.append(index_to_time)
# alpha_bpp.append(alpha)
continue
if (np.max(V_check) < MAX_voltage-10):
print('THIS IS NOT SUPPOSED TO HAPPEN CHECK VOLTAGE SWEEP AND ALLOWED RANGE')
Ti.append(-1000)
Ti_err.append(-1000)
Vfl.append(-1000)
R.append(-1000)
I_sat.append(-1000)
R_err.append(-1000)
Vfl_err.append(-1000)
I_sat_err.append(-1000)
index.append(iv_start)
recognition.append('unfitable V_check 2 ')
iv_time.append(index_to_time)
continue
try:
rrpopt, rrpcov = curve_fit(ivchar_fit, V_short, I_short, bounds=(lobound, upbound))
except (RuntimeError,ValueError): # for the graph is dark and full of errors
Ti.append(-1000)
Ti_err.append(-1000)
Vfl.append(-1000)
R.append(-1000)
I_sat.append(-1000)
R_err.append(-1000)
Vfl_err.append(-1000)
I_sat_err.append(-1000)
index.append(iv_start)
recognition.append('unfitable FF')
iv_time.append(index_to_time)
# alpha_bpp.append(alpha)
continue
rcut_V = ((rrpopt[2] - potential_shift*rrpopt[2]) < V) & (V < MAX_voltage)
rI_short = I[cut_dat & cut_I & rcut_V]
rV_short = V[cut_dat & cut_I & rcut_V]
try:
cpopt, cpcov = curve_fit(Afterburner(rrpopt[2]), rV_short, rI_short, bounds=([Ti_lobound, R_lobound, Isat_lobound], [Ti_upboud, R_upbound, Isat_upbound]))
except (RuntimeError,ValueError): # for the graph is dark and full of errors
Ti.append(-1000)
Ti_err.append(-1000)
Vfl.append(-1000)
R.append(-1000)
I_sat.append(-1000)
R_err.append(-1000)
Vfl_err.append(-1000)
I_sat_err.append(-1000)
index.append(iv_start)
iv_time.append(index_to_time)
# alpha_bpp.append(alpha)
recognition.append('unfitable cpopt')
continue #### must be continue otherwise two values can be written instead of one
try:
condition = [(MAX_voltage-rrpopt[2]) < 2*cpopt[0] ,
abs(np.sqrt(cpcov[0][0])/cpopt[0]) < 0.6]
condition = all(condition)
if condition :
boundss = ([Ti_lobound, Isat_lobound], [Ti_upboud, Isat_upbound])
else:
boundss = ([Ti_lobound, R_lobound, Isat_lobound], [Ti_upboud, R_upbound, Isat_upbound])
Ccrrpopt, Ccrrpcov = q_curve_fit(boundss, Rzero(condition)(rrpopt[2]), rV_short, rI_short, bounds=boundss)
except (RuntimeError,ValueError): # for the graph is dark and full of errors
Ti.append(-1000)
Ti_err.append(-1000)
Vfl.append(-1000)
R.append(-1000)
I_sat.append(-1000)
R_err.append(-1000)
Vfl_err.append(-1000)
I_sat_err.append(-1000)
index.append(iv_start)
recognition.append('unfitable errcode 1')
iv_time.append(index_to_time)
pass
rel_err = 0.6
cond1 = [abs(rrpopt[2]) < 5, np.sqrt(rrpcov[2][2]) < 3] #### From ISTTOK we have fluct.level of Te is 0.25. Thus 0.25*Te = 3 (\Te cca 13/);
if all(cond1):
# print(1)
typical_filtering_cond = [np.divide(Ccrrpcov[0],Ccrrpopt[0]) < rel_err ,
(MAX_voltage-rrpopt[2])/Ccrrpopt[0] >1,
Ccrrpopt[0] > 0.1]
else:
# print(2)
typical_filtering_cond = [np.divide(Ccrrpcov[0],Ccrrpopt[0]) < rel_err ,
abs(np.divide(np.sqrt(rrpcov[2][2]),rrpopt[2])) < rel_err,
(MAX_voltage-rrpopt[2])/Ccrrpopt[0] >1,
Ccrrpopt[0] > 0.1]
if all(typical_filtering_cond):
Ti.append(Ccrrpopt[0])
Ti_err.append(Ccrrpcov[0])
R.append(Ccrrpopt[1])
Vfl.append(rrpopt[2])
I_sat.append(Ccrrpopt[2])
R_err.append(Ccrrpcov[1])
Vfl_err.append(np.sqrt(rrpcov[2][2]))
I_sat_err.append(Ccrrpcov[2])
index.append(iv_start)
recognition.append('standard fit')
iv_time.append(index_to_time)
else:
Ti.append(-1000)
Ti_err.append(-1000)
Vfl.append(-1000)
R.append(-1000)
I_sat.append(-1000)
R_err.append(-1000)
Vfl_err.append(-1000)
I_sat_err.append(-1000)
index.append(iv_start)
recognition.append('unfitable')
iv_time.append(index_to_time)
# alpha_bpp.append(alpha)
rel_err= 0.6
Ti_filt = []
Ti_err_filt = []
Vp_probe_filt = []
Vp_err_filt = []
time_filter = []
Isat_filter = []
Index_filter = []
V_BPP_floating_pot_filter = []
statement = []
B_IV_filt =[]
# alpha_bpp_filt =[]
i = 0
for i in range(0,len(Ti)):
rules = [np.divide(Ti_err[i],Ti[i]) < rel_err ,
#abs(np.divide(Vfl_err[i],Vfl[i])) < rel_err, ### we do not need this. Vp err handeled already
Ti[i]<(MAX_voltage - np.array(Vfl[i]))/1,
Ti[i]<Ti_upboud-1,
Ti_err[i]>0.0,
]
if all(rules):
Ti_filt.append(Ti[i])
Ti_err_filt.append(Ti_err[i])
Vp_probe_filt.append(Vfl[i])
Vp_err_filt.append(Vfl_err[i])
time_filter.append(iv_time[i])
Isat_filter.append(I_sat[i])
Index_filter.append(index[i])
statement.append(recognition[i])
# alpha_bpp_filt.append(alpha_bpp[i])
V_BPP_floating_pot_filter.append(V_BPP_floating_pot[i])
B_IV_filt.append(B_IV[i])
print('nuber of succesful fits = '+ str(len(Ti_filt))+' / ' +str(len(Ti)))
nuber of succesful fits = 102 / 302
match_standard =[]
match_cutoff =[]
for i in range(len(statement)):
if 'standard fit' in statement[i]:
match_standard.append(i)
if 'cutoff_fit; precision: '+str(cutoff_resolution) in statement[i]:
match_cutoff.append(i)
plt.figure(figsize = (10,6))
plt.title('Histeresis checker')
for i in range(len(Index_filter)):
if Index_filter[i] % 2 == 0:
plt.plot(Index_filter[i], Vp_probe_filt[i], color='red', marker ='o', label = 'párne')
else:
plt.plot(Index_filter[i], Vp_probe_filt[i], color='orange', marker ='o',label='nepárne')
# plt.plot(Index_filter, V_BPP_floating_pot_filter, 'D-', label = 'BPP floating potential')
plt.plot(Index_filter, Vp_probe_filt, color = 'gray',ls = '--')
plt.xlabel('fit indentification number')
plt.ylabel('Plasma potential')
plt.savefig('Results/Plasma_potential.png')
fig,ax = plt.subplots(4, figsize = (8,10))
ax[0].plot(Uloop.t,Uloop.V,linewidth = 3 ,)
ax[0].set_ylabel('Uloop [V]')
ax[0].set_xlim(fromm,untill)
ax[0].grid()
ax[1].plot(Ip.t,Ip.I,linewidth = 3 ,color = 'red')
ax[1].set_ylabel('Ip [kA]')
ax[1].set_xlim(fromm,untill)
ax[1].grid()
ax[2].plot(Bt.t,Bt.B,linewidth = 3 , color = 'green')
ax[2].set_ylabel('B [T]')
ax[2].set_xlim(fromm,untill)
ax[2].grid()
ax[3].plot(time_filter, Ti_filt, '*',markersize = 10 , color = 'red' , label ='ion temperature ' + str(shot_number))
ax[3].errorbar(time_filter, Ti_filt, yerr=Ti_err_filt, xerr=None, fmt='None', ecolor='gray', elinewidth=None, capsize=None, barsabove=True, lolims=False, uplims=False, xlolims=False, xuplims=False, errorevery=1, capthick=None, )
ax[3].set_ylim(0,80)
ax[3].set_xlim(fromm,untill)
ax[3].legend(loc= 'upper right', fontsize = 10)
ax[3].set_xlabel('time [ms]', fontsize = 12)
ax[3].set_ylabel('$T_i$ [eV]', fontsize = 12)
#ax[3].set_yticks([0,5,10,15,20,30,40,50])
ax[3].grid()
fig.savefig('Results/ALL_params_' + str(radial_probe_position ) +'mm_'+str(shot_number)+'.png' )
fig, ax = subplots(figsize= (8,6))
# ax.plot(times, Tes, 'o',markersize = 1.5,label = 'electron temperature #'+ str(shot_number_float))
# ax.plot(times, Te_smooth,color = 'darkorange', label = 'electron temperature smoothed #'+ str(shot_number_float))
ax.plot(time_filter, Ti_filt, 'o',markersize = 8 , color = 'red' , label ='ion temperature #' + str(shot_number))
plt.errorbar(time_filter, Ti_filt, yerr=Ti_err_filt, xerr=None, fmt='None', ecolor='red', elinewidth=0.5, capsize=3, barsabove=True, lolims=False, uplims=False, xlolims=False, xuplims=False, errorevery=1, capthick=None, )
#ax.set_xlim(-0.005,0.02)
y_ticks = np.arange(0, 110, 10)
ax.set_yticks(y_ticks)
ax.set_ylim(0,80)
ax.legend(fontsize = 10, loc = 'upper right')
ax.set_xlabel('time [ms]', fontsize = 12)
ax.set_ylabel('Ion and electron temperature [eV]', fontsize = 12)
ax.grid()
plt.savefig('icon-fig.png' )
fig.savefig('Results/Ti_' + str(radial_probe_position ) +'mm_'+str(shot_number)+'.png' )
shot_number_float = shot_number
lp = pd.read_csv(urlopen('http://golem.fjfi.cvut.cz/shots/' + str(shot_number) +'/Diagnostics/PetiProbe/U_LP_fl.csv'),names = ['t','V'])
Vfl_LP = signal_cleaner(lp.V)
Vfl_LP_t = signal_cleaner(lp.t)*1e3
lo = nearest(Vfl_LP_t, fromm)-100
hi = nearest(Vfl_LP_t, untill)+100
Vfl_LP_smooth = my_lowpass_filter(Vfl_LP[lo:hi], 1 ,fs = 1/(lp.t[1]*1e3-lp.t[0]*1e3), order = 2)
# Vfl_LP_smooth = Vfl_LP[lo:hi]
xxxx = scipy.interpolate.interp1d(Vfl_LP_t[lo:hi],Vfl_LP_smooth)
V_LP_floating_pot = xxxx(iv_time)
V_BPP_floating_pot_smooth = running_mean(V_BPP_floating_pot, 100 )
# Tes =(V_BPP_floating_pot_smooth-V_LP_floating_pot)/2.5
Tes = []
Te_raw =[]
for i in range(len(V_BPP_floating_pot_smooth)):
Tes.append((V_BPP_floating_pot_smooth[i]-V_LP_floating_pot[i])/alpha_for_Te[i])
Te_raw.append((V_BPP_floating_pot[i]-V_LP_floating_pot[i])/alpha_for_Te[i])
Te_smooth = lowpasss(iv_time, Tes, 100)
V_p =[]
V_p_raw =[]
for i in range(len(Te_smooth)):
V_p.append(V_BPP_floating_pot_smooth[i] + alpha* Te_smooth[i])
V_p_raw.append(V_BPP_floating_pot[i] + alpha* Te_raw[i])
cut_lob = 0.
cut_hib = 0.
times = iv_time[nearest(iv_time, fromm+cut_lob) :nearest(iv_time, untill-cut_hib)]
Te_smooth = Te_smooth[nearest(iv_time, fromm+cut_lob) :nearest(iv_time, untill-cut_hib)]
V_p = V_p[nearest(iv_time, fromm+cut_lob) :nearest(iv_time, untill-cut_hib)]
Tes = Tes[nearest(iv_time, fromm+cut_lob) :nearest(iv_time, untill-cut_hib)]
inf = 10441 inf = 10442 inf = 10443 inf = 10444 inf = 10487 inf = 10488 inf = 10489 inf = 10490 inf = 10491 inf = 10493 inf = 10494 inf = 10495 inf = 10496 inf = 10497 inf = 10533 inf = 10534 inf = 10535 inf = 10536 inf = 10539 inf = 10542 inf = 10544 inf = 10545 inf = 10546 inf = 10549 inf = 10552 inf = 10579 inf = 10582 inf = 10584 inf = 10585 inf = 10586 inf = 10590 inf = 10591 inf = 10592 inf = 10593 inf = 10594 inf = 10595 inf = 10596 inf = 10598 inf = 10599 inf = 10601 inf = 10602 inf = 10603 inf = 10604 inf = 10627 inf = 10628 inf = 10630 inf = 10631 inf = 10632 inf = 10634 inf = 10636 inf = 10639 inf = 10645 inf = 10647 inf = 10648 inf = 10649 inf = 10650 inf = 10652 inf = 10654 inf = 10655 inf = 10656 inf = 10658 inf = 10676 inf = 10677 inf = 10678 inf = 10679 inf = 10680 inf = 10682 inf = 10683 inf = 10684 inf = 10685 10686 inf = 10687 inf = 10688 inf = 10689 inf = 10690 inf = 10691 inf = 10692 inf = 10693 inf = 10695 inf = 10699 inf = 10702 inf = 10704 inf = 10707 inf = 10708 inf = 10726 inf = 10727 inf = 10729 inf = 10730 inf = 10732 inf = 10733 inf = 10734 inf = 10735 inf = 10736 inf = 10738 inf = 10739 inf = 10740 inf = 10743 inf = 10744 inf = 10745 inf = 10747 inf = 10748 inf = 10749 inf = 10751 inf = 10752 inf = 10755 inf = 10756 inf = 10758 10775 inf = 10776 inf = 10777 inf = 10778 inf = 10779 inf = 10780 inf = 10782 inf = 10785 inf = 10786 inf = 10789 inf = 10790 inf = 10792 inf = 10793 inf = 10796 inf = 10797 inf = 10798 inf = 10799 inf = 10801 inf = 10803 inf = 10807 inf = 10824 inf = 10827 inf = 10828 inf = 10830 inf = 10831 inf = 10832 inf = 10833 inf = 10835 inf = 10837 inf = 10838 inf = 10839 inf = 10840 inf = 10844 inf = 10845 inf = 10848 inf = 10849 inf = 10850 inf = 10851 inf = 10852 inf = 10853 inf = 10855 inf = 10856 inf = 10858 inf = 10877 inf = 10879 inf = 10880 inf = 10881 inf = 10883 inf = 10887 inf = 10888 inf = 10889 inf = 10891 inf = 10892 inf = 10895 inf = 10896 inf = 10898 10899 inf = 10900 inf = 10901 inf = 10902 inf = 10927 inf = 10928 inf = 10931 10932 inf = 10934 inf = 10935 inf = 10936 inf = 10937 inf = 10941 inf = 10942 inf = 10943 inf = 10945 inf = 10946 10947 inf = 10950 inf = 10951 10952 inf = 10953 inf = 10978 inf = 10982 inf = 10983 inf = 10984 inf = 10987 inf = 10989 inf = 10990 inf = 10993 inf = 10994 10995 inf = 10996 inf = 10997 inf = 10998 inf = 10999 inf = 11000 inf = 11001 inf = 11002 inf = 11030 11035 inf = 11037 inf = 11038 11039 11040 inf = 11041 inf = 11044 inf = 11045 inf = 11046 inf = 11047 inf = 11082 inf = 11083 inf = 11086 inf = 11088 inf = 11090 11091 inf = 11093 inf = 11095 inf = 11096 inf = 11097 11098 inf = 11099 inf = 11134 inf = 11135 inf = 11136 inf = 11137 inf = 11138 11140 inf = 11141 inf = 11142 inf = 11143 inf = 11144 inf = 11145 inf = 11146 inf = 11187 inf = 11188 inf = 11190 inf = 11191 inf = 11192 inf = 11238 inf = 11242 11288 inf = 11289 inf = 11290 inf = 11292 inf = 11293 inf = 23666 inf = 23667
plt.figure(figsize = (8,6))
# plt.plot(iv_time,V_BPP_floating_pot, 'blue', marker ='o',alpha = .4, label = 'Vfl BPP raw')
# plt.plot(iv_time,V_BPP_floating_pot_smooth, 'blue', marker ='o', label = 'Vfl BPP ')
plt.plot(iv_time,V_p_raw, marker ='o',markersize = 4, color='green',alpha = .3, label = '$\Phi$ from BPP & LP')
plt.plot(time_filter,Vp_probe_filt, 'red', marker ='o', markersize = 7, ls ='--',alpha = .5, label = '$\Phi$ 4-p fit ')
plt.plot(times,V_p, linewidth = 3, color='lime', label = '$\Phi$ smooth from BPP & LP')
plt.plot(time_filter,running_mean(Vp_probe_filt,50), color = 'red',linewidth = 3,label = '$\Phi$ smooth from fits ')
plt.ylim(-60,60)
plt.xlabel('time [ms]', fontsize = 12)
plt.ylabel('Potential [V]', fontsize = 12)
plt.legend()
plt.savefig('Results/Potential_check' + str(radial_probe_position ) +'mm_'+str(shot_number)+'.png' )
BPP should not fluctuate more than lp. probably the noise from coax --> histeresis
plt.figure()
plt.plot(Vfl_LP_t[lo:hi], Vfl_LP[lo:hi], linewidth = 3, color = 'blue',alpha = .4, label = '$V_{fl}^{LP}$ [V]')
plt.plot(iv_time,V_BPP_floating_pot,'-o', linewidth = 1.5, color = 'red',alpha = .4, label = '$\Phi_{fl}^{BPP}$ [V]')
plt.plot(iv_time,V_LP_floating_pot, linewidth = 3, color = 'blue', label = 'smooth $V_{fl}^{LP}$ [V]')
plt.plot(iv_time,V_BPP_floating_pot_smooth,'-', linewidth = 3, color = 'orange', label = 'smooth $\Phi_{fl}^{BPP}$ [V]')
plt.xlabel('time [ms]', fontsize = 12)
plt.ylabel('Potential [V]', fontsize = 12)
plt.legend()
plt.savefig('Results/Vfl' + str(radial_probe_position ) +'mm_'+str(shot_number)+'.png' )
plt.figure()
plt.plot(times ,Tes, '-',label = 'Electron temperature [eV]' )
plt.plot(times,Te_smooth, linewidth = 3, color = 'orange', label = 'Electron temperature smooth [eV]')
plt.ylim(0,40)
plt.legend()
<matplotlib.legend.Legend at 0x7f07ebb1f9a0>
fig, ax = subplots(figsize= (8,6))
ax.plot([],[],ls='None', label = 'OH-regime; $R = $'+str(radial_probe_position) + ' mm')
# ax.plot(times, Tes, 'o',markersize = 1.5,label = 'electron temperature #'+ str(shot_number_float))
ax.plot(times, Te_smooth,color = 'blue', linewidth =3, label = '$T_e$ smoothed #'+ str(shot_number_float))
ax.plot(time_filter, Ti_filt, 'o',markersize = 5 , color = 'red' , label ='$T_i$ 4-p fit #' + str(shot_number))
plt.errorbar(time_filter, Ti_filt, yerr=Ti_err_filt, xerr=None, fmt='None', ecolor='red', elinewidth=0.5, capsize=3, barsabove=True, lolims=False, uplims=False, xlolims=False, xuplims=False, errorevery=1, capthick=None, )
#ax.set_xlim(-0.005,0.02)
ax.set_ylim(0,70)
ax.legend(fontsize = 10, loc = 'upper right')
ax.set_xlabel('time [ms]', fontsize = 12)
ax.set_ylabel('Ion and electron temperature [eV]', fontsize = 12)
plt.savefig('icon-fig.png' )
plt.savefig('Results/Ti_' + str(radial_probe_position ) +'mm_'+str(shot_number)+'.png' )
rel_p_for = np.ones(len(iv_time))
rel_p_forTi = rel_p_for-1 + radial_probe_position
My_Temperatures_full = Table([iv_time, Ti,Ti_err, Vfl,Vfl_err, R,R_err, I_sat,I_sat_err,rel_p_forTi, B_IV], names=[' t[ms] ', 'Ti[eV]', 'Ti_err[eV]', 'Vp_cut[V]','Vp_err[V]','R','R_err','I_sat','I_sat_err','p_position[m]','B_IV'])
ascii.write(My_Temperatures_full, 'Results/Ti_BPP_profile_Fast_full_'+str(shot_number)+'.txt')
#$$$$$$$$$ ULOZENIE FILTERED DAT $$$$$$$###############
p_fil = np.ones(len(time_filter))
p_filt = p_fil-1 + radial_probe_position
from astropy.table import QTable, Table, Column
My_Temperatures = Table([time_filter, Ti_filt,Ti_err_filt, Vp_probe_filt,Vp_err_filt, Isat_filter ,p_filt, B_IV_filt], names=[' t[ms] ', 'Ti[eV]', 'Ti_err[eV]', 'Vp_cut[V]','Vp_err[V]','Isat','Ti_norm_positions', 'B_IV_filt'])
ascii.write(My_Temperatures, 'Results/Ti_BPP_profile_filtered'+str(shot_number)+'.txt')
#$$$$$$$$$ ULOZENIE Te DAT $$$$$$$###############
Te_p_fil = np.ones(len(times))
Te_p_filt = Te_p_fil-1 + radial_probe_position
from astropy.table import QTable, Table, Column
Te_table = Table([times, Tes,Te_smooth,V_p, Te_p_filt , B_IV[:-1],V_LP_floating_pot[:-1] , V_BPP_floating_pot[:-1]], names=['t[ms] ', 'Te[eV]', 'Te_smooth','Plasma_potential', 'pos', 'B_IV', 'V_LP_floating_pot' , 'V_BPP_floating_pot'])
ascii.write(Te_table, 'Results/Te_'+str(shot_number)+'.txt')