%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.299998194
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.1756341972290209
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.60205858538632
iv_start = 1158
# iv_start = Index_filter[2]
index_to_time = interval_pivots[iv_start]
# print('For current B_tor = '+str('{0:4.2f}'.format(B_interpolated[nearest(time_ax[st:ed] , index_to_time)]))+ ' T ' + ' ln(alpha_BPP) is '+str('{0:4.2f}'.format(alpha)))
time_bin = 15
Ti_lobound= 0
Ti_upboud = 50
Vp_lobound = -20 ### The Vp can also be neg on golem
Vp_upbound = 40
R_lobound = 0
R_upbound = 1
Isat_lobound = -10## mA
Isat_upbound = 0
MAX_voltage = sweep_amplitude
potential_shift = 0.1 ### shifting the obtained Vp before cutting to the left by a factor: potetial_shift*rrpopt[2]
lobound = [Ti_lobound, R_lobound, Vp_lobound, Isat_lobound] # fit lower bound for Ti, R, Vf, Isat
upbound = [Ti_upboud, R_upbound, Vp_upbound, Isat_upbound] # fit upper bound for Ti, R, Vf, Isat
binsize = 5
# ------------------------------------------------------------------------------------------------------------
iv_stop = iv_start + 1 # +1 perioda # LOL
cut_err = (interval_pivots[iv_start-time_bin//2] < t) & (t < interval_pivots[iv_start+time_bin//2])
cut_dat = (interval_pivots[iv_start] < t) & (t < interval_pivots[iv_stop])
I_for_float_potential_index = nearest(I[cut_dat],0.0)
V_for_float_potential = V[cut_dat][I_for_float_potential_index]
cut_V = (-20 < V) & (V < MAX_voltage)
cut_I = True
cut_I = (I < 0.0)
I_short = I[cut_dat & cut_I & cut_V]
V_short = V[cut_dat & cut_I & cut_V]
I_long = I[cut_err & cut_I & cut_V]
V_long = V[cut_err & cut_I & cut_V]
if -np.nanmin(I_short) < abs(amplitude_of_cap_curr_lobound)*5:
print('WARNING: EFFECTIVE CURRENT LESS THAN 5 TIMES THE CAPACITIVE')
binned_data = bin_average(I_long, V_long, binsize) # Compute the standard deviations for each bin
verrs = binned_data['V_err'] # voltage stds
ierrs = binned_data['I_err'] # current stds
vbins = binned_data['V'] # un-center the bins
ibins = binned_data['I']
current_err = np.ones_like(I_short) # create an array of ones with the same size as I
voltage_err = np.ones_like(V_short) # which is also the size of V, so both are the same size
for i, vbin in enumerate(vbins): # iterate on the bins
current_err[V_short >= vbin] = ierrs[i]
voltage_err[V_short >= vbin] = verrs[i]
# This is the unweighted fit
popt, pcov = curve_fit(ivchar_fit, V_short, I_short, bounds=(lobound, upbound))
# This is the weighted fit
popt2, pcov2 = curve_fit(ivchar_fit, V_short, I_short, bounds=(lobound, upbound), sigma=current_err, absolute_sigma=False)
################# repeat the fit for better Vfl evaluation ##########################################################
# This is the unweighted fit
rpopt, rpcov = curve_fit(ivchar_fit, V_short, I_short, bounds=([Ti_lobound, R_lobound, popt[2]-10, Isat_lobound], [Ti_upboud, R_upbound, popt[2]+10, Isat_upbound]))
# This is the weighted fit
rpopt2, rpcov2 = curve_fit(ivchar_fit, V_short, I_short, bounds=([Ti_lobound, R_lobound, popt2[2]-10, Isat_lobound], [Ti_upboud, R_upbound, popt2[2]+10, Isat_upbound]), sigma=current_err, absolute_sigma=False)
################# 2x repeat the fit for better Vfl evaluation ######################################################
# This is the unweighted fit
rrpopt, rrpcov = curve_fit(ivchar_fit, V_short, I_short, bounds=([Ti_lobound, R_lobound, rpopt[2]-5, Isat_lobound], [Ti_upboud, R_upbound, rpopt[2]+5, Isat_upbound]))
# This is the weighted fit
rrpopt2, rrpcov2 = curve_fit(ivchar_fit, V_short, I_short, bounds=([Ti_lobound, R_lobound, rpopt2[2]-5, Isat_lobound], [Ti_upboud, R_upbound, rpopt2[2]+5, Isat_upbound]), sigma=current_err, absolute_sigma=False)
rcut_V = ((rrpopt[2] - potential_shift*rrpopt[2]) < V) & (V < MAX_voltage)
rI_short = I[cut_dat & cut_I & rcut_V]
rV_short = V[cut_dat & cut_I & rcut_V]
rI_long = I[cut_err & cut_I & rcut_V]
rV_long = V[cut_err & cut_I & rcut_V]
binned_data = bin_average(rI_long, rV_long, binsize) # Compute the standard deviations for each bin
rverrs = binned_data['V_err'] # voltage stds
rierrs = binned_data['I_err'] # current stds
rvbins = binned_data['V']
ribins = binned_data['I']
r_current_err = np.ones_like(rI_short) # create an array of ones with the same size as I
r_voltage_err = np.ones_like(rV_short) # which is also the size of V, so both are the same size
for i, vbin in enumerate(rvbins): # iterate on the bins
r_current_err[rV_short >= vbin] = rierrs[i]
r_voltage_err[rV_short >= vbin] = rverrs[i]
# This is the unweighted fit
crrpopt, crrpcov = curve_fit(Afterburner(rrpopt[2]), rV_short, rI_short, bounds=([Ti_lobound, R_lobound, Isat_lobound], [Ti_upboud, R_upbound, Isat_upbound]))
#This is the weighted fit
crrpopt2, crrpcov2 = curve_fit(Afterburner(rrpopt2[2]), rV_short, rI_short, bounds=([Ti_lobound, R_lobound, Isat_lobound], [Ti_upboud, R_upbound, Isat_upbound]), sigma=r_current_err, absolute_sigma=False)
conditions = [(MAX_voltage-rrpopt[2]) < 3*crrpopt[0] , np.sqrt(crrpcov[0][0])/crrpopt[0] < 0.6]
print(conditions)
condition = all(conditions)
if condition :
boundss = ([Ti_lobound, Isat_lobound], [Ti_upboud, Isat_upbound])
else:
boundss = ([Ti_lobound, R_lobound, Isat_lobound], [Ti_upboud, R_upbound, Isat_upbound])
#()() means call chaining... FUN()(xyz) means calling the second function returned by the first one
Ccrrpopt, Ccrrpcov = q_curve_fit(boundss, Rzero(condition)(rrpopt[2]), rV_short, rI_short, bounds=boundss)
Ccrrpopt2, Ccrrpcov2 = q_curve_fit(boundss, Rzero(condition)(rrpopt2[2]), rV_short, rI_short, bounds=boundss, sigma= r_current_err, absolute_sigma=False)
fig, ax = plt.subplots(figsize = (8,6))
ax.plot([],[],ls ='None', label ='#' +str(shot_number)+', time: '+str('{0:4.3f}'.format(index_to_time))+' s; '+ ' B_tor = '+'{0:4.2f}'.format(B_interpolated[nearest(time_ax[st:ed] , index_to_time)])+ ' T; ' +' alpha_BPP = '+str('{0:4.2f}'.format(alpha)))
ax.plot(V_short, -I_short, color='k', Marker = 'o',alpha = 0.3, linestyle='None')
ax.plot(rV_short, -rI_short, color='k', label='data', Marker = 'o', linestyle='None')
# ax.plot(vbins , -ibins , '*',color ='gray', label = 'bins through ' +str(time_bin)+' IVs')
#ax.errorbar(vbins , -ibins ,ierrs)
xx = np.linspace(rrpopt[2]- 0.1*rrpopt[2], np.nanmax(V_short), num=100) ## cutted IV
#ax.errorbar(V_short, -I_short, yerr=current_err, color='grey', label='sigma', Marker = '', linestyle='None')
##### FIRST Fits ###########
x = np.linspace(np.nanmin(V_short) , np.nanmax(rV_short), num=100) ## noncutted IV
ax.plot(x, -ivchar_fit(x, *rrpopt), color='gray', linewidth=4, alpha=.53 ,label='UNweighted_FF Ti = {0:4.2f} ± {1:4.2f} eV, Vfl = {2:4.2f} V, Isat = {3:4.2f} A'.format(rrpopt[0], np.sqrt(rrpcov[0][0]),rrpopt[2], -rrpopt[3]))
# ax.plot(x, -ivchar_fit(x, *rrpopt2), color='green', linewidth=4, alpha=.3, label='Weighted_FF Ti = {0:4.2f} ± {1:4.2f} eV, Vfl = {2:4.2f} V, Isat = {3:4.2f} A'.format(rrpopt2[0], np.sqrt(rrpcov2[0][0]),rrpopt2[2], -rrpopt2[3]))
###################### Results of normal fits ########
# ax.plot(xx, -ivchar_fit_Weighted(xx, *Ccrrpopt2), color='red', linewidth=4, alpha=.5,
# label='Weighted Ti Rzero= {0:4.2f} ± {1:4.2f} eV, Vfl = {2:4.2f} ± {3:4.2f} V, Isat = {4:4.2f} A'.format(Ccrrpopt2[0], Ccrrpcov2[0], rrpopt2[2], np.sqrt(rrpcov2[2][2]) ,-Ccrrpopt2[2]))
ax.plot(xx, -ivchar_fit_Unweighted(xx, *Ccrrpopt), color='blue', linewidth=4, alpha=.5,
label='fit: Ti= {0:4.1f} ± {1:4.1f} eV, V_pl = {2:4.1f} ± {3:4.1f} V, Isat = {4:4.3f} A'.format(Ccrrpopt[0], Ccrrpcov[0], rrpopt[2],np.sqrt(rrpcov[2][2]) , -Ccrrpopt[2]))
# ax.set_xlim(0,)
# ax.set_ylim(0,)
ax.axvline(x=0, color = 'black')
ax.axvline(x=V_for_float_potential, color = 'gray', ls = '--', label = 'BPP floating potential')
ax.axhline(y=0, color = 'black')
ax.axhline(y=-amplitude_of_cap_curr_lobound, label='stray current max. amplitude')
ax.axhline(y=-amplitude_of_cap_curr_upbound)
ax.set_xlabel('Voltage [V]', fontsize=12)
ax.set_ylabel('Current [mA]', fontsize=12)
ax.legend(loc='lower right', fontsize=10)
ax.legend(loc='upper left', fontsize=10)
plt.savefig('Results/Example_'+str(index_to_time)+'_shot'+str(shot_number)+'.png')
WARNING: EFFECTIVE CURRENT LESS THAN 5 TIMES THE CAPACITIVE [False, True]
<ipython-input-1-d607b0275828>:123: MatplotlibDeprecationWarning: Case-insensitive properties were deprecated in 3.3 and support will be removed two minor releases later ax.plot(V_short, -I_short, color='k', Marker = 'o',alpha = 0.3, linestyle='None') <ipython-input-1-d607b0275828>:124: MatplotlibDeprecationWarning: Case-insensitive properties were deprecated in 3.3 and support will be removed two minor releases later ax.plot(rV_short, -rI_short, color='k', label='data', Marker = 'o', linestyle='None')
time_bin = 15
Ti_lobound= 0
Ti_upboud = 50
Vp_lobound = -20 ### The Vp can also be neg on golem
Vp_upbound = 40
R_lobound = 0
R_upbound = 1
Isat_lobound = -10## mA
Isat_upbound = 0
MAX_voltage = sweep_amplitude
potential_shift = 0.1 ### shifting the obtained Vp before cutting to the left by a factor: potetial_shift*rrpopt[2]
min_float_potential = -50
max_float_potential = 40
time_start = fromm
time_end = untill
iv_start = nearest(interval_pivots,time_start)
iv_terminate = nearest(interval_pivots,time_end)
# iv_start = 666
# iv_terminate = 667
j = iv_start
cutoff_resolution = 5
resolution = 1 # radial temperature resolution. It is a multiplication constant! #############################1
###############################################################################################################
Ti = []
Vfl = []
R = []
I_sat = []
Ti_err = []
R_err = []
Vfl_err = []
I_sat_err= []
index = []
recognition = []
iv_time = []
V_BPP_floating_pot = []
B_IV = []
alpha_for_Te = []
for j in range(iv_start,iv_terminate, 1):
iv_start = j
lobound = [Ti_lobound, R_lobound, Vp_lobound, Isat_lobound] # fit lower bound for Ti, R, Vf, Isat
upbound = [Ti_upboud, R_upbound, Vp_upbound, Isat_upbound] # fit upper bound for Ti, R, Vf, Isat
# ------------------------------------------------------------------------------------------------------------
iv_stop = iv_start + resolution #+ 1 perioda # LOL
cut_dat = (interval_pivots[iv_start] < t) & (t < interval_pivots[iv_stop])
####### estimation of BPP floating potential##############
cut_V_for_float = (min_float_potential < V) & (V < max_float_potential)
I_for_float_potential_index = nearest(I[cut_dat & cut_V_for_float],0.0)
V_for_float_potential = V[cut_dat & cut_V_for_float][I_for_float_potential_index]
V_BPP_floating_pot.append(V_for_float_potential)
#############################################################
index_to_time = interval_pivots[iv_start]
alpha_for_Te.append(1.89 * B_interpolated[nearest(time_ax[st:ed] , index_to_time)] + 1.85)
B_IV.append(B_interpolated[nearest(time_ax[st:ed] , index_to_time)])
# cut_V = (V_for_float_potential < V) & (V < MAX_voltage)
# cut_V = (-30 < V) & (V < MAX_voltage)
cut_V = (-20 < V) & (V < MAX_voltage)
cut_I = (I < 0.0)
I_short = I[cut_dat & cut_I & cut_V]
V_short = V[cut_dat & cut_I & cut_V]
V_check = V[cut_dat & cut_V ]
if ( len(V_check) < 1 ):
Ti.append(-1000)
Ti_err.append(-1000)
Vfl.append(-1000)
R.append(-1000)
I_sat.append(-1000)
R_err.append(-1000)
Vfl_err.append(-1000)
I_sat_err.append(-1000)
index.append(iv_start)
recognition.append('unfitable V_check')
iv_time.append(index_to_time)
# alpha_bpp.append(alpha)
continue
if (np.max(V_check) < MAX_voltage-10):
print('THIS IS NOT SUPPOSED TO HAPPEN CHECK VOLTAGE SWEEP AND ALLOWED RANGE')
Ti.append(-1000)
Ti_err.append(-1000)
Vfl.append(-1000)
R.append(-1000)
I_sat.append(-1000)
R_err.append(-1000)
Vfl_err.append(-1000)
I_sat_err.append(-1000)
index.append(iv_start)
recognition.append('unfitable V_check 2 ')
iv_time.append(index_to_time)
continue
try:
rrpopt, rrpcov = curve_fit(ivchar_fit, V_short, I_short, bounds=(lobound, upbound))
except (RuntimeError,ValueError): # for the graph is dark and full of errors
Ti.append(-1000)
Ti_err.append(-1000)
Vfl.append(-1000)
R.append(-1000)
I_sat.append(-1000)
R_err.append(-1000)
Vfl_err.append(-1000)
I_sat_err.append(-1000)
index.append(iv_start)
recognition.append('unfitable FF')
iv_time.append(index_to_time)
# alpha_bpp.append(alpha)
continue
rcut_V = ((rrpopt[2] - potential_shift*rrpopt[2]) < V) & (V < MAX_voltage)
rI_short = I[cut_dat & cut_I & rcut_V]
rV_short = V[cut_dat & cut_I & rcut_V]
try:
cpopt, cpcov = curve_fit(Afterburner(rrpopt[2]), rV_short, rI_short, bounds=([Ti_lobound, R_lobound, Isat_lobound], [Ti_upboud, R_upbound, Isat_upbound]))
except (RuntimeError,ValueError): # for the graph is dark and full of errors
Ti.append(-1000)
Ti_err.append(-1000)
Vfl.append(-1000)
R.append(-1000)
I_sat.append(-1000)
R_err.append(-1000)
Vfl_err.append(-1000)
I_sat_err.append(-1000)
index.append(iv_start)
iv_time.append(index_to_time)
# alpha_bpp.append(alpha)
recognition.append('unfitable cpopt')
continue #### must be continue otherwise two values can be written instead of one
try:
condition = [(MAX_voltage-rrpopt[2]) < 2*cpopt[0] ,
abs(np.sqrt(cpcov[0][0])/cpopt[0]) < 0.6]
condition = all(condition)
if condition :
boundss = ([Ti_lobound, Isat_lobound], [Ti_upboud, Isat_upbound])
else:
boundss = ([Ti_lobound, R_lobound, Isat_lobound], [Ti_upboud, R_upbound, Isat_upbound])
Ccrrpopt, Ccrrpcov = q_curve_fit(boundss, Rzero(condition)(rrpopt[2]), rV_short, rI_short, bounds=boundss)
except (RuntimeError,ValueError): # for the graph is dark and full of errors
Ti.append(-1000)
Ti_err.append(-1000)
Vfl.append(-1000)
R.append(-1000)
I_sat.append(-1000)
R_err.append(-1000)
Vfl_err.append(-1000)
I_sat_err.append(-1000)
index.append(iv_start)
recognition.append('unfitable errcode 1')
iv_time.append(index_to_time)
pass
rel_err = 0.6
cond1 = [abs(rrpopt[2]) < 5, np.sqrt(rrpcov[2][2]) < 3] #### From ISTTOK we have fluct.level of Te is 0.25. Thus 0.25*Te = 3 (\Te cca 13/);
if all(cond1):
# print(1)
typical_filtering_cond = [np.divide(Ccrrpcov[0],Ccrrpopt[0]) < rel_err ,
(MAX_voltage-rrpopt[2])/Ccrrpopt[0] >1,
Ccrrpopt[0] > 0.1]
else:
# print(2)
typical_filtering_cond = [np.divide(Ccrrpcov[0],Ccrrpopt[0]) < rel_err ,
abs(np.divide(np.sqrt(rrpcov[2][2]),rrpopt[2])) < rel_err,
(MAX_voltage-rrpopt[2])/Ccrrpopt[0] >1,
Ccrrpopt[0] > 0.1]
if all(typical_filtering_cond):
Ti.append(Ccrrpopt[0])
Ti_err.append(Ccrrpcov[0])
R.append(Ccrrpopt[1])
Vfl.append(rrpopt[2])
I_sat.append(Ccrrpopt[2])
R_err.append(Ccrrpcov[1])
Vfl_err.append(np.sqrt(rrpcov[2][2]))
I_sat_err.append(Ccrrpcov[2])
index.append(iv_start)
recognition.append('standard fit')
iv_time.append(index_to_time)
else:
Ti.append(-1000)
Ti_err.append(-1000)
Vfl.append(-1000)
R.append(-1000)
I_sat.append(-1000)
R_err.append(-1000)
Vfl_err.append(-1000)
I_sat_err.append(-1000)
index.append(iv_start)
recognition.append('unfitable')
iv_time.append(index_to_time)
# alpha_bpp.append(alpha)
rel_err= 0.6
Ti_filt = []
Ti_err_filt = []
Vp_probe_filt = []
Vp_err_filt = []
time_filter = []
Isat_filter = []
Index_filter = []
V_BPP_floating_pot_filter = []
statement = []
B_IV_filt =[]
# alpha_bpp_filt =[]
i = 0
for i in range(0,len(Ti)):
rules = [np.divide(Ti_err[i],Ti[i]) < rel_err ,
#abs(np.divide(Vfl_err[i],Vfl[i])) < rel_err, ### we do not need this. Vp err handeled already
Ti[i]<(MAX_voltage - np.array(Vfl[i]))/1,
Ti[i]<Ti_upboud-1,
Ti_err[i]>0.0,
]
if all(rules):
Ti_filt.append(Ti[i])
Ti_err_filt.append(Ti_err[i])
Vp_probe_filt.append(Vfl[i])
Vp_err_filt.append(Vfl_err[i])
time_filter.append(iv_time[i])
Isat_filter.append(I_sat[i])
Index_filter.append(index[i])
statement.append(recognition[i])
# alpha_bpp_filt.append(alpha_bpp[i])
V_BPP_floating_pot_filter.append(V_BPP_floating_pot[i])
B_IV_filt.append(B_IV[i])
print('nuber of succesful fits = '+ str(len(Ti_filt))+' / ' +str(len(Ti)))
nuber of succesful fits = 118 / 301
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 = 10605 inf = 10606 inf = 10607 inf = 10652 inf = 10653 10655 inf = 10656 inf = 10657 inf = 10658 inf = 10659 inf = 10662 inf = 10698 inf = 10699 inf = 10700 inf = 10702 10703 inf = 10704 inf = 10705 inf = 10706 inf = 10707 inf = 10710 inf = 10711 inf = 10713 inf = 10714 inf = 10715 inf = 10746 inf = 10748 inf = 10749 inf = 10750 inf = 10751 inf = 10755 inf = 10756 inf = 10758 inf = 10761 inf = 10762 inf = 10763 inf = 10766 inf = 10793 inf = 10794 10797 inf = 10798 inf = 10801 inf = 10802 inf = 10804 inf = 10809 inf = 10810 inf = 10811 inf = 10812 inf = 10813 inf = 10814 inf = 10815 inf = 10816 inf = 10818 inf = 10819 inf = 10841 inf = 10842 inf = 10844 inf = 10847 inf = 10850 inf = 10852 inf = 10853 inf = 10855 inf = 10856 inf = 10861 inf = 10863 inf = 10864 inf = 10866 inf = 10867 inf = 10868 inf = 10869 inf = 10870 inf = 10871 inf = 10889 inf = 10892 inf = 10893 inf = 10894 inf = 10896 inf = 10897 inf = 10898 inf = 10899 inf = 10901 inf = 10903 inf = 10904 inf = 10906 inf = 10911 inf = 10912 inf = 10915 inf = 10916 inf = 10917 inf = 10918 inf = 10919 inf = 10920 inf = 10922 inf = 10940 inf = 10941 inf = 10942 inf = 10943 10944 inf = 10945 inf = 10946 inf = 10947 inf = 10948 inf = 10949 inf = 10950 inf = 10952 inf = 10953 inf = 10954 inf = 10955 inf = 10958 inf = 10962 inf = 10963 inf = 10964 inf = 10965 inf = 10966 inf = 10969 inf = 10988 inf = 10989 inf = 10991 inf = 10992 inf = 10993 inf = 10994 inf = 10996 inf = 11002 inf = 11004 inf = 11005 inf = 11008 inf = 11011 inf = 11012 inf = 11013 inf = 11014 inf = 11015 inf = 11016 11017 11018 inf = 11020 inf = 11021 inf = 11044 inf = 11046 inf = 11050 inf = 11053 inf = 11054 inf = 11056 inf = 11057 inf = 11058 inf = 11059 inf = 11065 inf = 11067 inf = 11068 inf = 11069 inf = 11094 inf = 11098 11099 inf = 11104 inf = 11105 inf = 11106 inf = 11108 inf = 11109 inf = 11110 inf = 11112 inf = 11114 inf = 11145 inf = 11148 inf = 11149 inf = 11150 inf = 11151 inf = 11154 inf = 11155 inf = 11156 inf = 11157 inf = 11158 inf = 11159 inf = 11162 inf = 11200 inf = 11201 inf = 11204 inf = 11205 inf = 11206 inf = 11207 inf = 11208 inf = 11211 inf = 11248 inf = 11249 inf = 11251 inf = 11252 inf = 11253 inf = 11254 inf = 11255 11256 inf = 11257 inf = 11258 inf = 11260 inf = 11261 inf = 11302 inf = 11303 inf = 11304 inf = 11305 inf = 11307 inf = 11308 inf = 11309 inf = 11354 inf = 11356 inf = 11357 inf = 11358 inf = 11359 inf = 11404 11405 inf = 11406 11407 11451 inf = 11452 inf = 11453 inf = 11454 inf = 11455 11457 inf = 11458
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 0x7f2fd7e56100>
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')