%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.300002366
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.1645926531392217
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.18201722642695
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')
[True, False]
<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 = 109 / 303
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 = 10622 inf = 10625 inf = 10626 inf = 10670 inf = 10671 inf = 10672 inf = 10673 inf = 10674 inf = 10675 inf = 10676 inf = 10677 inf = 10679 inf = 10680 inf = 10719 inf = 10720 inf = 10722 inf = 10723 inf = 10728 inf = 10730 inf = 10731 inf = 10732 inf = 10766 inf = 10768 inf = 10771 inf = 10772 inf = 10774 inf = 10776 inf = 10778 inf = 10782 inf = 10783 inf = 10784 inf = 10813 inf = 10814 inf = 10815 inf = 10816 inf = 10817 inf = 10819 inf = 10821 inf = 10824 inf = 10827 inf = 10830 inf = 10831 inf = 10832 inf = 10833 inf = 10837 inf = 10859 inf = 10861 inf = 10864 inf = 10865 inf = 10868 inf = 10869 inf = 10870 inf = 10872 inf = 10874 inf = 10875 inf = 10876 inf = 10877 inf = 10879 inf = 10881 inf = 10882 inf = 10883 inf = 10884 inf = 10886 inf = 10887 inf = 10888 inf = 10908 inf = 10910 inf = 10911 inf = 10912 inf = 10913 inf = 10914 inf = 10916 inf = 10919 inf = 10922 inf = 10923 inf = 10925 10927 inf = 10928 inf = 10930 inf = 10931 inf = 10932 inf = 10934 inf = 10938 10939 inf = 10957 inf = 10958 inf = 10959 inf = 10960 inf = 10961 inf = 10962 inf = 10965 inf = 10966 inf = 10967 inf = 10968 inf = 10970 inf = 10971 inf = 10972 inf = 10973 inf = 10974 inf = 10978 inf = 10979 inf = 10981 inf = 10982 inf = 10985 inf = 10986 inf = 10989 inf = 11008 inf = 11009 inf = 11010 inf = 11013 inf = 11014 inf = 11015 inf = 11017 inf = 11018 inf = 11019 inf = 11020 inf = 11023 inf = 11025 inf = 11026 11027 inf = 11028 inf = 11029 inf = 11030 inf = 11031 inf = 11032 inf = 11033 inf = 11034 inf = 11035 inf = 11036 inf = 11062 11063 inf = 11065 inf = 11067 inf = 11068 inf = 11071 inf = 11073 inf = 11074 inf = 11075 inf = 11077 inf = 11078 inf = 11079 inf = 11080 inf = 11081 inf = 11083 inf = 11084 inf = 11085 inf = 11112 inf = 11113 inf = 11119 inf = 11122 inf = 11123 inf = 11125 inf = 11126 inf = 11128 inf = 11131 inf = 11133 inf = 11161 inf = 11162 inf = 11163 inf = 11165 inf = 11166 inf = 11168 inf = 11169 inf = 11171 inf = 11176 inf = 11177 inf = 11178 inf = 11179 inf = 11181 inf = 11183 inf = 11215 inf = 11216 inf = 11217 inf = 11219 inf = 11220 inf = 11222 inf = 11224 inf = 11225 inf = 11227 11229 inf = 11230 inf = 11269 inf = 11270 11271 inf = 11272 inf = 11274 inf = 11276 inf = 11278 inf = 11323 inf = 11324 11368 inf = 11371 inf = 11372 inf = 11374 inf = 11378 inf = 11420 inf = 11421 inf = 11422 inf = 11423 inf = 11424 inf = 11471 inf = 11473 inf = 11474 inf = 11520 inf = 11521 inf = 11523 inf = 52441
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 0x7feff015efa0>
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')