%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.300013691
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.20299338052956145
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 = 132.04059437551686
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')
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-1-d607b0275828> in <module> 39 I_long = I[cut_err & cut_I & cut_V] 40 V_long = V[cut_err & cut_I & cut_V] ---> 41 if -np.nanmin(I_short) < abs(amplitude_of_cap_curr_lobound)*5: 42 print('WARNING: EFFECTIVE CURRENT LESS THAN 5 TIMES THE CAPACITIVE') 43 <__array_function__ internals> in nanmin(*args, **kwargs) /opt/anaconda3/lib/python3.8/site-packages/numpy/lib/nanfunctions.py in nanmin(a, axis, out, keepdims) 317 # Fast, but not safe for subclasses of ndarray, or object arrays, 318 # which do not implement isnan (gh-9009), or fmin correctly (gh-8975) --> 319 res = np.fmin.reduce(a, axis=axis, out=out, **kwargs) 320 if np.isnan(res).any(): 321 warnings.warn("All-NaN slice encountered", RuntimeWarning, ValueError: zero-size array to reduction operation fmin which has no identity
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)
/opt/anaconda3/lib/python3.8/site-packages/scipy/optimize/minpack.py:828: OptimizeWarning: Covariance of the parameters could not be estimated warnings.warn('Covariance of the parameters could not be estimated',
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 = 0 / 300
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 = 10806 inf = 10855 inf = 10857 inf = 10858 inf = 10903 inf = 10904 inf = 10907 10911 inf = 10912 inf = 10913 inf = 10915 inf = 10916 inf = 10917 inf = 10918 inf = 10920 inf = 10922 inf = 10925 inf = 10926 inf = 10928 inf = 10929 inf = 10930 inf = 10931 inf = 10932 inf = 10933 inf = 10934 inf = 10935 10938 inf = 10939 inf = 10940 inf = 10942 inf = 10945 inf = 10947 inf = 10948 inf = 10949 inf = 10951 inf = 10953 inf = 10954 inf = 10955 inf = 10956 inf = 10958 inf = 10959 inf = 10960 inf = 10962 inf = 10963 inf = 10964 inf = 10965 inf = 10967 inf = 10969 inf = 10970 inf = 10971 inf = 10973 inf = 10974 10975 inf = 10976 inf = 10978 inf = 10979 inf = 10981 inf = 10983 inf = 10984 inf = 10988 inf = 10989 inf = 10990 inf = 10992 inf = 10993 inf = 10996 inf = 10998 inf = 11000 inf = 11001 inf = 11002 inf = 11003 inf = 11004 inf = 11005 inf = 11006 11007 inf = 11008 inf = 11011 inf = 11013 inf = 11014 inf = 11015 inf = 11016 inf = 11020 inf = 11021 inf = 11023 inf = 11025 inf = 11027 inf = 11028 inf = 11029 inf = 11031 inf = 11033 inf = 11034 11039 inf = 11041 inf = 11043 inf = 11046 inf = 11047 inf = 11050 inf = 11051 inf = 11052 inf = 11053 inf = 11055 inf = 11056 inf = 11058 inf = 11059 inf = 11060 inf = 11062 inf = 11066 inf = 11067 inf = 11068 inf = 11069 inf = 11070 inf = 11073 inf = 11074 inf = 11075 inf = 11078 inf = 11080 inf = 11081 inf = 11088 inf = 11089 inf = 11090 inf = 11093 inf = 11094 inf = 11097 inf = 11098 inf = 11099 inf = 11101 inf = 11102 inf = 11103 inf = 11105 inf = 11110 inf = 11111 inf = 54142
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' )
--------------------------------------------------------------------------- ValueError Traceback (most recent call last) <ipython-input-1-8d177e5f52a9> in <module> 6 7 plt.plot(times,V_p, linewidth = 3, color='lime', label = '$\Phi$ smooth from BPP & LP') ----> 8 plt.plot(time_filter,running_mean(Vp_probe_filt,50), color = 'red',linewidth = 3,label = '$\Phi$ smooth from fits ') 9 plt.ylim(-60,60) 10 plt.xlabel('time [ms]', fontsize = 12) <ipython-input-1-5dbf927b7065> in running_mean(l, N) 197 back = np.zeros(N//2) 198 for i in range(1, (N//2)*2, 2): --> 199 front[i//2] = np.convolve(l[:i], np.ones((i,))/i, mode = 'valid') 200 for i in range(1, (N//2)*2, 2): 201 back[i//2] = np.convolve(l[-i:], np.ones((i,))/i, mode = 'valid') <__array_function__ internals> in convolve(*args, **kwargs) /opt/anaconda3/lib/python3.8/site-packages/numpy/core/numeric.py in convolve(a, v, mode) 813 raise ValueError('a cannot be empty') 814 if len(v) == 0: --> 815 raise ValueError('v cannot be empty') 816 mode = _mode_from_name(mode) 817 return multiarray.correlate(a, v[::-1], mode) ValueError: v cannot be empty
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 0x7f931ce53eb0>
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')