# # 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']
