#%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
shot_no = 36657
shot_number = shot_no
radial_probe_position =70
sweep_frequency = 50 # [kHz]
sweep_amplitude = 120 # [V]
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)
### ICE reading pickles ####################################################
# current_BPP = pd.read_pickle("current_BPP" +str(shot_number) + ".pkl")
# voltage_BPP = pd.read_pickle("voltage_BPP" +str(shot_number) + ".pkl")
# Bt= pd.read_pickle( "Bt"+str(shot_number) + ".pkl")
# time_ax = 1e3*current_BPP.t
# Vfl_LP = = pd.read_pickle( "Vfl_LP"+str(shot_number) + ".pkl")
############################################################################
# lp = pd.read_csv(urlopen('http://golem.fjfi.cvut.cz/shots/' + shot_number +'/Diagnostics/PetiProbe/DAS_raw_data_dir/ch4.csv'),names = ['t','V'])
# lp.t = 1e3*lp.t
# time_ax = lp.t
# lp_V = lp.V *1
# ## #offsetremoval
# tt=1
# time_ax =np.array(time_ax)
# OffsetLP = np.mean(lp_V[10:nearest(time_ax, tt)])
# print('ofset_lp: '+str(OffsetLP))
# Vfl_LP = lp_V - OffsetLP
#BPP
current_BPP = data_file['CH2']
bpp_time = 1e3*data_file['TIME']
current = (current_BPP/47)*1e3 ##odpor bol 47 *1000 to [mA]
voltage_BPP = data_file['CH1']
voltage_time = 1e3*data_file['TIME']
voltage = -1*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
################# PICKLE-ing of the df's "for the database is dark and full of errors" ##################
# current_BPP.to_pickle("current_BPP" +str(shot_number) + ".pkl")
# voltage_BPP.to_pickle("voltage_BPP" +str(shot_number) + ".pkl")
# Bt.to_pickle("Bt" +str(shot_number) + ".pkl")
# Ip.to_pickle("Ip" +str(shot_number) + ".pkl")
# Uloop.to_pickle("Uloop" +str(shot_number) + ".pkl")
# Vfl_LP.to_pickle("Vfl_LP" +str(shot_number) + ".pkl")
############ ICN ICN ICN ICN ICN ICN ICN ICN ICN ICN [In Case of NaNs] ##########################
#Vfl_LP.replace([np.inf, -np.inf, np.nan, NaN], value = 0)
DAS freqency = 12500.000000004846 MHz
fig,ax = plt.subplots(3)
fromm = 8 ### usefull signal start est.
untill = 16 ##16
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)
# print('alpha average calculated for present average magetic field= ' + str(alpha_calc_avg))
plt.savefig(str(shot_number)+'plasma_params.png')
alpha used = 0.25 B min calculated = 0.29240784
fig,ax = plt.subplots(3 , 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(str(shot_number)+'raw_probes.png')
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-1-f80f9aaf94b9> in <module> 18 19 #ax[1].plot(bpp_time, lp_I , color= 'Darkkhaki', label = 'lp Current original') ---> 20 ax[2].plot(time_ax, Vfl_LP , color= 'green', label = 'Lp_floating') 21 ax[2].set_ylabel('U [V]') 22 #ax[2].set_ylim(-6,6) NameError: name 'Vfl_LP' is not defined
# ------------------------------------------------------------------------------------------------------------
# ---------------------------------- 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*7.4 # 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 = 5))
real = my_lowpass_filter(current[nearest(bpp_time, 0.1):nearest(bpp_time, tcap)][:-1], lowpass, FS, order = 5)
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(str(shot_number)+'cleaning1A.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 = int(N_points)
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 = 5))
real = my_lowpass_filter(current_corrected_1[nearest(bpp_time, 0.1):nearest(bpp_time, tcap)], lowpass, FS, order = 5)[:-1]
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=artificial, shift=shift1), real)
ax[0].plot(artificial, linear_fit(artificial, popt[0], popt[1]), linewidth=5, color='orange', label = 'linear fit')
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))
volt_derived = np.diff(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
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, '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, '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(str(shot_number)+'cleaning1B.png')
<ipython-input-1-c35bf443bacc>: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-c35bf443bacc>: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-c35bf443bacc>:87: 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-c35bf443bacc>:88: 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')
# plt.figure()
# plt.plot( bpp_time[nearest(bpp_time, 0.1):nearest(bpp_time, tcap)],my_lowpass_filter(abc[nearest(bpp_time, 0.1):nearest(bpp_time, tcap)], lowpass, FS, order = 5), label = 'real current' )
# plt.plot(bpp_time[nearest(bpp_time, 0.1):nearest(bpp_time, tcap)][:-1] , np.roll(a=artificial, shift=shift1)* popt[0] , label = 'shifted artificial')
# plt.legend()
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(str(shot_number)+'cleaning2.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-796e1c0e631f>:41: 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.2543927302942374
VAchar_voltage_smoothed=running_mean(V,5)
VAchar_current_smoothed=running_mean(I,5)
ylim_min=np.min(VAchar_current_smoothed)
ylim_max=np.max(VAchar_current_smoothed)
xlim_min=np.min(VAchar_voltage_smoothed)
xlim_max=np.max(VAchar_voltage_smoothed)
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])
#np.savetxt('Data/extrems',extrems, fmt="%s")
#np.savetxt('Data/minima',minima, fmt="%s")
#np.savetxt('Data/maxima',maxima, fmt="%s")
# 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])
# e = match_standard[11]
# print('Ti '+str(Ti_filt[e]))
# print('Ti_err '+str(Ti_err_filt[e]))
# print('statement: '+str(statement[e]))
# print('index '+str(Index_filter[e]))
# print(alpha_bpp_filt[e])
iv_start = 1053
# 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 = 150
Vp_lobound = 0
Vp_upbound = 150
R_lobound = 0
R_upbound = 1
Isat_lobound = -50
Isat_upbound = 0
MAX_voltage = sweep_amplitude
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 = (0 < 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] - 0.0*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 , '*', label = 'bins through ' +str(time_bin)+' IVs')
#ax.errorbar(vbins , -ibins ,ierrs)
xx = np.linspace(rrpopt[2]- 0.0*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(str(iv_start)+'Example_'+str(index_to_time)+'_shot'+str(shot_number)+'.png')
WARNING: EFFECTIVE CURRENT LESS THAN 5 TIMES THE CAPACITIVE [True, False]
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-1-65d770ae3126> in <module> 119 120 fig, ax = plt.subplots(figsize = (8,6)) --> 121 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))) 122 ax.plot(V_short, -I_short, color='k', Marker = 'o',alpha = 0.3, linestyle='None') 123 ax.plot(rV_short, -rI_short, color='k', label='data', Marker = 'o', linestyle='None') NameError: name 'B_interpolated' is not defined
version 5.0 \ unweighted due to non-constant plasma parameters
# e = 13
# print('Ti '+str(Ti_filt[e]))
# print('Ti_err '+str(Ti_err_filt[e]))
# print('statement: '+str(statement[e]))
# print('index '+str(Index_filter[e]))
# iv_start = 1466 1100 NICE
# # iv_start = Index_filter[e]
# # iv_start = 1466
# iv_start = 697
# index_to_time = interval_pivots[iv_start]
# time_bin = 15
# Ti_lobound= 0
# Ti_upboud = 100
# Vp_lobound = -50 ### The Vp can also be neg on golem
# Vp_upbound = 100
# R_lobound = 0
# R_upbound = 1
# Isat_lobound = -100 ## mA
# Isat_upbound = 0
# MAX_voltage = 120
# potential_shift = 0.0 ### 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 = (V_for_float_potential < 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]
# I_long = I[cut_err & cut_I & cut_V]
# V_long = V[cut_err & cut_I & cut_V]
# 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]
# #ax.errorbar(V_short, -I_short, yerr=current_err, color='grey', label='sigma', Marker = '', linestyle='None')
# ##### FIRST Fits ###########
# #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]))
# def cutoff_fitting(current, voltage,r_current_err, n_segments, plasmapotential = rrpopt[2], additional_graph = False):
# if abs((np.sqrt(rrpcov[2][2])/plasmapotential))>0.6:
# print('Plasma potential could not be estimated with precision required')
# Final_Ti = NaN
# Final_Ti_err = NaN
# else:
# fig, ax = plt.subplots(2,figsize = (10,10))
# from matplotlib import gridspec
# gs = gridspec.GridSpec(2, 1, height_ratios=[3/2, 1])
# ax[0] = plt.subplot(gs[0])
# ax[1] = plt.subplot(gs[1])
# x = np.linspace(np.nanmin(V_short) , np.nanmax(rV_short), num=100)
# ax[0].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; ' +' $α_{{BPP}}$ = '+str('{0:4.2f}'.format(alpha)))
# ax[0].plot(V_short, -I_short, color='k', Marker = 'o',alpha = 0.3,markersize = 8, linestyle='None')
# ax[0].plot(rV_short, -rI_short, color='k', label='data', Marker = 'o',markersize = 8, linestyle='None')
# if voltage[0]>voltage[len(voltage)-1]:
# voltage = list(voltage) ### reversing into ascending order & checking for asc/desc if conditions change
# current = list(current)
# r_current_err = list(r_current_err)
# voltage.reverse()
# current.reverse()
# r_current_err.reverse()
# min_voltage = np.nanmin(voltage)
# max_voltage = np.nanmax(voltage)
# step_indexes = len(voltage)//n_segments
# stops= []
# V_fits = np.zeros(n_segments, dtype=object)
# I_fits = np.zeros(n_segments, dtype=object)
# crrpopt2 = np.zeros(n_segments, dtype=object)
# crrpcov2 = np.zeros(n_segments, dtype=object)
# crrpopt = np.zeros(n_segments, dtype=object)
# crrpcov = np.zeros(n_segments, dtype=object)
# Ccrrpopt = np.zeros(n_segments, dtype=object)
# Ccrrpcov = np.zeros(n_segments, dtype=object)
# Ccrrpopt2 = np.zeros(n_segments, dtype=object)
# Ccrrpcov2 = np.zeros(n_segments, dtype=object)
# Ti_UN_for_plateau = np.zeros(n_segments, dtype=object)
# Voltage_for_plateau = np.zeros(n_segments, dtype=object)
# Ti_UN_err_for_plateau = np.zeros(n_segments, dtype=object)
# voltage_cutoff = 0
# for i in range(n_segments):
# voltage_cutoff += step_indexes
# if i == n_segments-1: ### dealing with edge to Make the fit til the end. (-1 cus range is 'missing' last value thus this cond would never apply)
# stops.append(len(voltage))
# stop = len(voltage)
# else:
# stops.append(voltage_cutoff)
# stop = voltage_cutoff
# V_fits[i] = voltage[0:stop]
# I_fits[i] = current[0:stop]
# Voltage_for_plateau[i] = voltage[stop-1]
# # WEIGHTED
# try:
# params, covs = curve_fit(Afterburner(plasmapotential), V_fits[i], I_fits[i] , bounds=([Ti_lobound, R_lobound, Isat_lobound], [Ti_upboud, R_upbound, Isat_upbound]), sigma=r_current_err[0:stop], absolute_sigma=False)
# crrpopt2[i]= params
# crrpcov2[i]= covs
# except(RuntimeError,ValueError):
# crrpopt2[i] = [NaN,NaN,NaN]
# crrpcov2[i] = [NaN,NaN,NaN]
# Ccrrpopt2[i]= [NaN,NaN,NaN]
# Ccrrpcov2[i]= [NaN,NaN,NaN]
# pass ##### ATTENTION use continue for the leading fit (has to be pass otherwise a missing index wold be created if secondary method would not converge)
# #UNWEIGHTED
# try:
# params_UN, covs_UN = curve_fit(Afterburner(plasmapotential), V_fits[i], I_fits[i], bounds=([Ti_lobound, R_lobound, Isat_lobound], [Ti_upboud, R_upbound, Isat_upbound]))
# crrpopt[i]= params_UN
# crrpcov[i]= covs_UN
# except (RuntimeError,ValueError):
# crrpopt[i] = [NaN,NaN,NaN]
# crrpcov[i] = [NaN,NaN,NaN]
# Ccrrpopt[i]= [NaN,NaN,NaN]
# Ccrrpcov[i]= [NaN,NaN,NaN]
# continue
# conditions = [(Voltage_for_plateau[i]-plasmapotential) < 3*params_UN[0] , (np.sqrt(covs_UN[0][0])/params_UN[0]) < 0.6] ####ATTENTION CHOOSE Vp and TI FROM LEADING FIT + Voltage_for_plateau[i] is max voltage for particular fit#########
# condition = all(conditions)
# # print(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])
# try:
# params2_UN, covs2_UN = q_curve_fit(boundss, Rzero(condition)(rrpopt[2]), V_fits[i], I_fits[i], bounds=boundss)
# Ccrrpopt[i]= params2_UN
# Ccrrpcov[i]= covs2_UN
# Ti_UN_for_plateau[i]= params2_UN[0] #####CHANGE THIS################################################
# Ti_UN_err_for_plateau[i]= covs2_UN[0]
# except (RuntimeError,ValueError): # change this for the leading fit
# Ccrrpopt[i] = NaN
# Ccrrpcov[i] = NaN
# pass
# try:
# params3, covs3 = q_curve_fit(boundss, Rzero(condition)(rrpopt2[2]), V_fits[i], I_fits[i], bounds=boundss, sigma= r_current_err[0:stop], absolute_sigma=False)
# Ccrrpopt2[i]= params3
# Ccrrpcov2[i]= covs3
# except (RuntimeError,ValueError): # change this for the leading fit
# Ccrrpopt2[i] = NaN
# Ccrrpcov2[i] = NaN
# pass
# colors = ['black','gray','red','orange','yellow','gold','green','darkblue','steelblue','blue', 'indigo','magenta', 'red','orange','yellow','gray','red','orange','yellow','gold','green','darkblue','steelblue','blue', 'indigo','magenta', 'red','orange','yellow','gray','red','orange','yellow','gold','green','darkblue','steelblue','blue', 'indigo','magenta', 'red','orange','yellow','black','gray','red','orange','yellow','gold','green','darkblue','steelblue','blue', 'indigo','magenta', 'red','orange','yellow','gray','red','orange','yellow','gold','green','darkblue','steelblue','blue', 'indigo','magenta', 'red','orange','yellow','gray','red','orange','yellow','gold','green','darkblue','steelblue','blue', 'indigo','magenta', 'red','orange','yellow']
# alphas =[0.9,0.9,0.8,0.7,0.6,0.5,0.4,0.4,0.4,0.4,0.3, 0.2,0.2,0.3, 0.2,0.2,0.9,0.8,0.7,0.6,0.5,0.4,0.4,0.4,0.4,0.3, 0.2,0.2,0.3, 0.2,0.2,0.9,0.8,0.7,0.6,0.5,0.4,0.4,0.4,0.4,0.3, 0.2,0.2,0.3, 0.2,0.2,0.9,0.9,0.9,0.8,0.7,0.6,0.5,0.4,0.4,0.4,0.4,0.3, 0.2,0.2,0.3, 0.2,0.2,0.9,0.8,0.7,0.6,0.5,0.4,0.4,0.4,0.4,0.3, 0.2,0.2,0.3, 0.2,0.2,0.9,0.8,0.7,0.6,0.5,0.4,0.4,0.4,0.4,0.3, 0.2,0.2,0.3, 0.2,0.2,0.9]
# xx = np.linspace(rrpopt[2]- potential_shift*rrpopt[2], np.nanmax(V_fits[i]), num=100) ## cutted IV
# ax[0].plot(xx, -ivchar_fit_Unweighted(xx, *Ccrrpopt[i]), color=colors[i], linewidth=4, alpha=alphas[i],
# label='fit: Ti= {0:4.1f} ± {1:4.1f} eV, Vfl = {2:4.1f} ± {3:4.1f} V, Isat = {4:4.3f} A, $R$ = {5:.1e}'.format(Ccrrpopt[i][0], Ccrrpcov[i][0], rrpopt[2],np.sqrt(rrpcov[2][2]) , -Ccrrpopt[i][2], Ccrrpopt[i][1]))
# ######### Final result
# difference = []
# abs_difference=[]
# for j in range(0,len(Ccrrpopt)-1,1): #####change for weighted method
# if (Ccrrpopt[j][0] < Ti_upboud-1) and ((Ccrrpcov[j][0]/Ccrrpopt[j][0]) < 0.6) and (0.9*Ccrrpopt[j][0]) < (Voltage_for_plateau[j]-plasmapotential):
# #difference.append( abs(((Ccrrpopt[j+1][0]/Ccrrpopt[j][0])-1)*100))
# deriv = (Ccrrpopt[j+1][0] - Ccrrpopt[j][0]) / (Voltage_for_plateau[j+1] - Voltage_for_plateau[j])
# if abs(deriv) < 0.5:
# difference.append(deriv)
# abs_difference.append(abs(deriv))
# else:
# difference.append(NaN)
# abs_difference.append(NaN)
# else:
# difference.append(NaN)
# abs_difference.append(NaN)
# print('difference: '+str( difference )) ### What if there is no plaetau?
# min_value = np.nanmin(abs_difference) ### finds deriv with smallest difference
# try:
# index = abs_difference.index(min_value)
# total_index = abs_difference.index(min_value) ### finds an index of a min derivative
# # print('first index = '+str(total_index))
# except(ValueError):
# total_index = len(Ccrrpopt)-2 ### if there are only nan values of derivs. we take the last couple of fits and decide which one has lower covariance
# index = len(Ccrrpopt2)-2
# pass
# if Ccrrpcov[total_index][0] < Ccrrpcov[total_index+1][0]: ####change this for weighted
# Final_Ti =Ccrrpopt[total_index]
# Final_Ti_err = Ccrrpcov[total_index]
# else: ####change this for weighted
# total_index = total_index+1
# Final_Ti =Ccrrpopt[total_index]
# Final_Ti_err = Ccrrpcov[total_index]
# print('total index = '+str(total_index))
# print('Final_Ti: '+str(Final_Ti))
# print('Final_Ti_err: '+str(Final_Ti_err))
# ax[0].set_xlim(-20,)
# ax[0].set_ylim(0,)
# ax[0].axvline(0, color='black')
# ax[0].set_xlabel('Voltage [V]', fontsize=12)
# ax[0].set_ylabel('Current [mA]', fontsize=12)
# ax[0].legend(loc='lower right', fontsize=8)
# Voltage_for_diff = []
# faulty_fits = []
# for e in range(len(Voltage_for_plateau)-1): ### differences shifted between
# Vol = (Voltage_for_plateau[e] + Voltage_for_plateau[e+1])/2
# Voltage_for_diff.append(Vol)
# for e in range(len(Voltage_for_plateau)):
# if Ccrrpcov[e][0] == inf:
# faulty_fits.append(Ccrrpopt[e][0])
# else:
# faulty_fits.append(NaN)
# if additional_graph:
# # fig, ax = plt.subplots(figsize = (8,6))
# ax[1].errorbar(Voltage_for_plateau, Ti_UN_for_plateau, yerr=Ti_UN_err_for_plateau, color='red', Marker = 'o', markersize = 8, linestyle='None', label='Ti [eV]', zorder=0)
# ax[1].plot(Voltage_for_plateau[total_index] , Ti_UN_for_plateau[total_index], 'o', color = 'green', markersize = 10)
# ax[1].plot(Voltage_for_plateau , faulty_fits, 'o', color = 'white', markersize = 10, alpha= 1)
# ax[1].plot(Voltage_for_diff , difference, 'D', color = 'magenta', markersize = 8, label= 'derivatives')
# ax[1].plot(Voltage_for_diff[index] , difference[index], 'D', color = 'green', markersize = 8, label= 'minimal derivative')
# ax[1].axhline(y=0, color='black')
# ax[1].set_xlim(0,)
# ax[1].set_ylim(-2,80)
# ax[1].legend(fontsize = 10)
# ax[1].set_xlabel('Cutoff voltage [V]', fontsize=12)
# ax[1].set_ylabel('Ti [eV]', fontsize=12)
# ax[1].grid()
# fig.tight_layout()
# fig.savefig(str(iv_start)+'Example_B_'+str(index_to_time)+'_shot'+str(shot_number)+'.png')
# return {'stops':stops,'V_fits':V_fits,'I_fits':I_fits,'voltage':voltage, 'crrpopt':crrpopt,'crrpcov':crrpcov, 'crrpopt2':crrpopt2,'crrpcov2':crrpcov2, 'Ccrrpopt':Ccrrpopt,'Ccrrpcov':Ccrrpcov, 'Ccrrpopt2':Ccrrpopt2,'Ccrrpcov2':Ccrrpcov2, 'Ti_UN_for_plateau':Ti_UN_for_plateau,'Voltage_for_plateau':Voltage_for_plateau }
# cutoff_fitting(rI_short,rV_short,r_current_err,len(rV_short)//1, additional_graph= True)['stops']
# e = match_cutoff[5]
# print('Ti '+str(Ti_filt[e]))
# print('Ti_err '+str(Ti_err_filt[e]))
# print('statement: '+str(statement[e]))
# print('index '+str(Index_filter[e]))
# # print(alpha_bpp_filt[e])
time_bin = 15
Ti_lobound= 0
Ti_upboud = 100
Vp_lobound = -50 ### The Vp can also be neg on golem
Vp_upbound = 100
R_lobound = 0
R_upbound = 1
Isat_lobound = -100 ## mA
Isat_upbound = 0
MAX_voltage = 150
potential_shift = 0.0 ### shifting the obtained Vp before cutting to the left by a factor: potetial_shift*rrpopt[2]
min_float_potential = -100
max_float_potential = 50
# version 5.0
def cutoff_fitting_for_profile_UNweighted(current, voltage, n_segments, additional_graph = False):
plasmapotential = rrpopt[2]
if voltage[0]>voltage[len(voltage)-1]:
voltage = list(voltage) ### reversing into ascending order & checking for asc/desc if conditions change
current = list(current)
voltage.reverse()
current.reverse()
min_voltage = np.nanmin(voltage)
max_voltage = np.nanmax(voltage)
step_indexes = len(voltage)//n_segments
stops= []
V_fits = np.zeros(n_segments, dtype=object)
I_fits = np.zeros(n_segments, dtype=object)
crrpopt = np.zeros(n_segments, dtype=object)
crrpcov = np.zeros(n_segments, dtype=object)
Ccrrpopt = np.zeros(n_segments, dtype=object)
Ccrrpcov = np.zeros(n_segments, dtype=object)
Ti_UN_for_plateau = np.zeros(n_segments, dtype=object)
Voltage_for_plateau = np.zeros(n_segments, dtype=object)
Ti_UN_err_for_plateau = np.zeros(n_segments, dtype=object)
voltage_cutoff = 0
for i in range(n_segments):
voltage_cutoff += step_indexes
if i == n_segments-1: ### dealing with edge to Make the fit til the end. (-1 cus range is 'missing' last value thus this cond would never apply)
stops.append(len(voltage))
stop = len(voltage)
else:
stops.append(voltage_cutoff)
stop = voltage_cutoff
# print(stop)
V_fits[i] = voltage[0:stop]
I_fits[i] = current[0:stop]
Voltage_for_plateau[i] = voltage[stop-1]
#UNWEIGHTED
try:
params_UN, covs_UN = curve_fit(Afterburner(plasmapotential), V_fits[i], I_fits[i], bounds=([Ti_lobound, R_lobound, Isat_lobound], [Ti_upboud, R_upbound, Isat_upbound]))
crrpopt[i]= params_UN
crrpcov[i]= covs_UN
except (RuntimeError,ValueError):
crrpopt[i] = [NaN,NaN,NaN]
crrpcov[i] = [NaN,NaN,NaN]
Ccrrpopt[i]= [NaN,NaN,NaN]
Ccrrpcov[i]= [NaN,NaN,NaN]
continue
conditions = [(Voltage_for_plateau[i]-plasmapotential) < 3*params_UN[0] , (np.sqrt(covs_UN[0][0])/params_UN[0]) < 0.6] ####ATTENTION CHOOSE Vp and TI FROM LEADING FIT + Voltage_for_plateau[i] is max voltage for particular fit#########
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])
try:
params2_UN, covs2_UN = q_curve_fit(boundss, Rzero(condition)(rrpopt[2]), V_fits[i], I_fits[i], bounds=boundss)
Ccrrpopt[i]= params2_UN
Ccrrpcov[i]= covs2_UN
Ti_UN_for_plateau[i]= params2_UN[0] #####CHANGE THIS################################################
Ti_UN_err_for_plateau[i]= covs2_UN[0]
except (RuntimeError,ValueError): # change this for the leading fit
Ccrrpopt[i] = NaN
Ccrrpcov[i] = NaN
pass
######### Final result
difference = []
abs_difference=[]
for j in range(0,len(Ccrrpopt)-1,1): #####change for weighted method
if (Ccrrpopt[j][0] < Ti_upboud-1) and ((Ccrrpcov[j][0]/Ccrrpopt[j][0]) < 0.6) and (0.9*Ccrrpopt[j][0]) < (Voltage_for_plateau[j]-plasmapotential):
deriv = (Ccrrpopt[j+1][0] - Ccrrpopt[j][0]) / (Voltage_for_plateau[j+1] - Voltage_for_plateau[j])
if abs(deriv) < 0.5:
difference.append(deriv)
abs_difference.append(abs(deriv))
else:
difference.append(NaN)
abs_difference.append(NaN)
else:
difference.append(NaN)
abs_difference.append(NaN)
# print('difference: '+str( difference )) ### What if there is no plaetau?
min_value = np.nanmin(abs_difference) ### finds deriv with smallest difference
try:
index = abs_difference.index(min_value)
total_index = abs_difference.index(min_value) ### finds an index of a min derivative
# print('first index = '+str(total_index))
except(ValueError):
total_index = len(Ccrrpopt)-2 ### if there are only nan values of derivs. we take the last couple of fits and decide which one has lower covariance
index = len(Ccrrpopt)-2
pass
if Ccrrpcov[total_index][0] < Ccrrpcov[total_index+1][0]: ####change this for weighted
Final_Ti =Ccrrpopt[total_index]
Final_Ti_err = Ccrrpcov[total_index]
else: ####change this for weighted
total_index = total_index+1
Final_Ti =Ccrrpopt[total_index]
Final_Ti_err = Ccrrpcov[total_index]
# print('total_index = '+str(total_index))
# print(Voltage_for_plateau[total_index])
if additional_graph:
fig, ax = plt.subplots(figsize = (8,6))
ax.errorbar(Voltage_for_plateau, Ti_UN_for_plateau, yerr=Ti_UN_err_for_plateau, color='red', Marker = 'o', markersize = 8, linestyle='None', label='Ti [eV]', zorder=0)
ax.plot(Voltage_for_plateau[total_index] , Ti_UN_for_plateau[total_index], 'o', color = 'green', markersize = 10)
# ax.plot(Voltage_for_plateau , faulty_fits, 'o', color = 'white', markersize = 10, alpha= 1)
# ax.plot(Voltage_for_diff , difference, 'D', color = 'magenta', markersize = 8, label= 'derivatives')
ax.plot(Voltage_for_diff[index] , difference[index], 'D', color = 'green', markersize = 8, label= 'minimal derivative')
ax.axhline(y=0, color='black')
ax.set_xlim(0,)
ax.set_ylim(-5,)
ax.legend(fontsize = 10)
return { 'Final_Ti':Final_Ti, 'Final_Ti_err':Final_Ti_err }
version 5.2
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 = 4
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 = (-150 < 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) - np.min(V_check) ) < MAX_voltage-10) or (np.max(V_check) < MAX_voltage-10) ):
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]) < 3*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 ,
(Ti_upboud-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,
(Ti_upboud-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)
elif (all(cond1) or abs(np.divide(np.sqrt(rrpcov[2][2]),rrpopt[2])) < rel_err): ### potential can not be changed by cuttoff thus no need to try
print('cutoff is trying to improve: error ' + str(Ccrrpcov[0]) + ' Ti estimate='+str(Ccrrpopt[0]))
try:
result = cutoff_fitting_for_profile_UNweighted(rI_short,rV_short,len(rV_short)//cutoff_resolution, additional_graph=False)
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_by_cuttoff')
iv_time.append(index_to_time)
continue
else:
Ti.append(result['Final_Ti'][0])
Ti_err.append(result['Final_Ti_err'][0])
R.append(result['Final_Ti'][1])
Vfl.append(rrpopt[2])
I_sat.append(result['Final_Ti'][2])
R_err.append(result['Final_Ti_err'][1])
Vfl_err.append(np.sqrt(rrpcov[2][2]))
I_sat_err.append(result['Final_Ti_err'][2])
index.append(iv_start)
recognition.append('cutoff_fit; precision: '+str(cutoff_resolution))
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)
# fig, ax = plt.subplots(figsize = (10,8))
# ax.set_title(str(shot_number)+' index: '+str(iv_start))
# 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 , '*', label = 'bins through ' +str(time_bin)+' IVs')
# #ax.errorbar(vbins , -ibins ,ierrs)
# x = np.linspace(np.nanmin(V_short) , np.nanmax(rV_short), num=len(rV_short)) ## noncutted IV
# xx = np.linspace(rrpopt2[2]- 0.1*rrpopt2[2], np.nanmax(V_short), num=len(V_short)) ## cutted IV
# ax.errorbar(V_short, -I_short, yerr=current_err, color='grey', label='sigma', Marker = '', linestyle='None')
# ##### FIRST Fits ###########
# #ax.plot(x, -ivchar_fit(x, *rrpopt), color='blue', 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.set_xlabel('Voltage [V]', fontsize=12)
# ax.set_ylabel('Current [A]', fontsize=12)
# ax.legend(loc='upper left', fontsize=8)
# print(R2_value)
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-1-266cf2fdc889> in <module> 43 ############################################################# 44 index_to_time = interval_pivots[iv_start] ---> 45 alpha_for_Te.append(1.89 * B_interpolated[nearest(time_ax[st:ed] , index_to_time)] + 1.85) 46 B_IV.append(B_interpolated[nearest(time_ax[st:ed] , index_to_time)]) 47 cut_V = (V_for_float_potential < V) & (V < MAX_voltage) NameError: name 'B_interpolated' is not defined
rel_err= 0.6
Ti_filt = []
Ti_err_filt = []
Vp_probe_filt = []
Vp_err_filt = []
time_filter = []
Isat_filter = []
Index_filter = []
V_BPP_floating_pot_filter = []
statement = []
B_IV_filt =[]
# alpha_bpp_filt =[]
i = 0
for i in range(0,len(Ti)):
rules = [np.divide(Ti_err[i],Ti[i]) < rel_err ,
#abs(np.divide(Vfl_err[i],Vfl[i])) < rel_err, ### we do not need this. Vp err handeled already
Ti[i]<(MAX_voltage - np.array(Vfl[i]))/1,
Ti[i]<Ti_upboud-1,
Ti_err[i]>0.0,
]
if all(rules):
Ti_filt.append(Ti[i])
Ti_err_filt.append(Ti_err[i])
Vp_probe_filt.append(Vfl[i])
Vp_err_filt.append(Vfl_err[i])
time_filter.append(iv_time[i])
Isat_filter.append(I_sat[i])
Index_filter.append(index[i])
statement.append(recognition[i])
# alpha_bpp_filt.append(alpha_bpp[i])
V_BPP_floating_pot_filter.append(V_BPP_floating_pot[i])
B_IV_filt.append(B_IV[i])
print('nuber of succesful fits = '+ str(len(Ti_filt))+' / ' +str(len(Ti)))
nuber of succesful fits = 0 / 0
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('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('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('Ti_' + str(radial_probe_position ) +'mm_'+str(shot_number)+'.png' )
2.7-0.25
2.45
shot_number_float = shot_number
lp = pd.read_csv(urlopen('http://golem.fjfi.cvut.cz/shots/' + shot_number +'/Diagnostics/PetiProbe/DAS_raw_data_dir/ch4.csv'),names = ['t','V'])
Vfl_LP = signal_cleaner(lp.V)
lo = nearest(time_ax, fromm)-100
hi = nearest(time_ax, untill)+100
Vfl_LP_smooth = my_lowpass_filter(Vfl_LP[lo:hi], 1 ,fs = FS, order = 2)
# Vfl_LP_smooth = Vfl_LP[lo:hi]
xxxx = scipy.interpolate.interp1d(time_ax[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)]
--------------------------------------------------------------------------- TypeError Traceback (most recent call last) <ipython-input-1-8d48a1c5e7d5> in <module> 1 shot_number_float = shot_number ----> 2 lp = pd.read_csv(urlopen('http://golem.fjfi.cvut.cz/shots/' + shot_number +'/Diagnostics/PetiProbe/DAS_raw_data_dir/ch4.csv'),names = ['t','V']) 3 Vfl_LP = signal_cleaner(lp.V) 4 lo = nearest(time_ax, fromm)-100 5 hi = nearest(time_ax, untill)+100 TypeError: can only concatenate str (not "int") to str
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('Potential_check' + str(radial_probe_position ) +'mm_'+str(shot_number)+'.png' )
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-1-e9a77cabbbb6> in <module> 2 # plt.plot(iv_time,V_BPP_floating_pot, 'blue', marker ='o',alpha = .4, label = 'Vfl BPP raw') 3 # plt.plot(iv_time,V_BPP_floating_pot_smooth, 'blue', marker ='o', label = 'Vfl BPP ') ----> 4 plt.plot(iv_time,V_p_raw, marker ='o',markersize = 4, color='green',alpha = .3, label = '$\Phi$ from BPP & LP') 5 plt.plot(time_filter,Vp_probe_filt, 'red', marker ='o', markersize = 7, ls ='--',alpha = .5, label = '$\Phi$ 4-p fit ') 6 NameError: name 'V_p_raw' is not defined
<Figure size 576x432 with 0 Axes>
BPP should not fluctuate more than lp. probably the noise from coax --> histeresis
plt.figure()
plt.plot(time_ax[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('Vfl' + str(radial_probe_position ) +'mm_'+str(shot_number)+'.png' )
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-1-8cfd2673b85d> in <module> 1 plt.figure() ----> 2 plt.plot(time_ax[lo:hi], Vfl_LP[lo:hi], linewidth = 3, color = 'blue',alpha = .4, label = '$V_{fl}^{LP}$ [V]') 3 plt.plot(iv_time,V_BPP_floating_pot,'-o', linewidth = 1.5, color = 'red',alpha = .4, label = '$\Phi_{fl}^{BPP}$ [V]') 4 plt.plot(iv_time,V_LP_floating_pot, linewidth = 3, color = 'blue', label = 'smooth $V_{fl}^{LP}$ [V]') 5 plt.plot(iv_time,V_BPP_floating_pot_smooth,'-', linewidth = 3, color = 'orange', label = 'smooth $\Phi_{fl}^{BPP}$ [V]') NameError: name 'lo' is not defined
<Figure size 432x288 with 0 Axes>
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()
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-1-06e5b6d35e1e> in <module> 1 plt.figure() ----> 2 plt.plot(times ,Tes, '-',label = 'Electron temperature [eV]' ) 3 plt.plot(times,Te_smooth, linewidth = 3, color = 'orange', label = 'Electron temperature smooth [eV]') 4 plt.ylim(0,40) 5 plt.legend() NameError: name 'times' is not defined
<Figure size 432x288 with 0 Axes>
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,110)
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('Ti_' + str(radial_probe_position ) +'mm_'+str(shot_number)+'.png' )
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-1-02966dae2b4f> in <module> 2 ax.plot([],[],ls='None', label = 'OH-regime; $R = $'+str(radial_probe_position) + ' mm') 3 # ax.plot(times, Tes, 'o',markersize = 1.5,label = 'electron temperature #'+ str(shot_number_float)) ----> 4 ax.plot(times, Te_smooth,color = 'blue', linewidth =3, label = '$T_e$ smoothed #'+ str(shot_number_float)) 5 ax.plot(time_filter, Ti_filt, 'o',markersize = 5 , color = 'red' , label ='$T_i$ 4-p fit #' + str(shot_number)) 6 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, ) NameError: name 'times' is not defined
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, '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, 'Ti_BPP_profile_2ms_bin_averaged_Weighted_NewLP'+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, 'Te_'+str(shot_number)+'.txt')
--------------------------------------------------------------------------- NameError Traceback (most recent call last) <ipython-input-1-398f37b3fb11> in <module> 12 13 #$$$$$$$$$ ULOZENIE Te DAT $$$$$$$############### ---> 14 Te_p_fil = np.ones(len(times)) 15 Te_p_filt = Te_p_fil-1 + radial_probe_position 16 from astropy.table import QTable, Table, Column NameError: name 'times' is not defined
# time_start = 9.5
# time_end = 10
# iv_start = nearest(interval_pivots,time_start)
# iv_terminate = nearest(interval_pivots,time_end)
# j = iv_start
# 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 = []
# #R2arr = []
# for j in range(iv_start,iv_terminate, 1):
# iv_start = j
# time_bin = 10 # periods for deviation estimation
# lobound = [0., 0., 0., -1.] # fit lower bound for Ti, R, Vf, Isat
# upbound = [150., 1, 80., 0.] # fit upper bound for Ti, R, Vf, Isat
# # ------------------------------------------------------------------------------------------------------------
# iv_stop = iv_start + resolution #+ 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])
# cut_V = (0 < V) & (V < 100)
# 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]
# ##### Calculate the error bars from binned data and assign to raw data
# try:
# binsize = 5. # size of the voltage intervals for the calculation of the standard deviation
# 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
# 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]
# except (RuntimeError,ValueError): # for the graph is dark and full of errors
# Ti.append(-10)
# Ti_err.append(-10)
# Vfl.append(-10)
# R.append(-10)
# I_sat.append(-10)
# R_err.append(-10)
# Vfl_err.append(-10)
# I_sat_err.append(-10)
# index.append(-10)
# continue
# try:
# popt2, pcov2 = curve_fit(ivchar_fit, V_short, I_short, bounds=(lobound, upbound), sigma=current_err, absolute_sigma=False)
# rpopt2, rpcov2 = curve_fit(ivchar_fit, V_short, I_short, bounds=([0., 0., popt2[2]-10, -10.], [150., 1, popt2[2]+10, 10.]), sigma=current_err, absolute_sigma=False)
# rrpopt2, rrpcov2 = curve_fit(ivchar_fit, V_short, I_short, bounds=([0., 0., rpopt2[2]-5, -10.], [150., 1, rpopt2[2]+5, 10.]), sigma=current_err, absolute_sigma=False)
# except (RuntimeError,ValueError): # for the graph is dark and full of errors
# Ti.append(-10)
# Ti_err.append(-10)
# Vfl.append(-10)
# R.append(-10)
# I_sat.append(-10)
# R_err.append(-10)
# Vfl_err.append(-10)
# I_sat_err.append(-10)
# index.append(-10)
# # R2arr.append(-10)
# continue
# rcut_V = ((rrpopt2[2] - 0.1*rrpopt2[2]) < V) & (V < 100)
# 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]
# ##### Calculate the error bars from binned data and assign to raw data
# try:
# rbinned_data = bin_average(rI_long, rV_long, binsize) # Compute the standard deviations for each bin
# rverrs = rbinned_data['V_err'] # voltage stds
# rierrs = rbinned_data['I_err'] # current stds
# rvbins = rbinned_data['V'] # un-center the bins
# 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, rvbin in enumerate(rvbins): # iterate on the bins
# r_current_err[rV_short >= rvbin] = rierrs[i]
# r_voltage_err[rV_short >= rvbin] = rverrs[i]
# except (RuntimeError,ValueError): # for the graph is dark and full of errors
# Ti.append(-11)
# Ti_err.append(-11)
# Vfl.append(-11)
# R.append(-11)
# I_sat.append(-11)
# R_err.append(-11)
# Vfl_err.append(-11)
# I_sat_err.append(-11)
# index.append(-11)
# continue
# try:
# cpopt2, cpcov2 = curve_fit(Afterburner(rrpopt2[2]), rV_short, rI_short, bounds=([0., 0., -1.], [150., 1, 0.]), sigma=r_current_err, absolute_sigma=False)
# mean_of_data=(1/len(rI_short))*np.sum(-rI_short) #######np.mean(-rI_short)########
# SS_tott= np.sum((-rI_short-mean_of_data)**2)
# SS_ress=np.sum((-rI_short+ivchar_fit_Weighted(rV_short, *cpopt2))**2)
# R2_value = 1-(SS_ress/SS_tott)
# except (RuntimeError,ValueError): # for the graph is dark and full of errors
# Ti.append(-11)
# Ti_err.append(-11)
# Vfl.append(-11)
# R.append(-11)
# I_sat.append(-11)
# R_err.append(-11)
# Vfl_err.append(-11)
# I_sat_err.append(-11)
# index.append(-11)
# #R2arr.append(-11)
# continue
# try:
# condition = (85-rrpopt2[2]) < 3*cpopt2[0]
# if condition :
# boundss = ([0., -1.], [150., 0.])
# else:
# boundss = ([0., 0., -1.], [150., 1, 0.])
# Ccrrpopt2, Ccrrpcov2 = q_curve_fit(boundss, Rzero(condition)(rrpopt2[2]), rV_short, rI_short, bounds=boundss, sigma= r_current_err, absolute_sigma=False)
# except (RuntimeError,ValueError): # for the graph is dark and full of errors
# Ti.append(-12)
# Ti_err.append(-12)
# Vfl.append(-12)
# R.append(-12)
# I_sat.append(-12)
# R_err.append(-12)
# Vfl_err.append(-12)
# I_sat_err.append(-12)
# index.append(-11)
# # R2arr.append(-12)
# else:
# Ti.append(Ccrrpopt2[0])
# Ti_err.append(Ccrrpcov2[0])
# R.append(Ccrrpopt2[1])
# Vfl.append(rrpopt2[2])
# I_sat.append(Ccrrpopt2[2])
# R_err.append(Ccrrpcov2[1])
# Vfl_err.append(np.sqrt(rrpcov2[2][2]))
# I_sat_err.append(Ccrrpcov2[2])
# index.append(iv_start)
# # R2arr.append(R2_value)
# fig, ax = plt.subplots(figsize = (8,6))
# ax.set_title(str(shot_number)+' index: '+str(iv_start))
# 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 , '*', label = 'bins through ' +str(time_bin)+' IVs')
# #ax.errorbar(vbins , -ibins ,ierrs)
# x = np.linspace(np.nanmin(V_short) , np.nanmax(rV_short), num=len(rV_short)) ## noncutted IV
# xx = np.linspace(rrpopt2[2]- 0.1*rrpopt2[2], np.nanmax(V_short), num=len(V_short)) ## cutted IV
# ax.errorbar(V_short, -I_short, yerr=current_err, color='grey', label='sigma', Marker = '', linestyle='None')
# ##### FIRST Fits ###########
# #ax.plot(x, -ivchar_fit(x, *rrpopt), color='blue', 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.set_xlabel('Voltage [V]', fontsize=12)
# ax.set_ylabel('Current [A]', fontsize=12)
# ax.set_xlim(0,80)
# ax.legend(loc='upper left', fontsize=8)